%-----------------------------%
% Exercises unit 2 --- part 2 %
%-----------------------------%

%%
%-----------------------------%
% Ex 1.0.1
clear all
close all
clc

% the half-wave rectifier
x=linspace(-2*pi,2*pi,100);
f=sin(x);
g=f
for i=1:length(x)
    if f(i)<0
        g(i)=0;
    end
end

hold on
title('the half-wave rectifier')
plot(x,f,'g--',x,g,'ro-')
xlabel('x')
%ylabel('f,g')
legend('sin(x)','rectified sin(x)')
hold off
%%
% example use of meshgrid
clear all
close all
clc
%fix number of points on the grid

% generate x,y mesh
xmesh=linspace(-1,1,3)
ymesh=linspace(1,2,5)
[x,y]=meshgrid(xmesh,ymesh)

r=sqrt(x.*x+y.*y);
f=r.^2;
mesh(x,y,f)

%%
%-----------------------------%
% Ex 1.0.2
clear all
close all
clc
%fix number of points on the grid
N=100;

% generate x,y mesh
xmesh=linspace(-1,1,N);
ymesh=linspace(-1,1,N);
[x,y]=meshgrid(xmesh,ymesh);

r=sqrt(x.*x+y.*y);
f=r.^2-r.^4/2;

%2D plot of f
%contour(x,y,f)
%surf(x,y,f)
%surfl(x,y,f)
mesh(x,y,f)
colorbar

%Find the GLOBAL maximum value of f

%compare with the most negative number of the computer
fmax=-Inf;
imax=1; jmax=1;
for i=1:N
    for j=1:N
        if f(i,j)>fmax
           fmax=f(i,j);
           imax=i;
           jmax=j;
        end
    end
end
fmax

%compare with
max(max(f))

%the last position (in the loop) (i,j) 
%for which f(i,j) is equal to fmax is give by:
[imax,jmax]


%%
%-----------------------------%
% Ex 2.0.3
clear all
close all
clc
% type the integer number you want to know
% if it is a prime number or not
n=input('type an integer number');
div=2;
% define a a variable 'prime' which can either be 1 or 0
% if 'prime' is 1 the number is prime otherwise is not
prime=1;
while div<=n-1
    if mod(n,div)==0
        %if the rest of the division is zero n cannot be a prime number
         prime=0;
         %terminate the execution of while
         break
    end
    div=div+1;
end
if prime==1
      disp('n is a prime number')
   else
      disp('n is not a prime number')
      disp('and the smaller number it can be divided by is')
      div
end
%%
%-----------------------------%
% Ex 2.0.3 with for instead than while
clear all

% type the integer number you want to know
% if it is a prime number or not
n=input('type an integer number');

% define a a variable 'prime' which can either be 1 or 0
% if 'prime' is 1 the number is prime otherwise is not
prime=1;

for div=2:n-1
    if mod(n,div)==0
    %if the rest of the division is zero n cannot be a prime number
    prime=0;
    %terminate the execution of while
    break
    end
    div=div+1;
end

if prime
      disp('n is a prime number')
   else
      disp('n is not a prime number')
      disp('and the smaller number it can be divided by is')
      div
end
%%
%-----------------------------%
% Ex 2.0.3 solved in class by Eduardo in a very nice and neat way

clear all 

f=input('El siguiente numero');

n=2;
while mod(f,n)>0
    n=n+1;
end
if n==f
    disp('Es un numero primo')
else
    disp('No es un numero primo')
    disp('y el numero mas pequeno por el cual se puede dividir es')
    n
end

%%
%------------------------------------------%
% Additional excercises of Unit 2 - part 2 %
%------------------------------------------%
% Ex 3.0.4
clear all
close all
clc

% define a column vector
r1=[1.5; 0.5];

% angle
theta=pi/4;

% rotation matrix
A=[cos(theta) -sin(theta); sin(theta) cos(theta)];

% rotated vector (of pi/4 clock-wise)
r2=A*r1

% check: angle between vectors
costheta=r1'*r2/(norm(r1)*norm(r2));
acos(costheta)
% in degrees
acos(costheta)*180/pi

% rotated vector (of pi/4 anticlock-wise)
r3=A'*r1

% inverse matrix of A
A'*A
A*A'
inv(A)

%%
% Ex 3.0.5
% Rotation and translation of an ellipse
clear all
close all
clc
a=2; b=1;
N=100;
theta=linspace(0,2*pi,N);

x=a*cos(theta);
y=b*sin(theta);

%plot before rotating
figure(1)
hold on
plot(x,y,'b-')
Legend('ellipse')
xlabel('x')
ylabel('y')
axis('equal')
hold off

% rotate
theta0=pi/8;
% rotation matrix
A=[cos(theta0) -sin(theta0); sin(theta0) cos(theta0)];

%rotate and shift
xnew=2+A(1,1)*x+A(1,2)*y;
ynew=3+A(2,1)*x+A(2,2)*y;

%plot after rotating and shifting
figure(2)
hold on
plot(xnew,ynew,'r-')
Legend('ellipse rotated and shifted')
xlabel('x')
ylabel('y')
axis('equal')
hold off



%%
%-----------------------------%
% Ex 3.0.6
clear all
clc
close all

k_B=1.38e-23;
m=67e-27;
T=300;

