(* MECANICA TEORICA / FC - 
            UABC / PROF .  HECTOR ACEVES *)
(* PENDULO ROTANTE *)

(* Cargar primero el paquete de Calculo Variacional *)
Needs["Calculus`VariationalMethods`"]
  

(* Limpiar todas las variables *)
ClearAll["Global`*"]
Off[General::spell]

(* La energia cinetica,  
          potencial y el Lagrangiano .  
              La energia potencial es solo gravitacional  *)
T=(1/2) m (((x'[t])^2)+((y'[t])^2));
Vgr = -m g y[t];
polarrule = {x -> ((a  Cos[Ω  #]   +  b * Sin[θ[#]]) &),  
              y -> ((-a  Sin[Ω  #]   +  
                        b * Cos[θ[#]]) &)};
Vx = x'[t] /.polarrule
Vy =  y'[t] /.polarrule
Kinetic  =  
          T  /. polarrule  // Simplify
Potential  =  Vgr  /. polarrule
L = T - Vgr  /. 
              polarrule  //Simplify

(* Encontramos la ecuacion para θ[t] *)
EQN=Simplify[ EulerEquations[(1/2) m (2 b g Cos[θ[t]]+a (a (Ω^2) - 
                          2 g Sin[t Ω])-2 a b Ω Sin[t Ω - θ[t]] θ'[t]+(b^2) ((θ'[t])^2)) ,
θ[t],t] ]

eqMotion = Solve[EQN, θ''[t]]  /. 
                Rule ->  Equal  // Simplify  // Flatten // 
          ExpandAll  // ColumnForm

(* Graficamos una solucion,  
          con condiciones iniciales (initialc) tales que el pendulo  no suba mas que la posicion de su punto de suspencion .  Se obtienen las ecuaciones de movimiento con las condiciones iniciales (neqc) y luevo se resuelven numericamente (nSolc) .  Despues se define funcion interpolante θc para poder graficarla posteriormente*) 
valuesc = {a -> 1, b -> 2,  
          g -> 9.8, 
          m -> 
            0.5,  Ω  ->  (1.0*
                Sqrt[9.8/2]) }

      τ  =  1/Sqrt[g/b]     /. valuesc

      initialc = {θ[0] == π/2, θ'[0] ==  
                0.0} /. valuesc;

      neqc = {eqMotion, initialc} /. Rule -> Equal  /. 
            valuesc//Flatten

nSolc = 
        NDSolve[{θ''[t] == 
                2.45 Cos[2.21359 t  -  θ[t]] - 
                  4.9 Sin[θ[t]],  θ[0] == 
                1.5708,  θ'[0] == 0.0}, θ[t], {t, 0, 
              100*τ}, MaxSteps -> 3000 ] // 
          Flatten;
θc[t_] = θ[t] /. nSolc;

(* Graficamos la variacion temporal de la posicion y luego de la coordenada angular .  
              Y luego las combinamos para ver el movimiento de la particula masa en el plano XY *)

Plot[ θc[t]   // Evaluate, {t, 0, 
            50*τ}, 
PlotStyle -> {Hue[ .7]}  , 
          TextStyle -> {FontFamily -> "Times-Italic", 
              FontSize -> 14},  
          AxesLabel ->  {"t", "θ[t]" }    ];

grafica1 = 
        ParametricPlot[ {  1*Cos[2.21359* t]  +  2*Sin[θc[t]], 
              1*Sin[2.21359*t] - 2*Cos[θc[t]]}
  // 
            Evaluate , {t, 0, 
            50*τ}, 
PlotStyle -> Hue[0.7],  
          TextStyle -> {FontFamily -> "Times-Italic", 
              FontSize -> 12},  
AspectRatio -> 
            1, AxesLabel ->  {"x", "y"},  
          PlotRange -> {{-3.5, 3.5}, {-4, 3}},  
          Epilog -> {Circle[{0, 0}, #] & /@ {1, 3},   PointSize[
                0.04], Hue[1], Point[{3, 0}],   PointSize[0.03], Hue[0.1], 
              Point[{1, 0}]}
];

(* Trayectoria en el espacio fase :   (θ, θ') *)
pt1 = 
          ParametricPlot[ {θc[t], θc'[t]  } // 
              Evaluate, {t, 0, 
              50*τ}, 
PlotStyle -> {Hue[ .7]}, 
            TextStyle -> {FontFamily -> "Times-Italic", 
                FontSize -> 12},  
            AxesLabel -> {"θ", "θ'"}, 
            PlotRange -> {{-2, 2}, {-4, 4}}];

(* Animacion para la primer condicion inicial,  
          con X numero de frames .  
              Para sacar ya una movie mas decente cambiar,  
          por ejemplo frames = 200,  
          y el nombre del archivo y su direccion a donde se quiera .  El rango tambien puede cambiarse a 0, 50 tau, 
          Aqui solo se ponen dos cuadros para ilustrar. *)
\(line[
          t_] := 
Show[\ 
Graphics[
{\ \ \ {\ \ \ \ \ Line[{\ \ {0, 
                      0}, {1\ *Cos[\ 2.21359  t], 
                      1*Sin[2.21359  
                            t\ ]}\ \ }]\ , 
{\ \ \ \ PointSize[0.02], Hue[0.1], 
                  Point[\ {\ 1*Cos[2.21359\ t]\ , 
                      Sin[2.21359\ t]\ }\ ]\ \ \ \ }, 
\ 
                Line[{\ \ {1\ *Cos[\ 2.21359  t], 
                      1*Sin[2.21359  t\ ]}, 
{\ \ 1*
                          Cos[2.21359*\ t]\  + \ 2*Sin[θc[t]], 
                      1*Sin[2.21359*t] - 
                        2*Cos[θc[
                              t]]}\ \ }], 
{\ \ \ \ PointSize[0.05], Hue[0.6], 
                  Point[\ \ {\ \ 1*Cos[2.21359*\ t]\  + \ 
                        2*Sin[θc[t]], 
                      1*Sin[2.21359*t] - 
                        2*Cos[θc[
                              t]]}\ \ \ \ ]\ \ \ \ }
\ }
}, 
PlotRange -> {{\(-3.5\), 
                  3.5}, {\(-3.5\), 3.5}}, \ Frame -> False, 
            AspectRatio -> 1, 
            TextStyle -> {FontFamily -> "Times-Italic", 
                FontSize -> 12}, \ 
            Background\  -> \ GrayLevel[0.1], 
            Epilog -> {\(Circle[{0, 0}, #] &\) /@ {1}, \ \ PointSize[
                  0.02], Hue[0], 
                Point[{0, 
                    0}]}
]
]


      \(frames = 2;\)

      \(movie = \(line[#] &\)\  /@ \ 
            Range[0, 50*τ, 
              50*τ/\((\ frames - 1)\)];\)\ 

      Export["c:/test.gif", movie, 
        ConversionOptions -> {"AnimationDisplayTime" -> 
              0.5\ \ , \ "Loop" -> False}, 
        ImageSize -> {500, 500}]\)