Matlab code using ode45 difficulty

  • Thread starter zfolwick
  • Start date
  • #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.
 

Answers and Replies

  • #2
551
1
I've gotten nothing but errors.
Don't you think it would be helpful if you post the errors?
 
Last edited:
  • #3
1,196
0
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?
 
  • #4
551
1
It also wouldn't hurt to use
Code:
 tags to preserve the indentation, to help readability.
 

Related Threads on Matlab code using ode45 difficulty

  • Last Post
Replies
2
Views
3K
Replies
1
Views
887
  • Last Post
Replies
5
Views
2K
  • Last Post
Replies
1
Views
8K
  • Last Post
Replies
2
Views
2K
Replies
1
Views
2K
  • Last Post
Replies
14
Views
4K
  • Last Post
Replies
1
Views
1K
Replies
2
Views
10K
  • Last Post
Replies
3
Views
1K
Top