%% ENGN0040 Dynamics and Vibrations
%% Homework #2, 2018
%% Solutions to problem 3
%
% This is the simple version of the problem
% Define values for the constants

clear all
T1 = 90000;
T2 = 200000;
m = 7500;
c1 = 0.7;   % c1 is not used in the simple version
c2 = 1.1;

%
% Calculate the velocity during phase 1 by integrating the (constant) acceleration

% Phase (i)
syms v1(t) x1(t) t tend
a1 = T1/m;
v1(t) = int(a1,t)
x1(t) = int(v1,t)
tend = solve(v1(t)==150,t);    % Find the time when speed reaches 150 m/s
x1(tend);

%
% Calculate the velocity during phase 2 by integrating the (constant) acceleration
% We have written the integral as
% \int_{v1\left(tend\right)}^v \mathit{dv}=\int_{tend}^t \mathit{adt}

% Phase (ii)
syms v2(t) x2(t)
a2 = T2/m;
v2(t) = int(a2,t,[tend,t]) + v1(tend)
x2(t) = int(v2,t,[tend,t]) + x1(tend)
double(v2(tend+20))
double(x2(tend+20))

%
% Calculate the velocity during phase 3 by separating variables and integrating both sides
% We have written the integral as
% \int_{\text{v}2\left(\text{tend}+20\right)}^v \frac{\mathit{\text{dv}}}{v^2 }=\int_{\text{tend}+20}^t −\frac{c}{m}\mathit{\text{dt}}

% Phase (iii)
syms vv x3(t) V v3(t)
veq = int(1/V^2,[v2(tend+20),vv])==int(-c2/m,t,[tend+20,t]); % Separate variables and integrate
v3(t) = simplify(solve(veq,vv))
x3(t) = simplify(int(v3,t,[tend+20,t]) + x2(tend+20))

%
% The plots

fplot(v1(t),[0 double(tend)],'Color',[0.85 0.325 0.098],...
     'LineWidth',2);
hold on 
fplot(v2(t),[double(tend) double(tend)+20],'Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(v3(t),[double(tend)+20 120],'Color',[0.85 0.325 0.098],...
     'LineWidth',2);
axes1 = gca;
xlabel({'Time (s)'},'FontAngle','italic','FontName','Times New Roman');
title('Predicted velocity for Bloodhound test');
ylabel({'Speed (m/s)'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');

%
% This is the optional advanced version of the problem
% Define values for the constants

clear all
T1 = 90000.;
T2 = 200000.;
m = 7500.;
c1 = 0.7;
c2 = 1.1;

%

syms V vv v1(t) x1(t) t tend
assume(vv,'real')
assume(t,'real')
assume(T1/m-(c1/m)*vv^2>0);
a1 = T1/m - (c1/m)*V^2;
veq1 = int(1/a1,V,[0 vv]) == t;
v1(t) = solve(veq1,vv)
tend = solve(v1(t)==150,t)

%

syms v2(t) 
assume(T2/m-(c1/m)*vv^2>0)
a2 = T2/m - (c1/m)*V^2;
veq2 = int(1/a2,V,[v1(tend) vv]) == t-tend;
v2(t) = simplify(solve(veq2,vv))

%

% Phase (iii)
syms vv x3(t) V v3(t)
veq = int(1/V^2,[v2(tend+20),vv])==int(-c2/m,t,[tend+20,t]); % Separate variables and integrate
v3(t) = simplify(solve(veq,vv))

%

figure
fplot(v1(t),[0 double(tend)],'Color',[0.85 0.325 0.098],...
     'LineWidth',2);
hold on 
fplot(v2(t),[double(tend) double(tend)+20],'Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(v3(t),[double(tend)+20 120],'Color',[0.85 0.325 0.098],...
     'LineWidth',2);
axes1 = gca;
xlabel({'Time (s)'},'FontAngle','italic','FontName','Times New Roman');
title('Predicted velocity for Bloodhound test');
ylabel({'Speed (m/s)'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');