%% Constants
g = 1;
R = 1;

epsi_0 = 1*pi/180.;             % epsi_0 is in degrees--convert to radians

theta_0 = acos(g/(2*R));           % convert theta_0 in degrees->NO

%% define initial conditions:    theta = y(1) = theta_0 + epsi_0 
%                               theta'= y(2)=0

theta_ini=theta_0+epsi_0;

ICs = [theta_ini; 0]; % Initial Conditions

%% Perform ODE solve
[t, y] = ode45(@oscillator_eq, [0 30], ICs);

%% Plot solution
%subplot(1,2,1);
hold off;
plot(t,y(:,1),'b');

title('Exact solution of oscillator equation \epsilon_{0}= 1','FontSize',14);
xlabel('time t','FontSize',14);
ylabel('Exact solution \Theta (t)','FontSize',14);
%legend('\Theta (t)') 
hold on;

%% Calculate approximate solution

%subplot(1,2,2);
omega_prime = sqrt(2-(g/(2*R)));  % This is correct numerically, but it should really be sqrt(omega^2-(g/(omega*R))^2)
t1 = 0:0.1:30;
delta = 0 %t1.*pi;
approx_sol=theta_0+epsi_0*cos(omega_prime.*t1-delta);

%% Plot approximate solution
plot (t1,approx_sol, 'r');
title('Approximate solution of oscillator equation','FontSize',14);
xlabel('time t','FontSize',14);
ylabel('Approximate solution \Theta (t)','FontSize',14);
legend('\Theta (t)','Approx. sol. \Theta (t)')  
