%% examen 22/01/2018
%%
%--------------%
% Exercise 1   %
%--------------%
clear all
close all
clc

exact=pi^2/6

Nmax=500;
N=(5:1:Nmax);

for i=1:length(N)
    n=(0:1:N(i));
    f=1./(n.^2+2*n+1);
    
    a(i)=sum(f);
end

figure()
subplot(2,1,1)
hold on
plot(N,a,'r')
plot(N,exact*ones(length(N),1),'--k')
xlabel('N','FontSize',18)
ylabel('a(N)','FontSize',18)
hold off
subplot(2,1,2)
hold on
plot(log(1./N),a,'r')
plot(log(1./N),exact*ones(length(N),1),'--k')
xlabel('log(N^{-1})','FontSize',18)
ylabel('a(N)','FontSize',18)
hold off

% compare with
a(length(N))

% finite sum for different values of x:
x=linspace(-1+eps,1-eps,200);
n=(0:1:500);

for j=1:length(x)
    fx=1./(n.^2+2*n+1-x(j)^2);
    ax(j)=sum(fx);
end

% comparison with the analytical expression
fun=@(y)((1-pi*y./tan(pi*y))./(2*y.^2));
    
figure()
hold on
plot(x,fun(x),'b-')
plot(x,ax,'rx')
legend('exact','numerical')
xlabel('x')
axis([-1,1,1.4,10])
hold off

%%
%--------------%
% Exercise 2   %
%--------------%
clear all
close all
clc

f=@(t)( sin(2*t)./(t.*(t.^2+3)) );

N=300;
s=linspace(0+eps,2*pi,N);
il=(1:1:N-1);ir=(2:1:N);
ds=s(2)-s(1);
trap=ds*(f(s(il))+f(s(ir)))/2;

% integral indefinido
g(1)=0;
g(ir)=cumsum(trap);

const=pi/6*(1 - exp(-2*sqrt(3)));

hold on
plot(s,const*ones(N,1),'k--')
plot(s,g,'ro')
xlabel('x','FontSize',18)
ylabel('g(x)','FontSize',18)
hold off
%%
%--------------%
% Exercise 3   %
%--------------%
% Driven Damped Harmonic Oscillator: resonance
clear all
close all
clc


w0=1; beta=0.5; F0=0.2
% let's consider first a fixed value of w
w=0.1;
F0=@(v,x,t)(-w0^2*x-2*beta*v+F0*sin(w*t));

% and some initial conditions
x0=1.0; v0=1.0;

N=1000;
t=linspace(0,400,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)+F0(v(i),x(i),t(i))*dt/2;
    x(i+1)=x(i)+vc*dt;
    v(i+1)=v(i)+F0(vc,xc,t(i)+dt/2)*dt;
end

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

%%
% now let's change values of the driving frequency w


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

Nw=200;
w=linspace(0.1,6.0,Nw);
for n=1:Nw
    x0=1.0; v0=1.0;
    N=5000;
    t=linspace(0,400,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 + 4*beta^2*w.^2 );

figure()
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

%%
%--------------%
% Exercise 4   %
%--------------%
clear all
close all
clc

% (1)
% effective energy potential
theta=linspace(-pi+0.1,pi-0.1,1000);

a1=0.5;
[ Ueff1, th_min1, Ueff_min1] = Ueff_fun( theta,a1 );

a2=2;
[ Ueff2, th_min2, Ueff_min2] = Ueff_fun( theta,a2 );

figure()
hold on
plot(theta,Ueff1,'r','Linewidth',2)
plot(theta,Ueff2,'b','Linewidth',2)
plot(th_min1,Ueff_min1,'xk')
plot(th_min2,Ueff_min2,'sk')
legend('U_{eff}(\theta,a=0.5)','U_{eff}(\theta,a=2)','min a=0.5','min a=2',...
    'Location','SouthEast')
xlabel('\theta')
ylabel('U_{eff}(\theta)')
axis([-pi, pi, 1, 5])
hold off

% (2)
% phase diagram obtained from energy conservation

% a=0.5;
En=(Ueff_min1+0.1:0.2:Ueff_min1+3)

figure()
hold on
title('phase diagram a=0.5')
for j=1:length(En)
    [ theta_cut1, thetap_pl1, thetap_mn1 ] = phase_d_energy( En(j), Ueff1, theta );
    a=rand(1,1); b=rand(1,1); c=rand(1,1)
    plot(theta_cut1,thetap_pl1,'o','MarkerSize',2,'Color',[a b c])
    plot(theta_cut1,thetap_mn1,'o','MarkerSize',2,'Color',[a b c])
end

%legend()
s=ylabel('$\dot{\theta}$')
set(s,'interpreter','latex')
xlabel('\theta')
hold off

% a=2;
En=(Ueff_min2+0.1:0.7:Ueff_min2+5)

figure()
hold on
title('phase diagram a=2')
for j=1:length(En)
    [ theta_cut2, thetap_pl2, thetap_mn2 ] = phase_d_energy( En(j), Ueff2, theta );
    a=rand(1,1); b=rand(1,1); c=rand(1,1)
    plot(theta_cut2,thetap_pl2,'o','MarkerSize',2,'Color',[a b c])
    plot(theta_cut2,thetap_mn2,'o','MarkerSize',2,'Color',[a b c])
end
%legend()
s=ylabel('$\dot{\theta}$')
set(s,'interpreter','latex')
xlabel('\theta')
hold off

%%
clear all
close all
clc
% (3)
% phase diagram obtained from energy conservation
% equation of motion
% carefull with theta=0!!!
a=2;
f=@(q) ( -4*a*(1-cos(q)).^2./(sin(q)).^3+sin(q) )

tin=0; tfin=10; N=300;
theta0=pi/2; v0=0;

[ thetaRK, vRK, t, dt ] = ODE2_2RK( tin, tfin, N, theta0, v0, f );

figure()
subplot(2,1,1)
hold on
plot(t,thetaRK)
xlabel('t')
ylabel('\theta(t)')
hold off
subplot(2,1,2)
hold on
plot(t,vRK)
xlabel('t')
s=ylabel('$\dot{\theta}(t)$')
set(s,'interpreter','latex')
hold off

theta=linspace(-pi,pi,100);
[ Ueff, th_min, Ueff_min] = Ueff_fun( theta,a );
En=1 + 2*a*(1-cos(theta0)).^2./(sin(theta0)).^2 + cos(theta0);
[ theta_cut, thetap_pl, thetap_mn ] = phase_d_energy( En, Ueff, theta );

figure()
hold on
plot(thetaRK,vRK)
plot(theta_cut,thetap_pl,'xr')
plot(theta_cut,thetap_mn,'xr')
legend('numerical solution','from energy conservation')
s=ylabel('$\dot{\theta}$')
set(s,'interpreter','latex')
xlabel('\theta')
hold off

