% Compare Euler and Euler-Cromer

clear
g=9.8;
L=0.5;
m=0.005;
per=2*pi*sqrt(L/g);

a=0;
%b=5*per;
b=8;
N=500;
h=(b-a)/N;

% Initial Condition
t(1)=0;
theta(1)=pi/6; % Value used for book
theta(1) = 10*pi/180;
omega(1)=0;

% Euler's Method
for i=2:N+1
    omega(i)=omega(i-1)-g/L*h*sin(theta(i-1));
    theta(i)=theta(i-1)+h*omega(i-1);
    t(i)=t(i-1)+h;
end

figure(1)
plot(t,theta,'k-')
xlabel('t (s)')
ylabel('\theta (\pi rad)')
set(gca,'YTick',-3*pi/4:pi/4:3*pi/4)
set(gca,'YTickLabel',{'-3/4','-1/2','-1/4','0','1/4','1/2','3/4'})
hold on

figure(2)
E1 = 1/2*m*L^2*omega.^2+m*g*L*(1-cos(theta)); 
plot(t,E1,'k-')
xlabel('t (s)')
ylabel('E (J)')
hold on

% Initial Condition
t2(1)=0;
%theta2(1)=pi/6;
theta2(1)=theta(1);
omega2(1)=0;

% Euler-Cromer Method
for i=2:N+1
    omega2(i)=omega2(i-1)-g/L*h*sin(theta2(i-1));
    theta2(i)=theta2(i-1)+h*omega2(i);
    t2(i)=t2(i-1)+h;
end

figure(1)
plot(t2,theta2,'k--')
legend('Euler Method','Euler-Cromer','location','southwest')
axis([0,8,-3*pi/4,3*pi/4])
hold off

figure(2)
E2 = 1/2*m*L^2*omega2.^2+m*g*L*(1-cos(theta2)); 
plot(t2,E2,'k--')
xlabel('t (s)')
ylabel('E (J)')
axis([0,8,0,0.03])
legend('Euler Method','Euler-Cromer','location','northwest')
hold off

% ---- Added plots ----  

figure(3)
plot(t,theta,'k-')
xlabel('t (s)')
ylabel('\theta (\pi rad)')
set(gca,'YTick',-3*pi/4:pi/4:3*pi/4)
set(gca,'YTickLabel',{'-3/4','-1/2','-1/4','0','1/4','1/2','3/4'})
axis([0,8,-3*pi/4,3*pi/4])

figure(4)
plot(t,theta,'k-')
xlabel('t (s)')
ylabel('\theta (\pi rad)')
set(gca,'YTick',-3*pi/4:pi/4:3*pi/4)
set(gca,'YTickLabel',{'-3/4','-1/2','-1/4','0','1/4','1/2','3/4'})
axis([0,8,-3*pi/4,3*pi/4])
hold on

N=1000;
h=(b-a)/N;

% Initial Condition
t3(1)=0;
theta3(1)=theta(1);
omega3(1)=0;

% Euler's Method
for i=2:N+1
    omega3(i)=omega3(i-1)-g/L*h*sin(theta3(i-1));
    theta3(i)=theta3(i-1)+h*omega3(i-1);
    t3(i)=t3(i-1)+h;
end

plot(t3,theta3,'k--')

N=2000;
h=(b-a)/N;

% Initial Condition
t4(1)=0;
theta4(1)=theta(1);
omega4(1)=0;

% Euler's Method
for i=2:N+1
    omega4(i)=omega4(i-1)-g/L*h*sin(theta4(i-1));
    theta4(i)=theta4(i-1)+h*omega4(i-1);
    t4(i)=t4(i-1)+h;
end

plot(t4,theta4,'k-.')
legend('N = 500','N = 1000','N = 2000','location','southwest')
hold off

figure(5)
E3 = 1/2*m*L^2*omega3.^2+m*g*L*(1-cos(theta3)); 
E4 = 1/2*m*L^2*omega4.^2+m*g*L*(1-cos(theta4)); 

plot(t,E1,'k-',t3,E3,'k--',t4,E4,'k-.',t2,E2,'k:')
xlabel('t (s)')
ylabel('E (J)')
legend('Euler N = 500','Euler N = 1000','Euler N = 2000','Euler-Cromer N = 500','location','northwest')
axis([0,8,0,0.03])

% For solutions Manual
figure(6)
plot(t2,180*theta2/pi,'b-','linewidth',2)
xlabel('t (s)')
ylabel('\theta (degrees)')
%set(gca,'YTick',-90:10:90)
%set(gca,'YTickLabel',{'-270','-180','-90','0','90','180','270'})
hold on
plot(t2,180*theta2(1)/pi*cos(2*pi*t2/per),'r--','linewidth',2)
legend('Euler-Cromer Method','Linear Pendulum','location','southwest')
axis([0,8,-200*theta(1)/pi,200*theta(1)/pi])
grid on
grid(gca,'minor')
hold off