%% class 17/11/2017
%%
% vector fields: some examples
%%
% 2D vectors
clear all 
close all
clc

N=2;
range=1;
xmin=-range; xmax=range;
ymin=-range; ymax=range;

xmesh=linspace(xmin,xmax,N);
ymesh=linspace(ymin,ymax,N);
[x,y]=meshgrid(xmesh,ymesh);

% generate a 2D fixed vector field
vx=0.9*ones(N,N)
vy=2*ones(N,N)

S=0
quiver(x,y,vx,vy,S)
axis image
%%
% example figure 1 of rotatated vectors
clear all
close all
clc

r1=[1.5; 0.5];
theta=pi/4;
A=[cos(theta) -sin(theta); sin(theta) cos(theta)];
r2=A*r1

S=0
hold on
quiver(0,0,r1(1),r1(2),S)
quiver(0,0,r1(2),r2(2),S)
xlabel('x')
ylabel('y')
axis image
hold off
print -f -depsc 'vectors_fig1.eps'


%%
% 2D random vector field
clear all 
close all
clc
N=10;
range=2;

xmin=-range; xmax=range;
ymin=-range; ymax=range;

xmesh=linspace(xmin,xmax,N);
ymesh=linspace(ymin,ymax,N);
[x,y]=meshgrid(xmesh,ymesh);

% generate a 2D random vector field
vx=rand(N,N)-0.5;
vy=rand(N,N)-0.5;

S=1
quiver(x,y,vx,vy,S)
axis image
%%
close all
clear all
clc

xm=linspace(-3,3,10);
ym=linspace(-3,3,10);
[x,y]=meshgrid(xm,ym);

x0=-1;y0=0;
r=sqrt((x-x0).^2+(y-y0).^2);
Ex=(x-x0)./(r.^3);
Ey=(y-y0)./(r.^3);

quiver(x,y,Ex,Ey)

%%
close all
clear all
clc

xm=linspace(-3,4,20);
ym=linspace(-3,3,20);
[x,y]=meshgrid(xm,ym);

x01=-1;y01=0;
r1=sqrt((x-x01).^2+(y-y01).^2);
Ex1=(x-x01)./(r1.^3);
Ey1=(y-y01)./(r1.^3);

x02=+1;y02=0;
r2=sqrt((x-x02).^2+(y-y02).^2);
Ex2=-(x-x02)./(r2.^3);
Ey2=-(y-y02)./(r2.^3);

S=0.8
quiver(x,y,Ex1+Ex2,Ey1+Ey2,S)

%%
close all
clear all
clc

xm=linspace(-3,4,10);
ym=linspace(-3,3,10);
[x,y]=meshgrid(xm,ym);

x01=-1;y01=0;
r1=sqrt((x-x01).^2+(y-y01).^2);
Ex1=(x-x01)./(r1.^3);
Ey1=(y-y01)./(r1.^3);

x02=+1;y02=0;
r2=sqrt((x-x02).^2+(y-y02).^2);
Ex2=-(x-x02)./(r2.^3);
Ey2=-(y-y02)./(r2.^3);

S=0.8
quiver(x,y,Ex1+Ex2,Ey1+Ey2,S)

%% ex 8.0.1
clear all
close all
clc

% function defining the first derivative
% N.B. i define it as an anonymous function
% in general of t AND x
f=@(x,t)(-2*(t-5/2));

% initial and final conditions
t0=1;
x0=-9/4;
tN=4; N=5

% use Euler method for 1st order ODE
[ x, t ] = Euler_1ODE( t0, x0, tN, N, f )

x_exact=@(s)(-(s-5/2).^2);
% tangent at t2
xt2=x_exact(t(2));
x_tang=@(s)( xt2 +( -2*(t(2)-5/2) )*(s-t(2)) )

tplot=linspace(t0, tN, 100);

