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

Создание математической модели, позволяющей учесть эволюцию орбиты МКА. Исследование движения центра масс МКА под действием различных возмущающих ускорений. Разработка алгоритма коррекции, ликвидирующего ошибки выведения МКА, расчет массы топлива.

Рубрика Астрономия и космонавтика
Вид дипломная работа
Язык русский
Дата добавления 07.09.2010
Размер файла 245,7 K

Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже

Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.

  • 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;
  • }
  • 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
    • 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
    • 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. Файл 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
    • 7. Программа построения временных диаграмм
    • 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

    Подобные документы

    • Разработка ракет с широким применением унифицированных базовых конструкций и доступной элементной базой. Тактико-технические характеристики ракет-носителей "Виктория-К", "Волна", "Единство". Описание двигателей, определение центра масс в процессе полета.

      курсовая работа [2,2 M], добавлен 11.12.2014

    • Особенности и основные способы проектирования электрореактивной двигательной установки космического аппарата. Этапы разработки циклограммы энергопотребления, анализ чертежа движителя. Характеристика космических электроракетных двигательных установок.

      дипломная работа [496,1 K], добавлен 18.12.2012

    • Ограниченная круговая задача трех тел и уравнения движения. Типы ограниченных орбит в окрестности точек либрации и гравитационная задача. Затенённость орбит и моделирование движения космического аппарата. Проекция долгопериодической орбиты на плоскость.

      курсовая работа [3,6 M], добавлен 01.07.2017

    • Анализ баллистических характеристик космического аппарата. Расчет масс служебных систем, элементов топлива. Зона обзора на поверхности Земли и полоса обзора. Изучение системы электроснабжения, обеспечения теплового режима, бортового комплекса управления.

      курсовая работа [53,7 K], добавлен 10.07.2012

    • Характер и обоснование движения тел солнечной системы. Элементы эллиптической орбиты и их назначение. Особенности движения Земли и Луны. Феномен солнечного затмения, причины и условия его наступления. Специфика лунных затмений и их влияние на Землю.

      курсовая работа [4,0 M], добавлен 27.06.2010

    • История проблемы выхода на орбиту. Расчет возможности вывода тела на орбиту одним толчком. Признаки тела переменной массы. Моделирование обстоятельств наблюдения искусственных спутников земли. Математическое моделирование движения ракеты-носителя.

      реферат [120,6 K], добавлен 14.10.2015

    • Скорость вращения галактики как скорость вращения различных компонентов галактики вокруг её центра. Особенности движения газа и звёзд. Распределение звезд, анализ их поля скоростей как информация о движении в галактике, оценка вероятности столкновения.

      статья [34,3 K], добавлен 01.10.2010

    • Общая характеристика и направления деятельности организации. Общие сведения об энергоснабжении космических аппаратов, особенности использования солнечных батарей. Химические источники тока. Выбор параметров солнечных батарей и буферных накопителей.

      отчет по практике [195,1 K], добавлен 16.04.2016

    • Уравнения движения системы в инерциальной и неинерциальной системе отсчета. Оценка области местонахождения планет земного типа в тройной системе тел. Исследование устойчивости точек либрации. Группировка космических станций в окололунном пространстве.

      дипломная работа [1,3 M], добавлен 11.02.2013

    • Описание кометы как тела Солнечной системы, особенности ее строения. Траектория и характер движения этого космического объекта. История наблюдения астрономами движения кометы Галлея. Наиболее известные периодические кометы и специфика их орбиты.

      презентация [3,8 M], добавлен 20.05.2015

    Работы в архивах красиво оформлены согласно требованиям ВУЗов и содержат рисунки, диаграммы, формулы и т.д.
    PPT, PPTX и PDF-файлы представлены только в архивах.
    Рекомендуем скачать работу.