/* Maxima:     
Computeralgebra-System (OpenSource)
Webseite: maxima.sourceforge.net  
*/

/* - Am Ende eines Befehls ; (wenn man die Ausgabe sehen will) oder $ (wenn nicht)
   - Ausführen mit Shift+Enter
*/

load(draw);

/*  Parametrisierung der Kugel durch Kugelkoordinaten */
x: a*cos(theta[1])*cos(theta[2]);
y: b*cos(theta[1])*sin(theta[2]);
z: c*sin(theta[1]);
F:[x,y,z];

/* Zeichnen für a=1,b=2,c=3 
   s. auch http://www.austromath.at/daten/maxima/zusatz/Grafiken_mit_Maxima.pdf
*/
ellipsoid:parametric_surface(cos(theta)*cos(phi), 2*cos(theta)*sin(phi), 3*sin(theta), theta, -%pi, %pi, phi, 0,2*%pi);

draw3d(nticks=400,color=orange,ellipsoid);
/* nticks -- Anzahl der Stützstellen */

/* Berechnen der induzierten Metrik*/
g[i,j]:= diff(F, theta[i]).diff(F, theta[j]);

/* Gibt die zugehörige Matrix (Name: lg) der induzierten Metrik  aus
-- trigsimp performs simplifications that use the Pythagorean identity
-- genmatrix (a, i_2, j_2, i_1, j_1): Generiert eine Matrix aus einem Array a. Das erste Element der Matrix ist der Wert a[i_1,j_1] und das letzte Element der Matrix ist a[i_2,j_2]. 
*/
lg:trigsimp(genmatrix(g,2,2,1,1));

/* Setzt ug auf die inverse Matrix von lg*/
ug:invert(lg)$

/* Laden des ctensor-Pakets - das ist ein Paket zum Manipulieren von Tensoren (wie Metrik, Krümmungstensor...)
Anleitung/Definition (auch um die Index- und Vorzeichenkonventionen zu sehen)  
unter http://maxima.sourceforge.net/docs/manual/maxima_26.html
-- darunter uebergeben wir die Dimension und die Koordinaten*/
load(ctensor);

dim:2;
ct_coords:[theta[1], theta[2]];

/* Berechnet die Christoffelsymbole, Krümmungen ... - 
mcs berechnet \Gamma_{ij}^k 
lcs berechnet \Gamma_{ijl}:=\Gamma_{ij}^kg_{kl}
Es werden immer nur die Ausdrücke angezeigt, die nicht Null sind.
*/
christof(mcs)$

mcs[1,1,1];

/* Riemanntensor, Riccitensor, Skalarkrümmung 
    (Reihenfolge Indizes wird von der bei uns abweichen, s. http://maxima.sourceforge.net/docs/manual/maxima_26.html )*/
riemann(1)$
ricci(1)$
scal:scurvature()$

/*Skalarkrümmung für Kugel*/
subst([a=1,b=1,c=1], ''scal)$

/* trigsimp -- Vereinfachen trigonometrischer Ausdrücke 
   (Allg. Übersicht von Maximabefehlen s. http://www.siart.de/service/maxima-referenz.pdf )
   % bezieht sich auf den zuletzt berechneten Ausdruck */
trigsimp(%);

/* Wir wollen den Teil der Geodätengleichung aufstellen, der -Gamma_{ij}^k \dot x^i \dot x^j ist. Hier sei u[i]=\dot x^i. */
z1:-sum(sum(mcs[i,j,1]*u[i]*u[j], j,1,2),i,1,2);
z2:-sum(sum(mcs[i,j,2]*u[i]*u[j], j,1,2),i,1,2);

/* Wahl der Parameter */
b1:subst([a=1,b=2,c=3], ''z1);
b2:subst([a=1,b=2,c=3], ''z2);
/* " - steht hier für evaluate...*/

