%%
% class of 3/11/2017
%% ex 3 of week 3 homework
clear all
clc
close all

% initial conditions
x0=2.8; %m


t0=0; t1=18; %s
N=10000;
t=linspace(t0, t1, N);

% parameters
m=1.0 % kg
k=1.0 %N/m
c=0.1 %kg/s

[ x, gamma, omega0, sq ] = dho_position( m, c, k, x0, t );

figure()
subplot(2,1,2)
title('m=1, k=1, c=0.1')
hold on
xlabel('t')
ylabel('x(t)')
plot(t,x)
plot(t,x0/sq * exp(-gamma*t),'r--')
plot(t,-x0/sq * exp(-gamma*t),'r--')
hold off

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

[ x, gamma, omega0, sq ] = dho_position( m, c, k, x0, t );

subplot(2,1,1)
title('m=1.4, k=6.5, c=0.8')
hold on
xlabel('t')
ylabel('x(t)')
plot(t,x)
plot(t,x0/sq * exp(-gamma*t),'r--')
plot(t,-x0/sq * exp(-gamma*t),'r--')
hold off

% velocity = first derivative of position
[ v, il, ir, dt ] = first_derivative( t, x );
% N.B.!!!
% even if in the function i  use different names as in the following
%[ derivative_fun, il, ir, dvar ] = first_derivative( var, fun )
% i can call the function as above giving different names BUT not 
% different positions!!

% acceleration = first derivative of velocity
[ a, ila, ira, dt ] = first_derivative( t, v );

figure()
subplot(3,1,1)
title('m=1.4, k=6.5, c=0.8')
hold on
xlabel('t')
ylabel('x(t)')
plot(t,x)
plot(t,x0/sq * exp(-gamma*t),'r--')
plot(t,-x0/sq * exp(-gamma*t),'r--')
hold off

subplot(3,1,2)
hold on
xlabel('t')
ylabel('v(t)')
plot(t(il)+dt/2,v)
plot(t,0*ones(N,1),'k--')
hold off

subplot(3,1,3)
hold on
xlabel('t')
ylabel('a(t)')
plot(t(ila)+dt/2,a)
plot(t,0*ones(N,1),'k--')
hold off

% point 3 check numerically that the diff eq is satisfied
% N.B. x, v and a have a different number of points!!!
% let's use a different differentiation routine

[ v2, ic, dvar ] = first_derivative_vs2( t, x );

figure()
subplot(3,1,1)
title('m=1.4, k=6.5, c=0.8')
hold on
xlabel('t')
ylabel('x(t)')
plot(t,x)
plot(t,x0/sq * exp(-gamma*t),'r--')
plot(t,-x0/sq * exp(-gamma*t),'r--')
hold off

% check both routines give the same result!!
subplot(3,1,2)
hold on
xlabel('t')
ylabel('v(t)')
plot(t(il)+dt/2,v)
plot(t(ic),v2(ic),'or')
legend('derivative vs1', 'derivative vs2')
plot(t,0*ones(N,1),'k--')
hold off

% also for acceleration
[ a2, ic_a, dvar ] = first_derivative_vs2( t, v );

subplot(3,1,3)
hold on
xlabel('t')
ylabel('a(t)')
plot(t(ila)+dt/2,a)
plot(t(ic_a),a2(ic_a),'or')
legend('derivative vs1', 'derivative vs2')
plot(t,0*ones(N,1),'k--')
hold off

% so now we can use the reduced grid ic_a for checking that
% a(t) + 2*gamma* v(t) + omega0^2*x(t) = 0
% is numerically (approximatively) satisfied
figure()
hold on
plot(t(ic_a),a2(ic_a)+2*gamma*v2(ic_a)+omega0^2*x(ic_a))
axis([t0, t1, -1, 1])
hold off

%% ex 3 of week 4 homework
% finding zeros of a function
% or minimal distance between two vectors
clear all
close all
clc

x0=2; y0=-3; R=5;
N=500;
x=linspace(x0-R,x0+R,N);

c1=y0+sqrt(-(x-x0).^2+R^2);
c2=y0-sqrt(-(x-x0).^2+R^2);

f=exp(-x/2)-4;

% option 1
[min_val1, imin1]=min(abs(c1-f))
[min_val2, imin2]=min(abs(c2-f))

% option 2
g1=@(x) (y0+sqrt(-(x-x0).^2+R^2) - (exp(-x/2)-4))
g2=@(x) (y0-sqrt(-(x-x0).^2+R^2) - (exp(-x/2)-4))

x1min=fzero(g1, [x0-R,0])
x2min=fzero(g2, [4, x0+R])

figure()
hold on
plot(x,c1)
plot(x,c2)
plot(x,f)
plot(x(imin1),f(imin1),'sr')
plot(x(imin2),f(imin2),'sg')
plot(x1min,y0+sqrt(-(x1min-x0).^2+R^2),'*k')
plot(x2min,y0-sqrt(-(x2min-x0).^2+R^2),'*k')
axis('equal')
hold off

figure()
hold on
plot(x,abs(c1-f),'r')
plot(x,abs(c2-f),'b')
plot(x(imin1),abs(c1(imin1)-f(imin1)),'sr')
plot(x(imin2),abs(c2(imin2)-f(imin2)),'sg')
plot(x1min,abs(g1(x1min)),'*k')
plot(x2min,abs(g2(x2min)),'*k')
hold off

figure()
hold on
plot(x,c1-f,'r')
plot(x,c2-f,'b')
plot(x(imin1),c1(imin1)-f(imin1),'sr')
plot(x(imin2),c2(imin2)-f(imin2),'sg')
plot(x1min,g1(x1min),'*k')
plot(x2min,g2(x2min),'*k')
plot(x,0*ones(N,1),'k--')
hold off

%% ex 3 of week 4 homework
% option 3 with double loop 
% not the best and most efficient way!
clear all
close all
clc


x0=2; y0=-3; R=5;
N=1000;

phi=linspace(0,2*pi,N);
x1=x0+R*cos(phi);
y1=y0+R*sin(phi);

x2=linspace(-4,8,N);
y2=exp(-x2/2)-4;

for i=1:N
    for j=1:N
        d(i,j)=sqrt( (x1(i)-x2(j))^2 + (y1(i)-y2(j))^2 );
    end
end

dvect=d(:);
[min_value, index_min]=min(dvect)

dvect(index_min)
d(index_min)
% which value of i and j

for i=1:N
    for j=1:N
        if (index_min==(j-1)*N+i)
            i_min=i
            j_min=j
        end
    end
end

% in this way we find one only crossing point! find the next one
dvect_red=[dvect(1:index_min-1); dvect(index_min+1:N*N)];
[min_value2, index_min2]=min(dvect_red)

% which value of i and j
for i=1:N
    for j=1:N
        if (index_min2==(j-1)*N+i)
            i_min2=i
            j_min2=j
        end
    end
end

figure()
hold on
plot(x1,y1)
plot(x2,y2)
plot(x1(i_min),y1(i_min),'rs')
plot(x1(i_min2),y1(i_min2),'rs')
legend('circle=(x1,y1)','f(x)=(x2,y2)')
xlabel('x')
ylabel('y')
axis('equal')
hold off
