# Kapitel  3_f: Massenträgheitsmomente regelmäßiger Körper
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart:with(LinearAlgebra): with(plottools): with(plots):
# Es wird eine Maple-Prozedur zur Auswertung von Dreifachintegralen (s.h. TM 3, Kap. 3.1) bereitgestellt. 
> Proc_Calc_18:=proc(K1::list,K2::list,GR::listlist,rho::anything)
> #-------------------------------------------------------------------------------------
> #Eingabe: 
> #          K1: Kartesische Koordinaten [x,y,z]
> #          K2: Allgemeine  Koordinaten [u,v,w]
> #          GR: Grenzen der allgemeinen Koordinaten
> #         rho: Dichte des Materials
> #Ausgabe:  
> #          JM: Jacobimatrix
> #          DE: Determinante der Jacobimatrix
> #           V: Volumen des Rotationskörpers
> #           m: Masse des Rotationskörpers
> #          rS: Volumenmittelpunkt (Schwerpunkt)
> #      Thetaq: Trägheitsmatrix bezogen auf den Ursprung des Ausgangskoordinatensystems
> #       Theta: Trägheitsmatrix bezogen auf Zentralachsen (ZA)
> #          EW: Eigenwerte der Matrix Theta    (Hauptträgheitsmomente)
> #          EV: Eigenvektoren der Matrix Theta ( Hauptträgheitsachsen)
> #-------------------------------------------------------------------------------------
> local R1,ID,INQ,IN,VB;
> global JM,DE,V,m,rS,Thetaq,Theta,EW,EV;
> description "Berechnung von Massenträgheitsmomenten regelmäßiger homogener Körper";
> R1:=Vector(K1); rS:=Vector(3); ID:=IdentityMatrix(3);
> # Jacobimatrix und deren Determinante
> (JM,DE):=VectorCalculus[Jacobian](K1,K2,'determinant'); 
> # Volumen
> V:=factor(int(int(int(DE,K2[1]=GR[1,1]..GR[1,2]),K2[2]=GR[2,1]..GR[2,2]),K2[3]=GR[3,1]..GR[3,2]));
> #Masse
> m:=rho*V;
> #Volumenmittelpunkt
> rS :=simplify(map(int,map(int,map(int,DE*R1,K2[1]=GR[1,1]..GR[1,2]),K2[2]=GR[2,1]..GR[2,2]),K2[3]=GR[3,1]..GR[3,2])/V);
> #Integrand
> INQ :=DE*(LinearAlgebra[DotProduct](R1,R1,conjugate=false)*ID-R1.LinearAlgebra[Transpose](R1));
> #Trägheitsmomente bezogen auf das Ausgangskoordinatensystem
> Thetaq:=rho*map(int,map(int,map(int,INQ,K2[1]=GR[1,1]..GR[1,2]),K2[2]=GR[2,1]..GR[2,2]),K2[3]=GR[3,1]..GR[3,2]);
> #Hauptachsentransformation
> # 1.Parallelverschiebung des Koordinatensystems in den Volumenmittelpunkt (Schwerpunkt)
> Theta:=Thetaq-m*(LinearAlgebra[DotProduct](rS,rS,conjugate=false)*ID - rS.Transpose(rS));
> # 2. Anschließende Drehung des Koordinatensystems in das Zentralachsensystem
> (EW,EV):= LinearAlgebra[Eigenvectors](Theta);
> #Ausgabe der Ergebnisse
> print(`Jacobimatrix =       `,simplify(JM,trig));
> print(`Determinante der Jacobimatrix = `,simplify(DE));
> print(`Volumen = `,V);
> print(`Masse m = `,m);
> print(`Schwerpunktskoordinaten (rS) = `,rS);
> print(`Trägheitsmomente im Ausgangssystem (Thetaq) =  `, map(factor,Thetaq));
> print(`Trägheitsmomente im Zentralachsensystem (Theta)  = `, map(factor, Theta));
> print(`Hauptträgheitsmomente = `,map(factor,EW),`Eigenvektoren = `, map(factor,EV));
> end proc:
# Um die Prozedur für nachfolgende Arbeitsblätter anwendbar zu machen, speichern wir diese ab. Das geht wie folgt:
> save Proc_Calc_18, "Proc_Calc_18.m":
# Beispiel 3-12:
# Homogener Quader der Breite 0 d x1  d  a,  Höhe  0 d  x2 d b,  Tiefe 0 d  x3  d c
> display(cuboid([0,0,0],[4,2,3]),axes=normal,orientation=[-145,75,0],scaling=constrained,lightmodel=light4,shading=xyz,labels=["x1","x2","x3"]);
> Proc_Calc_18([x,y,z],[x,y,z],[[0,a],[0,b],[0,c]],rho);
# Beispiel 3-13:
# Torus  (R: Abstand vom Toruszentrum zum Mittelpunkt des Meridians, r: Radius des Meridians)
> display(torus([0,0,0], 1, 2), scaling=constrained, lightmodel=light2, shading=xyz, labels=["x1","x2","x3"]);
> x:=(R+r*sin(theta))*cos(phi); y:=(R+r*sin(theta))*sin(phi); z:=r*cos(theta);
> Proc_Calc_18([x,y,z],[r,theta,phi],[[0,a],[0,2*Pi],[0,2*Pi]],rho);
# Beispiel 3-14:  (Rotationskörper )
# Wir stellen zunächst eine Prozedur zur Animation des Rotationsvorganges zur Verfügung:
> PROC_GRAF_01:=proc(f,a,b)
> #---------------------------------------------------------
# #  Eingabe:
# #        f:  Erzeugende Funktion
# #      a,b:  Linke und rechte Grenze in z-Richtung
> #  Ausgabe:  Animation des Rotationsvorganges
> #---------------------------------------------------------
> local c:
> description "3D-Animation der Rotation einer Funktion f(x3) um die x3-Achse";
> c:=[f(z)*cos(t*phi),f(z)*sin(t*phi),z];
> plots[animate3d](c,z=a..b,phi=0..2*Pi,t=0..1,style=patchcontour,shading=zhue,frames=20,labels=["x1","x2","x3"],axes=boxed,orientation=[155,25,125]);
> end proc:
# Rotationskörper (Erzeugende Funktion) 
> r:='r': phi:='phi': z:='z':
> a:=unapply(-(1/8)*z^4+(19/12)*z^3-7*z^2+(295/24)*z-23/4,z);
> za:=2; ze:=5;  #Grenzen in x3-Richtung
;
# Wir animieren den  Rotationsvorgang:
> PROC_GRAF_01(a(z),za,ze);
> x:=r*cos(phi); y:=r*sin(phi);
> Proc_Calc_18([x,y,z],[r,phi,z],[[0,a(z)],[0,2*Pi],[za,ze]],rho);
# Gerader Kreiszylinder mit dem Radius a und der Höhe h
> display(cylinder([0,0,0],1,3,capped=true,strips=40),orientation=[45, 70],scaling=constrained,lightmodel=light2,shading=xyz, labels=["x1","x2","x3"]);
> r:='r': phi:='phi': z:='z': a:='a':
> x:=r*cos(phi); y:=r*sin(phi);
> Proc_Calc_18([x,y,z],[r,phi,z],[[0,a],[0,2*Pi],[0,h]],rho);
# Für einen Hohlzylinder mit dem Innenradius ri und dem Außenradius ra und der Höhe h sind die Trägheitsmomente
> Proc_Calc_18([x,y,z],[r,phi,z],[[ri,ra],[0,2*Pi],[0,h]],rho);
# Gerader Kreiskegel mit dem Radius a und der Höhe h
> display(cone([0,0,0],1,3),orientation=[45, 70],scaling=constrained,lightmodel=light2,shading=xyz,labels=["x1","x2","x3"]);
> r:='r': phi:='phi': z:='z': a:='a':
> x:=r*cos(phi); y:=r*sin(phi);
> Proc_Calc_18([x,y,z],[r,phi,z],[[0,a/h*z],[0,2*Pi],[0,h]],rho);
# Gerader Kegelstumpf:  r(z = 0) = a1, r(z = h) = a2.
> a:=unapply(a1+(a2-a1)/h*z,z);
> a1:=1: a2:=2: h:=5:
> plot3d(a(z),phi=0..2*Pi,z=0..h,coords=cylindrical,orientation=[75,15,120],scaling=constrained,labels=["x1","x2","x3"]);
> x:=r*cos(phi); y:=r*sin(phi);
> Proc_Calc_18([x,y,z],[r,phi,z],[[0,a(z)],[0,2*Pi],[0,h]],rho);
# Kugel  mit Radius a
> display(sphere([0, 0, 0], 3.0), scaling=constrained,axes=boxed,labels=["x1","x2","x3"]);
> x:=r*cos(phi)*sin(theta); y:=r*sin(phi)*sin(theta); z:=r*cos(theta);
> a:='a':
> Proc_Calc_18([x,y,z],[r,theta,phi],[[0,a],[0,Pi],[0,2*Pi]],rho);
# Halbkugel  mit Radius a
> display(hemisphere([0, 0, 0], -1),scaling=constrained,axes=boxed,orientation=[-145,75,0],color = green,labels=["x1","x2","x3"]);
> x:=r*cos(phi)*sin(theta); y:=r*sin(phi)*sin(theta); z:=r*cos(theta);
> Proc_Calc_18([x,y,z],[r,theta,phi],[[0,a],[0,Pi/2],[0,2*Pi]],rho);
# Im folgenden Beispiel (Spielkreisel) verläuft der Radius a des Rotationskörpers längs der Koordinate x3 unstetig. 
> r:='r': phi:='phi': z:='z': a:='a':
> gl:=piecewise(z < h1,r1/h1*z,z < h1+h2 ,r2,z < h1+h2+h3,r3);
> a:=unapply(gl,z);
# Alle Längenmaße werden in cm eingegeben.  Für die Erdbeschleunigung folgt g = 981 cm/s2.
> h1:=2: h2:=h1/2: h3:= 6/5*h1: r1:=h1: r2:=2*h1: r3:=h1/4:
> plot3d(a(z),phi=0..2*Pi,z=0..h1+h2+h3,coords=cylindrical,orientation=[160,85,-130],scaling=constrained,labels=["x1","x2","x3"]);
> x:=r*cos(phi); y:=r*sin(phi);
> Proc_Calc_18([x,y,z],[r,phi,z],[[0,a(z)],[0,2*Pi],[0,h1+h2+h3]],rho);
> 
;
# 