a1:subst(delta, u[2], subst( gamma, u[1], subst(beta, theta[2],subst(alpha, theta[1], b1))));
a2:subst(delta, u[2], subst( gamma, u[1], subst(beta, theta[2],subst(alpha, theta[1], b2))));
/* aus irgendeinem Grund will der nächste Befehl darunter u[1]... nicht als Variablen, 
    deshalb benenn ich hier um -- geht bestimmt besser/einfacher */

/* Geodätengleichung für ein c(t)=F(alpha(t),beta(t)) (in den Koordinaten (alpha=theta[1], beta=theta[2]) - als System erster Ordnung: 
    Setze alpha'(t)=gamma(t), beta'(t)=delta(t)*/
   
sol: rk([gamma, delta,  a1, a2],[alpha, beta, gamma,delta],[0.,0.,1.,0.],[s,0,6,0.01])$
    /* Dies löst die Geodätengleichung numerisch (rk=Runge-Kutta Verfahren) als erstes Ordnungssystem
     --   gamma=alpha′(=theta_1') und delta=beta'(=theta_2') sind die neuen Variablen
     --   [gamma, delta,  a1, a2],[alpha, beta, gamma,delta] bedeutet:
            alpha'=gamma; beta'=delta, gamma'=a1, delta'=a2
     --    [0.,0.,1.,0.] sind  die gewählten Anfangswerte für [alpha, beta, gamma,delta] 
     --   [s,0,6,0.01] bedeutet, dass die DGL für s in [0,6] und in diskreten Schritten der Länge 0.01 gelöst wird. 
     --   Das Ergebnis des Runge-Kutta-Verfahrenwird als Liste in sol gespeichert. 
          Diese Liste besteht aus Fünftupeln der Form[s,alpha(s),beta(s), gamma(s), delta(s)].
    */

pts:makelist([cos(p[2])*cos(p[3]),2*cos(p[2])*sin(p[3]), 3*sin(p[2])],p,sol)$
/* aus obigem soll wird hier eine Liste der F(alpha(s), beta(s)) gemacht*/
curve:points(pts)$

draw3d(nticks=200,surface_hide=true,color=orange,ellipsoid, color=blue, curve);

sol2: rk([gamma, delta,  a1, a2],[alpha, beta, gamma,delta],[0.,0.,0.,1.],[s,0,6,0.01])$
pts2:makelist([cos(p[2])*cos(p[3]),2*cos(p[2])*sin(p[3]), 3*sin(p[2])],p,sol2)$
curve2:points(pts2)$

draw3d(nticks=200,surface_hide=true,color=orange,ellipsoid, color=blue, curve, color=red, curve2);

sol3: rk([gamma, delta,  a1, a2],[alpha, beta, gamma,delta],[0.,-%pi/2.,1.,0.],[s,0,6,0.01])$
pts3:makelist([cos(p[2])*cos(p[3]),2*cos(p[2])*sin(p[3]), 3*sin(p[2])],p,sol3)$
curve3:points(pts3)$
draw3d(nticks=200,surface_hide=true,
    color=orange,ellipsoid,  
    color=blue, curve, 
    color=red, curve2, 
    color=green,  curve3)$

sol4: rk([gamma, delta,  a1, a2],[alpha, beta, gamma,delta],[0.,0,1.,-1.],[s,0,6,0.01])$
pts4:makelist([cos(p[2])*cos(p[3]),2*cos(p[2])*sin(p[3]), 3*sin(p[2])],p,sol4)$
curve4:points(pts4)$
draw3d(nticks=200,surface_hide=true,
    color=orange,ellipsoid,  
    color=blue, curve, 
    color=red, curve2, 
    color=green,  curve3,
    color=black,  curve4)$

pause([options]):=block([tsecs],
    tsecs:assoc('pausetime,options,0),
    if tsecs=0 then
       read("Execution Paused...enter any character then CTRL-ENTER")
    else(
        disp(sconcat("paused for ", tsecs," seconds")),
       ?sleep(tsecs)),
    return("")
);

pause(pausetime=5);