# Kapitel 2_f:
# Translation und Rotation, Homogene Koordinaten
# © 2015 Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> interface(rtablesize=infinity): restart: with(LinearAlgebra): with(plots): with(plottools): unprotect(D);
# Beispiel: 2-14:
# Gegeben sind die Eckpunkte P1(2,2,1), P2(6,2,1), P3(6,5,1), P4(2,5,1), P5(2,2,5), P6(6,2,5), P7(6,5,5), P8(2,5,5) eines starren Quaders, der um die Raumdiagonale  mit dem Winkel  Õ  gedreht und anschließend mit dem Vektor v verschoben wird. Die Drehachse verläuft durch den Punkt P1. Gesucht werden die Eckpunkte  des so transformierten Quaders. Es soll eine Maple-Prozedur bereitgestellt werden, mit der die oben beschriebene Transformation automatisiert durchgeführt werden kann. 
# Als Eingabe werden benötigt:
# 1.)  Die Lage und Orientierung der raumfesten Drehachse n  und der Drehwinkel Õ.
# 2.)  Die in Spalten zusammengefassten Koordinaten der Originalpunkte
# 
> Proc_Calc_13:=proc(rd,n,phi,v,OP,ED)
> #----------------------------------------------------------------------------------------
# #  Eingabe:
# #     rd  : Ein Punkt auf der Drehachse
# #     n   : Richtungsvektor der Drehachse (wird programmintern auf die Länge 1 normiert)
# #     phi : Drehwinkel im Bogenmaß um die positive Drehachse
> 
> #     v   : Verschiebungsvektor
> #     OP  : In Spalten zusammengefasste Koordinaten der Originalpunkte
> #     ED  : In Zeilen angeordnete Endknoten der Linienelemente: 
> #           ED[i,1]: Anfangspunkt des Linienelements i
> #           ED[i,2]: Endpunkt     des Linienelements i
> #  Ausgabe:
> #     TP  : Koordinaten der transformierten Originalpunkte
> #-----------------------------------------------------------------------------------------
> local np,OPH,i,nn,ndn,id3,nke,R,M1H,M2H,M3H,M4H,text1,text2,text3,linie1,linie2,linie3,SKNR1,SKNR2,pao,peo,pav,pev,lambda,da,de; 
> global TH,TPH,TP;
> description "Homogene Koordinaten, Translation und Rotation eines Quaders";
> TPH:=Matrix(4,4): OPH:=Matrix(4,4):
> np :=LinearAlgebra[ColumnDimension](OP);
> id3:= Matrix(3,3,shape=identity): 
> #---Aufbau der homogenisierten Matrix der Originalpunkte
> OPH:=ArrayTools[Concatenate](1,OP,Vector[row](np,1)):
> #Berechnung des Versors R
> nn := LinearAlgebra[Normalize](n,Euclidean):
> ndn:= nn.LinearAlgebra[Transpose](nn):          
> nke:= Matrix(3,3,shape=antisymmetric):          
> nke[1,2]:=-nn[3]: nke[1,3]:= nn[2]: nke[2,3]:=-nn[1]:
> R:=ndn + (id3-ndn)*cos(phi) + nke*sin(phi): print(`R = `,R,Determinant(R));
> #---Aufbau der homogenisierten Drehmatrix
> M2H:= Matrix([[R[1,1],R[1,2],R[1,3],0],
>               [R[2,1],R[2,2],R[2,3],0],
>               [R[3,1],R[3,2],R[3,3],0],
>               [0     ,0     ,0     ,1]]):
> #---Aufbau der homogenisierten Verschiebungsmatrizen
> M1H:=Matrix([[1,0,0,-rd[1]],[0,1,0,-rd[2]],[0,0,1,-rd[3]],[0,0,0,1]]);
> M3H:=Matrix([[1,0,0, rd[1]],[0,1,0, rd[2]],[0,0,1, rd[3]],[0,0,0,1]]);
> M4H:=Matrix([[1,0,0,  v[1]],[0,1,0,  v[2]],[0,0,1,  v[3]],[0,0,0,1]]);
> #Berechnung der homogenisierten Transformationsmatrix TH
> TH:= M4H.M3H.M2H.M1H;
> #Berechnung der Koordinaten der transformierten Punkte
> TPH:= TH.OPH;
> TP:=LinearAlgebra[SubMatrix](TPH,[1..3],[1..np]);
> #---Ausgabe der  Transformationsmatrix
> print(`Homogenisierte Transformationsmatrix TH = `,TH);
> print(`Homogenisierte Koordinaten der Originalpunkte = `,OPH);
> print(`Homogenisierte Koordinaten der transformierten Punkte = `,TPH);
> #------Grafische Ausgabe der Ergebnisse
> text1:={}: text2:={}: linie1:={}: linie2:={}:
> for i to np do
> text1:=text1 union {plots[textplot3d]([OP[1,i],OP[2,i],OP[3,i],i],font=[HELVETICA,BOLD,9],align={ABOVE,RIGHT},color=blue)};
> text2:=text2 union {plots[textplot3d]([TP[1,i],TP[2,i],TP[3,i],i],font=[HELVETICA,BOLD,9],align={ABOVE,RIGHT},color=black)};
> end do:
> for i to LinearAlgebra[RowDimension](ED) do
>   SKNR1:=ED[i,1]:SKNR2:=ED[i,2]:
>   pao:=[OP[1,SKNR1],OP[2,SKNR1],OP[3,SKNR1]]:
>   peo:=[OP[1,SKNR2],OP[2,SKNR2],OP[3,SKNR2]]:
>   pav:=[TP[1,SKNR1],TP[2,SKNR1],TP[3,SKNR1]]:
>   pev:=[TP[1,SKNR2],TP[2,SKNR2],TP[3,SKNR2]]:
>   linie1:=linie1 union {plottools[line](pao,peo,thickness=1,color=blue,linestyle=dash)}:
>   linie2:=linie2 union {plottools[line](pav,pev,thickness=1,color=black)}:
> end do:
> lambda:=5:
> da:=seq(rd[i]-lambda*nn[i],i=1..3):
> de:=seq(rd[i]+lambda*nn[i],i=1..3):
> linie3:=plottools[line]([da],[de],thickness=2,color=gold,linestyle=dash):
> text3 :=plots[textplot3d]([da,"Drehachse"],font=[HELVETICA,BOLD,9],align={ABOVE,RIGHT},color=gold);
> plots[display](linie1,linie2,linie3,text1,text2,text3,scaling=constrained,axes=boxed,labels=[x1,x2,x3],title="Ausgangslage (blau, gestrichelt)\nTransformierte Lage (schwarz)",orientation=[60,60,0]);
> end proc:
# Aufbau der Knotendatei für einen Quader
> P1 :=Vector([2,2,1]): P2 :=Vector([6,2,1]): P3 :=Vector([6,5,1]): P4 :=Vector([2,5,1]):
> P5 :=Vector([2,2,5]): P6 :=Vector([6,2,5]): P7 :=Vector([6,5,5]): P8 :=Vector([2,5,5]):
> OP:=Matrix([P1,P2,P3,P4,P5,P6,P7,P8]);
# Zugehörige Elementdatei 
> IZES:=Matrix([[1,2],[2,3],[3,4],[4,1],[1,5],[2,6],[3,7],[4,8],[5,6],[6,7],[7,8],[8,5]]);
# Die Drehachse n (hier die Raumdiagonale) mit dem Punkt rd verläuft durch den Punkt P1, der Drehwinkel beträgt Õ = -À/3 (-60°) und der Verschiebungsvektor ist v. Die Reihenfolge ist 
# 1. Drehung mit Õ
# 2. Verschiebung mit v
> rd:=P1: n:=Vector([1,1,1]): phi:=-Pi/3: v:=Vector([-8,-2,0]):
> Proc_Calc_13(rd,n,phi,v,OP,IZES);
> 
;
