:: perturbed.tm :: :: @(#)perturbed.tm 1.2 U of T Physics 10/25/92 :: :: This is the mathlink template file for the perturbed pendulum program. :: It specifies the pattern for the Mathematica function Pendulum[] and the :: way in which the Pendulum[] arguments will be translated into arguments :: to the C function pendulum(). :: :: The function can be called from Mathematica in these ways: :: :: Pendulum[{q0, p0}, {pts, dt, sample}, {h, gamma, shift, every}, double] :: Pendulum[{q0, p0}, {pts, dt, sample}] :: Pendulum[{q0, p0}, {pts, dt, sample}, {h, gamma, shift, every}] :: Pendulum[{q0, p0}, {pts, dt, sample}, double] :: :: {q0, p0} are the initial conditions; pts is the number of points to :: return; dt is the timestep length; sample is (in the unperturbed case) :: the sampling rate: every sampleth point will be returned, the rest :: discarded; h and gamma are the perturbation parameters; shift is the :: position of the poincare section; every is (in the perturbed case) whether :: to save every timestep or every period; double is either True or False :: depending on whether you want the phase portrait to cover -2Pi < q < 2Pi :: or just -Pi < q < Pi, respectively. :Begin: :: Specify what the name of the C function to be called is: :Function: pendulum :: Define what the arguments to Pendulum[] are and what patterns they must :: match. (The contortions with Im[] etc. are necessary to allow any real :: numbers. See page 238 of the Mathematica book.) :Pattern: Pendulum[ point0:{ (q0_ /; NumberQ[N[q0]] && Im[N[q0]]==0), (p0_ /; NumberQ[N[p0]] && Im[N[p0]]==0) }, time:{ (pts_Integer /; pts > 0), (dt_ /; NumberQ[N[dt]] && Im[N[dt]]==0), (sample_Integer /; sample > 0) }, perturb:{ (h_ /; NumberQ[N[h]] && Im[N[h]]==0), (gamma_ /; NumberQ[N[gamma]] && Im[N[gamma]]==0), (shift_ /; NumberQ[N[shift]] && Im[N[shift]]==0), (every_ /; every === TimeStep || every === Period) }, (double_ /; double === True || double === False) ] :: Then define what the arguments to the C function will be. All the reals :: are converted to C doubles, the integers are sent as is, every (which has :: a symbolic value) is converted to an integer with Switch, and double :: (which has a boolean value) is converted to an integer with If. :Arguments: { N[q0, 25], N[p0, 25], pts, N[dt, 25], sample, N[h, 25], N[gamma, 25], N[shift, 25], Switch[every, TimeStep, 1, Period, 0], If[double, 1, 0] } :ArgumentTypes: { Real, Real, Integer, Real, Integer, Real, Real, Real, Integer, Integer } :: The return value of Pendulum[] will not be the return value of the C :: function but will instead be explicitly sent down the link, in other :: words sent manually. This is because there is no way of returning a :: Mathematica list from a C function. :ReturnType: Manual :End: :: Set up some alternate calling sequences for the Mathematica funtion: :Evaluate: Pendulum[point0_, time_] := Pendulum[point0, time, { 0.0, 0.0, 0.0, TimeStep}, True] :Evaluate: Pendulum[point0_, time_, perturb:{_, _, _, _}] := Pendulum[point0, time, perturb, True] :Evaluate: Pendulum[point0_, time_, (double_ /; double === True || double === False)] := Pendulum[point0, time, { 0.0, 0.0, 0.0, TimeStep}, double] :: A help message so that you can type "?Pendulum" from within Mathematica :: to get information about the function: :Evaluate: Pendulum::usage = "Pendulum[{q0, p0}, {pts, dt, sample}, {h, gamma, shift, every}, double] generates a list of {q, p, E} triples. The initial position and momentum are denoted by \"q0\" and \"p0\"; the total number of points to return (which is not the same as the number of points to calculate) is given by \"pts\"; \"dt\" is the timestep size; \"sample\" (which is only used if there is no perturbation) is the sampling period (in units of dt of course), so if it is 1 then every point will be returned. If gravitational perturbation is desired then \"h\" and \"gamma\" will be non-zero; these are the amplitude and frequency of the perturbation. If you want to take a poincare section, then \"every\" should be Period, and the poincare section will be taken at gamma t == 2 m Pi + shift; otherwise it should be TimeStep. If you want to return every point twice (in order to plot a phase portrait with -2Pi < q < 2Pi instead of -Pi < q < Pi) then \"double\" should be True. The arguments \"{h, gamma, shift, every}\" and \"double\" default to no perturbation (i.e. {0, 0, 0, TimeStep}) and True, respectively."