Code covered by the BSD License

# Using Analytical Tools to Gain Insight and Speed-up Num. Analysis in MATLAB & Symbolic Math Toolbox

### Deepak Ramaswamy (view profile)

28 Jan 2013 (Updated )

files from the webinar

simulateRun(x0)
```%
% Copyright 2012 The MathWorks, Inc.
function [time, vel, wheelBrakeTime, stopTime, acel] = simulateRun(x0)
% The dynamical eqations which determine the motion of the car are given by
% Newton's Second Law: The change in momentum is equal to the resultant
% force.
%
%  dP/dt = F_r
%
% where P = M*v. Substituting this in,
%
%  dM/dt * v + dv/dt * M = F_r
%
% which is equivalent to
%
%  dv/dt = [F_r - v*dM/dt]/M
%
rocketStartTime = x0(1);
engineOffTime = x0(2);
chute1Time = x0(3);
chute2Time = x0(4);

% Drag Force
DF = @(t,v) surfaceDrag(t, rocketStartTime, engineOffTime) + aerodynamicDrag(v);

% Mass terms
M = @(t) Mass(t, rocketStartTime, engineOffTime);

%Stage1 - Jet Only
dvdt_1 = @(t,v) (jetThrust(t) - DF(t,v))/M(t);
[t1, v1] = ode45(dvdt_1, [0, rocketStartTime], 0);

% Stage2 - Jet + Rocket
dvdt_2 = @(t,v) (jetThrust(t) + rocketThrust(t-rocketStartTime) - DF(t,v))/M(t);
[t2, v2] = ode45(dvdt_2, [rocketStartTime, engineOffTime], v1(end));

% Stage3 - Coasting
dvdt_3 = @(t,v) -(DF(t,v))/M(t);
[t3, v3] = ode45(dvdt_3, [engineOffTime, chute1Time], v2(end));

% Stage4 - Chute 1 deployed
dvdt_4 = @(t,v) -(DF(t,v) + chute1Drag(v))/M(t);
[t4, v4] = ode45(dvdt_4, [chute1Time,chute2Time], v3(end));

% Stage5 - Chute 2 deployed
dvdt_5 = @(t,v) -(DF(t,v) + chute2Drag(v))/M(t);
[t5, v5] = ode45(dvdt_5, chute2Time+[0, 50], v4(end));

wheelBrakeVel = mph2mps(200);
wheelBrakeTime = spline(v5, t5, wheelBrakeVel);

% Stage 6 - Wheel Brakes only
dvdt_6 = @(t,v) -(DF(t,v))/M(t) - wheelBreakAccel;
[t6, v6] = ode45(dvdt_6, wheelBrakeTime+[0, 50], wheelBrakeVel);

stopTime = spline(v6, t6, 0);

% Prepare time and velocity vectors for concatination
t1(end) =[];
v1(end) = [];

t2(end) = [];
v2(end) = [];

t3(end) = [];
v3(end) = [];

t4(end) = [];
v4(end) = [];

t5 = t5( t5 < wheelBrakeTime );
v5 = v5( t5 < wheelBrakeTime );

t6 = [t6( t6 < stopTime ); stopTime];
v6 = [v6( t6 < stopTime ); 0];

time = vertcat(t1,t2,t3,t4,t5,t6);
vel = vertcat(v1,v2,v3,v4,v5,v6);

acel = [];

% a1 = arrayfun(dvdt_1, t1, v1);
% a2 = arrayfun(dvdt_2, t2, v2);
% a3 = arrayfun(dvdt_3, t3, v3);
% a4 = arrayfun(dvdt_4, t4, v4);
% a5 = arrayfun(dvdt_5, t5, v5);
% a6 = arrayfun(dvdt_6, t6, v6);
%
% acel = vertcat(a1,a2,a3,a4,a5,a6);
```