: Исследование движения центра масс межпланетных космических аппаратов

zs = A*sin(Tet_s)*sin(e_0);

real Tet_l = 0+w_l*t;

real Om_l = 0-ww_l*t;

real i_l = acos(cos(e_0)*cos(5.15*g_r)-sin(e_0)*sin(5.15*g_r)*cos(Om_l));

real rsr_l = 3.8448e8;

xl = rsr_l*(cos(Tet_l)*cos(Om_l)-cos(i_l)*sin(Tet_l)*sin(Om_l));

yl = rsr_l*(cos(Tet_l)*sin(Om_l)+cos(i_l)*sin(Tet_l)*cos(Om_l));

zl = rsr_l*sin(i_l)*sin(Tet_l);

real R_ka = sqrt(x*x+y*y+z*z);

real Fz_x = -mu_z*x/pow(R_ka,3.);

real Fz_y = -mu_z*y/pow(R_ka,3.);

real Fz_z = -mu_z*z/pow(R_ka,3.);

real mu_sd = p*sm_s*A*A/m;

real R_s = sqrt((x-xs)*(x-xs)+(y-ys)*(y-ys)+(z-zs)*(z-zs));

real Fs_x = -(mu_s-mu_sd)*x/pow(R_s,3.);

real Fs_y = -(mu_s-mu_sd)*y/pow(R_s,3.);

real Fs_z = -(mu_s-mu_sd)*z/pow(R_s,3.);

real R_l = sqrt((x-xl)*(x-xl)+(y-yl)*(y-yl)+(z-zl)*(z-zl));

real Fl_x = -mu_l*x/pow(R_l,3.);

real Fl_y = -mu_l*y/pow(R_l,3.);

real Fl_z = -mu_l*z/pow(R_l,3.);

real V_ka = sqrt(Vx*Vx+Vy*Vy+Vz*Vz);

real Fa_x = (-Cx*sm_a/(2*m))*ro*V_ka*Vx;

real Fa_y = (-Cx*sm_a/(2*m))*ro*V_ka*Vy;

real Fa_z = (-Cx*sm_a/(2*m))*ro*V_ka*Vz;

const real c20 = -1.09808e-3;

const real c22 = 5.74e-6;

const real d22 = -1.58e-6;

const real r_e = 6378137.;

real cr = mu_z*r_e*r_e/pow(R_ka,5);

real lr = 2*atan(y/x);

real mr = 3*(c22*cos(lr)+d22*sin(lr));

real U20_x = cr*x*(c20*(1.5-7.5*z*z/pow(R_ka,2))+mr*(5*z*z/pow(R_ka,2)-3));

real U20_y = cr*y*(c20*(1.5-7.5*z*z/pow(R_ka,2))+mr*(5*z*z/pow(R_ka,2)-3));

real U20_z = cr*z*(c20*(4.5-7.5*z*z/pow(R_ka,2))+5*mr*(z*z/pow(R_ka,2)-1));

dery[0] = Vx;

dery[1] = Vy;

dery[2] = Vz;

dery[3] = (Fz_x+U20_x+Fs_x+Fl_x+Fa_x+akor[0]);

dery[4] = (Fz_y+U20_y+Fs_y+Fl_y+Fa_y+akor[1]);

dery[5] = (Fz_z+U20_z+Fs_z+Fl_z+Fa_z+akor[2]);

Fz = sqrt(Fz_x*Fz_x+Fz_y*Fz_y+Fz_z*Fz_z);

Fs = sqrt(Fs_x*Fs_x+Fs_y*Fs_y+Fs_z*Fs_z);

Fl = sqrt(Fl_x*Fl_x+Fl_y*Fl_y+Fl_z*Fl_z);

Fa = sqrt(Fa_x*Fa_x+Fa_y*Fa_y+Fa_z*Fa_z);

U20 = sqrt(U20_x*U20_x+U20_y*U20_y+U20_z*U20_z);

parn[3] = parn[3]+w_s*t;

par_or(f,par);

korr(t,f,dery);

if ((u_last-par[7]) > 300*g_r)

Fl_u = 1;

u_last = par[7];

}

