Lab 10: Variable length pendulae
This lab will follow up from class.

In class we wrote the procedure DPenddT, and visualizePend 
function dPdT = dPenddT(t,p)
u = p(1);
v = p(2);

g = -9.8;
L = 1;
uPrime = v;
vPrime = (g/L) * sin(u);  % no friction
dPdT = [uPrime vPrime]';
................ 
function visualizePend(p)
%loop over each time point
for ix = 1:size(p,1);
    %plot the top of the pendulum
    plot(0,0,'+');hold on;
    %plot the bottom of the pendulum
    plot(sin(p(ix,1)),-cos(p(ix,1)),'o');hold off;
    % force it to draw
    axis([-1 1 -1 1]);
    drawnow;
end
Now, suppose that the pendulum is attached to a helicopter, and the weight at the bottom is slowly being winched up at a rate alpha. That is, the length of the pendulum L decreases by "alpha" during each time unit. A rate of 0.1 works well, so you can use that as a constant in the below equations and calculate L as 1-alpha*T. 

(old version)   y'' = -(g / L) * sin(y), we have:

(new version)   y'' = -(g / L) * sin(y) - 2*(L'/L) * y';
Where the new term is the "centripetal acceleration" of pulling something into the center of a circle.

Task 1 Adjust both the dPenddT and the visualizePend functions to be able to simulate and visualize the behavior of the pendulum as it becomes shorter. Look up the "xlabel" and "ylabel" commands and make a plot of the phase space with the axes labelled. When you're done, run and visualize the command: 

[t p] = ode23('dPenddT',[0:0.01:8],[pi/2 -1]);
Your phase space plot should look something like the one pictured below. I've drawn an extra 'o' at the end of the plot (to tell it from the beginning.

Show the TA both your visualization of the pendulum swinging and your phase space plot for your first check-off.

Task 2 Instead of including alpha as a "hard coded" value that we've set to 1, we can include L as another parameter in our differential equation. This requires that you make dPenddT(t,p) take a vector p of 3 arguments, where the 3rd argument is the "length" of the pendulum. Then, dPenddT must return a 3 element vector, dPdt, that you might set as:

dPdT = [uPrime vPrime alphaPrime]';
Where the alphaPrime is set the to the constant -0.1 to recreate your previous simulation. When you call the ode23 solver on this, you starting condition include the length of the pendulum (the "1" at the end of the ode23 command below). 
[t p] = ode23('dPenddT',[0:0.01:10],[pi/2 -1 1]);

This will give you a result "p" which is is an array with the shape:

"number of time steps" x 3
where the first element is the theta angle, the second is the angular velocity and the third is the length of the pendulum. Now, fix your visualization code to use the third column of "p" as the length of the pendulum and see that you get the same visualization as before. Unless you are have trouble, you need not have the TA check this off

Task 3 Now, we are going to make a schedule of how to winch up the pendulum, without it going faster and faster. The goal is to winch the pendulum all the way up to be zero length, and to never let the "theta" angle exceed pi/6.

To do this we are going to allow the "alpha" term in you dPenddT function to vary. Instead of setting alphaPrime = -0.1, you must figure out "when" to winch the pendulum in. (Imagine holding a pendulum by a string. When does shortening the string make it go faster? When does it make it go slower?).

This might include terms like: 

if abs(u)>0.1,    % if the winch isn't at the bottom of its swing, pull up
  alphaPrime = -0.1
else              % otherwise, let it down (which may slow down the motion)
  alphaPrime = 0.1
end

--- or --- 

if abs(u)<0.1,    % if the winch is going slowly, pull up
  alphaPrime = -0.1
else              % otherwise, keep string the same length
  alphaPrime = 0.0
end

--- or --- 

Many other possibilities
At the end, you need to show the TA:
  1. A visualization of a run with:
    [t p] = ode23('dPenddT',[0:0.01:10],[pi/6 0 1]);
    which sets an initial string angle of pi/6 and length of 1, *and* where (as the equation continues) the pendulum is eventually winched to no more that 10% its original length. You may let the similution go for more that 10 seconds.
  2. You dPenddT code, with an explanation of "when you are winching".