Динамика полета летательного аппарата

function [h]=pol()
Rz=6378.136;
t0=70000;
ra=240+Rz; %radius apogeya
rp=190+Rz; %radius perigeya
m=398602; %grav parametr zemli
i=63; %naklonenie
W=15; %ugol voshod uzla
v=0; %ist anomaliya
w=0; %argument perigeya
u=w+v; 
e=(ra-rp)/(ra+rp); %ekscentrisitet
p=ra*(1-e); %fokal parmetr
Vr=sqrt(m/p)*e*sind(v);
Vn=sqrt(m/p)*(1+e*cosd(v));
R=p/(1+e*cosd(v));  %radius vektor
%perehod v koordinaty
x0=R*(cosd(W)*cosd(u)-sind(W)*sind(u)*cosd(i));
y0=R*(sind(W)*cosd(u)+cosd(W)*sind(u)*cosd(i));
z0=R*sind(u)*sind(i);
Vx0=Vr*(cosd(W)*cosd(u)-sind(W)*sind(u)*cosd(i))-Vn*(cosd(W)*sind(u)+sind(W)*cosd(u)*cosd(i));
Vy0=Vr*(sind(W)*cosd(u)+cosd(W)*sind(u)*cosd(i))-Vn*(sind(W)*sind(u)-cosd(W)*cosd(u)*cosd(i));
Vz0=Vr*sind(u)*sind(i)+Vn*cosd(u)*sind(i);
f0=[x0 y0 z0 Vx0 Vy0 Vz0];
options=odeset('RelTol',1e-6);
[t h]=ode23(@model,[0 70000],f0,options);
r=sqrt(h(:,1).^2+h(:,2).^2+h(:,3).^2)-Rz;
Vv=sqrt(h(:,4).^2+h(:,5).^2+h(:,6).^2);
r1=sqrt(h(:,1).^2+h(:,2).^2+h(:,3).^2);
a=(m*r1(1,1))/(2*m-r1(1,1)*Vv(1,1).^2);
p=((h(1,1)*h(1,5)-h(1,2)*h(1,4)).^2+(h(1,2)*h(1,6)-h(1,3)*h(1,5)).^2+(h(1,1)*h(1,6)-h(1,3)*h(1,4)).^2)/m;
ra=a*(1+sqrt(1-p/a));
rp=p/(1+sqrt(1-p/a));
ee=(ra-rp)/(ra+rp);
for i=1:1:10
    if t=t(:,1)/(2*i)
        v=2*acosd((r1(:,1)-p)./(r1(:,1)*ee));
    else
        v=2*acosd(-(r1(:,1)-p)./(r1(:,1)*ee));
    end
end
figure 
plot(h(:,2),h(:,1))
%comet3(h(:, 1), h(:, 2), h(:, 3));
figure 
plot(t,r)
%figure 
%plot(t,v)
end
 
 
function [f]=model(t,A)
f=zeros(6,1);
Rz=6378.136;
C20=-1082.6274*10^(-6);
C40=-2.371*10^(-8);
m=398602;
f(1)=A(4);
f(2)=A(5);
f(3)=A(6);
r=sqrt(A(1).^2+A(2).^2+A(3).^2);
f(4)=(-m/r^3)*A(1)+0*((3/2)*C20*(m/r^2)*(Rz^2/r^2)*((5*A(3)^2/r^2)-1)+((5/8)*C40*(m/r^2)*(Rz/r)^4*(3+(63*(A(3)^2)/(r^2)-42)*(A(3)^2)/r^2)));
f(5)=(-m/r^3)*A(2)+0*((3/2)*C20*(m/r^2)*(Rz^2/r^2)*((5*A(3)^2/r^2)-1)+((5/8)*C40*(m/r^2)*(Rz/r)^4*(3+(63*(A(3)^2)/(r^2)-42)*(A(3)^2)/r^2)));
f(6)=(-m/r^3)*A(3)+0*((3/2)*C20*(m/r^2)*(Rz^2/r^2)*((5*A(3)^2/r^2)-3)+((5/8)*0.5*C40*(m/r^2)*(Rz/r)^4*(30+(63*(A(3)^2)/(r^2)-70)*(A(3)^2)/r^2)));
end