#                                    MAT 331        Project 2
#                                        Fun on a Swing
#                                            part 1
# 
# 
> with(share): readshare(ODE,plots):
--------------------------------------------------------------------------------
# We tell maple about the system of differential equations which describe the swing, and give values to all of the 
# constants. Here we have a weight of mass  m   on rigid pendulum of length  l .   
# This model accounts for both friction (given by the  
# 

                                    epsilon nu
                                  - ----------
                                         m
# term)  and for "pumping"

                              arctan(f/l) arctan(5 nu)
                          - 4 ------------------------
                                         Pi
#  achieved by moving the mass  f  units from the end of the rod (the above term gives the angular offset achieved by 
# this motion). 
# 
> thdot :=  (t,theta,nu) -> nu;\
nudot:=  (t,theta,nu) -> -g/l * sin(theta -2/Pi*arctan(f/l)*arctan(5*nu)) - epsilon/m * nu;\
\
m:=150; \
g:= 32;\
epsilon := 4;\
f := 1;\
\
l := 10;

                           thdot := (t,theta,nu) -> nu

  nudot :=

                                        arctan(f/l) arctan(5 nu)
                        g sin(theta - 2 ------------------------)
                                                   Pi               epsilon nu
      (t,theta,nu) -> - ----------------------------------------- - ----------
                                            l                            m

                                    m := 150

                                     g := 32

                                  epsilon := 4

                                     f := 1

                                     l := 10
# 
--------------------------------------------------------------------------------
# Now we are set to make a picture or two of the motion of the swing. 
# 
# First, lets make a phase plot of the case with no friction, and no pumping.  To remove those forces, we set epsilon 
# and f to zero.
# 
> epsilon:=0: \
f:=0:\
phaseplot([thdot, nudot],-2*Pi..2*Pi,-5..5,{seq([0,0,0.5*nu], nu=-10..10)},\
                 stepsize=0.1,numsteps=200);\


   ** Maple V Graphics **

# Here we can see  a number of solutions for which the swing oscillates (the elliptical solutions in the middle),  3 
# solutions corresponding to continuous clockwise motion (the "running solutions" at the top), and 3  with continuous 
# counterclockwise solutions at the bottom. 
# 
# Here is a plot of  the angle as a function of time for two such swings starting with an initial angle of 0.  One has an 
# initial angular velocity of 4 (which is a running solution, so the angle continues to increase),  and the other with an 
# initial angular velocity of 2 (which oscillates). 
> plot({makelist(rungekuttahf([thdot,nudot],[0,0,4],0.1,100),1,2),\
         makelist(rungekuttahf([thdot,nudot],[0,0,2],0.1,100),1,2)}, t=0..10, theta=-3..20);\


   ** Maple V Graphics **

# 
--------------------------------------------------------------------------------
# We can write a maple procedure to animate these solutions, so we can see the swing actually moving.  This takes as 
# input the result of a call to rungekuttahf, and shows a series of plots in sequence.  This can be rather memory and 
# time consuming to run, because of the large number of plots it computes.
# 
> SwingMovie := proc (theTrajectory) \
  local numframes, i, theTitle;\
\
  numframes := op(2,op(2,eval(theTrajectory)));\
  theTitle := `swing with l=` .l;\
\
  plots[display]([\
    seq(\
      plot([0,0,  l*sin(theTrajectory[i][2]), -l*cos(theTrajectory[i][2])],\
            -l-1..l+1, -l-1..l+1, scaling=constrained, axes=boxed,\
            thickness=2,title=theTitle),\
       i=0..numframes)],\
    insequence=true);  \
\
end:
# 
# Here are two samples, namely, a  movie for each of the two swings above:
# 
> SwingMovie(rungekuttahf([thdot,nudot],[0,0,4],0.1,100));
> SwingMovie(rungekuttahf([thdot,nudot],[0,0,2],0.1,100));
--------------------------------------------------------------------------------
# 
# Now, lets add a little friction and see how things change. 
# 
# 
> epsilon:=4: \
f:=0:\
phaseplot([thdot, nudot],-2*Pi..2*Pi,-5..5,{[0,0,4],[0,0,2]},\
                 stepsize=0.1,numsteps=500);\


   ** Maple V Graphics **

# Here, we show only the two solutions investigated above, namely [0,0,4] and [0,0,2], and we run them for a longer 
# time (50 time steps). Note that the oscilatory solution winds down toward 0.   What happens to the other one?
# You could use either 
> plot({makelist(rungekuttahf([thdot,nudot],[0,0,4],0.1,500),1,2);
#  or
> SwingMovie(rungekuttahf([thdot,nudot],[0,0,4],.1,500));
# (or both) to find out, but I'll leave that to you.  Here is a graph of theta for the inner solution:
> nopump := rungekuttahf([thdot,nudot],[0,0,2],.1,500):\
plot(makelist(nopump,1,2),t=0..50,theta=-2..2);

   ** Maple V Graphics **

--------------------------------------------------------------------------------
# Now, lets allow the swinger to start pumping.  Once again, we take our equations, but set f to be nonzero.
# 
> epsilon:=4: \
f:=1:\
phaseplot([thdot, nudot],-2*Pi..2*Pi,-5..5,{[0,0,4],[0,0,2]},\
                 stepsize=0.1,numsteps=500);\


   ** Maple V Graphics **

# What is different here?  At first glance, it looks like the solution starting at [0,0,2] goes to 0 even quicker, but a 
# closer look shows that the spiral moves OUTWARDS.  This is because phaseplot is showing us the past history, as 
# well.  Lets make a graph of the angle as a function of time, to illustrate this.
# 
> pump := rungekuttahf([thdot,nudot],[0,0,2],.1,500):\
plot(makelist(pump,1,2),t=0..50,theta=-3..3);\


   ** Maple V Graphics **

# 
# You could compare the two  graphs using
> plot(pump, nopump);
# if you like. 
--------------------------------------------------------------------------------
>