# Kapitel 2_d: Kardanwinkel (±,²,³)
# © 2015 Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision=5): restart: with(LinearAlgebra): _local(gamma):
# Beispiel: 2-10:
# Wir beschaffen uns eine Prozedur, mittels derer wir bei Vorgabe der Kardanwinkel (±, ², ³) die Drehmatrix TK(±,²,³) = T(±)·T(²)·T(³)·automatisiert  berechnen können.
> Proc_Calc_09:=proc(KW::list,flag::integer)
> #-------------------------------------------------------
> # Eingabe:
> #        KW:     Kardanwinkel [alpha,beta,gamma]
> #      flag:     flag = 0 : Winkel im Gradmaß
> #                flag = 1 : Winkel in Radiant
> # Ausgabe:
> #       TK:      Drehmatrix TK(alpha,beta,gamma)
> #-------------------------------------------------------
> local w,fak,ca,cb,cg,sa,sb,sg;
> global Ta,Tb,Tg,TK,TKlin;
> description "Berechnung der Drehmatrix TK mittels der Kardanwinkel";
> if   flag = 0 then fak:=Pi/180;
> elif flag = 1 then fak:=1;
> else
>   print(`Wählen Sie\nflag = 0 (Winkel im Gradmaß)\nflag = 1 (Winkel im Bogenmaß)`); return;
> end if;
> w    := [seq(fak*KW[i],i=1..3)];
> ca:=cos(w[1]): cb:=cos(w[2]): cg:=cos(w[3]): 
> sa:=sin(w[1]): sb:=sin(w[2]): sg:=sin(w[3]):
> Ta:= Matrix(3,3,[[1 ,  0,  0],[0,ca,-sa],[  0,sa,ca]]):
> Tb:= Matrix(3,3,[[cb,  0, sb],[0, 1,  0],[-sb, 0,cb]]):
> Tg:= Matrix(3,3,[[cg,-sg,  0],[sg,cg, 0],[  0, 0, 1]]):
> TK:= simplify(Ta.Tb.Tg):
> TKlin:=Matrix(3,3,[[1,-w[3],  w[2]],[w[3],1,-w[1]],[-w[2],w[1], 1]]):
> print(`Drehmatrix TK = `,TK);
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab.
> save Proc_Calc_09, "Proc_Calc_09.m":
# Wir rufen die Prozedur auf:
> Proc_Calc_09([alpha,beta,gamma],1);
> print(Ta,Tb,Tg,TKlin);
# Beispiel: 2-11:
# Es wird eine Maple-Prozedur zur Verfügung gestellt, die bei Vorgabe einer Drehmatrix  TK(±,²,³)  die Kardanwinkel (±, ², ³) berechnet. 
> Proc_Calc_10:=proc(R::Matrix)
> #-------------------------------------------------------------------
# #  Eingabe:
# #      R   : Drehmatrix 
> 
> #  Ausgabe:
> #      Kardanwinkel (alpha, beta, gamma)
> #-------------------------------------------------------------------
> local cd,Z,zn,eps,R1,alpha,beta,gamma,ca,cb,cg,sa,sb,sg,i,j,a,b,g;
> global Ra,Rb,Rg,RK,gl,loe;
> description "Berechnung der Kardanwinkel aus der Drehmatrix";
> #Es wird zunächst geprüft, ob R orthogonal 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(`Drehmatrix ist nicht orthogonal`); return;
> end if;
> #Es wird geprüft, ob R eine Umklappung mit phi=(+,-)Pi/2 erzeugt 
> R1:=evalf(LinearAlgebra[Trace](R));
> if abs(R1 + 1) < eps then
>   print(`Es liegt mit der Drehmatrix eine Umklappung vor`); return;
> end if;
> #Aufbau der Drehmatrix als Funktion von (alpa,beta,gamma) 
> ca:=cos(alpha): cb:=cos(beta): cg:=cos(gamma): cg:=cos(gamma): 
> sa:=sin(alpha): sb:=sin(beta): sg:=sin(gamma): sg:=sin(gamma):
> Ra:= Matrix(3,3,[[1 ,  0,  0], [0,ca,-sa],[  0,sa,ca]]):
> Rb:= Matrix(3,3,[[cb,  0, sb], [0, 1,  0],[-sb, 0,cb]]):
> Rg:= Matrix(3,3,[[cg,-sg,  0], [sg,cg, 0],[  0, 0, 1]]):
> RK := Ra.Rb.Rg:
> #Aufstellung des nichtlinearen Gleichungssystems und Berechnung der Kardanwinkel
> gl:=NULL:
> for i to 3 do 
>   for j to 3 do
>     gl:=gl,R[i,j]-RK[i,j]:
>   end do:
> end do:
> loe:=solve([gl]);
> #Ausgabe der Kardanwinkel
> print(`Winkel in Radiant = `,allvalues(loe));
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab.
> save Proc_Calc_10, "Proc_Calc_10.m":
# Entsprechend der Aufgabenstellung ist
> w:=sqrt(3):
> TK:=1/8*Matrix([[6+w,3-2*w,2],[3-2*w,2+3*w,2*w],[-2,-2*w,4*w]]):
> print(TK,evalf(TK));
# Wir rufen die Prozedur auf:
> Proc_Calc_10(TK);
# Maple liefert uns in diesem Fall zwei Lösungen:
> loe1:=convert(allvalues(loe)[1],list);loe2:=convert(allvalues(loe)[2],list);
> winkel1rad:=seq(evalf(rhs(loe1[i])),i=1..3);
> winkel2rad:=seq(evalf(rhs(loe2[i])),i=1..3);
# Wir prüfen beide Ergebnisse:
> Proc_Calc_09([winkel1rad],1);
> Proc_Calc_09([winkel2rad],1);
# Die in Radiant gegebene Lösung wird in Gradmaß umgewandelt:
> winkel1grad:=seq(evalf(convert(rhs(loe1[i]),degrees)),i=1..3);
> winkel2grad:=seq(evalf(convert(rhs(loe2[i]),degrees)),i=1..3);
# Keine eindeutige Lösung liefert die Drehmatrix
> TK:=Matrix([[0,0,1],[1,0,0],[0,1,0]]);
> Proc_Calc_10(TK);
# Der Winkel ± kann hier beliebig gewählt werden. Wir prüfen das obige Ergebnis:
> Proc_Calc_09([alpha,Pi/2,Pi/2-alpha],1);
> 
;
