
(* : Package Version : 1.0*)
(* : Copyright : Copyright 2007, Jesse Welch*)
(* : Name : Graphics`Pendulum`*)
(* : Author : Jesse Welch*)


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";

ThetaPlot::usage = "ThetaPlot[{drivforce, drivfreq, drivphase, dampen}, {initt, fint, steps}, {initth, initw}]
 This gives graphical output of theta vs T.";

OmegaPlot::usage = "OmegaPlot[{drivforce, drivfreq, drivphase, dampen}, {initt, fint, steps}, {initth, initw}]
 This gives graphical output of omega vs T.";

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[driveforce_, drivefreq_, dampen_, phase_] :=
Module[{q = driveforce, p = phase, b = dampen, w = drivefreq},
   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]};
];

ThetaPlot[{driveforc_, drivefre_, dampe_, phas_}, {int_, fint_, step_}, {intx_, inty_}] :=
Module[{ddforce = driveforc, ddrivfre = drivefre, ddampe = dampe, pphase = phas, intt = int,
fintt = fint, intxx = intx, intyy = inty, stepp = step,finishh = finish},
   SetPendParam[ddforce, ddrivfre, ddampe, pphase];
   solu = Runge[intt, fintt, intxx, intyy, stepp];
   newsolu2 = Table[{solu[[i]][[3]], Mod[solu[[i]][[1]] - \[Pi], 2\[Pi]]}, {i, 1, stepp}];
   ListPlot[newsolu2];
];

OmegaPlot[{driveforc_, drivefre_, dampe_, phas_}, {int_, fint_, step_}, {intx_, inty_}] :=
Module[{ddforce = driveforc, ddrivfre = drivefre, ddampe = dampe, pphase = phas, intt = int,
fintt = fint, intxx = intx, intyy = inty, stepp = step,finishh = finish},
   SetPendParam[ddforce, ddrivfre, ddampe, pphase];
   solu = Runge[intt, fintt, intxx, intyy, stepp];
   newsolu2 = Table[{solu[[i]][[3]], Mod[solu[[i]][[2]] - \[Pi], 2\[Pi]]}, {i, 1, stepp}];
   ListPlot[newsolu2];
];

End[]
EndPackage[]
