%pendulum, not SHO
%the de's are
% theta'' + w02*sin(theta) = 0
% z' + w02*sin(theta) = 0
% theta' = z
% need initial conditions of theta0 = pi/2
% and theta'(0) = z(0) = 0

clear all

%% 
theta0 = pi;
z0 = 1;
w02 = 10;

h = 0.000001;
t = 0:h:8;

thetaNumerical = zeros(size(t));
zNumerical = zeros(size(t));
zPrime = zeros(size(t));
thetaPrime = zPrime;

%% this is the first time, initial conditions

thetaNumerical(1) = theta0;
zNumerical(1) = z0;
thetaPrime(1) = zNumerical(1);
zPrime(1) = -w02*sin(thetaNumerical(1));

%%

for i=2:(length(t))
    thetaNumerical(i) = thetaNumerical(i-1) + thetaPrime(i-1)*h;
    zNumerical(i) = zNumerical(i-1) + zPrime(i-1)*h;
    thetaPrime(i) = zNumerical(i);
    zPrime(i) = -w02*sin(thetaNumerical(i));
end;

thetaScaled = thetaNumerical./theta0;

plot(t,thetaScaled)



%%
hold;
simpleHarmonic = (pi/4)*cos(sqrt(w02)*t);
plot(t, simpleHarmonic);

%%

% %% this is time 2
% 
% thetaNumerical(2) = thetaNumerical(1) + thetaPrime(1)*h;
% 
% zNumerical(2) = zNumerical(1) + zPrime(1)*h;
% thetaPrime(2) = zNumerical(2);
% zPrime(2) = -w02*sin(thetaNumerical(2));
% %where thetaPrime(1) = z(1)
% 
% %% this is time 3
% 
% thetaNumerical(3) = thetaNumerical(2) + thetaPrime(2)*h;
% 
% zNumerical(3) = zNumerical(2) + zPrime(2)*h;
% thetaPrime(3) = zNumerical(3);
% zPrime(3) = -w02*sin(thetaNumerical(3));
% 
% %% this is time 4
% 
% thetaNumerical(4) = thetaNumerical(3) + thetaPrime(3)*h;
% 
% zNumerical(4) = zNumerical(3) + zPrime(3)*h;
% thetaPrime(4) = zNumerical(4);
% zPrime(4) = -w02*sin(thetaNumerical(4));
