function HW3_2017

%% ENGN0040 Dynamics and Vibrations
%% Homework #3, 2017: Solutions to MATLAB questions


 close all;
 
%% Problem 2
g = 9.81;      % Gravitational acceleration, m/s^2
R0 = 0.1;      % Outer radius of separator, m
U0 = 10;       % Flow speed, m/s
eta = 1.8*10^(-5);   % Air viscosity
rho = 500;     % Density, kg/m^3
b = 5*10^(-6); % Particle radius, m
c = 9*eta/(4*rho*b^2);  % Drag constant.

initial_w = [0.01;0;0;0;0;0];
options = odeset('RelTol',0.00001);
[times,sols] = ode45(@(t,w) separator(t,w,g,U0,R0,c),[0,50],initial_w,options);

x = sols(:,1).*cos(sols(:,2));
y = sols(:,1).*sin(sols(:,2));
z = sols(:,3);
figure1 = figure;
axes1 = axes('Parent',figure1);
plot3(x,y,z);
xlabel({'x(m)'});
zlabel({'z (m)'});
title({'Trajectory for 5 \mu m particle'});
ylabel({'y(m)'});
view(axes1,[-37.5 30]);
grid(axes1,'on');
set(axes1,'FontSize',14);


b = 10*10^(-6); % Particle radius, m
c = 9*eta/(4*rho*b^2);  % Drag constant.

initial_w = [0.01;0;0;0;0;0];
options = odeset('RelTol',0.00001);
[times,sols] = ode45(@(t,w) separator(t,w,g,U0,R0,c),[0,50],initial_w,options);

x = sols(:,1).*cos(sols(:,2));
y = sols(:,1).*sin(sols(:,2));
z = sols(:,3);
figure1 = figure;
axes1 = axes('Parent',figure1);
plot3(x,y,z);
xlabel({'x(m)'});
zlabel({'z (m)'});
title({'Trajectory for 10 \mu m particle'});
ylabel({'y(m)'});
view(axes1,[-37.5 30]);
grid(axes1,'on');
set(axes1,'FontSize',14); 

b = 20*10^(-6); % Particle radius, m
c = 9*eta/(4*rho*b^2);  % Drag constant.

initial_w = [0.01;0;0;0;0;0];
options = odeset('RelTol',0.00001);
[times,sols] = ode45(@(t,w) separator(t,w,g,U0,R0,c),[0,50],initial_w,options);

x = sols(:,1).*cos(sols(:,2));
y = sols(:,1).*sin(sols(:,2));
z = sols(:,3);
figure1 = figure;
axes1 = axes('Parent',figure1);
plot3(x,y,z);
xlabel({'x(m)'});
zlabel({'z (m)'});
title({'Trajectory for 20 \mu m particle'});
ylabel({'y(m)'});
view(axes1,[-37.5 30]);
grid(axes1,'on');
set(axes1,'FontSize',14);


%% Problem 3
Re = 6378000;   % Earths radius, m
m = 1.673*10^(-27);  % Proton mass, kg
q = 1.6022*10^(-19); % Proton charge, Coulombs
B0 = 3.07*10^(-5);   % Magnetic field strength, Tesla
c = q*B0*Re^3/m;

initial_w = [4*Re;0;0;0;7.7e06;13.3e06];

 options = odeset('RelTol',0.00001);
[times,sols] = ode45(@(t,w) proton(t,w,c),[0,500],initial_w,options);


% Load the globe
   load('topo.mat','topo','topomap1')
   props.AmbientStrength = 0.1;
   props.DiffuseStrength = 1;
   props.SpecularColorReflectance = .5;
   props.SpecularExponent = 20;
   props.SpecularStrength = 1;
   props.FaceColor= 'texture';
   props.EdgeColor = 'none';
   props.FaceLighting = 'phong';
   props.Cdata = topo;

   [XX,YY,ZZ] = sphere;
   XX = XX*Re;
   YY = YY*Re;
   ZZ = ZZ*Re;
   
