m_t = m*(1-exp(-dV_ss/W));
cout << "Потребный импульс: \n - перицентра dV_ps="
<< dV_ps << "\n апоцентра dV_as=" << dV_as
<< "\n Суммарный импульс=" << dV_ss << "Масса топлива=" << m_t
<< '\n';
dV_ps = 0;
dV_as = 0;
dV_is = 0;
dV_ss = 0;
m_t = 0;
}
}
}
void par_or(real *f, real *par)
{
real x = f[0];
real y = f[1];
real z = f[2];
real Vx = f[3];
real Vy = f[4];
real Vz = f[5];
real c1 = (y*Vz-z*Vy);
real c2 = (z*Vx-x*Vz);
real c3 = (x*Vy-y*Vx);
real C = sqrt(c1*c1+c2*c2+c3*c3);
par[0] = (C/mu_z)*C;
real R_ka = sqrt(x*x+y*y+z*z);
real V_ka = sqrt(Vx*Vx+Vy*Vy+Vz*Vz);
real f1 = (Vy*c3-Vz*c2)-(mu_z*x/R_ka);
real f2 = (Vz*c1-Vx*c3)-(mu_z*y/R_ka);
real f3 = (Vx*c2-Vy*c1)-(mu_z*z/R_ka);
real F = sqrt(f1*f1+f2*f2+f3*f3);
real h = V_ka*V_ka-(2*mu_z/R_ka);
if ((1+h*C*C/(mu_z*mu_z)) < 0)
{
cout << " Ошибка! \n";
getch();
}
par[1] = F/mu_z;
if ((1-par[1]*par[1]) < 0)
{
cout << " (1-e*e) < 0 Ошибка! \n";
getch();
}
par[2] = par[0]/(1-par[1]*par[1]);
par[4] = acos(c3/C);
real s_Om = c1/(C*sin(par[4]));
real c_Om = -c2/(C*sin(par[4]));
if (s_Om >= 0)
par[3] = acos(c_Om);
else
par[3] = 2*M_PI-acos(c_Om);
real c_om = (f1*cos(par[3])+f2*sin(par[3]))/F;
real s_om = f3/(F*sin(par[4]));
if (s_om > 0)
par[5] = acos(c_om);
else
par[5] = 2*M_PI - acos(c_om);
if (par[2] < 0)
{
cout << " Ошибка! \n";
getch();
}
par[6] = 2*M_PI*sqrt((par[2]/mu_z)*par[2]*par[2]);
real c_u = (x*cos(par[3])+y*sin(par[3]))/R_ka;
real s_u = z/(R_ka*sin(par[4]));
if (s_u > 0)
par[7] = acos(c_u);
else
par[7] = 2*M_PI - acos(c_u);
}
#include "rk5.h"
#include <iostream.h>
void Drkgs(real *prmt,real *y,real *dery,int ndim,int& ihlf,
void (*fct)(real &,real*,real*),
void (*out_p)(real,real*,real*,int,int,real*))
{
static real a[] = { 0.5, 0.292893218811345248, 1.70710678118665475,
0.16666666666666667 };
static real b[] = { 2.0, 1.0, 1.0, 2.0 };
static real c[] = { 0.5, 0.292893218811345248, 1.70710678118665475, 0.5 };
real *aux[8];
int i,j,imod,itest,irec,istep,iend;
real delt,aj,bj,cj,r,r1,r2,x,xend,h;
for (i=0; i<8; i++) aux[i] = new real[ndim];
for (i=0; i<ndim; i++) aux[7][i] = (1./15.)*dery[i];
x = prmt[0];
xend = prmt[1];
h = prmt[2];
prmt[4] = 0.0;
fct(x,y,dery);
r = h*(xend-x);
if (r <= 0.0)
{
ihlf = 13;
if (r == 0.0) ihlf = 12;
goto l39;
}
for(i=0; i<ndim; i++)
{
aux[0][i] = y[i];
aux[1][i] = dery[i];
aux[2][i] = 0.0;
aux[5][i] = 0.0;
}
irec = 0;
h = h+h;
ihlf = -1;
istep = 0;
iend = 0;
l4: r = (x+h-xend)*h;
if (r >= 0.0)
{
iend = 1;
if (r > 0.0) h = xend-x;
}
out_p(x,y,dery,irec,ndim,prmt);
if (prmt[4] != 0.0) goto l40;
itest = 0;
l9: istep++;
j = 0;
l10: aj = a[j];
bj =b[j];
cj = c[j];
for (i=0; i<ndim; i++)
{
r1 = h*dery[i];
r2 = aj*(r1-bj*aux[5][i]);
y[i] = y[i]+r2;
r2 = r2+r2+r2;
aux[5][i] += r2-cj*r1;
}
if (j-3 < 0)
{
j++;
if (j-2 != 0) x = x+0.5*h;
fct(x,y,dery);
goto l10;
}
if (itest <= 0)
{
for (i=0; i<ndim; i++) aux[3][i] = y[i];
itest = 1;
istep = istep+istep-2;
l18: ihlf++;
x = x-h;
h = 0.5*h;
for (i=0; i<ndim; i++)
{
y[i] = aux[0][i];
dery[i] = aux[1][i];
aux[5][i] = aux[2][i];
}
goto l9;
}
imod = istep/2;
if (istep-imod-imod != 0)
{
fct(x,y,dery);
for (i=0; i<ndim; i++)
{
aux[4][i] = y[i];
aux[6][i] = dery[i];
}
goto l9;
}
delt = 0.0;
for (i=0; i<ndim; i++)
delt += aux[7][i]*fabs(aux[3][i]-y[i]);
if (delt-prmt[3] > 0.0)
{
if (ihlf-10 >= 0)
{
ihlf = 11;
fct(x,y,dery);
goto l39;
}
for (i=0; i<ndim; i++) aux[3][i] = aux[4][i];
istep = istep+istep-4;
x = x-h;
iend = 0;
goto l18;
}
fct(x,y,dery);
for (i=0; i<ndim; i++)
{
aux[0][i] = y[i];
aux[1][i] = dery[i];
aux[2][i] = aux[5][i];
y[i] = aux[4][i];
dery[i] = aux[6][i];
}
out_p(x-h,y,dery,ihlf,ndim,prmt);
if (prmt[4] != 0) goto l40;
for (i=0; i<ndim; i++)
{
y[i] = aux[0][i];
dery[i] = aux[1][i];
}
irec = ihlf;
if (iend > 0) goto l39;
ihlf--;
istep = istep/2;
h = h+h;
if (ihlf < 0) goto l4;
imod = istep/2;
if ((istep-2*imod != 0) || (delt-0.02*prmt[3] > 0.0)) goto l4;
ihlf--;
istep = istep/2;
h = h+h;
goto l4;
l39: out_p(x,y,dery,ihlf,ndim,prmt);
l40: for (i=0; i<ndim; i++) delete aux[i];
return;
}
6.3. ФАЙЛ НАЧАЛЬНОЙ ИНИЦИАЛИЗАЦИИ INIT.H
ifndef _INIT
#define _INIT
#include "def.h"
#include <stdlib.h>
#include <fstream.h>
ifstream if_init;
void nex_ln (void);
void init_m()
{
Np = 150;
t_beg = 0;
t_end = 8000000;
dt = 2;
toler = .05;
dTp = (t_end-t_beg)/float(Np);
Curp = 0;
J1 = 532;
J2 = 563;
J3 = 697;
m = 597.;
W = 2200;
mu_z = 3.9858e14;
mu_s = 1.3249e20;
mu_l = 4.9027e12;
w_s = 2*M_PI/(365.2422*24*3600);
w_z = 2*M_PI/(24*3600);
w_l = 2*M_PI/(27.32*24*3600);
ww_l = 2*M_PI/(18.6*365.2422*24*3600);
parn[0] = 6952137.;
parn[1] = 0;
parn[2] = 6952137;
parn[3] = 28.1*g_r;
parn[4] = 97.6*g_r;
parn[5] = 63.1968*g_r;
parn[6] = 5769.;
parn[7] = 5.751*g_r;
Fl_u = 1;
u_last = parn[7];
Fl_ka = 0;
Fl_kp = 0;
Fl_ki = 0;
Fl_p = 0;
Fl_a = 0;
Fl_i = 0;
Fl_pkT = 0;
Tkor = 0;
T_vd = 0;
akor[0] = 0;
akor[1] = 0;
akor[2] = 0;
dV_ps = 0;
dV_as = 0;
dV_is = 0;
dV_ss = 0;
Fl_l0 = 0;
Fl_l1 = 0;
Fl_pki = 0;
real x0 = 6137262.9+7000;
real y0 = 3171846.1+7000;
real z0 = 689506.95+7000;
real Vx0 = -201.288+5;
real Vy0 = -1247.027+5;
real Vz0 = 7472.65+5;
prmt[0] = t_beg;
prmt[1] = t_end;
prmt[2] = dt;
prmt[3] = toler;
prmt[4] = 0.0;
y_main[0] = x0;
y_main[1] = y0;
y_main[2] = z0;
y_main[3] = Vx0;
y_main[4] = Vy0;
y_main[5] = Vz0;
}
void nex_ln (void)
{
char ch;
if_init.get(ch);
while (ch != '\n')
if_init.get(ch);
}
#endif
6.4 ФАЙЛ ОПИСАНИЯ ПЕРЕМЕННЫХ DEF.H
#ifndef _DEFH
#define _DEFH
#include <math.h>
typedef long double real;
extern const float g_r;
extern const float r_g;
extern int Np;
extern int Curp;
extern real dTp;
extern real t_beg;
extern real t_end;
extern real dt;
extern real toler;
extern real J1,J2,J3;
extern real mu_z;
extern real mu_s;
extern real mu_l;
extern real m;
extern real m_t;
extern real W;
extern real w_s;
extern real w_z;
extern real w_l;
extern real ww_l;
extern real xs,ys,zs;
extern real xl,yl,zl;
extern real Fz,Fs,Fl,Fa,U20;
extern int nomin;
extern real par[8];
extern real parn[8];
extern real a_p,e_p,p_p,Om_p,i_p,om_p,Rp_p,Ra_p;
extern real y_main[6];
extern real prmt[5];
extern int Fl_u;
extern real u_last;
extern int Fl_ka;
extern int Fl_kp;
extern int Fl_ki;
extern int Fl_i;
extern int Fl_p;
extern int Fl_a;
extern int Fl_lu;
extern int Fl_pkT;
extern real dl;
extern real T_vd;
extern real dRa;
extern real dRp;
extern int Sig;
extern int Sig_a;
extern real Vkor[3];
extern real akor[3];
extern real Tkor;
extern real Tkore;
extern real dV_ps;
extern real dV_as;
extern real dV_is;
extern real dV_ss;
extern int Fl_l0;
extern int Fl_l1;
extern int Fl_pki;
#endif
6.5 ФАЙЛ SFUN.H
#ifndef _SFUN
#define _SFUN
#include "def.h"
#include <iostream.h>
#include <conio.h>
#include <math.h>
void out_p(real x,real *y,real*,int,int,real *);
real interpl(real*,real*,int,real);
void fct(real& ,real *y,real *dery);
void par_or(real *,real *);
#endif
6.5 ФАЙЛ RK5.H
#ifndef _RK5
#define _RK5
#include "def.h"
#include <iostream.h>
#include <conio.h>
#include "sfun.h"
void Drkgs(real *prmt,real *y,real *dery,int ndim,int& ihlf,
void (*fct)(real&,real*,real*),
void (*out_p)(real,real*,real*,int,int,real*));
#endif
6.6 ПРОГРАММА ПОСТРОЕНИЯ ВРЕМЕННЫХ ДИАГРАММ
clc
g_r = pi/180;
r_g = 180/pi;
load m_y.dat
t = m_y(:,1);
x = m_y(:,2);
y = m_y(:,3);
z = m_y(:,4);
Vx = m_y(:,5);
Vy = m_y(:,6);
Vz = m_y(:,7);
clear m_y;
s_tmp = size(t);
s_m = s_tmp(1);
clear s_tmp;
load m_f.dat
Fz = m_f(:,2);
Fs = m_f(:,3);
Fl = m_f(:,4);
Fa = m_f(:,5);
U20 = m_f(:,6);
clear m_f;
load m_s.dat
xs = m_s(:,2);
ys = m_s(:,3);
zs = m_s(:,4);
clear m_s;
load m_par.dat
p = m_par(:,2);
e = m_par(:,3);
a = m_par(:,4);
Om = m_par(:,5);
i = m_par(:,6);
omg = m_par(:,7);
T = m_par(:,8);
u = m_par(:,9);
clear m_par;
p_n = 6952137.;
e_n = 0;
a_n = 6952137.;
Om_n0 = 28.1*g_r;
i_n = 97.6*g_r;
omg_n = 346.725*g_r;
T_n = 5765;
ws = 2*pi/(365.2422*24*3600);
for j = 1:s_m, tmp(j) = Om_n0+ws*t(j);
end
Om_n = tmp';
clear tmp;
map = [1,1,1];
colormap(map);
plot(t,p,'y-',[min(t) max(t)],[p_n p_n],'r-'), title (' Фокальный параметр '), grid on;
print -dwin;
pause;
plot(t,p-p_n,'y-'), title (' dp '), grid on;
print -dwin;
pause;
plot(t,e,'y-',[min(t) max(t)],[e_n e_n],'r-'), title (' Эксцентриситет '), grid on;
print -dwin;
pause;
plot(t,e-e_n,'y-'), title (' de '), grid on;
print -dwin;
pause;
plot(t,a,'y-',[min(t) max(t)],[a_n a_n],'r-'), title (' Большая полуось орбиты '), grid on;
print -dwin;
pause;
plot(t,a-a_n,'y-'), title (' da '), grid on;
print -dwin;
pause;
plot(t,Om*r_g,'y-',t,Om_n*r_g,'r-'), title (' Долгота восходящего узла '), grid on;
print -dwin;
pause;
plot(t,Om*r_g-Om_n*r_g,'y-'), title (' dOm '), grid on;
print -dwin;
pause;
plot(t,i*r_g,'y-',[min(t) max(t)],[i_n*r_g i_n*r_g],'r-'), title (' Наклонение '), grid on;
print -dwin;
pause;
plot(t,i*r_g-i_n*r_g,'y-'), title (' di '), grid on;
print -dwin;
pause;
plot(t,T,'y-',[min(t) max(t)],[T_n T_n], 'r-'), title (' Период '), grid on;
print -dwin;
pause;
plot(t,T-T_n,'y-'), title (' dT '), grid on;
print -dwin;
pause;
plot3(x,y,z,'b')
axis([min(x) max(x) min(y) max(y) min(z) max(z)])
set(gca,'box','on')
title (' Положение МКА ')
hold on
plt = plot3(0,0,0,'.','erasemode','xor','markersize',24);
dk = ceil(length(y)/2500);
for k = 1:dk:length(y)
set(plt,'xdata',x(k),'ydata',y(k),'zdata',z(k))
drawnow
end
hold off
pause;
plot(t,Fz,'y-'), title (' Гравитация Земли ' ), grid on;
print -dwin;
pause;
plot(t,Fs,'y-'), title (' Гравитация Солнца и солнечное давление '), grid on;
print -dwin;
pause;
plot(t,Fl,'y-'), title (' Гравитация Луны '), grid on;
print -dwin;
pause;
plot(t,Fa,'y-'), title (' Сопротивление атмосферы '), grid on;
print -dwin;
pause;
plot(t,U20,'y-'), title (' Нецентральность гравитационного поля Земли '), grid on;
print -dwin;
pause;
plot(t,Fz+Fs+Fl+Fa+U20,'y-'), title (' Суммарное возмущающее ускорение '), grid on;
print -dwin;
pause;
clear all
clc
g_r = pi/180;
r_g = 180/pi;
p_n = 6952137.;
e_n = 0;
a_n = 6952137.;
Om_n0 = 28.1*g_r;
i_n = 97.6*g_r;
omg_n = 346.725*g_r;
T_n = 5765;
load u_par.dat
t_u = u_par(:,1);
p_u = u_par(:,2);
e_u = u_par(:,3);
a_u = u_par(:,4);
Om_u = u_par(:,5);
i_u = u_par(:,6);
omg_u = u_par(:,7);
T_u = u_par(:,8);
u_u = u_par(:,9);
clear u_par;
load u_f.dat;
Fz_u = u_f(:,2);
Fs_u = u_f(:,3);
Fl_u = u_f(:,4);
Fa_u = u_f(:,5);
U20_u = u_f(:,6);
clear u_f;
s_tmp = size(t_u);
s_m_u = s_tmp(1);
clear s_tmp;
ws = 2*pi/(365.2422*24*3600);
for j = 1:s_m_u, tmp(j) = Om_n0+ws*t_u(j);
end
Om_n_u = tmp';
clear tmp;
plot(t_u,p_u,'y-',[min(t_u) max(t_u)],[p_n p_n],'r-'), title (' Фокальный параметр '), grid on;
print -dwin;
pause;
plot(t_u,p_u-p_n,'y-'), title (' dp '), grid on;
print -dwin;
pause;
plot(t_u,e_u,'y-',[min(t_u) max(t_u)],[e_n e_n],'r-'), title (' Эксцентриситет '), grid on;
print -dwin;
pause;
plot(t_u,e_u-e_n,'y-'), title (' de '), grid on;
print -dwin;
pause;
plot(t_u,a_u,'y-',[min(t_u) max(t_u)],[a_n a_n],'r-'), title (' Большая полуось орбиты '), grid on;
print -dwin;
pause;
plot(t_u,a_u-a_n,'y-'), title (' da '), grid on;
print -dwin;
pause;
plot(t_u,Om_u*r_g,'y-',t_u,Om_n_u*r_g,'r-'), title (' Долгота восходящего узла '), grid on;
print -dwin;
pause;
plot(t_u,Om_u*r_g-Om_n_u*r_g,'y-'), title (' dOm '), grid on;
print -dwin;
pause;
plot(t_u,i_u*r_g,'y-',[min(t_u) max(t_u)],[i_n*r_g i_n*r_g],'r-'), title (' Наклонение '), grid on;
print -dwin;
pause;
plot(t_u,i_u*r_g-i_n*r_g,'y-'), title (' di '), grid on;
print -dwin;
pause;
plot(t_u,T_u,'y-',[min(t_u) max(t_u)],[T_n T_n], 'r-'), title (' Период '), grid on;
print -dwin;
pause;
plot(t_u,T_u-T_n,'y-'), title (' dT '), grid on;
print -dwin;
pause;
clear all
Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14
