%% ENGN0040 Homework 3, 2021
%

syms t
r = 1+t/sqrt(2)
q = log(1+t/sqrt(2))
x = r*cos(q)
y = r*sin(q)
figure
fplot(x,y,[0,8],'LineWidth',2,'Color',[1,0,0])
axis square
axis equal
ylim([-1,7])
xlim([-4,4])
xlabel('x')
ylabel('y')
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');

%

clear all
syms t real
assume(t>0)
H=2; % Height
eta = 18.1e-06; %; Viscosity
rho = 1000; % Density
V0 = 4; % Initial speed
g = 9.81 ; % Gravitational accel
U0vals = [-1,0,1];
R = linspace(0.025,0.3,100)/1000;
figure
hold on
for j=1:3
    U0 = U0vals(j);
    for i=1:length(R)
      C = 2*rho*R(i)^2/(9*eta);
      z = H - C*g*t + (C^2)*g*(1-exp(-t/C));
      tval = vpasolve(z==0,t,20);
      x(i) = U0*tval + C*(V0-U0)*(1-exp(-tval/C));
    end
    plot(R*1000,x,'LineWidth',2,'DisplayName',['U_0=',num2str(U0),'m/s'])
end
xlabel('Droplet Radius (mm)')
ylabel('Distance traveled (m)')
grid on
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.6 0.2 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);