function project_anglefunction
close all
clc

% using ODE45 to solve for theta
ops=odeset('RelTol',1e-10);
a=.1; b=.1; theta_0=pi/2; thetap_0=0; g=9.8; tinit=0; tfinal=50;
increment=500;
test=linspace(tinit,tfinal,increment);
[t,u]=ode45(@anim,test,[theta_0,thetap_0],ops,a,b);

% the animation of the pendulum
parabstart=-.5;
parabstop=.5;
tparab=linspace(parabstart,parabstop);

tback=[];
uback=[];
figure

for time = 1:increment
    plot((a+b*t).*sin(u(:,1)),(a+b*t).*cos(u(:,1)),'k')
    hold on
    tback=[tback,t(time)];
    uback=[uback,u(time,1)];
    plot((a+b*tback).*sin(uback),(a+b*tback).*cos(uback),'g')
    plot(-tparab+(a+b*t(time)).*sin(u(time,1)),-tparab.^2+(a+b*t(time)).*cos(u(time,1)),'r')
    plot([0,(a+b*t(time)).*sin(u(time,1))],[0,(a+b*t(time)).*cos(u(time,1))])
    set(gca,'ydir','reverse')
    hold off
    axis([-sin(max(u(:,1)))+2*parabstart,sin(max(u(:,1)))+2*parabstop,0,(a+b*tfinal)])
    M(time) = getframe;
    
end

% play movie
%movie(M,0)


function uprime=anim(t,u,a,b)
uprime=zeros(2,1);
uprime(1)=u(2);
uprime(2)=-2*b*u(2)./(a+b.*t)-9.8*sin(u(1))./(a+b.*t);