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++)