- #1

- 36

- 0

## Homework Statement

I've been trying to get this code to work:

ode45(@(t,y) fallode(t,y,B), [0, tmax], [31330, 0])

## Homework Equations

function dy = fallode(t, y, B)

% ODE function to model a fall over a large range of

% altitudes, reaching up to high subsonic Mach numbers.

% y(1) is altitude, y(2) is vertical velocity (positive up).

% Variation of gravity with altitude is ignored.

g = 9.8; % Earth gravity, m/s^2

R = 287; % Specific gas constant of air, J/kg*K

gamma = 1.4; % Ratio of specific heats of air, dimensionless

[T, rho] = atmos(y(1));

dy(1,1) = y(2);

dy(2,1) = -g + rho*B*y(2)^2/(2*sqrt(1 - y(2)^2/(gamma*R*T)));

--------------

function [T, rho] = atmos(ha)

% Calculate the temperature and density at a given absolute

% altitude ha in the US standard atmosphere model (lowest

% three layers only)

re = 6378.e3; % Radius of Earth, mean equatorial, meters

% Geopotential altitude

h = re*ha/(re + ha);

R = 287; % Specific gas constant for air, J/(kg*K)

g0 = 9.8; % Earth surface gravity, m/s^2

hbreak1 = 11e3; % Altitude of break between 1st and 2nd layers in model

hbreak2 = 25e3; % Altitude of break between 2nd and 3rd layers in model

for ii = 1:length(ha)

if (h <= hbreak1)

a = -6.5e-3; % Temperature lapse rate in troposphere, K/m

rho1 = 1.225; % Surface air density in standard atm. model, kg/m^3

T1 = 288.16; % Surface temperature in standard atm. model, K

T(ii) = T1 + a*h(ii); % Temperature at altitude, K

rho(ii) = rho1*(T(ii)/T1)^(-g0/(a*R)-1); % Density at altitude, kg/m^3

elseif (h <= hbreak2)

T(ii) = 216.66; % Temperature between 11km and 25km altitude, K

rhos = 0.3642; % Density at 11km altitude, kg/m^3

rho(ii) = rhos*exp(-(h(ii) - hbreak1)*g0/(R*T)); % Density at altitude, kg/m^3

else

4a = 3e-3; % Temperature lapse rate in stratosphere, K/m

rho1 = 0.0401; % Density at 25km altitude, kg/m^3

T1 = 216.66; % Temperature at 25 km altitude, K

T(ii) = T1 + a*(h(ii) - hbreak2); % Temperature at altitude, K

rho(ii) = rho1*(T(ii)/T1)^(-g0/(a*R)-1); % Density at altitude, kg/m^3

end

end

## The Attempt at a Solution

I've gotten nothing but errors. Any help troubleshooting would be appreciated.