Homework Help: Matlab code using ode45 difficulty

1. Aug 13, 2011

zfolwick

1. The problem statement, all variables and given/known data

I've been trying to get this code to work:
ode45(@(t,y) fallode(t,y,B), [0, tmax], [31330, 0])

2. Relevant 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

3. The attempt at a solution

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

2. Aug 17, 2011

Nylex

Don't you think it would be helpful if you post the errors?

Last edited: Aug 17, 2011
3. Aug 17, 2011

GreenPrint

I concur, also why don't you just use fprintf to output to the display screen your results with units instead of using notes that don't get outputted to the display screen?

4. Aug 20, 2011

Nylex

It also wouldn't hurt to use
Code (Text):
tags to preserve the indentation, to help readability.