figure1 = figure;
axes1 = axes('Parent',figure1);
   plot3(sols(:,1),sols(:,2),sols(:,3))
   hold on
   axis square
   axis equal
   surface(XX,YY,ZZ,props);
xlabel({'x(m)'});
zlabel({'z (m)'});
title({'Trajectory of a proton in the earths magnetic field'});
ylabel({'y(m)'});
view(axes1,[-37.5 30]);
grid(axes1,'on');
set(axes1,'FontSize',14);



%% Problem 4: 
 k = 40;   % Cantilever stiffness, N/m
 m = 8.27e-12;  % AFM tip mass, kg
 y0 = 1*10^(-9);  % Vibration ampltude, m
 o = 2.4e06;  % Vibration frequency, rad/s
 d = 5*10^(-9); % Force range, m
 F0= 1.e-07; % Force max, N
 c = 2.0e-06;   % damping, Ns/m

 Delta = 10.e-09;
 initial_w = [Delta;0];
 [times,sols] = ode45(@(t,w) afm(t,w,k,m,y0,o,F0,d,Delta,c),[0,1.e-04],initial_w);
  
  figure1 = figure;

% Create axes
  axes1 = axes('Parent',figure1);
  hold(axes1,'on');
  plot(times,10^9*(sols(:,1)-Delta),'DisplayName','\Delta = 20nm');

  Delta = 40.e-09;
  initial_w = [Delta;0];
  [times,sols] = ode45(@(t,w) afm(t,w,k,m,y0,o,F0,d,Delta,c),[0,1.e-04],initial_w);
  plot(times,10^9*(sols(:,1)-Delta),'DisplayName','\Delta = 40nm')
  xlabel({'Time (s)'},'FontSize',20);
  title({'Vibration of an AFM cantilever'});
  ylabel({'Deflection (nm)'},'FontSize',20);
  box(axes1,'on');
  set(axes1,'XGrid','on','YGrid','on');
  legend1 = legend(axes1,'show');
  set(legend1,...
    'Position',[0.784670009941435 0.804277645569833 0.119047617030921 0.0883777216618056]);
  
  
  

end
function dwdt = separator(t,w,g,U0,R0,c)
    r = w(1); q = w(2); z = w(3); drdt = w(4); dqdt = w(5); dzdt = w(6);
    
    u = U0*( (r/R0)^2 - (r/R0)^20 )/(0.04+(r/R0)^3);
    if (r>R0) u=0; end
    
    d2qdt2 = (2*c*(u-r*dqdt) - 2*drdt*dqdt)/r;
    d2rdt2 = -2*c*drdt + r*dqdt^2;
    d2zdt2 = -g - 2*c*dzdt;
    
    dwdt = [drdt;dqdt;dzdt;d2rdt2;d2qdt2;d2zdt2];

end
function dwdt = afm(t,w,k,m,y0,o,F0,d,Delta,c)

   x = w(1); dxdt = w(2);
   
   y=y0*sin(o*t)+Delta;
   
   Ftip = F0*(x/d)*exp(-x/d);
   F = (k*(y-x) - Ftip-c*dxdt)/m;

   dwdt = [dxdt;F];

end
function dwdt = proton(t,w,c)
   x = w(1); y=w(2); z=w(3);
   vx = w(4); vy=w(5); vz=w(6);
   
   R = sqrt(x^2+y^2+z^2);
   dvxdt = c*(-vy*(2*z^2-x^2-y^2) + 3*vz*y*z)/R^5;
   dvydt = c*(vx*(2*z^2-x^2-y^2) - 3*vz*x*z)/R^5;
   dvzdt = -3*c*(vx*y*z-vy*x*z)/R^5;

   dwdt = [vx;vy;vz;dvxdt;dvydt;dvzdt];
end