%%
%-----------------------------%
% Exercise 7.0.1              %
%-----------------------------%
clear all
close all
clc

t1=0; tN=4*pi; N=50;
t=linspace(t1,tN,N); dt=t(2)-t(1);
x0=1; v0=0; x_exact=@(t)( cos(t) );

x(1)=x0; v(1)=v0;
for i=1:N-1
    xm(i)=x(i) + (dt/2)*v(i);
    vm(i)=v(i) + (dt/2)*(-x(i));
    
    x(i+1)=x(i) + dt*vm(i);
    v(i+1)=v(i) + dt*(-xm(i)); 
end

figure(1)
hold on
plot(t,x_exact(t),'r-','LineWidth',2)
plot(t,x,'g-')
legend('exact solution',...
    'Runge-Kutta approximation',...
    'Location','NorthWest')
hold off

%%
% Harmonic Oscillator -- ex. 7.0.1
clear all
close all
clc

m=1; w0=1;
x_ex=@(t,x0)(x0*sin(t+pi/2));
v_ex=@(t,x0)(x0*cos(t+pi/2));

f=@(x,v,t)(-w0^2*x);
x0=1; v0=0;

tin=0; tfin=4*pi; N=500;
t=linspace(tin, tfin, N); dt=t(2)-t(1);
x(1)=x0; v(1)=v0;
for i=1:N-1
    xc=x(i)+v(i)*dt/2;
    vc=v(i)+f(x(i),v(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+f(xc,vc,t(i)+dt/2)*dt;
end

figure(1)
hold on
plot(t,x_ex(t,x0),'r')
plot(t,x,'ro')
plot(t,v_ex(t,x0),'b')
plot(t,v,'bo')
legend('x(t) exact','x(t) RK', 'v(t) exact', 'v(t) RK',...
    'Location','SouthEast')
xlabel('t [s]','FontSize',18)
hold off

figure(2)
hold on
plot(x_ex(t,x0),v_ex(t,x0),'m')
plot(x,v,'mo')
legend('exact','RK')
axis('equal')
title('phase space')
xlabel('x [m]','FontSize',18)
ylabel('v [m/s]','FontSize',18)
hold off

% Energy

T=m*v_ex(t,x0).^2/2;
U=m*w0^2*x_ex(t,x0).^2/2;
E=T+U;
TRK=m*v.^2/2;
URK=m*w0^2*x.^2/2;
ERK=TRK+URK;

figure(3)
hold on
plot(t,E,'c')
plot(t,ERK,'co')
legend('exact','RK')
xlabel('t [s]','FontSize',18)
ylabel('Energy [J]','FontSize',18)
axis([tin tfin 0 2])
hold off


%%
% Damped Harmonic Oscillator -- ex. 7.0.2
clear all
close all
clc

m=1.4; k=6.5; c=0.8;
w0=sqrt(k/m); zeta=c/(2*sqrt(m*k)); gamma=c/(2*m);
phi=acos(sqrt(1-zeta^2));

x_ex=@(t,x0)(x0*exp(-gamma*t).*...
    cos(w0*sqrt(1-zeta^2)*t-phi)/sqrt(1-zeta^2));
v_ex=@(t,x0)(-gamma*x0*exp(-gamma*t).*...
    cos(w0*sqrt(1-zeta^2)*t-phi)/sqrt(1-zeta^2)-...
    w0*sqrt(1-zeta^2)*x0*exp(-gamma*t).*...
    sin(w0*sqrt(1-zeta^2)*t-phi)/sqrt(1-zeta^2));

f=@(x,v,t)(-w0^2*x-2*zeta*w0*v);
x0=2.8; v0=0;

tin=0; tfin=18; N=500;
t=linspace(tin, tfin, N); dt=t(2)-t(1);
x(1)=x0; v(1)=v0;
for i=1:N-1
    xc=x(i)+v(i)*dt/2;
    vc=v(i)+f(x(i),v(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+f(xc,vc,t(i)+dt/2)*dt;
end

figure(1)
hold on
plot(t,x_ex(t,x0),'r')
plot(t,x0*exp(-gamma*t)/sqrt(1-zeta^2),'r--')
plot(t,-x0*exp(-gamma*t)/sqrt(1-zeta^2),'r--')
plot(t,x,'ro')
plot(t,v_ex(t,x0),'b')
%plot(t,-gamma*x0*exp(-gamma*t)/sqrt(1-zeta^2)-...
%    w0*sqrt(1-zeta^2)*x0*exp(-gamma*t)/sqrt(1-zeta^2),'b--')
%plot(t,gamma*x0*exp(-gamma*t)/sqrt(1-zeta^2)+...
%    w0*sqrt(1-zeta^2)*x0*exp(-gamma*t)/sqrt(1-zeta^2),'b--')
plot(t,v,'bo')
legend('x(t) [m] exact','exp','exp',...
    'x(t) [m] RK', 'v(t) [m/s] exact', 'v(t) [m/s] RK',...
    'Location','SouthEast')
xlabel('t [s]','FontSize',18)
hold off


figure(2)
hold on
plot(x_ex(t,x0),v_ex(t,x0),'m')
plot(x,v,'mo')
legend('exact','RK')
axis('equal')
title('phase space')
xlabel('x [m]','FontSize',18)
ylabel('v [m/s]','FontSize',18)
hold off

% Energy

T=m*v_ex(t,x0).^2/2;
U=m*w0^2*x_ex(t,x0).^2/2;
E=T+U;
TRK=m*v.^2/2;
URK=m*w0^2*x.^2/2;
ERK=TRK+URK;

figure(3)
hold on
plot(t,E,'c')
plot(t,ERK,'co')
plot(t,T,'b')
plot(t,U,'r')
legend('exact','RK','T [J]','U [J]')
xlabel('t [s]','FontSize',18)
ylabel('Energy [J]','FontSize',18)
%axis([tin tfin 0 2])
hold off

%%
% Driven Damped Harmonic Oscillator -- ex. 7.0.3
clear all
close all
clc

w0=1; xi=0.2; F0=0.1; w= 1.2;
F=@(v,x,t)(-w0^2*x-xi*v+F0*sin(w*t));

x0=0.2; v0=0.8;

N=2000;
t=linspace(0,80,N); dt=t(2)-t(1);
x(1)=x0; v(1)=v0;
for i=1:N-1
    xc=x(i)+v(i)*dt/2;
    vc=v(i)+F(v(i),x(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+F(vc,xc,t(i)+dt/2)*dt;
end

figure(1)
subplot(1,2,1)
hold on
plot(t,x,'r')
xlabel('t','FontSize',18)
ylabel('x(t) [m]','FontSize',18)
hold off
subplot(1,2,2)
hold on
plot(t,v,'b')
xlabel('t','FontSize',18)
ylabel('v(t) [m/s]','FontSize',18)
hold off

x0an=F0./sqrt( (w0^2-w.^2).^2 + xi^2*w.^2 )


%%
% Driven Damped Harmonic Oscillator -- ex. 7.0.3
clear all
close all
clc

w0=1; xi=0.2; F0=0.1; w= 1.2
F=@(v,x,t)(-w0^2*x-xi*v+F0*sin(w*t));

x0=0.2; v0=0.8;

N=1000;
t=linspace(0,80,N); dt=t(2)-t(1);
x(1)=x0; v(1)=v0;
for i=1:N-1
    xc=x(i)+v(i)*dt/2;
    vc=v(i)+F(v(i),x(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+F(vc,xc,t(i)+dt/2)*dt;
end

figure(1)
hold on
plot(t,x,'r')
plot(t,v,'b')
legend('x(t) [m]','v(t) [m/s]')
xlabel('t','FontSize',18)
hold off

figure(2)
hold on
plot(x,v,'r')
xlabel('x [m]','FontSize',18)
ylabel('v [m/s]','FontSize',18)
hold off

T=1/2*v.^2;
U=1/2*w0^2*x.^2;
E=T+U;
figure(3)
hold on
plot(t,E,'r')
plot(t,T,'g')
plot(t,U,'b')
legend('total en', 'kinetic en','potential en')
xlabel('t [s]','FontSize',18)
ylabel('En [J]','FontSize',18)

%%
% Driven Damped Harmonic Oscillator -- ex. 7.0.3
clear all
close all
clc

w0=1; xi=0.2; F0=0.1; w= 1.2
F=@(v,x,t)(-w0^2*x-xi*v+F0*sin(w*t));

x0=0.2; v0=0.8;

N=1000;
t=linspace(0,100,N); dt=t(2)-t(1);
x(1)=x0; v(1)=v0;
for i=1:N-1
    xc=x(i)+v(i)*dt/2;
    vc=v(i)+F(v(i),x(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+F(vc,xc,t(i)+dt/2)*dt;
end

figure(1)
hold on
plot(t,x,'r')
plot(t,v,'b')
legend('x(t) [m]','v(t) [m/s]')
xlabel('t','FontSize',18)
hold off

%%
% Driven Damped Harmonic Oscillator -- ex. 7.0.3
clear all
close all
clc

w0=1; xi=0.2; F0=0.1; w= 1.5
F=@(v,x,t)(-w0^2*x-xi*v+F0*sin(w*t));

x0=0.2; v0=0.8;

N=1000;
t=linspace(0,200,N); dt=t(2)-t(1);
x(1)=x0; v(1)=v0;
for i=1:N-1
    xc=x(i)+v(i)*dt/2;
    vc=v(i)+F(v(i),x(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+F(vc,xc,t(i)+dt/2)*dt;
end

figure(1)
hold on
plot(t,x,'r')
plot(t,v,'b')
legend('x(t) [m]','v(t) [m/s]')
xlabel('t','FontSize',18)
hold off

%%
% Driven Damped Harmonic Oscillator -- new exercise:
% we want to see a resonance
clear all
close all
clc

w0=1; xi=0.2; F0=0.1
F=@(v,x,t,w)(-w0^2*x-xi*v+F0*sin(w*t));

Nw=100;
w=linspace(0.2,1.5,Nw);
for n=1:Nw
    x0=0.2; v0=0.8;
    N=5000;
    t=linspace(0,200,N); dt=t(2)-t(1);
    x(1)=x0; v(1)=v0;
    for i=1:N-1
        xc=x(i)+v(i)*dt/2;
        vc=v(i)+F(v(i),x(i),t(i),w(n))*dt/2;
        x(i+1)=x(i)+vc*dt;
        v(i+1)=v(i)+F(vc,xc,t(i)+dt/2,w(n))*dt;
    end
    
    i=(round(3*N/4):1:N);
    xmax(n)=max(x(i))
end

x0an=F0./sqrt( (w0^2-w.^2).^2 + xi^2*w.^2 );

hold on
plot(w,xmax,'b')
plot(w,x0an,'xr')
legend('xmax (RK)','analitico')
xlabel('w [s^{-1}]','FontSize',18)
ylabel('x0_{max} [m]','FontSize',18)
hold off