v=linspace(0,2.0e03,500);
f= 4*pi*(m/(2*pi*k_B*T))^(3/2)*v.^2.*exp(-m*v.^2/(2*k_B*T))

% evaluate the max
f_max=-1.0;

for i=1:length(f)
  if f(i)>f_max
      f_max=f(i);
      i_max=i;
  end
end
f_max
i_max
v_max=v(i_max)

figure(1)
hold on
plot(v,f,'-b','LineWidth',2)
plot(v(i_max),f_max,'ok','MarkerSize',12,'LineWidth',2)
title('Maxwell-Boltzmann distribution')
xlabel('v [m/s]')
ylabel('f_{MB}(v)')
hold off

% evaluate the max - quick way
[f_max,i_max]=max(f)
v(i_max)

%%
% plot of f_max and v_max when changing temeperature
T=linspace(1,1.0e03,50);

for i=1:length(T)
    f= 4*pi*(m/(2*pi*k_B*T(i)))^(3/2)*v.^2.*exp(-m*v.^2/(2*k_B*T(i)));
    [f_max,i_max]=max(f);
    fT_max(i)=f_max;
    vT_max(i)=v(i_max);
end

fT_max
vT_max

figure(2)
subplot(1,2,1)
hold on
plot(T,fT_max,'r-o')
xlabel('T [K]')
ylabel('fT_{max}')
hold off
subplot(1,2,2)
hold on
plot(T,vT_max,'b-x')
xlabel('T [K]')
ylabel('v_{max} [m/s]')
hold off

%%
%-----------------------------%
% Ex 3.0.7 - Minimal distance between two circles
clear all
close all
clc

th=linspace(0,2*pi,1000);
% circle 1
x01=-3;
y01=-0.5;
r1=1;
x1=x01+r1*cos(th);
y1=y01+r1*sin(th);

% circle 2
x02=3;
y02=2;
r2=2.5;
x2=x02+r2*cos(th);
y2=y02+r2*sin(th);

% exact value of the 
% minimal distance between the two circles
min_dis_exact=sqrt((x02-x01)^2+(y02-y01)^2)-(r1+r2)

% minimal distance between circles
min_dist=+100.0;

for i=1:length(th)
    for j=1:length(th)
        dist=sqrt((x2(j)-x1(i))^2+(y2(j)-y1(i))^2);
        if dist<min_dist
           min_dist=dist;
           xmin1=x1(i);
           xmin2=x2(j);
           ymin1=y1(i);
           ymin2=y2(j);
        end
    end
end

min_dist

x_seg=linspace(xmin1,xmin2,100);
slope=(ymin2-ymin1)/(xmin2-xmin1);
y_seg=ymin1+slope*(x_seg-xmin1);

% plot
hold on
patch (x1,y1,'yellow')
plot(x1,y1,'-b','LineWidth',2)
patch (x2,y2,'green')
plot(x2,y2,'-r','LineWidth',2)
plot(x_seg,y_seg,'k--')
plot(xmin1,ymin1,'ok','MarkerSize',12,'LineWidth',2)
plot(xmin2,ymin2,'ok','MarkerSize',12,'LineWidth',2)
xlabel('x')
ylabel('y')
set(gca,'FontSize',18)
%axis([-4.5 6 -3.5 7], 'square')
axis('equal')
hold off
print -f -depsc 'circle.eps'


%%
%-----------------------------%
% Ex 3.0.8 - Position, velocity, and acceleration

clear all
close all
clc
% parameters of the problem 
v0=2.3; Va=10; omega=2; gamma=0.2;
tmin=0; tmax=10;

N=100;
t=linspace(tmin, tmax, N);
[vel] = velo(v0, Va, omega, gamma, t);
% equivalent to:
%vel=v0 + A0*exp(-t/tau).*( cos(omega*t) - omega*tau*sin(omega*t) )

figure(1)
hold on
plot(t,vel,'-b','Linewidth',2)
Legend('particle velocity in interval [tmin, tmax]')
xlabel('t [s]')
ylabel('v [m/s]')
hold off

%%
%
il=(1:1:N-1); ir=(2:1:N);
dt=t(ir)-t(il);
dvel=vel(ir)-vel(il);
acc=dvel./dt;

figure(2)
hold on
plot(t(il),acc,'r-o')
plot(t(ir),acc,'g-*')
plot(t(il)+dt/2,acc,'ks')
Legend('appr 1', 'appr 2', 'appr 3 acceleration')
xlabel('t [s]')
ylabel('a [m/s^{2}]')
hold off

%%

% set initial condition to zero
x_l(1)=0;

for i=2:N
    t0=tmin; t1=t(i); dt_int=0.0001;
    t_int=(t0:dt_int:t1);
    [vel] = velo(v0, Va, omega, gamma, t_int);
    
    N_int=length(t_int);
    i1=(1:1:N_int-1);
    i2=(2:1:N_int);
    
    x_l(i)=sum( vel(i1)*dt_int );
end

%correct initial condition
x0=1;
x_l=x0+x_l

figure(3)
hold on
plot(t,x_l,'r-o')
Legend('appr 1 position')
xlabel('t [s]')
ylabel('x [m]')
hold off
%%
% 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,12);
ym=linspace(-3,3,12);
[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.7
quiver(x,y,Ex1+Ex2,Ey1+Ey2,S)