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 = 6.015*1e-06;
T = (T0 + 459.67)*exp(-b*c(3)) - 459.67;
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 = 0.0751265;
h = 3.158*1e-5;
p = po*exp(-h*c(3));
k = (p*pi*Cd)/(8*Cb)
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 = 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