akor[2] = 0;
cout 0)
Sig_a = -1;
else
Sig_a = 1;
cout R_n)
{
Rp_p = R_n;
Ra_p = R_t;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1-e_p);
}
else
{
Rp_p = R_t;
Ra_p = R_n;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1+e_p);
}
real dV = Vk-V_t;
real dVmax = 20*25./m;
cout dVmax)
{
akor[0] = Sig_a*(25./m)*f[3]/V_t;
akor[1] = Sig_a*(25./m)*f[4]/V_t;
akor[2] = Sig_a*(25./m)*f[5]/V_t;
cout fabs(dV))
{
dVmax = fabs(dV);
real Vk_r = Sig_a*dVmax+V_t;
real Ra_r = R_t;
real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;
real a_r = Ra_r/(1+e_r);
real p_r = a_r*(1-e_r*e_r);
real Rp_r = a_r*(1-e_r);
cout 0)
Sig_a = -1;
else
Sig_a = 1;
cout R_n)
{
Rp_p = R_n;
Ra_p = R_t;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1-e_p);
}
else
{
Rp_p = R_t;
Ra_p = R_n;
a_p = (Ra_p+Rp_p)/2.;
e_p = 1-Rp_p/a_p;
p_p = a_p*(1-e_p*e_p);
Vk = sqrt(mu_z/p_p)*(1+e_p);
}
real dV = Vk-V_t;
real dVmax = 20*25./m;
cout dVmax)
{
akor[0] = Sig_a*(25./m)*f[3]/V_t;
akor[1] = Sig_a*(25./m)*f[4]/V_t;
akor[2] = Sig_a*(25./m)*f[5]/V_t;
cout fabs(dV))
{
dVmax = fabs(dV);
real Vk_r = Sig_a*dVmax+V_t;
real Ra_r = R_t;
real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;
real a_r = Ra_r/(1+e_r);
real p_r = a_r*(1-e_r*e_r);
real Rp_r = a_r*(1-e_r);
cout par[5])
teta = 2*M_PI+par[7]-par[5];
else
teta = par[7]-par[5];
real Vr_i = sqrt(mu_z/par[0])*par[1]*sin(teta);
real Vn_i = sqrt(mu_z/par[0])*(1+par[1]*cos(teta));
V_t = sqrt(f[3]*f[3]+f[4]*f[4]+f[5]*f[5]);
Vk_x = -Vn_i*cos(parn[4])*sin(par[3])+Vr_i*cos(par[3]);
Vk_y = Vn_i*cos(parn[4])*cos(par[3])+Vr_i*sin(par[3]);
Vk_z = Vn_i*sin(parn[4]);
Vk = sqrt(Vk_x*Vk_x+Vk_y*Vk_y+Vk_z*Vk_z);
real dV_x = Vk_x-f[3];
real dV_y = Vk_y-f[4];
real dV_z = Vk_z-f[5];
real dV = sqrt(dV_x*dV_x+dV_y*dV_y+dV_z*dV_z);
real dVmax = 20*25./m;
cout dVmax)
{
akor[0] = (25./m)*dV_x/dV;
akor[1] = (25./m)*dV_y/dV;
akor[2] = (25./m)*dV_z/dV;
cout = 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 0)
par[7] = acos(c_u);
else
par[7] = 2*M_PI - acos(c_u);
}
#include "rk5.h"
#include
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= 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 0.0)
{
if (ihlf-10 >= 0)
{
ihlf = 11;
fct(x,y,dery);
goto l39;
}
for (i=0; i 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
#include
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
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
#include
#include
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
#include
#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
|