# Kapitel 3_l: 
# Der unsymmetrische schwere Pyramidenkreisel (Eulerwinkel)
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision=3): restart: with(LinearAlgebra): with(plottools): with(plots):
# Für die anstehenden  Aufgaben benötigen wir die folgenden Prozeduren:
> read "Proc_Calc_19.m"; #Berechnung der massengeometrischen Größen des Kreisels
;
> read "Proc_Calc_20.m"; #Berechnung der Zustandsgrößen, Animation des Bewegungsvorganges
;
> read "PROC_Graf_02.m"; #Erzeugung einer 3D-Plotstruktur des tetraederförmigen Kreisels
;
# Beispiel 3-23:  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 Kreisel eine Knoten- und eine Elementdatei.  Dazu beschreiben wir die Pyramide im raumfesten Koordinatensysterm durch ihre vier  Eckpunkte und fassen deren Koordinaten 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
# Wir schauen uns den Kreisel an:
> TITEL:="Der unsymmetrische schwere Pyramidenkreisel im körperfesten Bezugssystem": ABEZ:=["x1","x2","x3"]:
> Proc_Graf_02(COORD,IZES,TITEL,ABEZ);
# Wir beschaffen uns die Massenmomente des unsymmetrischen Kreisels.
> Proc_Calc_19(COORD,IZES,rho):
# Zur Lösung des nichtlinearen  Differenzialgleichungssystems benötigen wir die Hauptträgheitsmomente und die Hauptträgheitsachsen bezogen auf den Ursprung des körperfesten Punktes P1 die wir der folgenden Eigenwertberechnung entnehmen. Dabei ist zu beachten, dass die Eigenvektoren nur bis auf einen konstanten Faktor bestimmt sind, denn den Eigenvektoren sind  nur die Achsrichtungen und nicht die Achsorientierungen festgelegt, da neben ei auch -ei ein Eigenvektor ist. 
> (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: 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. 
# Der Aufbau eines rechtshändigen Basissystems kann übrigens auch manuell erfolgen. Dazu orientieren wir uns optisch an der Grafik im Buch TM 3 und wählen folgende Eigenvektoren:
> EVm:=Matrix([Column(EVq,2),Column(EVq,3),Column(EVq,1)]);
# Die zugehörigen Eigenwerte 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  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 schwere Pyramidenkreisel im raumfesten Bezugssystem mit Hauptachsendarstellung":
> Proc_Graf_01(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 3). Im gedrehten Hauptachsensystem sind dann die Eckpunkte des Tetraeders zeilenweise in der Matrix CH angeordnet. 
> CH:=COORD.EVm;
# Wir stellen das Tetraeder im gedrehten Hauptachsensystem grafisch dar.
> TITEL:="Der unsymmetrische schwere Pyramidenkreisel im Hauptachsensystem": ABEZ:=["x1H","x2H","x3H"]:
> Proc_Graf_02(CH,IZES,TITEL,ABEZ);
# Im nächsten Schritt wird der Bewegungsvorgang des unsymmetrischen schweren Kreisels in Eulerwinkeln berechnet werden. Die einzige äußere Kraft ist die im Schwerpunkt angreifende Gewichtskraft G = m g. Wir nehmen an, dass die Gewichtskraft mit G = -G e3  parallel zur raumfesten x3-Achse wirkt. 
> G:=MA*981;
# Im raumfesten 0-Basissystem ist dann
> G0:=Vector([0,0,-G]);
# Wir benötigen die Eulersche Drehmatrix
> RE:=Matrix(3,3,{(1,1)=cos(psi(t))*cos(phi(t))-sin(psi(t))*cos(theta(t))*sin(phi(t)),(1,2)=-cos(psi(t))*sin(phi(t))-sin(psi(t))*cos(theta(t))*cos(phi(t)),(1,3)=sin(psi(t))*sin(theta(t)),(2,1)=sin(psi(t))*cos(phi(t))+cos(psi(t))*cos(theta(t))*sin(phi(t)),(2,2)=-sin(psi(t))*sin(phi(t))+cos(psi(t))*cos(theta(t))*cos(phi(t)),(2,3)=-cos(psi(t))*sin(theta(t)),(3,1)=sin(theta(t))*sin(phi(t)),(3,2)=sin(theta(t))*cos(phi(t)),(3,3)=cos(theta(t))}):
> GH:=Transpose(RE).G0;
# Wir benötigen weiterhin den Ortsvektor des Volumenmittelpunktes (Schwerpunktes) im Hauptachsensystem
> RSH:=Transpose(EVm).RS;
# Das Moment aus Eigengewicht ist dann
> MH0:=CrossProduct(RSH,GH);
# Wir geben den Anfangszustand des Systems zum Zeitpunkt  t = 0 vor. Die Simulationszeit TD beträgt 1 Sekunde. Die Anzahl der Simulationsschritte ist N = 100.
> psi0:=0.: psi0p:=5.: theta0:=evalf(Pi/4): theta0p:=3.: phi0:=0.: phi0p:=5.: 
> AB:=[psi0,psi0p,theta0,theta0p,phi0,phi0p]:
> N:=100:TD:=1.:
> Proc_Calc_20(EWm,MH0,AB,N,TD,PL):
> ZG;
> AN;
> 
;
# 