void korr(real& t, real *f, real *)

{

if (t > (Tkor+172800.))

{

if ((fabs(dl) > 0.1*g_r) && (!Fl_ka) && (!Fl_kp) && (!Fl_ki))

{

Fl_kp = 1;

Fl_ka = 0;

Fl_ki = 0;

cout << "Результат измерений накоплен" << '\n';

cout << "Необходима коррекция периода. dl=" << dl*r_g << "град." << '\n';

cout << "Период ном.=" << parn [6] << "Период тек.=" << par[6] << '\n';

cout << "Параметры орбиты" << '\n';

cout << " Rp = " << par[2]*(1-par[1]) << '\n';

cout << " Ra = " << par[2]*(1+par[1]) << '\n';

cout << " p = " << par[0] << '\n';

cout << " a = " << par[2] << " e = " << par[1] << "\n T = "

<< par[6] << " w = " << par[5]*r_g << " u = " << par[7]*r_g

<< '\n';

clrscr();

}

}

Fl_a = 0;

Fl_p = 0;

Fl_lu = 0;

real da;

if (par[5] > par[7])

da = fabs(par[5]-par[7]-M_PI);

else

da = fabs(par[5]-par[7]+M_PI);

if (da < .1*g_r)

{

Fl_a = 1;

}

if (fabs(par[5] - par[7]) < .1*g_r)

{

Fl_p = 1;

}

if (par[7] < .1*g_r )

{

Fl_lu = 1;

}

real Vk;

if (T_vd)

if (t >= (T_vd +20))

{

T_vd = 0;

akor[0] = 0;

akor[1] = 0;

akor[2] = 0;

cout << "Выкл.дв. \n t = " << t;

}

if (((Fl_kp && Fl_a) || (Fl_ka && Fl_p) || (Fl_ki && Fl_lu)) && (!T_vd))

{

cout << " \n Коррекция \n";

cout << "\n Начало t=" << t << "сек \n";

int sim;

if ((t-Tkor) < 2500)

{

cout << "Не корректировать!";

return;

}

Tkor = t;

real R_t = sqrt(f[0]*f[0]+f[1]*f[1]+f[2]*f[2]);

real V_t = sqrt(f[3]*f[3]+f[4]*f[4]+f[5]*f[5]);

real R_n = parn[0];

if (Fl_a)

{

dRa = R_t-R_n;

dRp = par[2]*(1-par[1])-R_n;

cout << "Апоцентр dRp:" << dRp << "м \n";

cout << "dRa:" << dRa << "м \n";

cout << "w=" << par[5]*r_g << "u=" << par[7]*r_g << '\n';

real l,ln;

l = -(w_z-w_s)*par[6];

ln = -(w_z-w_s)*parn[6];

dl = -(w_z-w_s)*(par[6]-parn[6]);

cout << "T=" << par[6] << "Тном=" << parn[6] << " T-Tном="

<< par[6]-parn[6] << '\n' << "l=" << l*r_g << "lном="

<< ln*r_g << "l-lном=" << (l-ln)*r_g << "dl=" << dl

<< '\n';

if (dRp > 0)

Sig_a = -1;

else

Sig_a = 1;

cout << "Знак ускорения:" << Sig_a << '\n';

clrscr();

real Rp = par[2]*(1-par[1]);

real Ra_p = par[2]*(1+par[1]);

real Rp_p2 = Rp;

real Ra_p2 = R_t;

cout << "Rp=" << Rp_p2 << "Ra=" << Ra_p2 << '\n';

cout << "Ra_p=" << Ra_p << "\n Rt=" << R_t << '\n';

if (fabs(Rp - R_n) < 500)

{

Fl_kp = 0;

Fl_ka = 1;

cout << "Закончить коррекцию в апоцентре \n" << "dRp=" << Rp-R_n

<< "dRa=" << dRa << "t=" << t << '\n';

cout << "Параметры орбиты: \n" << "Rp=" << par[2]*(1-par[1])

<< "Ra=" << par[2]*(1+par[1]) << "\n p=" << par[0]

<< "a=" << par[2] << "e=" << par[1] << "\n T="

<< par[6] << "w=" << par[5]*r_g << "u=" << par[7]*r_g

<< '\n';

cout << "Суммарный импульс для коррекции перицентра=" << dV_ps << '\n';

clrscr();

}

else

{

if (R_t > 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 << "\n dVтреб=" << dV << "dVmax за 20 сек=" << dVmax;

if (fabs(dV) > 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 << "\n dV=" << dV << "dVmax=" << dVmax;

cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

<< '\t' << akor[2] << '\t' <<

sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';

dV_ps = dV_ps+dVmax;

cout << "Суммарный импульс=" << dV_ps << '\n';

}

else

{

akor[0] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[3]/V_t;

akor[1] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[4]/V_t;

akor[2] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[5]/V_t;

cout << "\n dV=" << dV << "dVmax=" << dVmax;

cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

<< '\t' << akor[2] << '\t' <<

sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';

dV_ps = dV_ps+fabs(dV);

cout << "Суммарный импульс=" << dV_ps << '\n';

}

if (dVmax > 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 << "Параметры орбиты: \n" << " Rp_r = " << Rp_r

<< " Ra_r = " << Ra_r << "\n p_r = " << p_r << " a_r = "

<< a_r << " e_r = " << e_r << '\n';

}

else

{

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 << "Параметры орбиты: \n" << " Rp_r = " << Rp_r

<< " Ra_r = " << Ra_r << "\n p_r = " << p_r << " a_r = "

<< a_r << " e_r = " << e_r << '\n';

}

T_vd = t;

cout << "Вкл.дв. t=" << T_vd << '\n';

}

}

if (Fl_p)

{

dRp = R_t-R_n;

dRa = par[2]*(1+par[1])-R_n;

cout << " Перицентра - dRp:" << dRp << "м \n";

cout << "dRa:" << dRa << "м. \n";

cout << "w=" << par[5]*r_g << "u=" << par[7]*r_g << '\n';

real l,ln;

l = -(w_z-w_s)*par[6];

ln = -(w_z-w_s)*parn[6];

dl = -(w_z-w_s)*(par[6]-parn[6]);

cout << "T=" << par[6] << "Tном=" << parn[6] << "T-Tном="

<< par[6]-parn[6] << '\n' << "l=" << l*r_g << "lном="

<< ln*r_g << "l-lном=" << (l-ln)*r_g << "dl=" << dl << '\n';

if (dRa > 0)

Sig_a = -1;

else

Sig_a = 1;

cout << "Знак ускорения:" << Sig_a << '\n';

clrscr();

real Ra = par[2]*(1+par[1]);

real Rp_p1 = R_t;

real Ra_p1 = Ra;

cout << "Rp=" << Rp_p1 << "Ra=" << Ra_p1 << '\n';

if ((fabs(Ra-R_n) < 500) || (fabs(dl*r_g) < .0001))

{

cout << "Закончить коррекцию в перицентре \n" << "dRa="

<< Ra-R_n << "dRp=" << dRp << "t=" << t << '\n';

cout << "Параметры орбиты: \n " << "Rp="

<< par[2]*(1-par[1]) << "Ra=" << par[2]*(1+par[1])

<< " \n p=" << par[0] << "a=" << par[2] << "e="

<< par[1] << " \n T=" << par[6] << "w=" << par[5]*r_g

<< "u=" << par[7]*r_g << '\n';

cout << "Суммарный импульс для коррекции перицентра=" << dV_as << '\n';

clrscr();

Fl_ka = 0;

Fl_ki = 1;

}

else

{

if (R_t > 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 << "\n dVнадо=" << dV << " dVmax за 20 сек=" << dVmax;

if (fabs(dV) > 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 << "\n dV=" << dV << "dVmax=" << dVmax;

cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

<< '\t' << akor[2] << '\t' <<

sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';

dV_as = dV_as+dVmax;

cout << "Суммарный импульс=" << dV_as << '\n';

}

else

{

akor[0] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[3]/V_t;

akor[1] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[4]/V_t;

akor[2] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[5]/V_t;

cout << "\n dV=" << dV << " dVmax=" << dVmax;

cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

<< '\t' << akor[2] << '\t' <<

sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';

dV_as = dV_as+fabs(dV);

cout << "Суммарный импульс=" << dV_as << '\n';

}

if (dVmax > 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 << "Параметры орбиты: \n" << "Rp_r=" << Rp_r

<< "Ra_r=" << Ra_r << "\n p_r=" << p_r << "a_r="

<< a_r << "e_r=" << e_r << '\n';

}

else

{

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 << "Параметры орбиты: \n" << "Rp_r=" << Rp_r

<< "Ra_r=" << Ra_r << "\n p_r=" << p_r << "a_r="

<< a_r << "e_r=" << e_r << '\n';

}

T_vd = t;

cout << "Вкл.дв. t=" << T_vd << '\n';

}

}

if (Fl_lu)

{

real di = par[4]-parn[4];

cout << "Линия узлов - di: " << di*r_g << "град \n";

cout << "w=" << par[5]*r_g << "u=" << par[7]*r_g << '\n';

real l,ln;

l = -(w_z-w_s)*par[6];

ln = -(w_z-w_s)*parn[6];

dl = -(w_z-w_s)*(par[6]-parn[6]);

cout << "T=" << par[6] << "Tном=" << parn[6] << "T-Tном="

<< par[6]-parn[6] << '\n' << "l=" << l*r_g << "lном="

<< ln*r_g << "l-lном=" << (l-ln)*r_g << "dl=" << dl

<< "\n i=" << par[4]*r_g << "iном=" << parn[4]*r_g << '\n';

cout << "Параметры орбиты: \n " << "Rp="

<< par[2]*(1-par[1]) << "Ra=" << par[2]*(1+par[1])

<< " \n p=" << par[0] << "a=" << par[2] << "e="

<< par[1] << " \n T=" << par[6] << "w=" << par[5]*r_g

<< "u=" << par[7]*r_g << " \n i=" << par[4]*r_g << '\n';

clrscr();

real Vk_x,Vk_y,Vk_z;

if (fabs(di) < .0001*g_r)

{

Fl_ki = 0;

cout << "Закончить коррекцию наклонения \n "

<< "di=" << (par[4]-parn[4])*r_g << "t=" << t << '\n';

cout << "Параметры орбиты: \n " << "Rp="

<< par[2]*(1-par[1]) << "Ra=" << par[2]*(1+par[1])

<< " \n p=" << par[0] << "a=" << par[2] << "e="

<< par[1] << " \n T=" << par[6] << "w=" << par[5]*r_g

<< "u=" << par[7]*r_g << " \n i=" << par[4]*r_g << '\n';

cout << "Суммарный импульс=" << dV_is

<< '\n';

clrscr();

}

else

{

real teta;

if (par[7] > 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 << "Vнач=" << V_t << "Vк=" << Vk << "teta=" << teta*r_g

<< '\n';

cout << "dV=" << dV << "dVmax за 20 сек=" << dVmax;

if (dV > dVmax)

{

akor[0] = (25./m)*dV_x/dV;

akor[1] = (25./m)*dV_y/dV;

akor[2] = (25./m)*dV_z/dV;

cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1] <<

'\t' << akor[2] << '\t' <<

sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';

dV_is = dV_is+dVmax;

cout << "Суммарный импульс=" << dV_is << '\n';

}

else

{

akor[0] = (fabs(dV)/dVmax)*(25./m)*dV_x/dV;

akor[1] = (fabs(dV)/dVmax)*(25./m)*dV_y/dV;

akor[2] = (fabs(dV)/dVmax)*(25./m)*dV_z/dV;

cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

<< '\t' << akor[2] << '\t'<<

sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';

dV_is = dV_is+fabs(dV);

cout << "Суммарный импульс=" << dV_is << '\n';

}

T_vd = t;

cout << "Вкл.дв. t=" << T_vd << '\n';

}

}

if ((!Fl_ka) && (!Fl_kp) && (!Fl_ki))

{

cout << "Коррекция окончена!" << '\n';

real m_t;

dV_ss = dV_ps+dV_as+dV_is;

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