# Kapitel 3_j: 
# Der unsymmetrische Pyramidenkreisel unter Momentenbelastung
# (Eulerwinkel)
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plottools): with(plots):
# Wir lesen die Prozedur zur Berechnung der massengeometrischen Größen des Kreisels ein:
> read "Proc_Calc_19.m"; 
# Beispiel 3-21:  Mit den Werten des Beispiels  ergeben sich folgende Systemwerte (Längen in [cm], Masse in [Gramm, g]).
> a:=6.; b:=3.; h:=9.; rho:=0.670; 
# Wir benötigen zur Beschreibung des Pyramidenkreisel eine Knoten- und eine Elementdatei.  Dazu beschreiben wir die Pyramide durch ihre vier  Eckpunkte und fassen diese in der Matrix COORD zusammen. 
> r1:=[0,0,0]: r2:=[a,0,h]:  r3:=[0,b,h]: r4:=[0,0,h]: 
> COORD:=Matrix([r1,r2,r3,r4]);
# Die  Elementmatrix (Koinzidenztafel) enthält die Zuordnung von jeweils drei Knoten zu einem Oberflächendreieck des Tetraeders . Es liegen insgesamt 4 Dreieckselemente vor.
> IZES:=Matrix([[1,3,2],[1,4,3],[1,3,2],[2,3,4]]);
# Mittels der folgenden Maple-Prozedur erzeugen wir für unseren Kreisel  eine global vereinbarte  3D-Plotstruktur mit dem Namen PL, die wir auch zur Animation des Bewegungsvorganges benötigen. 
> Proc_Graf_02:=proc(KNOTEN,ELEM,TEXT,LBEZ)
> #------------------------------------------------
> #Eingabe: 
> #          Knoten: Knotendatei  des Polyeders               
> #          ELEM  : Elementdatei des Polyeders
> #          TEXT  : Titel der Berechnung
> #          LBEZ  : Achsbezeichnungen
> #Ausgabe:  
> #              PL: 3D-Plotstruktur des Polyeders
> #------------------------------------------------
> local T,n,A,B,C,P;
> global PL;
> description "Grafische Darstellung eines Polyeders";
> T:=NULL:
> for n to RowDimension(ELEM) do
>   A:=[seq(KNOTEN[ELEM[n,1],i],i=1..3)];
>   B:=[seq(KNOTEN[ELEM[n,2],i],i=1..3)];
>   C:=[seq(KNOTEN[ELEM[n,3],i],i=1..3)];
>   P:=plots[polygonplot3d]([A,B,C]):
>   T:=T,P;
> end do;
> PL:=plots[display](T,scaling=constrained,axes=boxed,labels=LBEZ,title=TEXT,font=[title,"ROMAN",15],transparency=0.4);
> PL;
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab. Das geht wie folgt:
> save Proc_Graf_02, "Proc_Graf_02.m":
# Wir schauen uns den Kreisel an:
> TITEL:="Der unsymmetrische Pyramidenkreisel im körperfesten Bezugssystem": ABEZ:=["x1","x2","x3"]:
> Proc_Graf_02(COORD,IZES,TITEL,ABEZ);
# Wir berechnen die massengeometrischen Größen des Kreisels
> Proc_Calc_19(COORD,IZES,rho):
# Zur Lösung der nichtlinearen Bewegungsgleichungen benötigen wir die axialen Hauptträgheitsmomente und die Hauptträgheitsachsen bezogen auf den Ursprung des körperfesten Punktes P1, die wir der folgenden Eigenwertberechnung entnehmen. 
> (EWq,EVq):=LinearAlgebra[Eigenvectors](Matrix(Thetaq,shape=symmetric));
# Hinweis: Die Eigenvektoren sind nur bis auf einen konstanten Faktor bestimmt. Bilden die spaltenweise in der Eigenvektormatrix  EVq  erscheinenden und auf die Länge 1 normierten Eigenvektoren ein Rechtssystem, dann muss die Determinante der Eigenvektormatrix  den Wert +1 besitzen. Ergibt sich für die Determinante der Wert  -1, dann kann ein Rechtsystem gebildet werden, indem aus den ersten beiden Spaltenvektoren e1 und e2 ein dritter Vektor e3 = e1×e2  gebildet wird. Mit den Eigenvektoren sind nämlich nur die Achsrichtungen und nicht die Achsorientierungen festgelegt, da neben ei auch -ei (i = 1,2,3) ein Eigenvektor ist.  Wir überprüfen das und erhalten
> DET:=Determinant(EVq);
# Hinweis: 
# Zur Bildung der Hauptachsenbasis gehen wir hier etwas anders vor. Wir wählen  Column (EVq[1])  als Figurenachse und achten auf zyklische Vertauschung der beiden verbleibenden Achsen, damit das ursprünglich rechtshändige Basissystems erhalten bleibt, also
> EVm:=Matrix([Column(EVq,2),Column(EVq,3),Column(EVq,1)]);
# Die zugehörigen Eigenwerte (Hauptträgheitsmomente) sind:
> EWm:=Vector([EWq[2],EWq[3],EWq[1]]);
# Wir prüfen, ob es sich bei den obigen Einheitsvektoren um eine rechtshändige Basis handelt. Der resultierende Vektor muss dann verschwinden.
> CrossProduct(Column(EVm,1),Column(EVm,2))-Column(EVm,3);
# Wir beschaffen uns noch die Darstellung des Tetraeders mit Lage und Richtung der Hauptachsen. Der Ursprung des Hauptachsensystems ist identisch mit dem Ursprung des Bezugsachsensystems (Punkt P1 in der folgenden Grafik). 
> TITEL:="Der unsymmetrische Pyramidenkreisel im raumfesten Bezugssystem mit Hauptachsendarstellung":
> Proc_Graf_02(COORD,IZES,TITEL,ABEZ):
> lambda:=5.:li:={}:T1:={}:
> for i to 3 do
> evi:=Column(EVm,i):
> rai:=convert(-lambda*evi,list):
> rbi:=convert( lambda*evi,list):
> li:=li union {plottools[line](rai,rbi,color=blue,linestyle = dash)}:
> T1:=T1 union {[rai[1],rai[2],rai[3],i]}:
> end do:
> TP1:=plots[textplot3d](T1,'font'=["times","roman",14],'align'={'below','left'},color=blue):
> TP2:=plots[textplot3d]([0,0,0,"P1",'font'=["times","roman",16]],'align'={'below','left'},color=black):
> plots[display](PL,li,TP1,TP2,axes=boxed,scaling=constrained,orientation =[-75,55,-5],labels=["x1","x2","x3"]);
# Die rechtshändig orientierten Eigenvektoren in der Eigenvektormatrix EVm bilden die Matrixkomponenten der Drehmatrix (s.h. Buch TM III). Im gedrehten Hauptachsensystem sind dann die Eckpunkte des Tetraeders zeilenweise in der Matrix CH angeordnet. 
> CH:=COORD.EVm;
# Für die Koordinaten des Schwerpunktes im körperfesten Hauptachsensystem gilt
> RSH:=Transpose(EVm).RS;
# Wir stellen das Tetraeder im gedrehten Hauptachsensystem grafisch dar.
> TITEL:="Der unsymmetrische Pyramidenkreisel im Hauptachsensystem": ABEZ:=["x1H","x2H","x3H"]:
> Proc_Graf_02(CH,IZES,TITEL,ABEZ);
# Im nächsten Schritt soll der Bewegungsvorgang des unsymmetrischen Kreisels durch nummerische Integration der maßgebenden Differenzialgleichungen berechnet werden. Wir wählen zur Lagebeschreibung Eulerwinkel und stellen zur nummerischen Lösung des gekoppelten nichtlinearen  Differenzialgleichungssystems folgende Prozedur zur Verfügung. 
> Proc_Calc_20:=proc(TR::Vector,M::Vector,AW::list,N::integer,TE,POLY)
> #---------------------------------------------------------------------------
> #Eingabe: 
> #          TR: Massenträgheitsmomente 
> #              TR[1]: Trägheitsmoment (A)
> #              TR[2]: Trägheitsmoment (B)
> #              TR[3]: Trägheitsmoment (C, Figurenachse)
> #           M: Momente in Komponenten des körperfesten Hauptachsensystems
> #              M[1]: Moment um die 1-Achse
> #              M[2]: Moment um die 2-Achse
> #              M[3]: Moment um die 3-Achse            
> #          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 eta,zeta,gl1,gl2,gl3,AB,sys,i,val,p1,p2,p3,NT,j,w1,w2,w3,P,Q;
> global L,ZG,AN,T,q1,q2,q3,q4,q5,q6,q7,q8,q9;
> description "Lösung der Bewegungsgleichungen in Eulerwinkeln und Animation der Bewegung eines unsymmetrischen Kreisels";
> eta:=TR[2]/TR[1]; zeta:=TR[3]/TR[1];
> gl1:=(diff(psi(t),t,t))*sin(theta(t))*sin(phi(t))+(diff(psi(t),t))*(diff(theta(t),t))*cos(theta(t))*sin(phi(t))+(diff(psi(t), t))*sin(theta(t))*(diff(phi(t),t))*cos(phi(t))+(diff(theta(t),t,t))*cos(phi(t))-(diff(theta(t),t))*(diff(phi(t),t))*sin(phi(t))+(zeta-eta)*((diff(psi(t),t))*sin(theta(t))*cos(phi(t))-(diff(theta(t),t))*sin(phi(t)))*((diff(psi(t),t))*cos(theta(t))+diff(phi(t),t))-M[1]/TR[1];
> gl2:=eta*(diff(psi(t),t,t)*sin(theta(t))*cos(phi(t))+(diff(psi(t),t))*(diff(theta(t),t))*cos(theta(t))*cos(phi(t))-(diff(psi(t),t))*sin(theta(t))*(diff(phi(t),t))*sin(phi(t))-(diff(theta(t),t,t))*sin(phi(t))-(diff(theta(t),t))*(diff(phi(t),t))*cos(phi(t)))+(1-zeta)*((diff(psi(t),t))*sin(theta(t))*sin(phi(t))+(diff(theta(t),t))*cos(phi(t)))*((diff(psi(t),t))*cos(theta(t))+diff(phi(t),t))-M[2]/TR[1];
> gl3:=zeta*(diff(psi(t),t,t)*cos(theta(t))-(diff(psi(t),t))*(diff(theta(t),t))*sin(theta(t))+diff(phi(t),t,t))+(eta-1)*((diff(psi(t),t))*sin(theta(t))*sin(phi(t))+(diff(theta(t),t))*cos(phi(t)))*((diff(psi(t),t))*sin(theta(t))*cos(phi(t))-(diff(theta(t),t))*sin(phi(t)))-M[3]/TR[1];
> #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=green):
> #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=green):
> #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=green):
> #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])):
> #Animation des Bewegungsvorganges
> T:=NULL:
> NT:=LinearAlgebra[RowDimension](L[2,1]);
> for j to NT do  #Schleife über jeden Zeitschritt
>   w1:=L[2,1][j,2]; w2:=L[2,1][j,4]; w3:=L[2,1][j,6];
>   P:=plottools[rotate](plottools[rotate](plottools[rotate](POLY,0,0,-w1),-w3,0,0),0,0,-w2);
>   Q:=plots[display](P): T:= T,Q:
> end do:
> AN:=plots[display](T,insequence=true,scaling =constrained,orientation=[45, 45, 0],labels=[x1H,x2H,x3H]):
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab. Das geht wie folgt:
> save Proc_Calc_20, "Proc_Calc_20.m":
# Wir benötigen als Eingabe für die obige Prozedur noch den zeitabhängigen Momentenvektor in Komponenten des Hauptachsensystems. Entsprechend der Aufgabenstellung gilt für den Momentenvektor
> MHS:=Vector([0,0,piecewise(t > ta and t < te,M,0)]);
> ta:=0.40; te:=0.42; M:=5.*10^4;
# Der Bewegungsvorgang zwischen den Zeitpunkten 0  d  t  d TD  wird in  N  gleichgroße Zeitschritte zerlegt. 
> N:=100; TD:=1.;
# Wir geben den Anfangszustand des Systems zum Zeitpunkt  t = 0  vor:
> psi0:=0.: psi0p:=5.: theta0:=evalf(Pi/4): theta0p:=3.: phi0:=0.: phi0p:=5.: 
> AW:=[psi0,psi0p,theta0,theta0p,phi0,phi0p]:
> Proc_Calc_20(EWm,MHS,AW,N,TD,PL):
# Wir geben die Zustandsgrößen aus
> ZG;
# und animieren den Bewegungszustand
> AN; 
> 
;
