Matlab code using ode45 difficulty

[zfolwick]

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.

 

[Nylex]

I've gotten nothing but errors.

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

 

[GreenPrint]

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

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?

 

[Nylex]

It also wouldn't hurt to use

 tags to preserve the indentation, to help readability.

 

[DeeNos]

Based on the provided code, it seems like you are trying to model the fall of an object over a large range of altitudes, taking into account air density and temperature variations. The first thing you should do when troubleshooting a code is to check the error messages you are getting. This will give you a better understanding of where the problem lies.

In this case, it seems like the issue might be with the input arguments for the ode45 function. The first argument should be a function handle, but in your code, you are passing in a function call. Instead, it should be ode45(@fallode, [0, tmax], [31330, 0]).

Additionally, it seems like the function fallode is not defined in the code you provided. Make sure you have defined it correctly and have included it in your code.

Another potential issue could be with the values of tmax and B. Make sure they are defined correctly and are within the appropriate range for your problem.

Lastly, it is always a good idea to start with a simple test case and make sure your code is working before moving on to more complex scenarios. This will help you isolate any issues and make troubleshooting easier. Good luck!