(* :Package Version: 1.0 *)

(* :Copyright: Copyright 2007, Jesse Welch *)

(* :Name: Graphics`Pendulum` *)
 
(* :Author: Jesse Welch *)

(* Shadowing prevention *)

(*(* little utility that is really a holdover from an early version of
   the package; returns the second arg if the first arg is Automatic,
   otherwise return the first arg. *)
automatic[x_, value_] :=
    If[x === Automatic, value, x]*)


Unprotect["Graphics`Pendulum`*"];
ClearAll["Graphics`Pendulum`*"];

BeginPackage["Graphics`Pendulum`"]

Runge::usage = "Runge[initialtime, finaltime, initialx, initialy, numberofsteps] 
 Uses a fourth order runge kutta algorithm to integrate third order ODE's.  The output is in the form of a table of
 3 x steps values.  The first of the three values of the set is the x function, the second value is the y function, 
 and the third value is the time value, for each step";

SetPendParam::usage = "SetPendParam[drivforce,dampen,drivefreq] 
 This sets the values of the parameters of your derivative functions for the ODE solver.  This specifically puts in
 the derivative functions for the chaotic pendulum,though you could make your own function to set the functions and
 parameters of different derivative functions to solve a different problem";

PlotPendulum::usage = "PlotPendulum[string, drivforce, dampen, drivfreq, init\[Theta], init\[Omega], steps] 
 This gives graphical output of various graphs, depending on the input.  The output is a listplot, and the string
 variable input has the options ''th vs T'',''w vs T'',''th vs w'',''poin th vs w'' where the two single quotes 
 indicate the use of double quotes.  The next three inputs are used to define the parameters of the system to be
 solved, and the remaining inputs are those used by Runge to calculate the solution.  The output hides the first
 20% of the solutions to allow for the transient behavior to disappear.";
 
AnimatePendulum::usage = "AnimatePendulum[string, drivforce, dampen, drivfreq, init\[Theta], init\[Omega],
 steps, paraminit, paramfinal, paramdx] 
 This gives graphical output of various graphs, depending on the input.  The output is a series of listplots, and
 the string variable input has the options ''th vs T'',''w vs T'',''th vs w'',''poin th vs w'' where the two single
 quotes indicate the use of double quotes.  The next three inputs are used to define the parameters of the system to be
 solved, and the remaining inputs are those used by Runge to calculate the solution.  The output hides the first
 20% of the solutions to allow for the transient behavior to disappear.";

PeriodicPendulum::usage = "PeriodicPendulum[string,steps].  This will output whichever graph you tell it to input for 
however many steps you tell it, and output a listplot of periodic motion of the dampened driven pendulum.";

Begin["`Private`"]

Runge[initt_, finalt_, initx_, inity_, steps_] :=
Module[{it = initt, ft = finalt, j, m = steps}, 
   h = (ft - it) / m;
   T = Table[ a + (j - 1) h, {j, 1, m + 1}];
   R = Table[{ initx, inity, initt} , {j, 1, m + 1}];
   For[j = 1, j < m+1  , j++, 
      k1 = h F[T[[j]], R[[j]]];
      k2 = h F[T[[j]] + h/2, R[[j]] + k1 / 2];
      k3 = h F[T[[j]] + h/2, R[[j]] + k2 / 2];
      k4 = h F[T[[j]] + h, R[[j]] + k3];
      R[[j + 1]] = R[[j]] + 1/6 (k1 + 2 k2 + 2 k3 + k4)];
   Return[R]; 
];
      
SetPendParam[drivforce_, dampen_, drivfreq_,drivephase_] :=
Module[{q = drivforce, b = dampen, w = drivfreq, p = drivephase},
   f1[t_, x_, y_, z_] = y;
   f2[t_, x_, y_, z_] = q Cos[z + p] - Sin[x] - b y;
   f3[t_, x_, y_, z_] = w;
   F[t_, {x_, y_, z_}] = {f1[t, x, y, z], f2[t, x, y, z], f3[t, x, y, z]};
]; 

PlotPendulum[{axes1_, axes2_}, {drivforce_, dampen_, drivfreq_, drivephase_}, {initth_, initw_}, steps_, mfinalt_] :=
Module[{ax1 = axes1, ax2 = axes2, pa = drivforce, pb = dampen, pc = drivfreq, pd = drivephase, n = steps,tf = mfinalt},
   g = Evaluate[If[ax1 == "t", 1, 0,1]];
   G = Evaluate[If[ax2 == "t", 1, 0,1]];
   s = If[ax1 == "t", 3, If[ax1 == "w", 2, 1]];
   r = If[ax2 == "t", 3, If[ax2 == "w", 2, 1]];
   start = n * 0.2;
   SetPendParam[pa,pb,pc,pd];
   sol = Runge[0, tf, initth, initw, n];
   graph = Table[
      {(sol[[i]][[r]] * g) + ((Mod[sol[[i]][[r]] - Pi, 2 * Pi]) * (1 - g))
      ,(sol[[i]][[s]] * G) + ((Mod[sol[[i]][[s]] - Pi, 2 * Pi]) * (1 - G))}
   ,{i,1,n}];
   ListPlot[graph];
];

AnimatePendulum[string_, drivforce_, dampen_, drivfreq_, initt_, finalt_, initth_, initw_, steps_, start_, stop_, wid_] :=
Module[{new = string, pa = drivforce, pb = dampen, pc = drivfreq, tt = initt, tf = finalt, n2 = steps, init = start, final = stop},
   Do[PlotPendulum[new, o, pb, pc, tt, tf, initth, initw, n2], {o, init, final, wid}]
];

PeriodicPendulum[string_,steps_]:= 
Module[{},
PlotPendulum[string, 0.9, 0.5, 2/3, 0, 0, steps]
];

End[]
EndPackage[]