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_breakingVelsSpecified(rocketStartTime, engineOffTime)
```%
% Copyright 2012 The MathWorks, Inc.
function [time, vel, acel, chute1Time, chute2Time, wheelBrakeTime] = simulateRun_breakingVelsSpecified(rocketStartTime, engineOffTime)
% 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
%

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

% Mass terms
vdMdt = @(t,v) v * dMassDt(t, rocketStartTime, engineOffTime);
M = @(t) Mass(t, rocketStartTime, engineOffTime);

%Stage1 - Jet Only
dvdt_1 = @(t,v) (jetThrust(t) - DF(t,v) - vdMdt(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) - vdMdt(t,v))/M(t);
[t2, v2] = ode45(dvdt_2, [rocketStartTime, engineOffTime], v1(end));

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

chute1DeployVel =  mph2mps(700);
chute1Time = spline(v3, t3, chute1DeployVel);

% Stage4 - Chute 1 deployed
dvdt_4 = @(t,v) -(DF(t,v) + chute1Drag(v) + vdMdt(t,v))/M(t);
[t4, v4] = ode45(dvdt_4, chute1Time+[0, 30], chute1DeployVel);

chute2DeployVel = mph2mps(500);
chute2Time = spline(v4, t4, chute2DeployVel);

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

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

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

stopTime = spline(v6, t6, 0);

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

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

t3 = t3( t3 < chute1Time );
v3 = v3( t3 < chute1Time );

t4 = t4( t4 < chute2Time );
v4 = v4( t4 < chute2Time );

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);

% Acceleration (dv/dt)
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);
```