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

    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, 7, 8, 9, 10, 11, 12, 13, 14



Реклама
В соцсетях
рефераты скачать рефераты скачать рефераты скачать рефераты скачать рефераты скачать рефераты скачать рефераты скачать