%% class 01/12/2017
%%
% ex 15.0.5
clear all
close all
clc


f=@(x,t)(x.^2);

x_ex=@(s,x1)(x1./(1-x1*s));

% let's compare exact and numerical
% initial conditions
t0=0;
x0=0.2;
tN=7; N=100;

[ x_2RK, t ] = RK2ndO_1ODE( t0, x0, tN, N, f )

figure()
hold on
plot(t,x_ex(t,x0),'r--')
plot(t,x_2RK,'r')

x0=-2;
[ x_2RK, t ] = RK2ndO_1ODE( t0, x0, tN, N, f );

plot(t,x_ex(t,x0),'b--')
plot(t,x_2RK,'b')

legend('exact solution (x_0=0.2)', 'RK (x_0=0.2)',...
    'exact solution (x_0=-2)', 'RK (x_0=-2)')
plot(t,0*ones(N,1),'k--')
xlabel('t')
ylabel('x(t)')
axis([t0, tN, -2, 2])
hold off

% more than 1 initial condition
t0=0;
x0=[-1, -0.5, 0.5, 1];
tN=5; N=1000;

% for one initial condition we compare with exact solution
figure()
hold on
for i=1: length(x0)
    % 2nd order Runge-Kutta
    [ x_2RK, t ] = RK2ndO_1ODE( t0, x0(i), tN, N, f );
    plot(t,x_2RK)
end
legend('RK x1=-1','RK x1=-0.5','RK x1=0.5','RK x1=1')
plot(t, 0*ones(N,1),'k--')
xlabel('t')
ylabel('x(t)')
axis([t0, tN, -2, 2])
hold off

%% planar pendulum
clear all
close all
clc

g1=@(th, w, t)(w);
g2=@(th, w, t)(-sin(th));

% initial conditions
th0=pi/2;
w0=0.2;

% first we need to ckeck convergence w.r.t. the time step dt or N
figure()
hold on
title('study of the convergence w.r.t. N')

t0=0; tN=12*pi; N=200;
[ th, w, t ] = RK2ndO_2ODE( t0, tN, N, th0, w0, g1, g2 );
plot(t,th)

t0=0; tN=12*pi; N=300;
[ th, w, t ] = RK2ndO_2ODE( t0, tN, N, th0, w0, g1, g2 );
plot(t,th)

t0=0; tN=12*pi; N=500;
[ th, w, t ] = RK2ndO_2ODE( t0, tN, N, th0, w0, g1, g2 );
plot(t,th)

xlabel('t')
ylabel('\theta (t)')
hold off


% initial conditions
th0=pi/2;
w0=0.2;

t0=0; tN=12*pi; N=500;
[ th, w, t ] = RK2ndO_2ODE( t0, tN, N, th0, w0, g1, g2 );

figure()
subplot(3,1,1)
hold on
title('angle \theta(t)')
plot(t,th)
xlabel('t')
ylabel('\theta(t)')
hold off

subplot(3,1,2)
hold on
title('angular velocity \omega(t)')
plot(t,w)
xlabel('t')
ylabel('\omega(t)')
hold off

subplot(3,1,3)
hold on
title('phase diagram')
plot(th,w)
xlabel('\theta')
ylabel('\omega')
hold off

% study of the phase diagram:
% from energy conservation we know that the phase diagram
% can be derived from
% E = const = w^2/2 + U(th)
% where U(th)=-cos(th)

U=@(theta)(-cos(theta));

E=w0^2/2+U(th0)

w_plus=@(theta,e)( +sqrt(2*(e+cos(theta))) );
w_minus=@(theta,e)( -sqrt(2*(e+cos(theta))) );

phi=linspace(-acos(-E),acos(-E),500);

figure()
hold on
title('phase diagram: comparison with exact solution')
plot(th(1:10:N),w(1:10:N),'xr')
plot(phi,w_plus(phi,E),'b')
legend('numerical','exact')
plot(phi,w_minus(phi,E),'b')
xlabel('\theta')
ylabel('\omega')
axis('equal')
hold off

%%
% and now let's deduce the entire phase diagram
% by changing initial conditions
clear all
close all
clc

g1=@(th, w, t)(w);
g2=@(th, w, t)(-sin(th));

U=@(theta)(-cos(theta));
w_plus=@(theta,e)( +sqrt(2*(e+cos(theta))) );
w_minus=@(theta,e)( -sqrt(2*(e+cos(theta))) );

t0=-4*pi; tN=4*pi; N=250;

th0=0; w0=(-3:0.5:3);

figure()
hold on
title('entire phase diagram: comparison with exact solution')
for j=1:length(w0)
    [ th, w, t ] = RK2ndO_2ODE( t0, tN, N, th0, w0(j), g1, g2 );

    % phase diagram
    E=w0(j)^2/2+U(th0);
    
    if E<1
        phi=linspace(-acos(-E),acos(-E),500);
    else
        phi=linspace(-pi,pi,500);
    end
    a=rand(1,1); b=rand(1,1); c=rand(1,1)
    
    plot(th(1:3:N),w(1:3:N),'o','Color',[a b c])
    plot(phi,w_plus(phi,E),'Color',[a b c])
    plot(phi,w_minus(phi,E),'Color',[a b c])
    
end
axis([-pi,pi,w0(1),w0(length(w0))])
xlabel('\theta')
b=ylabel('$\dot{\theta}=\omega$')
set(b,'interpreter','latex')
hold off

% HOMEWORK: 
% 1. SHOW THAT FOR THETA <<1 YOU RECOVER THE HARMONIC OSCILLATOR SOLUTION
% 2. EVALUATE THE PERIOD OF OSCILLATIONS T AS A FUNCTION OF THE AMPLITUDE
% OF OSCILLATIONS DELTA THETA AND SHOW THAT THE PLANAR PENDULUM IS
% ISOCHRONOUS AND GIVEN BY 2PI ONLY FOR DELTA THETA SMALL

%% ex. 13.0.2 = damped harmonic oscillator
clear all
close all
clc

% parameters
m=1.4; % kg
k=6.5; % N/m
c=0.8; % kg/s

% get the intrinsic frequency of the harmonic oscillator
omega0=sqrt(k/m);
% friction frequency
gamma=c/(2*m);
% dimensionless friction parameter
zeta=c/(2*sqrt(m*k));
% because zeta<1 the oscillations are underdamped

g1=@(x, v, t)(v);
g2=@(x, v, t)(-omega0^2*x-2*gamma*v);

% initial condition
x0=2.8; % m
v0=0;

t0=0; tN=18; N=1000;
[x, v, t ] = RK2ndO_2ODE( t0, tN, N, x0, v0, g1, g2 );

% phase shift
phi=acos(sqrt(1-zeta^2));
x_ex=@(s)( x0/sqrt(1-zeta^2) * exp(-gamma*s) .* cos(omega0*s*sqrt(1-zeta^2) -phi) )

figure
hold on
title('underdamped harmonic oscillator')
plot(t,x_ex(t),'r')
plot(t(1:10:N),x(1:10:N),'ob')
legend('exact','numerical')
plot(t,x0/sqrt(1-zeta^2) * exp(-gamma*t),'r--')
plot(t,-x0/sqrt(1-zeta^2) * exp(-gamma*t),'r--')
xlabel('t [s^{-1}]')
ylabel('x(t) [m]')
hold off

