{

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;