g=9.8; % m/s^2
l=2;   % length of pendulum in m
theta0=17; %theta_0 in degrees
thetap0=0; %theta'_0 in degrees/s
tmax=20; %max time in seconds
h=1e-4; %h in seconds
damping=0; %Damping constant in 1/s

degrees_to_rad=2*pi/360; %Conversion of degrees to radians

%Convert to radians
theta0=theta0*degrees_to_rad;
thetap0=thetap0*degrees_to_rad;

%Calculate number of iterations
N=ceil(tmax/h);

%Preallocate for memory/speed vectors storing times and angles to be
%calcualted
time=1:N;
theta=time;

for i=1:N
thetapp0=-g/l*sin(theta0)-damping*thetap0;
thetap1=thetap0+h*thetapp0;
theta1=theta0+h*thetap1;

time(i)=h*i;
theta(i)=theta1;
[time(i),theta(i)]; %Reports values at each iteration--comment out for speed

theta0=theta1;
thetap0=thetap1;

end


%%
figure;
plot(time,theta);
xlabel('\omega_0 t')
ylabel('\theta(t) (rad)')


%%
%For comparison with analytic method
%we want to plot to omega_0*t=25, so
g=9.8; % m/s^2
l=2;   % length of pendulum in m
omega0=sqrt(g/l);
theta0=.1*180; %theta_0 in degrees
thetap0=0; %theta'_0 in degrees/s
tmax=25/omega0; %max time in seconds
h=1e-4/omega0; %h in seconds
damping=1; %Damping constant in 1/s


degrees_to_rad=2*pi/360; %Conversion of degrees to radians

%Convert to radians
theta0=theta0*degrees_to_rad;
thetap0=thetap0*degrees_to_rad;

%Calculate number of iterations
N=ceil(tmax/h);

%Preallocate for memory/speed vectors storing times and angles to be
%calcualted
time=1:N;
theta=time;

for i=1:N
thetapp0=-g/l*sin(theta0)-damping*thetap0;
thetap1=thetap0+h*thetapp0;
theta1=theta0+h*thetap1;

time(i)=h*i;
theta(i)=theta1;
[time(i),theta(i)]; %Reports values at each iteration--comment out for speed

theta0=theta1;
thetap0=thetap1;

end



figure;
plot(omega0*time,theta);
xlabel('\omega_0 t')
ylabel('\theta(t) (rad)')
