point mass model using ode45

this is the pointmass model programme for projectile using ode45 with event location, in this the drag coefficient (Cd) is varying with mach number. where mach number itself depends on variation of temperature with altitude. i want that at each altitude value (start from initial condition) calculated by odesolver(ode45) it would calculate the respective temperature and mach, and take the drag value from the conditions. which i have written in the programme. i do this same for density variation with altitude. i have written the formula for temperature and density variation with altitude in the programme. here is the programme

 defining function for equation

function [value,isterminal,dircn] = fnpointmassmodel(t,c,flag,g,w1,w2,w3,Cb)

T0 = 59;
%b = temperature -altitude decay factor;
b = 6.015*1e-06;

T = (T0 + 459.67)*exp(-b*c(3)) - 459.67;
%a = speed of sound;

a = 49.19*(sqrt(T + 459.67));

vel = sqrt((c(2) - w1)^2 + (c(4) - w2)^2 + (c(6) - w3)^2) ;

M = vel/a;

if M ==  0
    Cd = 0.25;
  elseif M>0 && M<0.8
    Cd = 0.25;
  elseif M>0.8 && M<0.9
    Cd = 0.252;
  elseif M>0.9 && M<0.95
    Cd = 0.288;
  elseif M>0.95 && M<1
    Cd = 0.362;
  elseif M>1 && M<1.05
    Cd = 0.443;
  end

% po = air density at firing site
po = 0.0751265;
% h = air density altitude decay factor
h = 3.158*1e-5;
p = po*exp(-h*c(3));
k = (p*pi*Cd)/(8*Cb)

% c(1) = x; c(2) = dx/dt = vx ; c(3) = y ; c(4) = dy/dt = vy ; c(5) = z
% c(6) = dc/dt = vz
vel = sqrt((c(2) - w1)^2 + (c(4) - w2)^2 + (c(6) - w3)^2) ;

if nargin<8 || isempty(flag)

       value = [c(2); -k*vel*(c(2) - w1) ; c(4); (-k*vel*(c(4) - w2)) - g ;c(6); -k*vel*(c(6) - w3) ];

else

      switch flag
          case 'events'

              value = c(3);
              isterminal = 1; 
              dircn = 0;

          otherwise
              error('function not programmed for this event');
      end
  end

solver programme

tspan  = [0 50];

% w1,w2,w3 are the wind components
w1 = 0;
w2 = 0;
w3 = 0;
Cb = 0.226;
v0 = 2800;
theta = 30*pi/180;
g = 32.174;

c0 = [0; v0*cos(theta); 0; v0*sin(theta); 0; 0];

options = odeset('events','on');

[t,c,te,ze,ie] = ode45('fnpointmassmodel',tspan,c0,options,g,w1,w2,w3,Cb);

i have doubt whether this programme will calculate the mach value at each step and take the respective dragcoefficient and give the accurate result or not.

reply soon,

thanks

pawankumar

[pawan kumar]

will ode solver give accurate result when i apply conditions in between steps of odesolver using formula example;

the drag coefficient(Cd) is varying with mach number.where mach number itself depends on variation of temperature with altitude,i want that at each altitude value (start from initial condition) calculated by odesolver(ode45)it would calculate the respective temperature and mach ,and take the drag value from the conditions .which i have written in the programme.

how can i do this

Answer 1

[Walter Roberson]

No, that program cannot work. Refer to the reference documentation here. There is no possibility with ode45() for passing additional variables after the "options" structure.