(* ::Package:: *)

(* Rechnernutzung WS23/24
"Upside-down" Pendel
>= Mathematica 8 notwendig 

Mathematisches Pendel, wobei sich der Aufhaengepunkt
in y-Richtung mit h(t) = a * sin(omega t) bewegt.

Ruhelage fuer t=0, theta=0: (x,y) = (0,l)
*)

(* siehe auch: https://www.youtube.com/watch?v=5oGYCxkgnHQ *)

(* Lagrange-Funktion -> Bewegungsgleichungen *)

Clear[x, y, h, v, U, T, L, dLdTheta, dLdThetaDot, ELEQ];

x[t_] := l*Sin[theta[t]]
y[t_] := h[t] + l*Cos[theta[t]]
h[t_] := a*Sin[omega*t]
v = D[{x[t], y[t]}, t];
T = 1/2*m*v.v;
U = m*g*y[t];
L = Simplify[T - U]
dLdTheta = D[L, theta[t]] // Simplify;
dLdThetaDot = D[L, theta'[t]] // Simplify;
ELEQ = dLdTheta == D[dLdThetaDot, t];
ELEQ = Simplify[ELEQ, {l != 0, m != 0}]


(* plotPendulum: zeichne Pendel zum Zeitpunkt t
   mit Auslenkung theta = theta(t)
   subs: Ersetzungsregel fuer Pendellaenge l, Auslenkung a und Kreisfrequenz omega *)
plotPendulum[theta_, t_, subs_List ] := Module[
  {p1, p2, le, ps, ra, lines, dots, ceil, text, text2, all, pendulum},
  (* Punkte: p1: Aufhaengepunkt; p2: Ende Pendel *)
  p1 = {0, h[t]} /. subs;
  p2 = {l*Sin[theta], h[t] + l*Cos[theta]} /. subs;
  (* Groessenskala *)
  le = (l + a) /. subs;
  (* dargestellter Bereich *)
  ra = le*{{-1.05, 1.05}, {-2.15, 2.15}};
  (* Aufhaengung, Linien und Punkte *)
  ceil  = {Brown, Thickness[0.05], Line[{{-l,0}/.subs, {l,0}/.subs}]};
  lines = {Black, Thickness[0.01], Line[{{0, 0}, p1}], Thickness[0.005], Line[{p1, p2}]};
  dots  = {Black, PointSize[0.04], Point[{{0, 0}, p1}], Red, PointSize[0.1], Point[{p2}]};
  (* Beschriftung *)
  text = {Black, Background -> White,
    Text[m, p2 + {1, 0}, {-1, 0}],
    Text[l, p1 + (p2 - p1)/2 - {0.05*le, 0}, {1,0}],
    Text[" t = " <> ToString[t] <> " s ",{-l*0.5,a}/.subs]};
  (* Fasse alles zusammen: *)
  all = Join[ceil, lines, dots, text];
  pendulum = Show[Graphics[all], PlotRange -> ra];
  Return[pendulum]
  ]


(* Beispiel zur Verwendung von plotPendulum[] *)
plotPendulum[2.001,Pi/4//N,{l->10,a->1.5,omega->2*Pi}]


(* Loese Differentialgleichung.
   ld: Pendellaenge
   ad: Auslenkung Aufhaengung
   omegad: Frequenz Aufhaengung
   theta0, thetaDot0: Anfangswerte
   Animiere Loesung. *)
upsideDownPendulum[{ld_, ad_, omegad_},{theta0_, thetaDot0_}] := Module[
  {gVal, tM, temp, sol, thetaSol, params, res},
  gVal = 9.81;
  tM = 60; (* timeMax *);
  (* Parameter einsetzen in Bewegungsgleichungen: *)
  params = {g -> gVal,l -> ld, a -> ad, omega -> omegad}; 
  temp = {ELEQ} /. params;
  (* Anfangsbedingungen hinzufuegen: *)
  temp = Join[temp, {theta[0] == theta0, theta'[0] == thetaDot0}];
  (* DGL loesen und Fkt. definieren: *)
  sol = NDSolve[temp, theta[t], {t, 0, tM}];
  thetaSol = theta[t] /. sol[[1]];
  (* Pendel animieren: *)
  res = Animate[plotPendulum[thetaSol /. {t -> tt}, tt, params],
   {tt, 0, tM, 0.01 (* timeInt *)}];
  Return[res]
  ]


(* "normale" Schwingung *)
upsideDownPendulum[{10, 1.5, 1.0}, {2.7, 0.01}]


(* upside-down Pendel: variiere Frequenz der Aufhaengung *)

n = 1;
upsideDownPendulum[{10, 1.5, 2*Pi*n}, {0.3, 0.01}]