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

Страница 14

ofstream m_f ("m_f.dat");

ofstream m_s ("m_s.dat");

ofstream m_l ("m_l.dat");

ofstream m_par ("m_par.dat");

ofstream u_f ("u_f.dat");

ofstream u_par ("u_par.dat");

ofstream k_par ("k_par.dat");

void out_p(real,real *,real*,int,int,real*);

void main()

{

clrscr();

init_m();

real dery[]={ .167, .167, .167, .167, .166, .166};

int ihlf;

int ndim = 6;

Drkgs(prmt,y_main,dery,ndim,ihlf,fct,out_p);

clrscr();

if (ihlf<11)

{

cout << '\n' << "Успешное завершение моделирования" << '\n';

cout << " t0 = " << t_beg << " tk = " << t_end << " dt = " << dt;

cout << '\n' << "Число делений шага=" << ihlf;

}

else

{

cout << '\n' << "Ненормальное завершение моделирования" << '\n';

cout << " t0 = " << t_beg << " tk = " << t_end << " dt = " << dt;

cout << '\n' << "Число делений шага=" << ihlf;

}

getch();

m_y.close();

m_f.close();

m_s.close();

m_l.close();

m_par.close();

u_f.close();

u_par.close();

k_par.close();

}

void out_p(real x,real *y,real*,int,int,real*)

{

if (x >= (dTp*Curp))

{

Curp++;

gotoxy(1,20);

cout << "Процесс выполнения:" << float(Curp)*100./Np << " % " << '\n';

cout.precision(7);

m_y << x << '\t' << y[0] << '\t' << y[1] << '\t' << y[2] << '\t'

<< y[3] << '\t' << y[4] << '\t' << y[5] << '\n';

m_f << x << '\t' << Fz << '\t' << Fs << '\t' << Fl << '\t' << Fa

<< '\t' << U20 << '\n';

m_s << x << '\t' << xs << '\t' << ys << '\t' << zs << '\n';

m_l << x << '\t' << xl << '\t' << yl << '\t' << zl << '\n';

m_par << x << '\t' << par[0] << '\t' << par[1] << '\t' << par[2]

<< '\t' << par[3] << '\t' << par[4] << '\t' << par[5]

<< '\t' << par[6] << '\t' << par[7] << '\n';

}

if (Fl_u && (par[7] > parn[7]))

{

Fl_u = 0;

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

u_par << x << '\t' << par[0] << '\t' << par[1] << '\t' << par[2]

<< '\t' << par[3] << '\t' << par[4] << '\t' << par[5]

<< '\t' << par[6] << '\t' << par[7] << '\n';

u_f << x << '\t' << Fz << '\t' << Fs << '\t' << Fl

<< '\t' << Fa << '\t' << U20 << '\n';

}

if ((x > 79000) && (x < 81000))

{

k_par << x << '\t' << par[5] << '\t' << par[7] << '\n';

}

}

6.2. ПОДПРОГРАММА РАСЧЕТА ВОЗМУЩАЮЩИХ УСКОРЕНИЙ, ПАРАМЕТРОВ ОРБИТЫ И КОРРЕКЦИИ SFUN.CPP

#include "sfun.h"

const real p = 4.64e-6;

const real sm_s = 8.;

const real A = 1.496e11;

const real Cx = 2.;

const real sm_a = 2.5;

const real ro = 5.098e-13;

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

void fct(real& t, real *f, real *dery)

{

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 Tet_s = (28.1+60*g_r)+w_s*t;

real e_0 = 23.45*g_r;

xs = A*cos(Tet_s);

ys = A*sin(Tet_s)*cos(e_0);

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";