# Kapitel 3_m: 
# Der symmetrische schwere Spielkreisel (Eulerwinkel)
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plottools): with(plots):
# Wir beschaffen uns zunächst die  Bewegungsgleichungen des  schweren Kreisels in Eulerkoordinaten. Dazu notieren wir die Winkelgeschwindigkeiten É im körperfesten Koordinatensystem, die dann in den Drallsatz eingesetzt werden. 
> omega[1](t):=diff(psi(t),t)*sin(theta(t))*sin(phi(t))+diff(theta(t),t)*cos(phi(t));
> omega[2](t):=diff(psi(t),t)*sin(theta(t))*cos(phi(t))-diff(theta(t),t)*sin(phi(t));
> omega[3](t):=diff(psi(t),t)*cos(theta(t))+diff(phi(t),t);
# Die Eulerschen Gleichungen für den schweren symmetrischen Kreisel sind (s.h. Buch TM 3, Kap. 4.1 ) 
> gl1:=expand(A*diff(omega[1](t),t)-(A-C)*omega[2](t)*omega[3](t)-G*s*sin(theta(t))*cos(phi(t)));
> gl2:=expand(A*diff(omega[2](t),t)-(C-A)*omega[1](t)*omega[3](t)+G*s*sin(theta(t))*sin(phi(t)));
> gl3:=diff(omega[3](t),t);
# Diese Gleichungen lassen sich noch vereinfachen. Dazu multiplizieren wir gl1 und gl2  wechselseitig mit cosÆ  bzw. sinÆ  und ziehen die dadurch entstandenden Gleichungen voneinander ab bzw. addieren diese. Das reduziert die neuen Bewegungsgleichungen erheblich. 
> C:=zeta*A: G:=sigma*A/s:
> g1:=gl1*cos(phi(t)):
> g2:=gl2*sin(phi(t)):
> simplify(g1-g2,trig)/A:
> gl1n:=collect(%,sin(theta(t)));
> g1:=gl1*sin(phi(t)):
> g2:=gl2*cos(phi(t)):
> simplify(g1+g2,trig)/A:
> gl2n:=collect(%,diff(theta(t), t));
# Zur nummerischen Lösung dieses nichtlinearen  Differenzialgleichungssystems stellen wir eine Prozedur zur Verfügung
> Proc_Calc_22:=proc(T1,T3,G,S,AW::list,N::integer,TE,POLY)
> #---------------------------------------------------------------------------
> #Eingabe: 
> #          T1: Massenträgheitsmoment (A = B)               in [kg*m^2]
> #          T3: Massenträgheitsmoment (C, Figurenachse)     in [kg*m^2]
> #          G : Gewicht des Kreisels                        in [kg*m/s^2]
> #          S : Schwerpunktskoordinate (x3S)                in [m]
> #          AW: Anfangswerte in Eulerkoordinaten und deren
> #              Ableitungen zum Zeitpunkt t = 0
> #              AW[1]: Präzessionswinkel psi0               in [rad]
> #              AW[2]: d(psi0)/dt                           in [rad/s]
> #              AW[3]: Nutationswinkel theta0               in [rad]
> #              AW[4]: d(theta0)/dt                         in [rad/s]
> #              AW[5]: Eigenrotationswinkel phi0 	           in [rad]                              
> #              AW[6]: d(phi0)/dt                           in [rad/s]  
> #          N : Anzahl der Simulationsschritte
> #          TE: Dauer der Simulationszeit	                   in [s]
> #        POLY: 3D-Geometrieobjekt des Kreisels
> #Ausgabe:  
> #          Eulerwinkel           
> #            q1: Präzessionswinkel    psi(t)               in [rad]
> #            q2: Nutationswinkel    theta(t)               in [rad]
> #            q3: Eigenrotationswinkel phi(t)               in [rad]
> #          Eulerwinkelgeschwindigkeiten
> #            q4: d( psi(t))/dt                             in [rad/s]
> #            q5: d(theta(t)/dt                             in [rad/s]
> #            q6: d( phi(t))/dt                             in [rad/s]
> #          Winkelgeschwindigkeiten omega                     
> #            q7: omega1                                    in [rad/s]
> #            q8: omega2                                    in [rad/s]
> #            q9: omega3                                    in [rad/s]
> #-------------------------------------------------------------------------
> local zeta,sigma,gl1,gl2,gl3,AB,sys,i,val,p1,p2,p3,NT,j,phi,psi,theta,P,Q;
> global L,ZG,AN,T,q1,q2,q3,q4,q5,q6,q7,q8,q9;
> description "Lösung der Bewegungsgleichungen und Animation der Bewegung eines symmetrischen schweren Kreisels";
> zeta:=T3/T1;sigma:=G*S/T1;
> gl1:=((diff(psi(t),t))^2*cos(theta(t))*zeta-(diff(psi(t),t))^2*cos(theta(t))+(diff(psi(t),t))*(diff(phi(t),t))*zeta-sigma)*sin(theta(t))+diff(theta(t),t,t);
> gl2:=(-(diff(psi(t),t))*cos(theta(t))*zeta+2*(diff(psi(t),t))*cos(theta(t))-(diff(phi(t),t))*zeta)*(diff(theta(t),t))+(diff(psi(t),t,t))*sin(theta(t));
> gl3:=(diff(psi(t),t,t))*cos(theta(t))-(diff(psi(t),t))*sin(theta(t))*(diff(theta(t),t))+diff(phi(t),t,t);
> #Anfangsbedingungen
> AB:=psi(0)= AW[1], D(psi)(0)= AW[2], theta(0)= AW[3], D(theta)(0)= AW[4], phi(0)= AW[5], D(phi)(0)= AW[6];
> sys:={gl1,gl2,gl3,AB};
> val:=Array(1..N+1); for i to N+1 do val[i]:=(i-1)*TE/N end do:
> L:=dsolve(sys,numeric,method=rkf45,output=val);
> #Die zeitabhängigen Eulerwinkel
> q1:=plots[odeplot](L,[t,  psi(t)],0..TE,color=   red):
> q2:=plots[odeplot](L,[t,theta(t)],0..TE,color=  blue):
> q3:=plots[odeplot](L,[t,  phi(t)],0..TE,color=yellow):
> #Die zeitabhängigen Eulerwinkelgeschwindigkeiten
> q4:=plots[odeplot](L,[t,diff(psi(t)  ,t)],0..TE,color=   red):
> q5:=plots[odeplot](L,[t,diff(theta(t),t)],0..TE,color=  blue):
> q6:=plots[odeplot](L,[t,diff(phi(t)  ,t)],0..TE,color=yellow):
> #Die Winkelgeschwindigkeiten omega im körperfesten Koordinatensystem
> q7:=plots[odeplot](L,[t,diff(psi(t),t)*sin(theta(t))*sin(phi(t))+diff(theta(t),t)*cos(phi(t))],0..TE,color=   red):
> q8:=plots[odeplot](L,[t,diff(psi(t),t)*sin(theta(t))*cos(phi(t))-diff(theta(t),t)*sin(phi(t))],0..TE,color=  blue):
> q9:=plots[odeplot](L,[t,diff(psi(t),t)*cos(theta(t))+diff(phi(t),t)]                          ,0..TE,color=yellow):
> #Grafische Darstellung der Ergebnisse 
> p1:=plots[display]([q1,q2,q3],axes=boxed,gridlines=true,title="\nEulerwinkel\npsi (rot), theta (blau), phi (gelb)\n",labels = ["t"," "]):
> p2:=plots[display]([q4,q5,q6],axes=boxed,gridlines=true,title="\nEulerwinkelgeschwindigkeiten\npsip (rot), thetap (blau), phip (gelb)\n",labels = ["t"," "]):
> p3:=plots[display]([q7,q8,q9],axes=boxed,gridlines=true,title="\nWinkelgeschwindigkeiten omega\nomega1 (rot), omega2 (blau), omega3 (gelb)\n",labels = ["t"," "]):
> ZG:=plots[display](Matrix(1,3,[p1,p2,p3]));ZG;
> #Animation des Bewegungsvorganges
> T:=NULL:
> NT:=LinearAlgebra[RowDimension](L[2,1]);
> for j to NT do  #Schleife über jeden Zeitschritt
>   phi:=L[2,1][j,2]; psi:=L[2,1][j,4]; theta:=L[2,1][j,6];
>   P:=plottools[rotate](plottools[rotate](plottools[rotate](POLY,0,0,-phi),-theta,0,0),0,0,-psi);
>   Q:=plots[display](P): T:= T,Q:
> end do:
> AN:=plots[display](T,insequence=true,scaling =constrained,orientation=[45, 45, 0],labels=[x1,x2,x3]);
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab. Das geht wie folgt:
> save Proc_Calc_22, "Proc_Calc_22.m":
# Zur Berechnung der massengeometrischen Größen des Kreisels benötigen wir folgende Prozedur:
> read "Proc_Calc_18.m";
# Beispiel 3-24:
# Wir ermitteln zunächst die massengeometrischen Größen des Spielkreisels:
> gl:=piecewise(x3 < h1,r1/h1*x3,x3 < h1+h2 ,r2,x3 < h1+h2+h3,r3);
> a:=unapply(gl,x3);
# Alle Längenmaße werden in cm eingegeben. 
> h1:=2: h2:=h1/2: h3:= 6/5*h1: r1:=h1: r2:=2*h1: r3:=h1/4:
> p:=plot3d(a(z),phi=0..2*Pi,z=0..h1+h2+h3,coords=cylindrical,orientation=[160,85,-130],scaling=constrained,labels=["x1","x2","x3"]): p;
> x1:=r*cos(phi): x2:=r*sin(phi):
> Proc_Calc_18([x1,x2,x3],[r,phi,x3],[[0,a(x3)],[0,2*Pi],[0,h1+h2+h3]],rho);
# Wir geben die nummerischen Werte für Á = 1 g/cm3  aus:
> rho:=1.:
> V:=evalf(V); m:=evalf(m); s:=evalf(rS[3]);  A:=evalf(Thetaq[1,1]); C:=evalf(Thetaq[3,3]); zeta:=evalf(C/A); G:=evalf(m*981); sigma:=evalf(G*s/A);
> TITEL:="Der symmetrische schwere Spielkreisel":
> x1:='x1':x2:='x2':
> PLS:=plots[display](p,scaling=constrained,axes=boxed,labels=["x1","x2","x3"],title=TITEL,font=[title,"ROMAN",15]):
# Wir wählen folgende Simulationsparameter
> N:=100:  #Anzahl der Simulationsschritte
;
> TE:=1.:  #Simulationszeit in [s]
;
# Der  allgemeine Fall: Es gelten entsprechend der Aufgabenstellung folgende Anfangswerte:
> psi0:=0.; psi0p:=0.; theta0:=evalf(Pi/6); theta0p:=0.; phi0:=0.; phi0p:=30.; 
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p];
# Wir lösen das allgemeine Anfangswertproblem und messen dazu die benötigte cpu-Zeit 
> settime:=time():                     #Beginn der Zeitmessung
;
> Proc_Calc_22(A,C,G,s,AW,N,TE,PLS):   #Das kann etwas dauern
;
> cpu_time:=(time()-settime)*seconds;  #Benötigte Rechenzeit
;
> ZG;AN;
# Wir prüfen, ob mit den gegebenen System- und Anfangswerten eine reguläre Präzession möglich ist.  Das erfordert folgende Grenzwerte der Präzessionswinkelgeschwindigkeiten  
> psi0p:='psi0p': 
> psip_gr:=solve((zeta-1)*psi0p^2*cos(theta0)+zeta*psi0p*phi0p-sigma,psi0p);
# Der waagerechte Kreisel
> theta0:=evalf(Pi/2); theta0p:=0.; phi0:=0.; phi0p:=30.; psi0:=0.; psi0p:=sigma/(zeta*phi0p);
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p];
> Proc_Calc_22(A,C,G,s,AW,N,TE,PLS): ZG; AN;
> psi0p:='psi0p':
# Wir setzen neue Anfangswerte, die die reguläre Präzession erzwingen.
# 1.) Die langsame Präzession
> psi0:=0.; psi0p:=psip_gr[1]; theta0:=evalf(Pi/6); theta0p:=0.; phi0:=0.; phi0p:=30.; 
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p];
> Proc_Calc_22(A,C,G,s,AW,N,TE,PLS): ZG; AN;
# 2.) Die schnelle Präzession
> psi0:=0.; psi0p:=psip_gr[2]; theta0:=evalf(Pi/6); theta0p:=0.; phi0:=0.; phi0p:=30.; 
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p];
> Proc_Calc_22(A,C,G,s,AW,N,TE,PLS): ZG; AN;
# 3.) Die mittlere Präzession
> psi0:=0.; psi0p:=(psip_gr[1]+psip_gr[2])/2.; theta0:=evalf(Pi/6); theta0p:=0.; phi0:=0.; phi0p:=30.; 
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p];
> Proc_Calc_22(A,C,G,s,AW,N,TE,PLS): ZG; AN;
#  
# Hinweis: Wir hätten die Lösung auch sofort mittels der bereits anderenorts bereitgestellten Prozedur  Proc_Calc_20  beschaffen können, die dazu eingelesen wird.
> read "Proc_Calc_20.m";
# Dazu benötigen wir das Moment aus dem Eigengewicht des Kreisels bezüglich des Drehpunktes (s.h. Buch TM 3, Gl. 3.70)
> M:=G*rS[3]*Vector([sin(theta(t))*cos(phi(t)),-sin(theta(t))*sin(phi(t)),0]);
# Wir wählen folgende Simulationsparameter
> N:=100:  #Anzahl der Simulationsschritte
;
> TE:=1.:  #Simulationszeit in [s]
;
# Wir behandeln hier nur den allgemeinen Fall.Es gelten entsprechend der Aufgabenstellung folgende Anfangswerte:
> psi0:=0.; psi0p:=0.; theta0:=evalf(Pi/6); theta0p:=0.; phi0:=0.; phi0p:=30.; 
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p];
# Wir lösen das allgemeine Anfangswertproblem und messen dazu die benötigte CPU-Zeit 
> settime:=time():                                 #Beginn der Zeitmessung
;
> Proc_Calc_20(Vector([A,A,C]),M,AW,N,TE,PLS):     #Das kann etwas dauern
;
> cpu_time:=(time()-settime)*seconds;              #Benötigte Rechenzeit
;
> ZG;AN;
> 
;
# Hinweis: Zur Berechnung des stehenden oder hängenden Kreisels sind Sonderbetrachtungen erforderlich (s.h. Buch TM 3), die in  Proc_Calc_22  nicht abgefangen werden.
