1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab code using ode45 difficulty

  1. Aug 13, 2011 #1
    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. jcsd
  3. Aug 17, 2011 #2
    Don't you think it would be helpful if you post the errors?
     
    Last edited: Aug 17, 2011
  4. Aug 17, 2011 #3
    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?
     
  5. Aug 20, 2011 #4
    It also wouldn't hurt to use
    Code (Text):
     tags to preserve the indentation, to help readability.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Matlab code using ode45 difficulty
  1. Matlab: ode45 help (Replies: 1)

  2. Matlab ode45 Problem (Replies: 5)

Loading...