// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ funcprot(0); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cv=%pi/180; // From degree to radian // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Latitude=52; // [°] L=30; // Length of pendulum M=1; // Initial data alpha=45; // [°] beta=0; // [°] Initial position of pendulum v0v=0.5; // Initial velocity x0=[ L*sin(alpha*cv)*cos(beta*cv);.. L*sin(alpha*cv)*sin(beta*cv);.. -L*cos(alpha*cv)]; // Initial pendulum position // disp(x0); v0=[ -v0v*cos(alpha*cv)*cos(beta*cv);.. -v0v*cos(alpha*cv)*sin(beta*cv);.. -v0v*sin(alpha*cv)]; // Initial pendulum velocity // disp(v0); t0=0; tend=1800; // 180 steps=1800; // 800 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // phi1 = geodedic latitude phi1=Latitude*cv; // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // ω (radians/s) = earth's angular spin rate omega=7.292115E-5; // omega=0; // ##################################################################### // Vector of coriolis acceleration // phi1 = geodedic (geographic) latitude function rbC=fbC(v) v1=v(1); v2=v(2); v3=v(3); sn=sin(phi1); cs=cos(phi1) rbC(1)= 2*omega*(-v2*sn+v3*cs); rbC(2)= 2*omega*( v1*sn); rbC(3)=-2*omega*(-v1*cs); endfunction // ##################################################################### function rbG=fbG() rbG(1)=0; rbG(2)=0; rbG(3)=-9.80616; // g45 endfunction // ##################################################################### // ##################################################################### function rvDot=fvDot(x,v) // return ph data x1=x(1); x2=x(2); x3=x(3); v1=v(1); v2=v(2); v3=v(3); r=sqrt(x1^2+x2^2+x3^2); vel2=v1^2+v2^2+v3^2; xu=x/r; bGC=fbG()+fbC(v); bGCDOTxu=bGC(1)*xu(1)+bGC(2)*xu(2)+bGC(3)*xu(3); bZWv=bGCDOTxu+(vel2/r); bZW=-bZWv*xu; b=bGC+bZW; rvDot(1:3)=v; rvDot(4:6)=b; endfunction // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function rPendulum=fPendulum(t,ph) // returns 'phDot data' x=ph(1:3); v=ph(4:6); rPendulum=fvDot(x,v); endfunction; // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ph0(1:3)=x0; ph0(4:6)=v0; tData=linspace(t0,tend,steps); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ phData=ode(ph0,t0,tData,fPendulum); // disp(phData); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ to_deg=180/%pi; // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ PhiData=[]; nbpoints=size(phData,"c"); for Ix=1:nbpoints PhiData(Ix)=atan(phData(2,Ix)/phData(1,Ix)); // disp(PhiData(Ix)); end RadData=[]; for Ix=1:nbpoints RadData(Ix)=sqrt(phData(1,Ix)^2+phData(2,Ix)^2+phData(3,Ix)^2); // disp(RadData(Ix)); end // disp(RadData); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ my_handle=scf(1); clf(my_handle,"reset"); a=gca(); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //drawlater(); // //xtitle ( " Radius der Pendelbewegung " ,.. // " Anzahl der Integrationsschritte", .. // " Radius " ); // //a.labels_font_size=3; //a.x_label.font_size = 3; //a.y_label.font_size = 3; //a.title.font_size = 3; // //plot2d(phData(1,:)', phData(3,:)'); // //drawnow(); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ drawlater(); xtitle ( " Radius der Pendelbewegung " ,.. " Anzahl der Integrationsschritte", .. " Radius " ); a.labels_font_size=3; a.x_label.font_size = 3; a.y_label.font_size = 3; a.title.font_size = 3; plot(tData,RadData); drawnow(); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //drawlater(); // //xtitle ( " Azimuthwinkel der Pendelbewegung " ,.. // " Anzahl der Integrationsschritte", .. // " Grad " ); // //a.labels_font_size=3; //a.x_label.font_size = 3; //a.y_label.font_size = 3; //a.title.font_size = 3; // //plot(tData,PhiData*to_deg); // // //drawnow(); // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++