% File:  pblm3_1.m  
%
% Computer Problem 3.1
%

% Part a.

%  ODE tolerance
tol=0.0001;
% Initial conditions
x0=[pi/2;0];
% Time interval to study
t0=0;tf=4;

% Solve with no approximations
[t,x] = ode23('func3_1a',t0,tf,x0,tol,0);
% Plot results
subplot(211);
plot(t,x(:,1));
hold on;
title('theta(t) for theta_initial=pi/2');
ylabel('theta in radians');
xlabel('time in sec');
subplot(212);
plot(t,x(:,2));
title('d/dt of theta(t) for d/dt of theta_initial=0');
ylabel('d/dt of theta in radians/sec');
xlabel('time in sec');
hold on;  % hold plot so can overplot sin(theta)=theta results

% solve with sin(theta)=theta approximation
[t,x] = ode23('func3_1a2',t0,tf,x0,tol,0);
% Plot results
subplot(211);
plot(t,x(:,1),'--');
subplot(212);
plot(t,x(:,2),'--');
legend('-','No approx.','--','sin(x)=x approx.');

% Part b.

% Now repeat the above calculation (no sin(x)=x approx.) for
%  theta_0=3.1 rad.

x0=[3.10;0];
[t,x] = ode23('func3_1a',t0,tf,x0,tol,0);
subplot(211);
plot(t,x(:,1),'-.');
subplot(212);
plot(t,x(:,2),'-.');
legend('-','No approx.','--','sin(x)=x approx.','-.','theta_0=3.1 rad');

% Part c.

figure(2);

% Calculate sqrt(g/L), the period for sin(x)=x approx.
T_perfect=2*pi/sqrt(9.8)

% Now calculate the period for various initial theta_0's from 0+ to 3.1 rad
tol=0.0001;

% First numerically calculate the period for sin(x)=x 
[t,x] = ode23('func3_1a2',t0,tf,x0,tol,0);  % func3_1a2 is for sin(x)=x
% Linearly interpolate about first zero crossing to get the period
ind=find(x(:,1) < 0);
y1=x(ind(1)-1,1);
y2=x(ind(1),1);
t1=t(ind(1)-1);
t2=t(ind(1));
a=(y2-y1)/(t2-t1);
T_perfect_num=4*(t1-y1/a)

% Now calculate the period for a range of theta_0's
% Note the step size of 0.1 is just a "guess"
k=0;
for i=0.1:0.1:3.1
   k=k+1;
   x0=[i,0];
   [t,x] = ode23('func3_1a',t0,tf,x0,tol,0);
   ind=find(x(:,1) < 0);
   y1=x(ind(1)-1,1);
   y2=x(ind(1),1);
   t1=t(ind(1)-1);
   t2=t(ind(1));
   a=(y2-y1)/(t2-t1);
   period(k)=4*(t1-y1/a);
   amplitude(k)=i;
end;

% Plot the period vs theta_0
plot(amplitude,period,'+');
title('Period of Oscillation vs Theta Initial');
ylabel('Period (sec)');
xlabel('Theta Initial Amplitude (rad)');
grid;
hold;
plot(amplitude,period);
% Locate the period for sin(x)=x
plot([t(1) t(length(t))],[T_perfect T_perfect],'r');
% Locate the period for sin(x)=x + 2%
plot([t(1) t(length(t))],[T_perfect*1.02 T_perfect*1.02],'g');
% We need to blowup the plot to see where the period is 2% larger than
%  sqrt(g/L)
axis([0 1 1.9 2.2]);
% Note the crossover is around 0.56
text(0.1,2.125,'Period increases by 2% near theta_0=0.56 radians');


