в этой функции все и расчитывается.
function dzdt = jump(t, w)
global m c s sw g v0 teta ro Vmax anglew azimyt Val Vmax_z Hmax_v n N wan flag MT MV SM0 TM VM
x = w(1);
y= w(2);
z = w(3);
vx = w(4);
vy = w(5);
vz= w(6);
xdot = vx;
ydot = vy;
zdot=vz;
ro=1.225*(1-w(3)/44300)^4.256; %плотность воздуха в зависимости от высоты
S=s;
Sw=sw;
Vxy=sqrt(vx^2+vy^2); % полная скорость в плоскости Х0У
if MV
if abs(Vxy) > VM
S=SM0;
Sw=SM0;
end;
end
if Val
[Vw,azimyt,Hnext]=getparamwind(w(3)); % функция которая возвращает скорость ветра и азимут в зависимости от высоты(табличные данные)
wx=Vw*sin(azimyt/180*pi); % состовляющая скорости ветра по оси Х
wy=Vw*cos(azimyt/180*pi); % состовляющая скорости ветра по оси У
if (wx & wy)
bw=(c*ro*Sw)/(2*m); % лобовое сопротивление ветра
vxwdot= -bw * wx * Vw; % ветровой напор по оси Х
vywdot= -bw * wy * Vw; % ветровой напор по оси У
else
vxwdot=0;
vywdot=0;
end
else
vxwdot=0;
vywdot=0;
end
b = (c*ro*S)/(2*m); % лобовое сопротивлени фрагмента
v = sqrt(w(4)^2 + w(5)^2 + w(6)^2); % общая скорость фрагмента
vxdot = -b * vx * v; % сила сопротивления по оси Х
vydot= -b * vy * v; % сила сопротивления по оси Y
vzdot = -g - b * vz * v; % сила сопротивления по оси Z
vxdot=vxdot+vxwdot; %сложение силы сопротивления фрагмента воздуху с силой сопротивления ветра по оси Х
vydot=vydot+vywdot; % сложение силы сопротивления фрагмента воздуху с силой сопротивления ветра по оси У
dzdt =[ xdot; ydot; zdot; vxdot; vydot; vzdot]; % данные отправляютя в решатель (внутренняя функция матлаб) ode45.
Отредактировано oper (2020-05-04 16:35:23)
