Files
Hayden-Project/code/T3dxdt.m
2026-03-28 16:33:32 -05:00

36 lines
1.4 KiB
Matlab

function TxDot = T3dxdt(t,x)
%% Threat is a 240 mm Rocket
% Compute Rocket xDot = [velocity; acceleration]'
g = 9.8066; % acceleration due to gravity
gravity = [0; 0; -g]; % gravity vector
diameter = 0.24; % missile diameter in m.
S = pi*(diameter/2)^2; % reference area in m^2
[acousticSpeed,rho] = getRho(x);
mach = norm(x(4:6))/acousticSpeed;
% Transcribed from the provided target-data plots.
timeCurve = [0 7.5 7.5001 45]'; %seconds
thrustCurve = [77000 77000 0 0]'; % Newtons
massCurve = [410 210 210 210]'; % kg
machCurve = [0 0.5 0.7 0.9 1.0 1.1 1.2 1.3 1.5 2.0 2.5 3.0 3.5 4.0];
dragCurve = [0.40 0.40 0.43 0.60 0.75 0.87 0.86 0.80 0.74 0.60 0.50 0.44 0.40 0.39];
machCurveB = [0 0.5 0.7 0.9 1.0 1.1 1.2 1.3 1.5 2.0 2.5 3.0 3.5 4.0];
dragCurveB = [0.31 0.31 0.34 0.47 0.62 0.78 0.77 0.72 0.64 0.50 0.40 0.34 0.30 0.27];
if(t<=timeCurve(end-1))
Cd = interp1(machCurveB,dragCurveB,mach,'linear','extrap');
else
Cd = interp1(machCurve,dragCurve,mach,'linear','extrap');
end
mass = interp1(timeCurve,massCurve,t,'linear','extrap');
thrust = interp1(timeCurve,thrustCurve,t,'linear','extrap');
%Compute acceleration --------------------------------------------
dragAccel = -0.5*rho*Cd*S*x(4:6)*norm(x(4:6))/mass;% drag acceleration in m/sec^2
thrustAccel = thrust*x(4:6)/(mass*(norm(x(4:6))));
accel = dragAccel + thrustAccel + gravity; % total acceleration
TxDot = [x(4:6);accel];
end