# Kapitel 2_a: 
# Berechnung des Versors bei Vorgabe der Drehachse und des Drehwinkels
# © 2015 Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision=4): restart: with(LinearAlgebra):
# Beispiel 2-2:
# Es wird eine Maple-Prozedur zur Berechnung der Matrix des Versors R bei Vorgabe des Drehachse n und des Drehwinkels Æ bereitgestellt. Die Drehachse hat die Orientierung des Vektors n, der programmintern auf die Länge 1 normiert wird. 
> Proc_Calc_02:=proc(nv::Vector,phi::anything,flag::integer)
> #----------------------------------------------------------------------------
# #  Eingabe:
# #     nv : Vektor der Drehachse(wird programmintern auf die Länge 1 normiert)
# #     phi: Drehwinkel im Bogenmaß um die positive Drehachse
> 
> #    flag: 
> #          flag = 0: Eigenwertberechnung wird nicht durchgeführt
> #          flag = 1: Eigenwertberechnung wird       durchgeführt
> #  Ausgabe:
> #       R:       Versor
> #   Lamba: Eigenwerte    von R (wenn flag = 1)
> #     Phi: Eigenvektoren von R (wenn flag = 1)
> #----------------------------------------------------------------------------
> local  ndn,id,ekn,n,i,j; 
> global R,Lambda,Phi;
> description "Berechnung des Versors R bei Vorgabe der Drehachse  n und des Drehwinkels phi, Lösung des Eigenwertproblems";
> n  := LinearAlgebra[Normalize](nv,Euclidean);
> ndn:= n.LinearAlgebra[Transpose](n):
> id := Matrix(3,3,shape=identity):
> ekn:= Matrix(3,3,shape=antisymmetric);
> ekn[1,2]:=-n[3];
> ekn[1,3]:= n[2];
> ekn[2,3]:=-n[1];
> R:=ndn + (id-ndn)*cos(phi) + ekn*sin(phi);
> if flag = 1 then
> #Berechnung der Eigenwerte und Eigenvektoren des Versors R
> (Lambda,Phi):=LinearAlgebra[Eigenvectors](R);
> for i to 3 do
>   Lambda[i]:=simplify(evalc(Lambda[i]),trig)    assuming phi::real;
>    for j to 3 do
>      Phi[i,j] :=simplify(evalc(Phi[i,j]), trig) assuming phi::real;
>    end do:
> end do;
> print(`Eigenwerte `, Lambda);
> print(`Eigenvektoren `, Phi);
> end if;
> end proc:
# Um die obige Prozedur für folgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab.
> save Proc_Calc_02, "Proc_Calc_02.m";
# Berechnung der Elementardrehungen
> n1 := UnitVector(1,3); Proc_Calc_02(n1,alpha,1): #Drehung um die raumfeste 1-Achse
;
> print(`Versor R =`, R);
> n2 := UnitVector(2,3); Proc_Calc_02(n2, beta,1): #Drehung um die raumfeste 2-Achse
;
> print(`Versor R =`, R);
> n3 := UnitVector(3,3); Proc_Calc_02(n3,gamma,1): #Drehung um die raumfeste 3-Achse
;
> print(`Versor R =`, R);
# Als Drehachse wird im nächsten Schritt die Raumdiagonale gewählt. Der Drehwinkel ist Æ:
> n:=Vector([1,1,1]): Proc_Calc_02(n, phi, 0): 
# Mit der Raumdiagonale als Drehachse und dem Drehwinkel  Æ = À/2 folgt:
> Proc_Calc_02(n, Pi/2, 1):
> print(`Versor R =`, R);
# Mit der Raumdiagonale als Drehachse und dem Drehwinkel  Æ = À (Umklappung) folgt:
> Proc_Calc_02(n, Pi, 1):
# Der linearisierte Versor (cosÆ H 1, sinÆ H Æ) für kleine Drehwinkel Æ ist
> Rlin:=subs(cos(phi)=1,sin(phi)=phi,R);
> 
;
