function HW3_2021

close all


%% Problem 4

H=2; % Height
eta = 18.1e-06; %; Viscosity
rho = 1000; % Density
V0 = 4; % Initial speed
g = 9.81 ; % Gravitational accel

R = linspace(0.025,0.3,100)/1000; % Range of droplet sizes

U0 = -0.5;
for i=1:length(R)
    C = 2*rho*R(i)^2/(9*eta);
    t = fzero(@(t) height(t,H,C,g),[0.05,30]);
    x(i) = U0*t + C*(V0-U0)*(1-exp(-t/C));
end
figure
plot(R,x)



%%Problem 5



L0 = 1;
g = 9.81;
initial_w = [0.;0.1];

dL = 0;
[time,sol] = ode45(@(t,w) pendulum(t,w,g,L0,dL), [0,7], initial_w);
figure
plot(time,sol(:,1),'LineWidth',2,'DisplayName','\Delta L = 0','Color',[1,0,0]);
hold on
dL = 0.1;
[time,sol] = ode45(@(t,w) pendulum(t,w,g,L0,dL), [0,7], initial_w);
plot(time,sol(:,1),'LineWidth',2,'DisplayName','\Delta L = 0.1','Color',[0,1,0]);
dL = 0.15;
[time,sol] = ode45(@(t,w) pendulum(t,w,g,L0,dL), [0,7], initial_w);
plot(time,sol(:,1),'LineWidth',2,'DisplayName','\Delta L = 0.15','Color',[0,0,1]);
xlabel('time (s)')
ylabel('Angle (radians)')
grid on
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.7 0.2 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);


dL = 0.1;
[time,sol] = ode45(@(t,w) pendulum(t,w,g,L0,dL), [0,30], initial_w);
figure
plot(time,sol(:,1),'LineWidth',2,'DisplayName','\Delta L = 0.15','Color',[0,0,1]);
xlabel('time (s)')
ylabel('Angle (radians)')
title('Original approach')
grid on
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');

dL = 0.1;
options = odeset('Events',@(t,w) stop_pendulum(t,w));
sols = [];
times = [];
for i=1:30
  [time_cycle,sol_cycle] = ode45(@(t,w) pendulum(t,w,g,L0,dL), [0,15], initial_w,options);
  initial_w = [0,sol_cycle(end,2)];
  if (isempty(times))
      times = time_cycle;
  else
     times = [times;times(end)+time_cycle];
  end
  sols = [sols;sol_cycle];
end

figure
plot(times,sols(:,1),'LineWidth',2,'Color',[1,0,0])
title('Modified approach')
xlabel('time (s)')
ylabel('Angle (radians)')
grid on
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');



%%




m = 1.e-06;
F0 = 1.e-06;
eta = 1.e-06;
L0 = 0.005;
V = 2.5/1000;
U0 = V/asin(V*eta/F0);
tstop = 2.5;

L = 0.25*L0;
initial_w = [0;L;0;-L;-V;0;V;0];
[times,sols ] = ode45(@(t,w) droplets(t,w,m,F0,U0,eta,L0),[0,tstop],initial_w);

figure
plot(sols(:,1),sols(:,2),'DisplayName','R=L_0/4','Color',[1,0,0])
hold on
plot(sols(:,3),sols(:,4),'DisplayName','R=L_0/4','LineStyle','--','Color',[1,0,0]);

tstop = 0.5;
L = 0.5*L0;
initial_w = [0;L;0;-L;-V;0;V;0];
[times,sols ] = ode45(@(t,w) droplets(t,w,m,F0,U0,eta,L0),[0,tstop],initial_w);

plot(sols(:,1),sols(:,2),'DisplayName','R=L_0/2','Color',[0,1,0])
plot(sols(:,3),sols(:,4),'DisplayName','R=L_0/2','LineStyle','--','Color',[0,1,0]);

tstop = 5;
L = 0.75*L0;
initial_w = [0;L;0;-L;-V;0;V;0];
[times,sols ] = ode45(@(t,w) droplets(t,w,m,F0,U0,eta,L0),[0,tstop],initial_w);

plot(sols(:,1),sols(:,2),'DisplayName','R=3L_0/4','Color',[0,0,1])
plot(sols(:,3),sols(:,4),'DisplayName','R=3L_0/4','LineStyle','--','Color',[0,0,1]);

tstop = 40;
L = 2.6*L0;
initial_w = [0;L;0;-L;-V;0;V;0];
[times,sols ] = ode45(@(t,w) droplets(t,w,m,F0,U0,eta,L0),[0,tstop],initial_w);

plot(sols(:,1),sols(:,2),'DisplayName','R=2.6 L_0','Color',[1,0,1])
plot(sols(:,3),sols(:,4),'DisplayName','R=2.6 L_0','LineStyle','--','Color',[1,0,1]);

xlim([-0.02,0.02])
ylim([-0.02,0.02])
axis equal

axes1 = gca
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.6 0.3 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

end
function z = height(t,H,C,g)


  z = H - C*g*t + (C^2)*g*(1-exp(-t/C));

end


function dwdt = pendulum(t,w,g,L0,dL)
   q = w(1); omega = w(2);
   
   dqdt = omega;
   Omega = sqrt(g/L0);
   L = L0 + dL*sin(2*Omega*t);
   domegadt = (-g*sin(q) - 4*dL*Omega*cos(2*Omega*t)*omega)/L;
   dwdt = [dqdt;domegadt];


end
function [ev,stop,dir] = stop_pendulum(t,w)
   q = w(1);
   ev = q;
   dir = 0;
   stop = 1;
   
end

function dwdt = droplets(t,w,m,F0,U0,eta,L0)

r1 = w(1:2); r2 = w(3:4); v1 = w(5:6); v2 = w(7:8);

V1 = norm(v1);
V2 = norm(v2);
d = norm(r2-r1);

F1 = F0*sin(V1/U0)*v1/V1 - eta*v1 + L0*F0*sin(2*pi*d/L0)/(d^2)*(r1-r2);
F2 = F0*sin(V2/U0)*v2/V2 - eta*v2 - L0*F0*sin(2*pi*d/L0)/(d^2)*(r1-r2);

dwdt = [v1;v2;F1/m;F2/m];

end