# Kapitel_02_b:
# Die Bewegung des starren Körpers
# © 2015  Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> interface(displayprecision=4): restart: with(LinearAlgebra): with(Physics): with(plots): with(plottools):
# Wir stellen zunächst eine Maple-Prozedur zur Verfügung, die das Kreuzprodukt eines Vektors n mit einem Tensor 2. Stufe A berechnet. Dabei ist aus der Tensorrechnung bekannt:
# `&x`(n, A);= -`&x`(A^T, n)^T*und;`&x`(A, n);= -`&x`(n, A^T)^T;. Insbesondere gilt für A = 1:  n*1 = 1*n;.
# Zur Bildung des Kreuzproduktes  `&x`(a, b) = `&x`(a[j]*e[j], b[k])*e[k] and `&x`(a[j]*e[j], b[k])*e[k] = a[j]*b[k]*varepsilon[jkl]*e[l];   der  beiden Vektoren a und b verwenden wir das dreifach indizierte Permutationssymbol  varepsilon[jkl] . ;
> Proc_Calc_03:=proc(n::Vector,A::Matrix)
> #-----------------------------------------------------------------
# #  Eingabe:
# #      n   : Vektor
# #      A   : Matrix
> 
> #  Ausgabe:
> #      nkA : Kreuzprodukt n x A (vektorielle  Linksmultiplikation)
> #      Akn : Kreuzprodukt A x n (vektorielle Rechtsmultiplikation)
> #-----------------------------------------------------------------
> local j,k,l,p;
> global nkA, Akn;
> description "Der Epsilon-Tensor, Berechnung der vektoriellen Links- und Rechtsmultiplikation eines Vektors mit einem Tensor 2.Stufe";
> nkA:= Matrix(3):Akn:= Matrix(3):
> for j to 3 do
>   for k to 3 do
>     for l to 3 do
>       for p to 3 do
>         nkA[p,l]:= nkA[p,l] + n[j]*A[k,l]*Physics[LeviCivita][j,k,p]:
>         Akn[k,p]:= Akn[k,p] + n[j]*A[k,l]*Physics[LeviCivita][l,j,p]:
>       end do:
>     end do:
>   end do:
> end do:
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab.
> save Proc_Calc_03, "Proc_Calc_03.m":
# Als Beispiel wird das Kreuzprodukt des Vektors  n mit der Einheitstensor 1 berechnet. 
> N:=Vector(3,symbol=n): A:=Matrix(3,shape=identity):
> Proc_Calc_03(N,A):
# Damit sind , und das gilt nur, wenn A der Einheitstensor ist:
> print(nkA,Akn);
# Ist A allgemein ein symmetrischer Tensor, dann ist auch B = n×A×n ein symmetrischer Tensor.
> A:=Matrix(3,shape=symmetric,symbol=a);
> Proc_Calc_03(N,Akn): B:=nkA;
# Beispiel 2-3: 
# Das in der Aufgabenstellung skizzierte starre Rechteck   soll durch Drehung um die Achse n mit dem Winkel Æ = 2/3À (120°) in eine benachbarte Lage gebracht werden. Es sind die Komponenten des Versors R und die Koordinaten der Eckpunkte des Rechtecks in der gedrehten Lage  ´ gesucht. Entsprechend der Aufgabenstellung sind
> n:=1/sqrt(3)*Vector([-1,-1,1]); phi:=2/3*Pi;
# Damit erhalten wir den Versor R
> read "Proc_Calc_02.m": Proc_Calc_02(n,phi,0);
# Die Ortsvektoren der  Eckpunkte A,B,C des starren Rechtecks werden spaltenweise in der Matrix rP gespeichert. 
> r0:=Vector([0,0,0]): rA:=Vector([2,0,0]): rB:=Vector([2,1,0]):rC:=Vector([0,1,0]): 
# Wir fassen die vier Ortsvektoren des Rechtecks in einer vierspaltigen Matrix zusammen:
> rP0:=Matrix([r0,rA,rB,rC]);
# Nach der Drehung sind
> rPs:=R.rP0;
# Wir animieren die Drehung des Rechtecks. Dazu verwenden wir den Befehl plottools[polygon].
> T:=NULL: NT:=25: phir:=0.: dphir:=phi/NT:
> for j to NT + 1 do   #Schleife über jeden Zeitschritt
>   Proc_Calc_02(n,phir,0);
>   P:=plottools[polygon](LinearAlgebra[Transpose](R.rP0), color = blue, thickness = 2):
>   Q:=plots[display](P): T:= T,Q:
>   phir:=phir + dphir;
> end do:
> AN1:=plots[display](T,insequence=true,axes = normal,scaling=constrained,orientation=[45, 45, 0],labels=[x1,x2,x3],transparency=0.3):
> AN1;
# Beispiel 2-4:
# Es wird eine Maple-Prozedur zur Berechnung des Einheitsvektors n der Drehachse, des  Drehwinkels Õ und des Drehvektors d = Õn  bei Vorgabe eines Versors R bereitgestellt.
> Proc_Calc_04:=proc(R::Matrix)
> #-----------------------------------------------------------------------
# #  Eingabe:
# #      R   : Versor
> 
> #  Ausgabe:
> #      Rx  : Vektor des Versors
> #      n   : Einheitsvektor in Richtung der Drehachse
> #      phir: Drehwinkel im Bogenmaß
> #      phid: Drehwinkel im Gradmaß
> #      d   : Drehvektor
> #-----------------------------------------------------------------------
> local cd,Z,zn,eps,R1,R2,R3,w,wb,c;
> global Rx,phir,phid,n,d;
> description "Berechnung des Drehvektors und Drehwinkels eines Versors R";
> #Es wird zunächst geprüft, ob R ein orthogonaler Tensor ist.
> cd := LinearAlgebra[ColumnDimension](R):
> Z  := R.LinearAlgebra[Transpose](R) - Matrix(cd,shape=identity);
> zn := LinearAlgebra[MatrixNorm](Z,Frobenius);
> eps:= 10.^(-Digits+2);
> if zn > eps then
>   print(`R ist kein orthogonaler Tensor`); return;
> end if;
> #Es wird geprüft, ob R ein uneigentlicher Versor ist. In diesem Fall 
> #liegt eine Klappung mit phi=(+,-)Pi/2 vor und die Berechnung wird abgebrochen. 
> R1:=evalf(LinearAlgebra[Trace](R));                 #1.Invariante
> if abs(R1 + 1) < eps then
>   print(`R ist ein uneigentlicher Versor`); return;
> end if;
> #Berechnung der Hauptinvarianten R2 und R3
> R2  :=0.5*(R1^2-LinearAlgebra[Trace](R.R)):         #2.Invariante
> R3  :=LinearAlgebra[Determinant](R):                #3.Invariante
> #Drehvektor und Drehwinkel
> Rx  := evalf(Vector([R[2,3]-R[3,2],R[3,1]-R[1,3],R[1,2]-R[2,1]]));
> w   :=-Rx/(1 + R1);
> wb  := evalf(LinearAlgebra[Norm](w,Euclidean));
> if wb < eps then
>   print(`Es liegt mit phi = 0 die identische Abbildung vor`);return;
> end if; 
> phir:= arctan(2*wb,(1-wb^2));
> phid:= evalf(convert(phir,degrees));
> n   := w/wb;
> d   := phir*n;
> # Ausgabe der Ergebnisse
> print(`Vektor des Versors =`,Rx):
> print(`Vektor n = `, n);
> print(`Drehwinkel im Bogenmaß = `,phir);
> print(`Drehwinkel im Gradmaß  = `,phid);
> print(`Drehvektor d = n . phi = `,d);
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab.
> save Proc_Calc_04, "Proc_Calc_04.m":
# Mit den Werten des Beispiels sind
> R1:=1/9*Matrix([[1,-4,8],[8,4,1],[-4,7,4]]); Proc_Calc_04(R1);
> (Lambda,Phi):=LinearAlgebra[Eigenvectors](R1);
> R2:=Matrix([[0,0,1],[1,0,0],[0,1,0]]); Proc_Calc_04(R2);
> (Lambda,Phi):=LinearAlgebra[Eigenvectors](R2);
# Beispiel 2-5:
# 1. Umklappung
> Proc_Calc_05:=proc(n::Vector)
> #--------------------------------------------------------------------
# #  Eingabe:
# #      n   : Vektor der Drehachse
> 
> #            (wird intern auf die Länge 1 normiert)
> #  Ausgabe:
> #      PSI  : Umklappungsmatrix R(n,Pi)
> #--------------------------------------------------------------------
> local nn,gl;
> global PSI;
> description "Berechnung der Umklappungsmatrix Psi";
> nn:=LinearAlgebra[Normalize](n,Euclidean);
> gl:={abs(n[1])^2+abs(n[2])^2+abs(n[3])^2=1};
> PSI:=simplify(2*nn.LinearAlgebra[Transpose](nn)-LinearAlgebra[IdentityMatrix](3),gl);
> print(`Umklappungsmatrix PSI =`,PSI):
> end proc:
# Es soll ein beliebiger Punkt r des Raumes um die Achse mit dem Richtungsvektor n mit dem Winkel Æ = À gedreht (umgeklappt) werden.
> n:='n': R:='R':
> N:=Vector(3,symbol=n);
> Proc_Calc_05(N);
# Für spätere Anwendungen wird ¨ gespeichert.
> PSIs:=PSI:
# Wir berechnen die Determinante der Umklappungsmatrix ¨
> gl:={2*(n[1]^2+n[2]^2+n[3]^2)=2};
> simplify(Determinant(PSI),gl);
# und die Eigenwerte und Eigenvektoren sind
> (Lambda,Phi):=LinearAlgebra[Eigenvectors](PSI);
# Wir zeigen, dass die symmetrische Umklappungsmatrix identisch ist mit ihrer Inversen:
> MatrixInverse(PSI)-PSI;
# Wir vereinfachen den obigen Ausdruck:
> simplify(MatrixInverse(PSI)-PSI,gl);
# Wir rechnen zur Umklappung ein Zahlenbeispiel. Der Vektor
> N:=Vector([1,1,1]);
# legt die Richtung der Drehachse (Raumdiagonale) fest, um die mit dem Winkel Æ = À gedreht (umgeklappt) wird . Damit ergibt sich die Umklappungsmatrix zu 
> Proc_Calc_05(N);
# Der planare Vektor
> r:=Vector(3,[1,1,0]);
# geht damit über in 
> rstr:=PSI.r;
# noch Beispiel 2-5: Wir wenden uns nun der Spiegelung eines Punktes P mit dem Ortsvektor r an einer Spiegelungsebene    zu. 
# 2. Spiegelung
> Proc_Calc_06:=proc(n::Vector)
> #--------------------------------------------------------------------
# #  Eingabe:
# #      n   : Vektor der Drehachse
> 
> #            (wird intern auf die Länge 1 normiert)
> #  Ausgabe:
> #      Q  :  Spiegelungsmatrix
> #--------------------------------------------------------------------
> local nn,gl;
> global Q;
> description "Berechnung der Spiegelungsmatrix Q";
> nn:=LinearAlgebra[Normalize](n,Euclidean);
> gl:={abs(n[1])^2+abs(n[2])^2+abs(n[3])^2=1};
> Q :=simplify(LinearAlgebra[IdentityMatrix](3)-2*nn.LinearAlgebra[Transpose](nn),gl);
> print(`Spiegelungsmatrix Q =`,Q):
> end proc:
# Wir berechnen die Spiegelungsmatrix
> N:=Vector(3,symbol=n);
> Proc_Calc_06(N);
# Für die Spiegelungsmatrix gilt: QÅQ = Q2 = 1 oder Q-1  = Q
> simplify(MatrixInverse(Q),gl);
# Es soll ein beliebiger Punkt r des Raumes an der Ebene   mit dem Normalenvektor n gespiegelt werden. 
> N:='N':
> r:=Vector(3,symbol=x); 
> rstr:=Q.r;
# Die Spiegelung kann immer als Hintereinderschaltung von Umklappung ¨ und einer Totalinversion mit (-1) dargestellt werden.
> PSI:='PSI':
> Q:=-IdentityMatrix(3).PSIs;
#  Maple ermöglicht uns mit der Prozedur plottools[reflect] die Spiegelung von 2-D oder 3-D Objekten.
> with(plots): with(plottools):
> poly := [[4, 1],[5,1],[6, 2], [6, 3], [4, 3] ]:
> p:=polygonplot(poly,scaling=constrained,color = "DarkGreen",transparency = .5):
> Punkte:=[0, 0], [6, 5]:
> gerade:=line(Punkte, color = red, linestyle = dash):
> q := reflect(p,[Punkte]):
> plots[display]([p,q,gerade],axes=boxed,gridlines=true);
# Beispiel 2-6:
# Es soll die Drehung eines Rechtecks mit dem Winkel Æ = À um die 3-Achse durch zwei Umklappungen realisiert werden. Wir wählen zunächst die durch den Nullpunkt verlaufenden Einheitsvektoren a und b, die voraussetzungsgemäß den Winkel Æ/2 = À/2 einschließen. Den Vektor a wählen wir so, dass er gegenüber der positiven x1-Achse um den (beliebigen) Winkel ±  geneigt ist. Die Lage des Vektors b ergibt sich dann aus der Forderung, dass die Drehrichtung von a nach b identisch sein muss mit der Drehrichtung von Æ um die 3-Achse.
> ID:=IdentityMatrix(3);
> a:=Vector([cos(alpha),sin(alpha),0]);
> Psia:=simplify(2*a.Transpose(a)-ID);
> b:=Vector([-sin(alpha),cos(alpha),0]);
> Psib:=2*b.Transpose(b)-ID;
# Unabhängig vom (beliebig) zu wählenden Winkel ± errechnen wir
> PSI:=simplify(Psib.Psia);
# 1. Schritt: Klappung des Rechtecks aus der Ausgangslage um die Achse a. Wir wählen ± = À/6. 
> alpha:=Pi/6:
# Die Achse a notieren wir in folgender Form
> aa:=[0, 0],2*[cos(alpha),sin(alpha)]:
> la:=line(aa,color = red, linestyle = dash):
# Die Achse b ist gegenüber der Achse a um À/2 gedreht:
> ab:=[0, 0],2*[cos(alpha+Pi/2),sin(alpha+Pi/2)]:
> lb:=line(ab,color = green, linestyle = dash):
# Wir fassen die Koordinaten der Eckpunkte des zu klappenden Rechtecks 
> r1:=Vector([0,0,0]): r2:=Vector([0,-1,0]): r3:=Vector([2,-1,0]): r4:=Vector([2,0,0]):
# in folgender Matrix zusammen:
> KM0:=Matrix([r1,r2,r3,r4]);
> knl0:=seq([seq(KM0[i,j],i=1..2)],j=1..ColumnDimension(KM0));
> p0:=polygon([knl0],color = blue, thickness=2,transparency=0.6):
> display(p0,title = "\nRechteck in der Ausgangslage", titlefont = ["ARIAL", 15], labels = ["x1", "x2"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],  axesfont = ["HELVETICA", "ROMAN", 8],scaling=constrained,gridlines=true,axes=boxed);
# Es erfolgt die Klappung des Rechtecks um die Achse a. 
> KMa:=Psia.KM0;
> knla:=seq([seq(KMa[i,j],i=1..2)],j=1..ColumnDimension(KMa));
> pa:=polygon([knla],color = red, thickness=2,transparency=0.6):
> display(p0,pa,la,title = "\nKlappung um die Achse a\n", titlefont = ["ARIAL", 15], labels = ["x1", "x2"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],  axesfont = ["HELVETICA", "ROMAN", 8],scaling=constrained,gridlines=true,axes=boxed);
# 2. Schritt: Umklappung um die Achse b
> KMb:=simplify(Psib.KMa);
> knlb:=seq([seq(KMb[i,j],i=1..2)],j=1..ColumnDimension(KMb));
> pb:=polygon([knlb],color=green,thickness=2,transparency=0.6):
> display(pa,pb,lb,title = "\nKlappung um die Achse b\n", titlefont = ["ARIAL", 15], labels = ["x1", "x2"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],  axesfont = ["HELVETICA", "ROMAN", 8],scaling=constrained,gridlines=true,axes=boxed);
# Beispiel 2-7:
# Wir fassen die Koordinaten der Eckpunkte 
> r1:=Vector([0,0,0]): r2:=Vector([2,0,0]): r3:=Vector([2,1,0]): r4:=Vector([0,1,0]):
# des Rechtecks in der Ausgangslage in der Matrix zusammen:
> KM0:=Matrix([r1,r2,r3,r4]);
# Die Einheitsvektoren der Drehachsen sind:
> Na:=Vector([1,0,0]);Nb:=Vector([0,1,0]);
# 1. Drehreihenfolge: 1. Drehung um die 2-Achse, anschließend Drehung um die 1-Achse, jeweils mit dem Drehwinkel Æ = À/2
> T:=NULL: Q:=NULL: NT:=25: alphar:=0: dalphar:=Pi/2/NT: betar:=0.: dbetar:=Pi/2/NT:
> for j from 1 to NT  do   #Schleife über jeden Zeitschritt
>   Proc_Calc_02(Nb,betar,0);
>   P:=plottools[polygon](LinearAlgebra[Transpose](R.KM0), color = blue, thickness = 2):
>   Q:=plots[display](P): T:= T,Q:
>   betar:=betar + dbetar;
> end do:
> Proc_Calc_02(Nb,betar,0); KMa:=R.KM0:
> for j from 1 to NT+1 do   #Schleife über jeden Zeitschritt
>   Proc_Calc_02(Na,alphar,0);
>   P:=plottools[polygon](LinearAlgebra[Transpose](R.KMa), color = red, thickness = 2):
>   Q:=plots[display](P): T:= T,Q:
>   alphar:=alphar + dalphar;
> end do:
> AN2:=plots[display](T,insequence=true,title = "\nDrehreihenfolge 2-1-Achse\n",axes = normal,scaling=constrained,orientation=[45, 45, 0],labels=[x1,x2,x3],transparency=0.3):
> AN2;
# 2. Drehreihenfolge: 1. Drehung um die 1-Achse, anschließend Drehung um die 2-Achse, jeweils mit dem Drehwinkel Æ = À/2
> T:=NULL: Q:=NULL: alphar:=0.: betar:=0.:
> for j from 1 to NT  do     #Schleife über jeden Zeitschritt
>   Proc_Calc_02(Na,alphar,0);
>   P:=plottools[polygon](LinearAlgebra[Transpose](R.KM0), color = blue, thickness = 2):
>   Q:=plots[display](P): T:= T,Q:
>   alphar:=alphar + dalphar; #print('dalphar = ',dalphar);
> end do:
> Proc_Calc_02(Na,alphar,0); KMb:=R.KM0:
> for j from 1 to NT + 1 do   #Schleife über jeden Zeitschritt
>   Proc_Calc_02(Nb,betar,0);
>   P:=plottools[polygon](LinearAlgebra[Transpose](R.KMb), color = red, thickness = 2):
>   Q:=plots[display](P): T:= T,Q:
>   betar:=betar + dbetar;
> end do:
> AN3:=plots[display](T,insequence=true,title="\nDrehreihenfolge 1-2-Achse\n",axes = normal,scaling=constrained,orientation=[45, 45, 0],labels=[x1,x2,x3],transparency=0.3):
> AN3;
> 
;