hold on
plot(tplot,x_exact(tplot),'r')
plot(t,x)
plot(tplot,x_tang(tplot),'r--')
legend('exact solution','numerical (Euler)','tangent of exact sol at t=t(2)','Location','South')
plot(t(2),x_exact(t(2)),'rs')
axis([t0, tN, -2.5, 1.5])
hold off

%% ex 8.0.2
clear all
close all
clc

f=@(x,t)(-3*x);

% initial and final conditions
t0=0;
x0=1;
tN=1; N=5

% use Euler method for 1st order ODE
[ x, t ] = Euler_1ODE( t0, x0, tN, N, f )

x_exact=@(s)( x0*exp(-3*s) );
tplot=linspace(t0, tN, 100);

% tangent at t2
xt2=x_exact(t(2));
x_tang=@(s)( xt2 +( -3*x0*exp(-3*t(2)) )*(s-t(2)) )

hold on
plot(tplot,x_exact(tplot),'r')
plot(t,x,'xb-')
plot(tplot,x_tang(tplot),'r--')
legend('exact solution','numerical (Euler)','tangent of exact sol at t=t(2)','Location','NorthEast')
plot(t(2),x_exact(t(2)),'rs')
axis([t0, tN, 0, 1])
hold off

%% ex 8.0.3
clear all
close all
clc

f=@(x,t)( -(3*t.^2-2*t+5).*(x-1) );

% initial and final conditions
t0=0;
x0_vector=[0, 1, 2];
tN=1; N=50;

hold on
for i=1:3
    x0=x0_vector(i);
    
    % use Euler method for 1st order ODE
    [ x, t ] = Euler_1ODE( t0, x0, tN, N, f );
    
    x_exact=@(s)( 1+(x0-1)*exp(-(s.^3-t0^3)+(s.^2-t0^2)-5*(s-t0)) );
    tplot=linspace(t0, tN, 100);

    plot(tplot,x_exact(tplot),'--')
    plot(t,x,'Linewidth',2)
end
legend('exact solution x_0=0','numerical (Euler) x_0=0',...
    'exact solution x_0=1','numerical (Euler) x_0=1',...
    'exact solution x_0=2','numerical (Euler) x_0=2',...
    'Location','SouthEast')
hold off

%% ex 9.0.4
clear all
close all
clc

% frequency of the harmonic oscillator
w0=1; %s^{-1}
g1=@(x,v,t)(v);
g2=@(x,v,t)(-w0^2*x);

%initial conditions
x0=1; v0=0;

% exact solution
delta=atan(x0/v0*w0);
x_exact=@(s)( sqrt(x0^2+v0^2/(w0^2))*sin(s*w0 + delta) )
v_exact=@(s)( w0*sqrt(x0^2+v0^2/(w0^2))*cos(s*w0 + delta) )

% interval
t0=0; tN=4*pi; N=1500;

[ x, v, t ] = Euler_2ODE( x0, v0, t0, tN, N, g1, g2 )

tplot=linspace(t0, tN, 200);

figure()
subplot(2,1,1)
hold on
plot(tplot,x_exact(tplot),'r','Linewidth',2)
plot(t,x)
legend('exact solution (w_0=1 s^{-1})', 'numerical (Euler)')
xlabel('t [s]')
ylabel('x(t) [m]')
hold off

subplot(2,1,2)
hold on
plot(tplot,v_exact(tplot),'r','Linewidth',2)
plot(t,v)
legend('exact solution (w_0=1 s^{-1})', 'numerical (Euler)')
xlabel('t [s]')
ylabel('v(t) [m/s]')
hold off

figure()
hold on
%plot(tplot,v_exact(tplot),'r','Linewidth',2)
plot(x,v)
% how to put the phase line directions? thanks Diego!
quiver(x(1:100:N),v(1:100:N),v(1:100:N),-w0^2*x(1:100:N))
%legend('exact solution (w_0=1 s^{-1})', 'numerical (Euler)')
xlabel('x [m]')
ylabel('v(t) [m/s]')
hold off