Исследование движения центра масс межпланетных космических аппаратов
Дипломная работа - Физика
Другие дипломы по предмету Физика
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
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
#include
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
typedef long double real;
extern const float g_r;
extern const float r_g;
extern int Np;
<