Rizzolo
27-11-2009, 17:30
C'è qualcuno disposto a compilarmi questo script?
Ho provato con free pascal ma non ne ho tirato fuori nulla.
Si tratta del propagatore SDP8 (magari avercelo in matlab).
Lo script potete scaricarlo da qui: http://canasta.altervista.org/sgp8sdp8.pas
unit SGP8SDP8;
{ Author: Dominik Brodowski }
{ Current Vision: 2001 Mar 03 }
{ Version: 1.00 }
{ uses the SGP4-PL2A TPU package by Dr TS Kelso }
{$N+}
INTERFACE
Uses Support,
SGP_Init,
SGP_Math,SGP_Time;
Procedure SGP8call(time : double; var pos,vel : vector);
Procedure SGP8(tsince : double; var iflag : integer; var pos,vel : vector);
Procedure SDP8(tsince : double; var iflag : integer; var pos,vel : vector);
IMPLEMENTATION
Uses SGP_Intf;
var
{dpinit}
eqsq,siniq,cosiq,rteqsq,ao,cosq2,sinomo,cosomo,
bsq,xlldot,omgdt,xnodot,xnodp : double;
{dpsec/dpper}
xll,omgasm,xnodes,_em,xinc,xn,t : double;
qoms2t : double;
Procedure Define_Derived_Constants;
begin
xke := Sqrt(3600*ge/Cube(xkmper)); {Sqrt(ge) ER^3/min^2}
qoms2t := Sqr(Sqr(qo-s)); {(qo-s)^4 ER^4}
end; {Define_Derived_Constants}
Procedure SGP8(tsince : double; var iflag : integer; var pos,vel : vector);
label
EndKepler;
const
a1 : double = 0; cosio : double = 0; sinio : double = 0;
delta1 : double = 0; a0 : double = 0; delta0 : double = 0;
B : double = 0; noii : double = 0; aoii : double = 0;
rho : double = (2.461*0.00001*XKMPER); beta : double = 0;
theta : double = 0; M1dot : double = 0; omega1dot: double= 0;
xnode1dot: double = 0; M2dot : double = 0; omega2dot: double= 0;
xnode2dot: double = 0; ldot : double = 0; omegadot: double = 0;
xnodedot : double = 0; tsi: double = 0; temp : double = 0;
eta : double = 0; phi : double = 0; alpha : double = 0;
C0 : double = 0; C1 : double = 0; D1 : double = 0;
D2 : double = 0; D3 : double = 0; D4 : double = 0;
D5 : double = 0; B1 : double = 0; B2 : double = 0;
B3 : double = 0; A30 : double = 0; C2 : double = 0;
C3 : double = 0; n0dot : double = 0; C4 : double = 0;
C5 : double = 0; D6 : double = 0; D7 : double = 0;
D8 : double = 0; e0dot : double = 0; etadot : double = 0;
alphadottoalpha : double = 0; C0dottoC0:double = 0;
tsidottotsi : double = 0; C1dottoC1:double = 0;
phidottophi : double = 0; C6 :double = 0;
D9 : double = 0; D10 : double = 0; D11 : double = 0;
D12 : double = 0; D13 : double = 0; D14 : double = 0;
D15 : double = 0; D1dot : double = 0; D2dot : double = 0;
D3dot : double = 0; D4dot : double = 0; D5dot : double = 0;
C2dot : double = 0; C3dot : double = 0; D16 : double = 0;
n0ddot : double = 0; e0ddot : double = 0; D17 : double = 0;
tsiddottotsi : double = 0; etaddot : double = 0;
D18 : double = 0; D19 : double = 0; D1ddot : double = 0;
n0tdot : double = 0; p : double = 0; gamma : double = 0;
nd : double = 0; q : double = 0; ed : double = 0;
n : double = 0; edot : double = 0; e : double = 0;
Z1 : double = 0; omega : double = 0; xnode : double = 0;
M : double = 0; Ew : double = 0; deltaEw : double = 0;
a : double = 0; sinf : double = 0; cosf : double = 0;
u : double = 0; sino : double = 0; coso : double = 0;
rii : double = 0; rdotii : double = 0; rfdotii : double = 0;
deltar : double = 0; deltardot:double= 0; deltau : double = 0;
deltaLambda:double= 0; r : double = 0; rdot : double = 0;
rfdot : double = 0; y4 : double = 0; y5 : double = 0;
lambda : double = 0; deltaI : double = 0; deltarfdot:double= 0;
f : double = 0; cosI2 : double = 0; sinu : double = 0;
cosu : double = 0; smallN0DDoT : double=0;
var
i : integer;
UV,VV : vector;
begin;
{constants I}
A1:=power((xke/xno),tothrd);
cosio := cos(xincl);
sinio := sin(xincl);
delta1 := 1.5*CK2*(3*cosio*cosio-1)/(sqr(a1)*power((1-sqr(eo)),1.5));
a0 := a1*(1-(1/3)*delta1-sqr(delta1)-(134/81)*sqr(delta1)*delta1);
delta0 := 1.5*CK2*(3*cosio*cosio-1)/(sqr(a0)*power((1-sqr(eo)),1.5));
noii := xno/(1+delta0);
aoii := a0/(1-delta0);
{B term}
B := 2*BStar/rho;
{constants II}
beta:=sqrt(1-sqr(eo));
theta:=cosio;
temp:=(noii*CK2)/(sqr(aoii)*sqr(beta)*beta);
M1dot:=-1.5*temp*(1-3*sqr(theta));
omega1dot:=-1.5*(temp/beta)*(1-5*sqr(theta));
xnode1dot:=-3*(temp/beta)*theta;
M2dot:=(3/16)*(temp*CK2/(sqr(aoii)*power(beta,4)))*(13-78*sqr(theta)+
137*power(theta,4));
omega2dot:=(3/16)*(temp*CK2/(sqr(aoii)*power(beta,5)))*(7-114*sqr(theta)+
395*power(theta,4))+1.25*(noii*CK4)/(power(aoii,4)*power(beta,8))*
(3-36*sqr(theta)+49*power(theta,4));
xnode2dot:=1.5*(temp*CK2/(sqr(aoii)*power(beta,5)))*theta*(4-19*sqr(theta))+
2.5*(noii*CK4)/(power(aoii,4)*power(beta,8))*theta*(3-7*sqr(theta));
ldot := noii + M1dot + M2dot;
omegadot := omega1dot+omega2dot;
xnodedot := xnode1dot + xnode2dot;
tsi := 1/(aoii*sqr(beta)-s);
eta := eo*s*tsi;
phi := sqrt(1-sqr(eta));
alpha := sqrt(1+sqr(eo));
C0:=0.5*B*rho*qoms2t*noii*aoii*power(tsi,4)/(alpha*power(phi,7));
C1:=1.5*noii*power(alpha,4)*C0;
D1:=(tsi*(1/sqr(phi)))/(aoii*sqr(beta));
D2:=12+36*sqr(eta)+4.5*power(eta,4);
D3:=15*sqr(eta)+2.5*power(eta,4);
D4:=5*eta+(15/4)*power(eta,3);
D5:=tsi/sqr(phi);
B1:=-CK2*(1-3*sqr(theta));
B2:=-CK2*(1-sqr(theta));
A30:=-XJ3*power(Ae,3)/CK2;
B3:=A30*sinio;
C2:=D1*D3*B2;
C3:=D4*D5*B3;
n0dot:=C1*((2+sqr(ETA)*(3+34*sqr(EO))+5*(Eo)*(ETA)*(4+sqr(ETA))+8.5*sqr(EO))+D1*D2*B1+ C2*COS(2*omegao)+C3*SIN(omegao));
D6:=30*eta+22.5*power(eta,3);
D7:=5*Eta+12.5*power(eta,3);
D8:=1+6.75*sqr(eta)+power(eta,4);
C4 := D1*D7*B2;
C5:=D5*D8*B3;
e0dot:=-C0*(ETA*(4+sqr(ETA)+sqr(EO)*(15.5+7*sqr(ETA)))+EO*(5+15*sqr(ETA))+D1*D6*B1+C4*COS(2*omegao)+C5*SIN(omegao));
alphadottoalpha:=eo*e0dot/(sqr(alpha));
C6:=(1/3)*(n0dot/noii);
tsidottotsi:=2*aoii*tsi*(C6*sqr(beta)+eo*e0dot);
etadot:=(e0dot+eo*tsidottotsi)*s*tsi;
phidottophi:=-eta*etadot/sqr(phi);
C0dottoC0:=C6+4*tsidottotsi-alphadottoalpha-7*phidottophi;
C1dottoC1:=n0dot/noii+4*alphadottoalpha+C0dottoC0;
D9:=6*eta+20*eo+15*eo*sqr(eta)+68*sqr(eo)*eta;
D10:=20*eta+5*power(eta,3)+17*eo+68*(eo)*sqr(eta);
D11:=eta*(72+18*sqr(eta));
D12:=eta*(30+10*sqr(eta));
D13:=5+11.25*sqr(eta);
D14:=tsidottotsi-2*phidottophi;
D15:=2*(C6+eo*e0dot/sqr(beta));
D1dot:=D1*(D14+D15);
D2dot:=etadot*D11;
D3dot:=etadot*D12;
D4dot:=etadot*D13;
D5dot:=D5*D14;
C2dot:=B2*(D1dot*D3+D1*D3dot);
C3dot:=B3*(D5dot*D4+D5*D4dot);
D16:=D9*etadot+D10*e0dot+B1*(D1dot*D2+D1*D2dot)+C2dot*cos(2*omegao)+
C3dot*sin(omegao)+omega1dot*(C3*cos(omegao)-2*C2*sin(2*omegao));
n0ddot := n0dot*C1dottoC1+C1*D16;
e0ddot := e0dot*C0dottoC0-C0*((4+3*sqr(eta)+30*eo*eta+15.5*sqr(eo)+21*sqr(eta)*sqr(eo))*etadot+
(5+15*sqr(eta)+31*eo*eta+14*eo*power(eta,3))*e0dot+
B1*(D1dot*D6+D1*etadot*(30+67.5*sqr(eta)))+
B2*(D1dot*D7+D1*etadot*(5+37.5*sqr(eta)))*cos(2*omegao)+
B3*(D5dot*D8+D5*eta*etadot*(13.5+4*sqr(eta)))*sin(omegao)+
omega1dot*(C5*cos(omegao)-2*C4*sin(2*omegao)));
D17:=n0ddot/noii-sqr(n0dot/noii);
tsiddottotsi:=2*(tsidottotsi-C6)*tsidottotsi+
2*aoii*tsi*((1/3)*D17*sqr(beta)-2*C6*eo*e0dot+sqr(e0dot)+eo*e0ddot);
etaddot:=(e0ddot+2*e0dot*tsidottotsi)*s*tsi+eta*tsiddottotsi;
D18:=tsiddottotsi-sqr(tsidottotsi);
D19:=-sqr(phidottophi)*(1+1/sqr(eta))-eta*etaddot/sqr(phi);
D1ddot:=D1dot*(D14+D15)+D1*(D18-2*D19+(2/3)*D17+2*sqr(alpha)*sqr(e0dot)/power(beta,4)+2*eo*e0ddot/sqr(beta));
N0tdot:=N0Dot*(2*TOTHRD*D17+3*(sqr(e0dot)+EO*E0DDOT)/sqr(ALPHA)-6.*sqr(EO*E0DOT/sqr(ALPHA)) +4.*D18-7.*D19 )+
C1DoTtoC1*N0DDoT;
N0tdot:=N0tdot+C1*(C1DoTtoC1*D16+D9*ETADDoT+D10*E0DDOT+sqr(etadot)*(6.+30.*(Eo)*ETA+68.*sqr(EO))+
ETaDoT*E0DOT*(40.+30.* sqr(ETA)+272.*Eo*ETA)+ sqr(e0dot)*(17.+68.*sqr(ETA)) + B1*(D1DDoT*D2+
2.*D1DoT*D2DoT+D1*(ETaDDoT*D11+sqr(etadot)*(72.+54.*sqr(ETA))))+
B2*(D1DDoT*D3+2.*D1DoT*D3DoT+D1*(ETaDDoT*D12+ sqr(etadot)*(30.+30.*sqr(ETA)))) *COS(2*omegao)+ B3*((D5DoT*D14+D5*
(D18- 2.*D19)) * D4+ 2.*D4DoT*D5DoT+ D5*(ETaDDoT*D13+22.5*ETA*sqr(etadot))) *SIN(omegao)+omega1dot* ((7.*(1/3)*
(n0dot/noii)+4.*EO*E0DOT/sqr(BETA))*(C3*COS(omegao)-2.*C2*SIN(2*omegao))+
((2*C3DoT*COS(omeGao)-4*C2DoT*SIn(2*omegao)-omega1dot*(C3*SIN(omegao)+4* C2*COS(2*omegao))))));
smallN0DDoT:=N0DDoT*1.E9;
TEMP:=sqr(smallN0DDoT)-N0DoT*1.E18*N0tdot;
p:=(TEMP+sqr(smallN0DDoT))/TEMP;
gamma:=-N0TDoT/(N0DDoT*(P-2));
nd:=n0dot/(p*gamma);
q:=1-e0ddot/(e0dot*gamma);
ed:=e0dot/(q*gamma);
{atmospheric drag and gravitation}
if (abs(n0dot/noii*xmnpda))<0.00216 then
begin;
n:=noii+n0dot*tsince;
edot:=-(2/3)*n0dot*(1-eo)/noii;
e:=eo+edot*tsince;
Z1:=0.5*n0dot*sqr(tsince);
end
else
begin;
n:=noii+nd*(1-power((1-gamma*tsince),p));
edot:=e0dot;
e:=eo+ed*(1-power((1-gamma*tsince),q));
z1:=gamma*tsince;
z1:=1-z1;
Z1:=power(Z1,(p+1));
Z1:=Z1-1;
Z1:=Z1/(gamma*(p+1));
Z1:=Z1+tsince;
Z1:=Z1*n0dot/(p*gamma);
end;
omega:=omegao+omegadot*tsince+(7*Z1)/(3*noii)*omega1dot;
xnode:=xnodeo+xnodedot*tsince+xnode1dot*(7*Z1)/(3*noii);
M:=xmo+tsince*ldot+Z1+M1dot*(7*Z1)/(3*noii);
{Kepler}
M:=fmod2p(M);
temp:=M+eo*sin(M)*(1+e*cos(M));
for i:=1 to 10 do
begin;
deltaEw:=(M+e*sin(temp)-temp)/(1-e*cos(temp));
Ew:=DeltaEw+temp;
if abs(deltaEw)<=e6a then goto EndKepler;
temp:=Ew;
end;
EndKepler:
{short-period periodics}
a:=power((xke/n),tothrd);
beta:=sqrt(1-sqr(e));
sinf:=beta*sin(Ew)*(1/(1-e*cos(Ew)));
cosf:=(cos(Ew)-e)*(1/(1-e*cos(Ew)));
f:=actan(sinf,cosf);
sinu:=sinf*cos(omega)+cosf*sin(omega);
cosu:=cosf*cos(omega)-sinf*sin(omega);
rii:=(a*(1-sqr(e)))/(1.+(e*cosf));
deltaR:=(0.5*CK2*(1/(a*(1-sqr(e)))))*((1-sqr(theta))*(2*sqr(cosu)-1)-3.*(3*sqr(theta)-1))-(0.25*A30*sinio)*sinu;
temp:=3.*(0.5*CK2*sqr(1/(a*(1-sqr(e)))))*sinio*(2*sqr(cosu)-1)-(0.25*A30*(1/(a*(1-sqr(e)))))*e*sin(omega);
deltaI:=temp*cosio;
deltaU:=sin(xincl/2)*((0.5*CK2*sqr(1/(a*(1-sqr(e)))))*(0.5*(1.-7.*sqr(theta))*(2*sinu*cosu)-
3.*(1-5*sqr(theta))*(f-M+e*sinf))-(0.25*A30*(1/(a*(1-sqr(e)))))*sinio*cosu*(2.+ (e*cosf)))-
0.5*(0.25*A30*(1/(a*(1-sqr(e)))))*sqr(theta)*e*cos(omega)/cos(xincl/2);
lambda:=f+omega+xnode+(0.5*CK2*sqr(1/(a*(1-sqr(e)))))*(0.5*(1.+6.*cosio-7.*sqr(theta))*(2*sinu*cosu)-3*((1-5*sqr(theta))+
2.*cosio)*(f-M+e*sinf))+(0.25*A30*(1/(a*(1-sqr(e)))))*sinio*(cosio*e*cos(omega)/(1.+cosio)-(2. +(e*cosf))*cosu);
Y4:=sin(xincl/2)*sinu+cosu*deltaU+0.5*sinu*cos(xincl/2)*deltaI;
Y5:=sin(xincl/2)*cosu-sinu*deltaU+0.5*cosu*cos(xincl/2)*deltaI;
r:=rii+deltaR;
rdot:=n*a*e*sinf/beta+(-n*sqr(a/rii))*(1*CK2*(1/(a*(1-sqr(e))))*(1-sqr(theta))*(2*sinu*cosu)+(0.25*A30*sinio)*cosu);
rfdot:=n*sqr(a)*beta/rii+(-n*sqr(a/rii))*deltaR+a*(n*(a/rii))*sinio*temp;
{Unit-orientation vector}
cosI2:=sqrt(1-sqr(y4)-sqr(y5));
UV[1]:=2*y4*(y5*sin(lambda)-y4*cos(lambda))+cos(lambda);
UV[2]:=-2*y4*(y5*cos(lambda)+y4*sin(lambda))+sin(lambda);
UV[3]:=2*y4*cosI2;
VV[1]:=2*y5*(y5*sin(lambda)-y4*cos(lambda))-sin(lambda);
VV[2]:=-2*y5*(y5*cos(lambda)+y4*sin(lambda))+cos(lambda);
VV[3]:=2*y5*cosI2;
for i:=1 to 3 do
begin;
pos[i]:=r*uv[i];
vel[i]:=rdot*uv[i]+rfdot*vv[i];
end;
end; {Procedure SGP8}
Procedure Deep(ideep : integer);
const
zns = 1.19459E-5; c1ss = 2.9864797E-6; zes = 0.01675;
znl = 1.5835218E-4; c1l = 4.7968065E-7; zel = 0.05490;
zcosis = 0.91744867; zsinis = 0.39785416; zsings = -0.98088458;
zcosgs = 0.1945905; zcoshs = 1; zsinhs = 0;
q22 = 1.7891679E-6; q31 = 2.1460748E-6; q33 = 2.2123015E-7;
g22 = 5.7686396; g32 = 0.95240898; g44 = 1.8014998;
g52 = 1.0508330; g54 = 4.4108898; root22 = 1.7891679E-6;
root32 = 3.7393792E-7; root44 = 7.3636953E-9; root52 = 1.1428639E-7;
root54 = 2.1765803E-9; thdt = 4.3752691E-3;
label
{dpinit}
5,10,20,30,40,45,50,55,60,65,70,80,
{dpsec}
90,100,105,110,120,125,130,140,150,152,154,160,165,170,175,180,
{dpper}
205,210,218,220,230;
const { Typed constants to retain values between ENTRY calls }
iresfl : integer = 0; isynfl : integer = 0;
iret : integer = 0; iretn : integer = 0;
ls : integer = 0;
a1 : double = 0; a2 : double = 0; a3 : double = 0;
a4 : double = 0; a5 : double = 0; a6 : double = 0;
a7 : double = 0; a8 : double = 0; a9 : double = 0;
a10 : double = 0; ainv2 : double = 0; alfdp : double = 0;
aqnv : double = 0; atime : double = 0; betdp : double = 0;
bfact : double = 0; c : double = 0; cc : double = 0;
cosis : double = 0; cosok : double = 0; cosq : double = 0;
ctem : double = 0; d2201 : double = 0; d2211 : double = 0;
d3210 : double = 0; d3222 : double = 0; d4410 : double = 0;
d4422 : double = 0; d5220 : double = 0; d5232 : double = 0;
d5421 : double = 0; d5433 : double = 0; dalf : double = 0;
day : double = 0; dbet : double = 0; del1 : double = 0;
del2 : double = 0; del3 : double = 0; delt : double = 0;
dls : double = 0; e3 : double = 0; ee2 : double = 0;
eoc : double = 0; eq : double = 0; f2 : double = 0;
f220 : double = 0; f221 : double = 0; f3 : double = 0;
f311 : double = 0; f321 : double = 0; f322 : double = 0;
f330 : double = 0; f441 : double = 0; f442 : double = 0;
f522 : double = 0; f523 : double = 0; f542 : double = 0;
f543 : double = 0; fasx2 : double = 0; fasx4 : double = 0;
fasx6 : double = 0; ft : double = 0; g200 : double = 0;
g201 : double = 0; g211 : double = 0; g300 : double = 0;
g310 : double = 0; g322 : double = 0; g410 : double = 0;
g422 : double = 0; g520 : double = 0; g521 : double = 0;
g532 : double = 0; g533 : double = 0; gam : double = 0;
omegaq : double = 0; pe : double = 0; pgh : double = 0;
ph : double = 0; pinc : double = 0; pl : double = 0;
preep : double = 0; s1 : double = 0; s2 : double = 0;
s3 : double = 0; s4 : double = 0; s5 : double = 0;
s6 : double = 0; s7 : double = 0; savtsn : double = 0;
se : double = 0; se2 : double = 0; se3 : double = 0;
sel : double = 0; ses : double = 0; sgh : double = 0;
sgh2 : double = 0; sgh3 : double = 0; sgh4 : double = 0;
sghl : double = 0; sghs : double = 0; sh : double = 0;
sh2 : double = 0; sh3 : double = 0; sh1 : double = 0;
shs : double = 0; si : double = 0; si2 : double = 0;
si3 : double = 0; sil : double = 0; sini2 : double = 0;
sinis : double = 0; sinok : double = 0; sinq : double = 0;
sinzf : double = 0; sis : double = 0; sl : double = 0;
sl2 : double = 0; sl3 : double = 0; sl4 : double = 0;
sll : double = 0; sls : double = 0; sse : double = 0;
ssg : double = 0; ssh : double = 0; ssi : double = 0;
ssl : double = 0; stem : double = 0; step2 : double = 0;
stepn : double = 0; stepp : double = 0; temp : double = 0;
temp1 : double = 0; thgr : double = 0; x1 : double = 0;
x2 : double = 0; x2li : double = 0; x2omi : double = 0;
x3 : double = 0; x4 : double = 0; x5 : double = 0;
x6 : double = 0; x7 : double = 0; x8 : double = 0;
xfact : double = 0; xgh2 : double = 0; xgh3 : double = 0;
xgh4 : double = 0; xh2 : double = 0; xh3 : double = 0;
xi2 : double = 0; xi3 : double = 0; xl : double = 0;
xl2 : double = 0; xl3 : double = 0; xl4 : double = 0;
xlamo : double = 0; xldot : double = 0; xli : double = 0;
xls : double = 0; xmao : double = 0; xnddt : double = 0;
xndot : double = 0; xni : double = 0; xno2 : double = 0;
xnodce : double = 0; xnoi : double = 0; xnq : double = 0;
xomi : double = 0; xpidot : double = 0; xqncl : double = 0;
z1 : double = 0; z11 : double = 0; z12 : double = 0;
z13 : double = 0; z2 : double = 0; z21 : double = 0;
z22 : double = 0; z23 : double = 0; z3 : double = 0;
z31 : double = 0; z32 : double = 0; z33 : double = 0;
zcosg : double = 0; zcosgl : double = 0; zcosh : double = 0;
zcoshl : double = 0; zcosi : double = 0; zcosil : double = 0;
ze : double = 0; zf : double = 0; zm : double = 0;
zmo : double = 0; zmol : double = 0; zmos : double = 0;
zn : double = 0; zsing : double = 0; zsingl : double = 0;
zsinh : double = 0; zsinhl : double = 0; zsini : double = 0;
zsinil : double = 0; zx : double = 0; zy : double = 0;
begin
case ideep of
dpinit : begin { Entrance for deep space initialization }
thgr := Thetag(epoch);
eq := eo;
xnq := xnodp;
aqnv := 1/ao;
xqncl := xincl;
xmao := xmo;
xpidot := omgdt + xnodot;
sinq := Sin(xnodeo);
cosq := Cos(xnodeo);
omegaq := omegao;
{ Initialize lunar solar terms }
5: day := ds50 + 18261.5; {Days since 1900 Jan 0.5}
if (day = preep) then
Goto 10;
preep := day;
xnodce := 4.5236020 - 9.2422029E-4*day;
stem := Sin(xnodce);
ctem := Cos(xnodce);
zcosil := 0.91375164 - 0.03568096*ctem;
zsinil := Sqrt(1 - zcosil*zcosil);
zsinhl := 0.089683511*stem/zsinil;
zcoshl := Sqrt(1 - zsinhl*zsinhl);
c := 4.7199672 + 0.22997150*day;
gam := 5.8351514 + 0.0019443680*day;
zmol := fmod2p(c - gam);
zx := 0.39785416*stem/zsinil;
zy := zcoshl*ctem + 0.91744867*zsinhl*stem;
zx := Actan(zx,zy);
zx := gam + zx - xnodce;
zcosgl := Cos(zx);
zsingl := Sin(zx);
zmos := 6.2565837 + 0.017201977*day;
zmos := fmod2p(zmos);
{ Do solar terms }
10: savtsn := 1E20;
zcosg := zcosgs;
zsing := zsings;
zcosi := zcosis;
zsini := zsinis;
zcosh := cosq;
zsinh := sinq;
cc := c1ss;
zn := zns;
ze := zes;
zmo := zmos;
xnoi := 1/xnq;
ls := 30; {assign 30 to ls}
20: a1 := zcosg*zcosh + zsing*zcosi*zsinh;
a3 := -zsing*zcosh + zcosg*zcosi*zsinh;
a7 := -zcosg*zsinh + zsing*zcosi*zcosh;
a8 := zsing*zsini;
a9 := zsing*zsinh + zcosg*zcosi*zcosh;
a10 := zcosg*zsini;
a2 := cosiq*a7 + siniq*a8;
a4 := cosiq*a9 + siniq*a10;
a5 := -siniq*a7 + cosiq*a8;
a6 := -siniq*a9 + cosiq*a10;
x1 := a1*cosomo + a2*sinomo;
x2 := a3*cosomo + a4*sinomo;
x3 := -a1*sinomo + a2*cosomo;
x4 := -a3*sinomo + a4*cosomo;
x5 := a5*sinomo;
x6 := a6*sinomo;
x7 := a5*cosomo;
x8 := a6*cosomo;
z31 := 12*x1*x1 - 3*x3*x3;
z32 := 24*x1*x2 - 6*x3*x4;
z33 := 12*x2*x2 - 3*x4*x4;
z1 := 3*(a1*a1 + a2*a2) + z31*eqsq;
z2 := 6*(a1*a3 + a2*a4) + z32*eqsq;
z3 := 3*(a3*a3 + a4*a4) + z33*eqsq;
z11 := -6*a1*a5 + eqsq*(-24*x1*x7 - 6*x3*x5);
z12 := -6*(a1*a6 + a3*a5)
+ eqsq*(-24*(x2*x7 + x1*x8) - 6*(x3*x6 + x4*x5));
z13 := -6*a3*a6 + eqsq*(-24*x2*x8 - 6*x4*x6);
z21 := 6*a2*a5 + eqsq*(24*x1*x5 - 6*x3*x7);
z22 := 6*(a4*a5 + a2*a6)
+ eqsq*(24*(x2*x5 + x1*x6) - 6*(x4*x7 + x3*x8));
z23 := 6*a4*a6 + eqsq*(24*x2*x6 - 6*x4*x8);
z1 := z1 + z1 + bsq*z31;
z2 := z2 + z2 + bsq*z32;
z3 := z3 + z3 + bsq*z33;
s3 := cc*xnoi;
s2 := -0.5*s3/rteqsq;
s4 := s3*rteqsq;
s1 := -15*eq*s4;
s5 := x1*x3 + x2*x4;
s6 := x2*x3 + x1*x4;
s7 := x2*x4 - x1*x3;
se := s1*zn*s5;
si := s2*zn*(z11 + z13);
sl := -zn*s3*(z1 + z3 - 14 - 6*eqsq);
sgh := s4*zn*(z31 + z33 - 6);
sh := -zn*s2*(z21 + z23);
if (xqncl < 5.2359877E-2) then
sh := 0;
ee2 := 2*s1*s6;
e3 := 2*s1*s7;
xi2 := 2*s2*z12;
xi3 := 2*s2*(z13 - z11);
xl2 := -2*s3*z2;
xl3 := -2*s3*(z3 - z1);
xl4 := -2*s3*(-21 - 9*eqsq)*ze;
xgh2 := 2*s4*z32;
xgh3 := 2*s4*(z33 - z31);
xgh4 := -18*s4*ze;
xh2 := -2*s2*z22;
xh3 := -2*s2*(z23 - z21);
case ls of
30 : Goto 30;
40 : Goto 40;
else
Halt;
end; {case}
{ Do lunar terms }
30: sse := se;
ssi := si;
ssl := sl;
ssh := sh/siniq;
ssg := sgh - cosiq*ssh;
se2 := ee2;
si2 := xi2;
sl2 := xl2;
sgh2 := xgh2;
sh2 := xh2;
se3 := e3;
si3 := xi3;
sl3 := xl3;
sgh3 := xgh3;
sh3 := xh3;
sl4 := xl4;
sgh4 := xgh4;
zcosg := zcosgl;
zsing := zsingl;
zcosi := zcosil;
zsini := zsinil;
zcosh := zcoshl*cosq + zsinhl*sinq;
zsinh := sinq*zcoshl - cosq*zsinhl;
zn := znl;
cc := c1l;
ze := zel;
zmo := zmol;
ls := 40; {assign 40 to ls}
Goto 20;
40: sse := sse + se;
ssi := ssi + si;
ssl := ssl + sl;
ssg := ssg + sgh - cosiq/siniq*sh;
ssh := ssh + sh/siniq;
{ Geopotential resonance initialization for 12 hour orbits }
iresfl := 0;
isynfl := 0;
if (xnq < 0.0052359877) and (xnq > 0.0034906585) then
Goto 70;
if (xnq < 8.26E-3) or (xnq > 9.24E-3) then
exit;
if (eq < 0.5) then
exit;
iresfl := 1;
eoc := eq*eqsq;
g201 := -0.306 - (eq - 0.64)*0.440;
if (eq > 0.65) then
Goto 45;
g211 := 3.616 - 13.247*eq + 16.290*eqsq;
g310 := -19.302 + 117.390*eq - 228.419*eqsq + 156.591*eoc;
g322 := -18.9068 + 109.7927*eq - 214.6334*eqsq + 146.5816*eoc;
g410 := -41.122 + 242.694*eq - 471.094*eqsq + 313.953*eoc;
g422 := -146.407 + 841.880*eq - 1629.014*eqsq + 1083.435*eoc;
g520 := -532.114 + 3017.977*eq - 5740*eqsq + 3708.276*eoc;
Goto 55;
45: g211 := -72.099 + 331.819*eq - 508.738*eqsq + 266.724*eoc;
g310 := -346.844 + 1582.851*eq - 2415.925*eqsq + 1246.113*eoc;
g322 := -342.585 + 1554.908*eq - 2366.899*eqsq + 1215.972*eoc;
g410 := -1052.797 + 4758.686*eq - 7193.992*eqsq + 3651.957*eoc;
g422 := -3581.69 + 16178.11*eq - 24462.77*eqsq + 12422.52*eoc;
if (eq > 0.715) then
Goto 50;
g520 := 1464.74 - 4664.75*eq + 3763.64*eqsq;
Goto 55;
50: g520 := -5149.66 + 29936.92*eq - 54087.36*eqsq + 31324.56*eoc;
55: if (eq >= (0.7)) then
Goto 60;
g533 := -919.2277 + 4988.61*eq - 9064.77*eqsq + 5542.21*eoc;
g521 := -822.71072 + 4568.6173*eq - 8491.4146*eqsq + 5337.524*eoc;
g532 := -853.666 + 4690.25*eq - 8624.77*eqsq + 5341.4*eoc;
Goto 65;
60: g533 := -37995.78 + 161616.52*eq - 229838.2*eqsq + 109377.94*eoc;
g521 := -51752.104 + 218913.95*eq - 309468.16*eqsq + 146349.42*eoc;
g532 := -40023.88 + 170470.89*eq - 242699.48*eqsq + 115605.82*eoc;
65: sini2 := siniq*siniq;
f220 := 0.75*(1 + 2*cosiq + cosq2);
f221 := 1.5*sini2;
f321 := 1.875*siniq*(1 - 2*cosiq - 3*cosq2);
f322 := -1.875*siniq*(1 + 2*cosiq - 3*cosq2);
f441 := 35*sini2*f220;
f442 := 39.3750*sini2*sini2;
f522 := 9.84375*siniq*(sini2*(1 - 2*cosiq - 5*cosq2)
+ 0.33333333*(-2 + 4*cosiq + 6*cosq2));
f523 := siniq*(4.92187512*sini2*(-2 - 4*cosiq + 10*cosq2)
+ 6.56250012*(1 + 2*cosiq - 3*cosq2));
f542 := 29.53125*siniq*(2 - 8*cosiq + cosq2*(-12 + 8*cosiq + 10*cosq2));
f543 := 29.53125*siniq*(-2 - 8*cosiq + cosq2*(12 + 8*cosiq - 10*cosq2));
xno2 := xnq*xnq;
ainv2 := aqnv*aqnv;
temp1 := 3*xno2*ainv2;
temp := temp1*root22;
d2201 := temp*f220*g201;
d2211 := temp*f221*g211;
temp1 := temp1*aqnv;
temp := temp1*root32;
d3210 := temp*f321*g310;
d3222 := temp*f322*g322;
temp1 := temp1*aqnv;
temp := 2*temp1*root44;
d4410 := temp*f441*g410;
d4422 := temp*f442*g422;
temp1 := temp1*aqnv;
temp := temp1*root52;
d5220 := temp*f522*g520;
d5232 := temp*f523*g532;
temp := 2*temp1*root54;
d5421 := temp*f542*g521;
d5433 := temp*f543*g533;
xlamo := xmao + xnodeo + xnodeo - thgr - thgr;
bfact := xlldot + xnodot + xnodot - thdt - thdt;
bfact := bfact + ssl + ssh + ssh;
Goto 80;
{ Synchronous resonance terms initialization }
70: iresfl := 1;
isynfl := 1;
g200 := 1 + eqsq*(-2.5 + 0.8125*eqsq);
g310 := 1 + 2*eqsq;
g300 := 1 + eqsq*(-6 + 6.60937*eqsq);
f220 := 0.75*(1 + cosiq)*(1 + cosiq);
f311 := 0.9375*siniq*siniq*(1 + 3*cosiq) - 0.75*(1 + cosiq);
f330 := 1 + cosiq;
f330 := 1.875*f330*f330*f330;
del1 := 3*xnq*xnq*aqnv*aqnv;
del2 := 2*del1*f220*g200*q22;
del3 := 3*del1*f330*g300*q33*aqnv;
del1 := del1*f311*g310*q31*aqnv;
fasx2 := 0.13130908;
fasx4 := 2.8843198;
fasx6 := 0.37448087;
xlamo := xmao + xnodeo + omegao - thgr;
bfact := xlldot + xpidot - thdt;
bfact := bfact + ssl + ssg + ssh;
80: xfact := bfact - xnq;
{ Initialize integrator }
xli := xlamo;
xni := xnq;
atime := 0;
stepp := 720;
stepn := -720;
step2 := 259200;
end; {dpinit}
dpsec : begin { Entrance for deep space secular effects }
xll := xll + ssl*t;
omgasm := omgasm + ssg*t;
xnodes := xnodes + ssh*t;
_em := eo + sse*t;
xinc := xincl + ssi*t;
if (xinc >= 0) then
Goto 90;
xinc := -xinc;
xnodes := xnodes + pi;
omgasm := omgasm - pi;
90: if (iresfl = 0) then
exit;
100: if (atime = 0) then
Goto 170;
if (t >= 0) and (atime < 0) then
Goto 170;
if (t < 0) and (atime >= 0) then
Goto 170;
105: if (Abs(t) >= Abs(atime)) then
Goto 120;
delt := stepp;
if (t >= 0) then
delt := stepn;
110: iret := 100; {assign 100 to iret}
Goto 160;
120: delt := stepn;
if (t > 0) then
delt := stepp;
125: if (Abs(t - atime) < stepp) then
Goto 130;
iret := 125; {assign 125 to iret}
Goto 160;
130: ft := t - atime;
iretn := 140; {assign 140 to iretn}
Goto 150;
140: xn := xni + xndot*ft + xnddt*ft*ft*0.5;
xl := xli + xldot*ft + xndot*ft*ft*0.5;
temp := -xnodes + thgr + t*thdt;
xll := xl - omgasm + temp;
if (isynfl = 0) then
xll := xl + temp + temp;
exit;
{ Dot terms calculated }
150: if (isynfl = 0) then
Goto 152;
xndot := del1*Sin(xli - fasx2) + del2*Sin(2*(xli - fasx4))
+ del3*Sin(3*(xli - fasx6));
xnddt := del1*Cos(xli - fasx2)
+ 2*del2*Cos(2*(xli - fasx4))
+ 3*del3*Cos(3*(xli - fasx6));
Goto 154;
152: xomi := omegaq + omgdt*atime;
x2omi := xomi + xomi;
x2li := xli + xli;
xndot := d2201*Sin(x2omi + xli - g22)
+ d2211*Sin(xli - g22)
+ d3210*Sin(xomi + xli - g32)
+ d3222*Sin(-xomi + xli - g32)
+ d4410*Sin(x2omi + x2li - g44)
+ d4422*Sin(x2li - g44)
+ d5220*Sin(xomi + xli - g52)
+ d5232*Sin(-xomi + xli - g52)
+ d5421*Sin(xomi + x2li - g54)
+ d5433*Sin(-xomi + x2li - g54);
xnddt := d2201*Cos(x2omi + xli - g22)
+ d2211*Cos(xli - g22)
+ d3210*Cos(xomi + xli - g32)
+ d3222*Cos(-xomi + xli - g32)
+ d5220*Cos(xomi + xli - g52)
+ d5232*Cos(-xomi + xli - g52)
+ 2*(d4410*Cos(x2omi + x2li - g44)
+ d4422*Cos(x2li - g44)
+ d5421*Cos(xomi + x2li - g54)
+ d5433*Cos(-xomi + x2li - g54));
154: xldot := xni + xfact;
xnddt := xnddt*xldot;
case iretn of
140 : Goto 140;
165 : Goto 165;
else
Halt;
end; {case}
{ Integrator }
160: iretn := 165; {assign 165 to iretn}
Goto 150;
165: xli := xli + xldot*delt + xndot*step2;
xni := xni + xndot*delt + xnddt*step2;
atime := atime + delt;
case iret of
100 : Goto 100;
125 : Goto 125;
else
Halt;
end; {case}
{ Epoch restart }
170: if (t >= 0) then
Goto 175;
delt := stepn;
Goto 180;
175: delt := stepp;
180: atime := 0;
xni := xnq;
xli := xlamo;
Goto 125;
end; {dpsec}
dpper : begin { Entrance for lunar-solar periodics }
sinis := Sin(xinc);
cosis := Cos(xinc);
if (Abs(savtsn - t) < 30) then
Goto 210;
savtsn := t;
zm := zmos + zns*t;
205: zf := zm + 2*zes*Sin(zm);
sinzf := Sin(zf);
f2 := 0.5*sinzf*sinzf - 0.25;
f3 := -0.5*sinzf*Cos(zf);
ses := se2*f2 + se3*f3;
sis := si2*f2 + si3*f3;
sls := sl2*f2 + sl3*f3 + sl4*sinzf;
sghs := sgh2*f2 + sgh3*f3 + sgh4*sinzf;
shs := sh2*f2 + sh3*f3;
zm := zmol + znl*t;
zf := zm + 2*zel*Sin(zm);
sinzf := Sin(zf);
f2 := 0.5*sinzf*sinzf - 0.25;
f3 := -0.5*sinzf*Cos(zf);
sel := ee2*f2 + e3*f3;
sil := xi2*f2 + xi3*f3;
sll := xl2*f2 + xl3*f3 + xl4*sinzf;
sghl := xgh2*f2 + xgh3*f3 + xgh4*sinzf;
sh1 := xh2*f2 + xh3*f3;
pe := ses + sel;
pinc := sis + sil;
pl := sls + sll;
210: pgh := sghs + sghl;
ph := shs + sh1;
xinc := xinc + pinc;
_em := _em + pe;
if (xqncl < 0.2) then
Goto 220;
Goto 218;
{ Apply periodics directly }
218: ph := ph/siniq;
pgh := pgh - cosiq*ph;
omgasm := omgasm + pgh;
xnodes := xnodes + ph;
xll := xll + pl;
Goto 230;
{ Apply periodics with Lyddane modification }
220: sinok := Sin(xnodes);
cosok := Cos(xnodes);
alfdp := sinis*sinok;
betdp := sinis*cosok;
dalf := ph*cosok + pinc*cosis*sinok;
dbet := -ph*sinok + pinc*cosis*cosok;
alfdp := alfdp + dalf;
betdp := betdp + dbet;
xls := xll + omgasm + cosis*xnodes;
dls := pl + pgh - pinc*xnodes*sinis;
xls := xls + dls;
xnodes := Actan(alfdp,betdp);
xll := xll + pl;
omgasm := xls - xll - Cos(xinc)*xnodes;
230: end; {dpper}
end; {case}
end; {Procedure Deep}
Procedure Call_dpinit(var eosq,sinio,cosio,betao,aodp,theta2,
sing,cosg,betao2,xmdot,omgdot,
xnodott,xnodpp : double);
begin
eqsq := eosq; siniq := sinio; cosiq := cosio; rteqsq := betao;
ao := aodp; cosq2 := theta2; sinomo := sing; cosomo := cosg;
bsq := betao2; xlldot := xmdot; omgdt := omgdot; xnodot := xnodott;
xnodp := xnodpp;
Deep(1);
eosq := eqsq; sinio := siniq; cosio := cosiq; betao := rteqsq;
aodp := ao; theta2 := cosq2; sing := sinomo; cosg := cosomo;
betao2 := bsq; xmdot := xlldot; omgdot := omgdt; xnodott := xnodot;
xnodpp := xnodp;
end; {Procedure Call_dpinit}
Procedure Call_dpsec(var xmdf,omgadf,xnode,emm,xincc,xnn,tsince : double);
begin
xll := xmdf; omgasm := omgadf; xnodes := xnode; {_em := emm;
xinc := xincc;} xn := xnn; t := tsince;
Deep(2);
xmdf := xll; omgadf := omgasm; xnode := xnodes; emm := _em;
xincc := xinc; xnn := xn; tsince := t;
end; {Procedure Call_dpsec}
Procedure Call_dpper(var e,xincc,omgadf,xnode,xmam : double);
begin
_em := e; xinc := xincc; omgasm := omgadf; xnodes := xnode;
xll := xmam;
Deep(3);
e := _em; xincc := xinc; omgadf := omgasm; xnode := xnodes;
xmam := xll;
end; {Procedure Call_dpper}
Procedure SDP8(tsince : double; var iflag : integer; var pos,vel : vector);
const rho: double=0.15696615;
var a1, cosi, theta2, tthmun, eosq, betao2, betao, del1, a0, delo, aodp, xnodp,b : double;
po, pom2, sini, sing, cosg, temp, sinio2, cosio2, theta4, unm5th, unmth2, a3cof : double;
pardt1, pardt2, pardt4, xmdt1, xgdt1, xhdt1, xlldot, omgdt,xnodot, tsi, eta, eta2,psim2, alpha2, eeta : double;
z1,z7,xmamdf,omgasm,xnodes, xn, em, xmam, zc2,sine,cose,zc5,cape : double;
cos2g, d1, d2, d3, d4, d5, b1, b2, b3, c0, c1, c4, c5, xndt, xndtn, edot : double;
am, beta2m,SINOS, cosos, AXNM, AYNM, PM, G1, G2,G3,G4, G5, BETA, SNF,CSF,FM : double;
dr, diwc, di, sni2du, xlamb, y4, y5, r, rdot, rvdot, snlamb, cslamb, ux, vx, uy, vy, uz, vz : double;
snfg, csfg, sn2f2g, cs2f2g, ecosf, g10, rm, aovr, g13, g14, sini2 : double;
i : integer;
label EndKepler;
begin;
A1:=power((XKE/XNO),TOTHRD);
COSI:=COS(XINCL);
THETA2:=COSI*COSI;
TTHMUN:=3*THETA2-1;
EOSQ:=EO*EO;
BETAO2:=1.-EOSQ;
BETAO:=SQRT(BETAO2);
DEL1:=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2);
AO:=A1*(1-DEL1*(.5*TOTHRD+DEL1*(1+134/81*DEL1)));
DELO:=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2);
AODP:=AO/(1.-DELO);
XNODP:=XNO/(1.+DELO);
B:=2.*BSTAR/RHO;
PO:=AODP*BETAO2;
POM2:=1./(PO*PO);
SINI:=SIN(XINCL);
SING:=SIN(OMEGAO);
COSG:=COS(OMEGAO);
TEMP:=0.5*XINCL;
SINIO2:=SIN(TEMP);
COSIO2:=COS(TEMP);
THETA4:=sqr(THETA2);
UNM5TH:=1-5*THETA2;
UNMTH2:=1-THETA2;
A3COF:=-XJ3/CK2*power(AE,3);
PARDT1:=3.*CK2*POM2*XNODP;
PARDT2:=PARDT1*CK2*POM2;
PARDT4:=1.25*CK4*POM2*POM2*XNODP;
XMDT1:=0.5*PARDT1*BETAO*TTHMUN;
XGDT1:=-0.5*PARDT1*UNM5TH;
XHDT1:=-PARDT1*COSI;
XLLDOT:=XNODP+XMDT1+0.0625*PARDT2*BETAO*(13-78*THETA2+137*THETA4);
OMGDT:=XGDT1+0.0625*PARDT2*(7-114*THETA2+395*THETA4)+PARDT4*(3-36*THETA2+49*THETA4);
XNODOT:=XHDT1+(0.5*PARDT2*(4-19*THETA2)+2*PARDT4*(3-7*THETA2))*COSI;
TSI:=1/(PO-S);
ETA:=EO*S*TSI;
ETA2:=sqr(ETA);
PSIM2:=ABS(1./(1.-ETA2));
ALPHA2:=1.+EOSQ;
EETA:=EO*ETA;
COS2G:=2*sqr(COSG)-1;
D5:=TSI*PSIM2;
D1:=D5/PO;
D2:=12.+ETA2*(36.+4.5*ETA2);
D3:=ETA2*(15.+2.5*ETA2);
D4:=ETA*(5.+3.75*ETA2);
B1:=CK2*TTHMUN;
B2:=-CK2*UNMTH2;
B3:=A3COF*SINI;
C0:=0.5*B*RHO*QOMS2T*XNODP*AODP*power(TSI,4)*power(PSIM2,3.5)/SQRT(ALPHA2);
C1:=1.5*XNODP*sqr(ALPHA2)*C0;
C4:=D1*D3*B2;
C5:=D5*D4*B3;
XNDT:=C1*((2+ETA2*(3+34*EOSQ)+5*EETA*(4+ETA2)+8.5*EOSQ)+D1*D2*B1+ C4*COS2G+C5*SING);
XNDTN:=XNDT/XNODP;
EDOT:=-TOTHRD*XNDTN*(1.-EO);
Call_DpINIT(EOSQ,SINI,COSI,BETAO,AODP,THETA2,SING,COSG,BETAO2,XLLDOT,OMGDT,XNODOT,XNODP);
Z1:=0.5*XNDT*TSINCE*TSINCE;
Z7:=3.5*TOTHRD*Z1/XNODP;
XMAMDF:=XMO+XLLDOT*TSINCE;
OMGASM:=OMEGAO+OMGDT*TSINCE+Z7*XGDT1;
XNODES:=XNODEO+XNODOT*TSINCE+Z7*XHDT1;
XN:=XNODP;
CALL_DPSEC(XMAMDF,OMGASM,XNODES,EM,XINC,XN,TSINCE);
XN:=XN+XNDT*TSINCE;
EM:=EM+EDOT*TSINCE;
XMAM:=XMAMDF+Z1+Z7*XMDT1;
CALL_DPPER(EM,XINC,OMGASM,XNODES,XMAM);
XMAM:=FMOD2P(XMAM);
ZC2:=XMAM+EM*SIN(XMAM)*(1+EM*COS(XMAM));
for I:=1 to 10 do
begin;
SINE:=SIN(ZC2);
COSE:=COS(ZC2);
ZC5:=1/(1-EM*COSE);
CAPE:=(XMAM+EM*SINE-ZC2)*ZC5+ZC2;
IF(ABS(CAPE-ZC2) < E6A) then GOTO Endkepler;
ZC2:=CAPE;
end;
EndKepler:
AM:=power((XKE/XN),TOTHRD);
BETA2M:=1.-EM*EM;
SINOS:=SIN(OMGASM);
COSOS:=COS(OMGASM);
AXNM:=EM*COSOS;
AYNM:=EM*SINOS;
PM:=AM*BETA2M;
G1:=1./PM;
G2:=0.5*CK2*G1;
G3:=G2*G1;
BETA:=SQRT(BETA2M);
G4:=0.25*A3COF*SINI;
G5:=0.25*A3COF*G1;
SNF:=BETA*SINE*ZC5;
CSF:=(COSE-EM)*ZC5;
FM:=ACTAN(SNF,CSF);
SNFG:=SNF*COSOS+CSF*SINOS;
CSFG:=CSF*COSOS-SNF*SINOS;
SN2F2G:=2.*SNFG*CSFG;
CS2F2G:=2.*sqr(CSFG)-1;
ECOSF:=EM*CSF;
G10:=FM-XMAM+EM*SNF;
RM:=PM/(1.+ECOSF);
AOVR:=AM/RM;
G13:=XN*AOVR;
G14:=-G13*AOVR;
DR:=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG;
DIWC:=3.*G3*SINI*CS2F2G-G5*AYNM;
DI:=DIWC*COSI;
SINI2:=SIN(0.5*XINC);
SNI2DU:=SINIO2*(G3*(0.5*(1-7*THETA2)*SN2F2G-3*UNM5TH*G10)-G5*SINI*CSFG*(2+ECOSF))-0.5*G5*THETA2*AXNM/COSIO2;
XLAMB:=FM+OMGASM+XNODES+G3*(0.5*(1+6*COSI-7*THETA2)*SN2F2G-3*(UNM5TH+2*COSI)*G10)+G5*SINI*(COSI*AXNM/(1+COSI)-
(2+ECOSF)*CSFG);
Y4:=SINI2*SNFG+CSFG*SNI2DU+0.5*SNFG*COSIO2*DI;
Y5:=SINI2*CSFG-SNFG*SNI2DU+0.5*CSFG*COSIO2*DI;
R:=RM+DR;
RDOT:=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG);
RVDOT:=XN*sqr(AM)*BETA/RM+G14*DR+AM*G13*SINI*DIWC;
SNLAMB:=SIN(XLAMB);
CSLAMB:=COS(XLAMB);
TEMP:=2.*(Y5*SNLAMB-Y4*CSLAMB);
UX:=Y4*TEMP+CSLAMB;
VX:=Y5*TEMP-SNLAMB;
TEMP:=2.*(Y5*CSLAMB+Y4*SNLAMB);
UY:=-Y4*TEMP+SNLAMB;
VY:=-Y5*TEMP+CSLAMB;
TEMP:=2.*SQRT(1.-Y4*Y4-Y5*Y5);
UZ:=Y4*TEMP;
VZ:=Y5*TEMP;
pos[1]:=R*UX;
pos[2]:=R*UY;
pos[3]:=R*UZ;
vel[1]:=RDOT*UX+RVDOT*VX;
vel[2]:=RDOT*UY+RVDOT*VY;
vel[3]:=RDOT*UZ+RVDOT*VZ;
end;
procedure SGP8call (time : double; var pos,vel : vector);
var
tsince : double;
begin
tsince := (time - julian_epoch) * xmnpda;
if ideep = 0 then
SGP8(tsince,iflag,pos,vel)
else
SDP8(tsince,iflag,pos,vel);
end; {Procedure SGP8call}
begin;
Define_Derived_Constants;
end.
Ho provato con free pascal ma non ne ho tirato fuori nulla.
Si tratta del propagatore SDP8 (magari avercelo in matlab).
Lo script potete scaricarlo da qui: http://canasta.altervista.org/sgp8sdp8.pas
unit SGP8SDP8;
{ Author: Dominik Brodowski }
{ Current Vision: 2001 Mar 03 }
{ Version: 1.00 }
{ uses the SGP4-PL2A TPU package by Dr TS Kelso }
{$N+}
INTERFACE
Uses Support,
SGP_Init,
SGP_Math,SGP_Time;
Procedure SGP8call(time : double; var pos,vel : vector);
Procedure SGP8(tsince : double; var iflag : integer; var pos,vel : vector);
Procedure SDP8(tsince : double; var iflag : integer; var pos,vel : vector);
IMPLEMENTATION
Uses SGP_Intf;
var
{dpinit}
eqsq,siniq,cosiq,rteqsq,ao,cosq2,sinomo,cosomo,
bsq,xlldot,omgdt,xnodot,xnodp : double;
{dpsec/dpper}
xll,omgasm,xnodes,_em,xinc,xn,t : double;
qoms2t : double;
Procedure Define_Derived_Constants;
begin
xke := Sqrt(3600*ge/Cube(xkmper)); {Sqrt(ge) ER^3/min^2}
qoms2t := Sqr(Sqr(qo-s)); {(qo-s)^4 ER^4}
end; {Define_Derived_Constants}
Procedure SGP8(tsince : double; var iflag : integer; var pos,vel : vector);
label
EndKepler;
const
a1 : double = 0; cosio : double = 0; sinio : double = 0;
delta1 : double = 0; a0 : double = 0; delta0 : double = 0;
B : double = 0; noii : double = 0; aoii : double = 0;
rho : double = (2.461*0.00001*XKMPER); beta : double = 0;
theta : double = 0; M1dot : double = 0; omega1dot: double= 0;
xnode1dot: double = 0; M2dot : double = 0; omega2dot: double= 0;
xnode2dot: double = 0; ldot : double = 0; omegadot: double = 0;
xnodedot : double = 0; tsi: double = 0; temp : double = 0;
eta : double = 0; phi : double = 0; alpha : double = 0;
C0 : double = 0; C1 : double = 0; D1 : double = 0;
D2 : double = 0; D3 : double = 0; D4 : double = 0;
D5 : double = 0; B1 : double = 0; B2 : double = 0;
B3 : double = 0; A30 : double = 0; C2 : double = 0;
C3 : double = 0; n0dot : double = 0; C4 : double = 0;
C5 : double = 0; D6 : double = 0; D7 : double = 0;
D8 : double = 0; e0dot : double = 0; etadot : double = 0;
alphadottoalpha : double = 0; C0dottoC0:double = 0;
tsidottotsi : double = 0; C1dottoC1:double = 0;
phidottophi : double = 0; C6 :double = 0;
D9 : double = 0; D10 : double = 0; D11 : double = 0;
D12 : double = 0; D13 : double = 0; D14 : double = 0;
D15 : double = 0; D1dot : double = 0; D2dot : double = 0;
D3dot : double = 0; D4dot : double = 0; D5dot : double = 0;
C2dot : double = 0; C3dot : double = 0; D16 : double = 0;
n0ddot : double = 0; e0ddot : double = 0; D17 : double = 0;
tsiddottotsi : double = 0; etaddot : double = 0;
D18 : double = 0; D19 : double = 0; D1ddot : double = 0;
n0tdot : double = 0; p : double = 0; gamma : double = 0;
nd : double = 0; q : double = 0; ed : double = 0;
n : double = 0; edot : double = 0; e : double = 0;
Z1 : double = 0; omega : double = 0; xnode : double = 0;
M : double = 0; Ew : double = 0; deltaEw : double = 0;
a : double = 0; sinf : double = 0; cosf : double = 0;
u : double = 0; sino : double = 0; coso : double = 0;
rii : double = 0; rdotii : double = 0; rfdotii : double = 0;
deltar : double = 0; deltardot:double= 0; deltau : double = 0;
deltaLambda:double= 0; r : double = 0; rdot : double = 0;
rfdot : double = 0; y4 : double = 0; y5 : double = 0;
lambda : double = 0; deltaI : double = 0; deltarfdot:double= 0;
f : double = 0; cosI2 : double = 0; sinu : double = 0;
cosu : double = 0; smallN0DDoT : double=0;
var
i : integer;
UV,VV : vector;
begin;
{constants I}
A1:=power((xke/xno),tothrd);
cosio := cos(xincl);
sinio := sin(xincl);
delta1 := 1.5*CK2*(3*cosio*cosio-1)/(sqr(a1)*power((1-sqr(eo)),1.5));
a0 := a1*(1-(1/3)*delta1-sqr(delta1)-(134/81)*sqr(delta1)*delta1);
delta0 := 1.5*CK2*(3*cosio*cosio-1)/(sqr(a0)*power((1-sqr(eo)),1.5));
noii := xno/(1+delta0);
aoii := a0/(1-delta0);
{B term}
B := 2*BStar/rho;
{constants II}
beta:=sqrt(1-sqr(eo));
theta:=cosio;
temp:=(noii*CK2)/(sqr(aoii)*sqr(beta)*beta);
M1dot:=-1.5*temp*(1-3*sqr(theta));
omega1dot:=-1.5*(temp/beta)*(1-5*sqr(theta));
xnode1dot:=-3*(temp/beta)*theta;
M2dot:=(3/16)*(temp*CK2/(sqr(aoii)*power(beta,4)))*(13-78*sqr(theta)+
137*power(theta,4));
omega2dot:=(3/16)*(temp*CK2/(sqr(aoii)*power(beta,5)))*(7-114*sqr(theta)+
395*power(theta,4))+1.25*(noii*CK4)/(power(aoii,4)*power(beta,8))*
(3-36*sqr(theta)+49*power(theta,4));
xnode2dot:=1.5*(temp*CK2/(sqr(aoii)*power(beta,5)))*theta*(4-19*sqr(theta))+
2.5*(noii*CK4)/(power(aoii,4)*power(beta,8))*theta*(3-7*sqr(theta));
ldot := noii + M1dot + M2dot;
omegadot := omega1dot+omega2dot;
xnodedot := xnode1dot + xnode2dot;
tsi := 1/(aoii*sqr(beta)-s);
eta := eo*s*tsi;
phi := sqrt(1-sqr(eta));
alpha := sqrt(1+sqr(eo));
C0:=0.5*B*rho*qoms2t*noii*aoii*power(tsi,4)/(alpha*power(phi,7));
C1:=1.5*noii*power(alpha,4)*C0;
D1:=(tsi*(1/sqr(phi)))/(aoii*sqr(beta));
D2:=12+36*sqr(eta)+4.5*power(eta,4);
D3:=15*sqr(eta)+2.5*power(eta,4);
D4:=5*eta+(15/4)*power(eta,3);
D5:=tsi/sqr(phi);
B1:=-CK2*(1-3*sqr(theta));
B2:=-CK2*(1-sqr(theta));
A30:=-XJ3*power(Ae,3)/CK2;
B3:=A30*sinio;
C2:=D1*D3*B2;
C3:=D4*D5*B3;
n0dot:=C1*((2+sqr(ETA)*(3+34*sqr(EO))+5*(Eo)*(ETA)*(4+sqr(ETA))+8.5*sqr(EO))+D1*D2*B1+ C2*COS(2*omegao)+C3*SIN(omegao));
D6:=30*eta+22.5*power(eta,3);
D7:=5*Eta+12.5*power(eta,3);
D8:=1+6.75*sqr(eta)+power(eta,4);
C4 := D1*D7*B2;
C5:=D5*D8*B3;
e0dot:=-C0*(ETA*(4+sqr(ETA)+sqr(EO)*(15.5+7*sqr(ETA)))+EO*(5+15*sqr(ETA))+D1*D6*B1+C4*COS(2*omegao)+C5*SIN(omegao));
alphadottoalpha:=eo*e0dot/(sqr(alpha));
C6:=(1/3)*(n0dot/noii);
tsidottotsi:=2*aoii*tsi*(C6*sqr(beta)+eo*e0dot);
etadot:=(e0dot+eo*tsidottotsi)*s*tsi;
phidottophi:=-eta*etadot/sqr(phi);
C0dottoC0:=C6+4*tsidottotsi-alphadottoalpha-7*phidottophi;
C1dottoC1:=n0dot/noii+4*alphadottoalpha+C0dottoC0;
D9:=6*eta+20*eo+15*eo*sqr(eta)+68*sqr(eo)*eta;
D10:=20*eta+5*power(eta,3)+17*eo+68*(eo)*sqr(eta);
D11:=eta*(72+18*sqr(eta));
D12:=eta*(30+10*sqr(eta));
D13:=5+11.25*sqr(eta);
D14:=tsidottotsi-2*phidottophi;
D15:=2*(C6+eo*e0dot/sqr(beta));
D1dot:=D1*(D14+D15);
D2dot:=etadot*D11;
D3dot:=etadot*D12;
D4dot:=etadot*D13;
D5dot:=D5*D14;
C2dot:=B2*(D1dot*D3+D1*D3dot);
C3dot:=B3*(D5dot*D4+D5*D4dot);
D16:=D9*etadot+D10*e0dot+B1*(D1dot*D2+D1*D2dot)+C2dot*cos(2*omegao)+
C3dot*sin(omegao)+omega1dot*(C3*cos(omegao)-2*C2*sin(2*omegao));
n0ddot := n0dot*C1dottoC1+C1*D16;
e0ddot := e0dot*C0dottoC0-C0*((4+3*sqr(eta)+30*eo*eta+15.5*sqr(eo)+21*sqr(eta)*sqr(eo))*etadot+
(5+15*sqr(eta)+31*eo*eta+14*eo*power(eta,3))*e0dot+
B1*(D1dot*D6+D1*etadot*(30+67.5*sqr(eta)))+
B2*(D1dot*D7+D1*etadot*(5+37.5*sqr(eta)))*cos(2*omegao)+
B3*(D5dot*D8+D5*eta*etadot*(13.5+4*sqr(eta)))*sin(omegao)+
omega1dot*(C5*cos(omegao)-2*C4*sin(2*omegao)));
D17:=n0ddot/noii-sqr(n0dot/noii);
tsiddottotsi:=2*(tsidottotsi-C6)*tsidottotsi+
2*aoii*tsi*((1/3)*D17*sqr(beta)-2*C6*eo*e0dot+sqr(e0dot)+eo*e0ddot);
etaddot:=(e0ddot+2*e0dot*tsidottotsi)*s*tsi+eta*tsiddottotsi;
D18:=tsiddottotsi-sqr(tsidottotsi);
D19:=-sqr(phidottophi)*(1+1/sqr(eta))-eta*etaddot/sqr(phi);
D1ddot:=D1dot*(D14+D15)+D1*(D18-2*D19+(2/3)*D17+2*sqr(alpha)*sqr(e0dot)/power(beta,4)+2*eo*e0ddot/sqr(beta));
N0tdot:=N0Dot*(2*TOTHRD*D17+3*(sqr(e0dot)+EO*E0DDOT)/sqr(ALPHA)-6.*sqr(EO*E0DOT/sqr(ALPHA)) +4.*D18-7.*D19 )+
C1DoTtoC1*N0DDoT;
N0tdot:=N0tdot+C1*(C1DoTtoC1*D16+D9*ETADDoT+D10*E0DDOT+sqr(etadot)*(6.+30.*(Eo)*ETA+68.*sqr(EO))+
ETaDoT*E0DOT*(40.+30.* sqr(ETA)+272.*Eo*ETA)+ sqr(e0dot)*(17.+68.*sqr(ETA)) + B1*(D1DDoT*D2+
2.*D1DoT*D2DoT+D1*(ETaDDoT*D11+sqr(etadot)*(72.+54.*sqr(ETA))))+
B2*(D1DDoT*D3+2.*D1DoT*D3DoT+D1*(ETaDDoT*D12+ sqr(etadot)*(30.+30.*sqr(ETA)))) *COS(2*omegao)+ B3*((D5DoT*D14+D5*
(D18- 2.*D19)) * D4+ 2.*D4DoT*D5DoT+ D5*(ETaDDoT*D13+22.5*ETA*sqr(etadot))) *SIN(omegao)+omega1dot* ((7.*(1/3)*
(n0dot/noii)+4.*EO*E0DOT/sqr(BETA))*(C3*COS(omegao)-2.*C2*SIN(2*omegao))+
((2*C3DoT*COS(omeGao)-4*C2DoT*SIn(2*omegao)-omega1dot*(C3*SIN(omegao)+4* C2*COS(2*omegao))))));
smallN0DDoT:=N0DDoT*1.E9;
TEMP:=sqr(smallN0DDoT)-N0DoT*1.E18*N0tdot;
p:=(TEMP+sqr(smallN0DDoT))/TEMP;
gamma:=-N0TDoT/(N0DDoT*(P-2));
nd:=n0dot/(p*gamma);
q:=1-e0ddot/(e0dot*gamma);
ed:=e0dot/(q*gamma);
{atmospheric drag and gravitation}
if (abs(n0dot/noii*xmnpda))<0.00216 then
begin;
n:=noii+n0dot*tsince;
edot:=-(2/3)*n0dot*(1-eo)/noii;
e:=eo+edot*tsince;
Z1:=0.5*n0dot*sqr(tsince);
end
else
begin;
n:=noii+nd*(1-power((1-gamma*tsince),p));
edot:=e0dot;
e:=eo+ed*(1-power((1-gamma*tsince),q));
z1:=gamma*tsince;
z1:=1-z1;
Z1:=power(Z1,(p+1));
Z1:=Z1-1;
Z1:=Z1/(gamma*(p+1));
Z1:=Z1+tsince;
Z1:=Z1*n0dot/(p*gamma);
end;
omega:=omegao+omegadot*tsince+(7*Z1)/(3*noii)*omega1dot;
xnode:=xnodeo+xnodedot*tsince+xnode1dot*(7*Z1)/(3*noii);
M:=xmo+tsince*ldot+Z1+M1dot*(7*Z1)/(3*noii);
{Kepler}
M:=fmod2p(M);
temp:=M+eo*sin(M)*(1+e*cos(M));
for i:=1 to 10 do
begin;
deltaEw:=(M+e*sin(temp)-temp)/(1-e*cos(temp));
Ew:=DeltaEw+temp;
if abs(deltaEw)<=e6a then goto EndKepler;
temp:=Ew;
end;
EndKepler:
{short-period periodics}
a:=power((xke/n),tothrd);
beta:=sqrt(1-sqr(e));
sinf:=beta*sin(Ew)*(1/(1-e*cos(Ew)));
cosf:=(cos(Ew)-e)*(1/(1-e*cos(Ew)));
f:=actan(sinf,cosf);
sinu:=sinf*cos(omega)+cosf*sin(omega);
cosu:=cosf*cos(omega)-sinf*sin(omega);
rii:=(a*(1-sqr(e)))/(1.+(e*cosf));
deltaR:=(0.5*CK2*(1/(a*(1-sqr(e)))))*((1-sqr(theta))*(2*sqr(cosu)-1)-3.*(3*sqr(theta)-1))-(0.25*A30*sinio)*sinu;
temp:=3.*(0.5*CK2*sqr(1/(a*(1-sqr(e)))))*sinio*(2*sqr(cosu)-1)-(0.25*A30*(1/(a*(1-sqr(e)))))*e*sin(omega);
deltaI:=temp*cosio;
deltaU:=sin(xincl/2)*((0.5*CK2*sqr(1/(a*(1-sqr(e)))))*(0.5*(1.-7.*sqr(theta))*(2*sinu*cosu)-
3.*(1-5*sqr(theta))*(f-M+e*sinf))-(0.25*A30*(1/(a*(1-sqr(e)))))*sinio*cosu*(2.+ (e*cosf)))-
0.5*(0.25*A30*(1/(a*(1-sqr(e)))))*sqr(theta)*e*cos(omega)/cos(xincl/2);
lambda:=f+omega+xnode+(0.5*CK2*sqr(1/(a*(1-sqr(e)))))*(0.5*(1.+6.*cosio-7.*sqr(theta))*(2*sinu*cosu)-3*((1-5*sqr(theta))+
2.*cosio)*(f-M+e*sinf))+(0.25*A30*(1/(a*(1-sqr(e)))))*sinio*(cosio*e*cos(omega)/(1.+cosio)-(2. +(e*cosf))*cosu);
Y4:=sin(xincl/2)*sinu+cosu*deltaU+0.5*sinu*cos(xincl/2)*deltaI;
Y5:=sin(xincl/2)*cosu-sinu*deltaU+0.5*cosu*cos(xincl/2)*deltaI;
r:=rii+deltaR;
rdot:=n*a*e*sinf/beta+(-n*sqr(a/rii))*(1*CK2*(1/(a*(1-sqr(e))))*(1-sqr(theta))*(2*sinu*cosu)+(0.25*A30*sinio)*cosu);
rfdot:=n*sqr(a)*beta/rii+(-n*sqr(a/rii))*deltaR+a*(n*(a/rii))*sinio*temp;
{Unit-orientation vector}
cosI2:=sqrt(1-sqr(y4)-sqr(y5));
UV[1]:=2*y4*(y5*sin(lambda)-y4*cos(lambda))+cos(lambda);
UV[2]:=-2*y4*(y5*cos(lambda)+y4*sin(lambda))+sin(lambda);
UV[3]:=2*y4*cosI2;
VV[1]:=2*y5*(y5*sin(lambda)-y4*cos(lambda))-sin(lambda);
VV[2]:=-2*y5*(y5*cos(lambda)+y4*sin(lambda))+cos(lambda);
VV[3]:=2*y5*cosI2;
for i:=1 to 3 do
begin;
pos[i]:=r*uv[i];
vel[i]:=rdot*uv[i]+rfdot*vv[i];
end;
end; {Procedure SGP8}
Procedure Deep(ideep : integer);
const
zns = 1.19459E-5; c1ss = 2.9864797E-6; zes = 0.01675;
znl = 1.5835218E-4; c1l = 4.7968065E-7; zel = 0.05490;
zcosis = 0.91744867; zsinis = 0.39785416; zsings = -0.98088458;
zcosgs = 0.1945905; zcoshs = 1; zsinhs = 0;
q22 = 1.7891679E-6; q31 = 2.1460748E-6; q33 = 2.2123015E-7;
g22 = 5.7686396; g32 = 0.95240898; g44 = 1.8014998;
g52 = 1.0508330; g54 = 4.4108898; root22 = 1.7891679E-6;
root32 = 3.7393792E-7; root44 = 7.3636953E-9; root52 = 1.1428639E-7;
root54 = 2.1765803E-9; thdt = 4.3752691E-3;
label
{dpinit}
5,10,20,30,40,45,50,55,60,65,70,80,
{dpsec}
90,100,105,110,120,125,130,140,150,152,154,160,165,170,175,180,
{dpper}
205,210,218,220,230;
const { Typed constants to retain values between ENTRY calls }
iresfl : integer = 0; isynfl : integer = 0;
iret : integer = 0; iretn : integer = 0;
ls : integer = 0;
a1 : double = 0; a2 : double = 0; a3 : double = 0;
a4 : double = 0; a5 : double = 0; a6 : double = 0;
a7 : double = 0; a8 : double = 0; a9 : double = 0;
a10 : double = 0; ainv2 : double = 0; alfdp : double = 0;
aqnv : double = 0; atime : double = 0; betdp : double = 0;
bfact : double = 0; c : double = 0; cc : double = 0;
cosis : double = 0; cosok : double = 0; cosq : double = 0;
ctem : double = 0; d2201 : double = 0; d2211 : double = 0;
d3210 : double = 0; d3222 : double = 0; d4410 : double = 0;
d4422 : double = 0; d5220 : double = 0; d5232 : double = 0;
d5421 : double = 0; d5433 : double = 0; dalf : double = 0;
day : double = 0; dbet : double = 0; del1 : double = 0;
del2 : double = 0; del3 : double = 0; delt : double = 0;
dls : double = 0; e3 : double = 0; ee2 : double = 0;
eoc : double = 0; eq : double = 0; f2 : double = 0;
f220 : double = 0; f221 : double = 0; f3 : double = 0;
f311 : double = 0; f321 : double = 0; f322 : double = 0;
f330 : double = 0; f441 : double = 0; f442 : double = 0;
f522 : double = 0; f523 : double = 0; f542 : double = 0;
f543 : double = 0; fasx2 : double = 0; fasx4 : double = 0;
fasx6 : double = 0; ft : double = 0; g200 : double = 0;
g201 : double = 0; g211 : double = 0; g300 : double = 0;
g310 : double = 0; g322 : double = 0; g410 : double = 0;
g422 : double = 0; g520 : double = 0; g521 : double = 0;
g532 : double = 0; g533 : double = 0; gam : double = 0;
omegaq : double = 0; pe : double = 0; pgh : double = 0;
ph : double = 0; pinc : double = 0; pl : double = 0;
preep : double = 0; s1 : double = 0; s2 : double = 0;
s3 : double = 0; s4 : double = 0; s5 : double = 0;
s6 : double = 0; s7 : double = 0; savtsn : double = 0;
se : double = 0; se2 : double = 0; se3 : double = 0;
sel : double = 0; ses : double = 0; sgh : double = 0;
sgh2 : double = 0; sgh3 : double = 0; sgh4 : double = 0;
sghl : double = 0; sghs : double = 0; sh : double = 0;
sh2 : double = 0; sh3 : double = 0; sh1 : double = 0;
shs : double = 0; si : double = 0; si2 : double = 0;
si3 : double = 0; sil : double = 0; sini2 : double = 0;
sinis : double = 0; sinok : double = 0; sinq : double = 0;
sinzf : double = 0; sis : double = 0; sl : double = 0;
sl2 : double = 0; sl3 : double = 0; sl4 : double = 0;
sll : double = 0; sls : double = 0; sse : double = 0;
ssg : double = 0; ssh : double = 0; ssi : double = 0;
ssl : double = 0; stem : double = 0; step2 : double = 0;
stepn : double = 0; stepp : double = 0; temp : double = 0;
temp1 : double = 0; thgr : double = 0; x1 : double = 0;
x2 : double = 0; x2li : double = 0; x2omi : double = 0;
x3 : double = 0; x4 : double = 0; x5 : double = 0;
x6 : double = 0; x7 : double = 0; x8 : double = 0;
xfact : double = 0; xgh2 : double = 0; xgh3 : double = 0;
xgh4 : double = 0; xh2 : double = 0; xh3 : double = 0;
xi2 : double = 0; xi3 : double = 0; xl : double = 0;
xl2 : double = 0; xl3 : double = 0; xl4 : double = 0;
xlamo : double = 0; xldot : double = 0; xli : double = 0;
xls : double = 0; xmao : double = 0; xnddt : double = 0;
xndot : double = 0; xni : double = 0; xno2 : double = 0;
xnodce : double = 0; xnoi : double = 0; xnq : double = 0;
xomi : double = 0; xpidot : double = 0; xqncl : double = 0;
z1 : double = 0; z11 : double = 0; z12 : double = 0;
z13 : double = 0; z2 : double = 0; z21 : double = 0;
z22 : double = 0; z23 : double = 0; z3 : double = 0;
z31 : double = 0; z32 : double = 0; z33 : double = 0;
zcosg : double = 0; zcosgl : double = 0; zcosh : double = 0;
zcoshl : double = 0; zcosi : double = 0; zcosil : double = 0;
ze : double = 0; zf : double = 0; zm : double = 0;
zmo : double = 0; zmol : double = 0; zmos : double = 0;
zn : double = 0; zsing : double = 0; zsingl : double = 0;
zsinh : double = 0; zsinhl : double = 0; zsini : double = 0;
zsinil : double = 0; zx : double = 0; zy : double = 0;
begin
case ideep of
dpinit : begin { Entrance for deep space initialization }
thgr := Thetag(epoch);
eq := eo;
xnq := xnodp;
aqnv := 1/ao;
xqncl := xincl;
xmao := xmo;
xpidot := omgdt + xnodot;
sinq := Sin(xnodeo);
cosq := Cos(xnodeo);
omegaq := omegao;
{ Initialize lunar solar terms }
5: day := ds50 + 18261.5; {Days since 1900 Jan 0.5}
if (day = preep) then
Goto 10;
preep := day;
xnodce := 4.5236020 - 9.2422029E-4*day;
stem := Sin(xnodce);
ctem := Cos(xnodce);
zcosil := 0.91375164 - 0.03568096*ctem;
zsinil := Sqrt(1 - zcosil*zcosil);
zsinhl := 0.089683511*stem/zsinil;
zcoshl := Sqrt(1 - zsinhl*zsinhl);
c := 4.7199672 + 0.22997150*day;
gam := 5.8351514 + 0.0019443680*day;
zmol := fmod2p(c - gam);
zx := 0.39785416*stem/zsinil;
zy := zcoshl*ctem + 0.91744867*zsinhl*stem;
zx := Actan(zx,zy);
zx := gam + zx - xnodce;
zcosgl := Cos(zx);
zsingl := Sin(zx);
zmos := 6.2565837 + 0.017201977*day;
zmos := fmod2p(zmos);
{ Do solar terms }
10: savtsn := 1E20;
zcosg := zcosgs;
zsing := zsings;
zcosi := zcosis;
zsini := zsinis;
zcosh := cosq;
zsinh := sinq;
cc := c1ss;
zn := zns;
ze := zes;
zmo := zmos;
xnoi := 1/xnq;
ls := 30; {assign 30 to ls}
20: a1 := zcosg*zcosh + zsing*zcosi*zsinh;
a3 := -zsing*zcosh + zcosg*zcosi*zsinh;
a7 := -zcosg*zsinh + zsing*zcosi*zcosh;
a8 := zsing*zsini;
a9 := zsing*zsinh + zcosg*zcosi*zcosh;
a10 := zcosg*zsini;
a2 := cosiq*a7 + siniq*a8;
a4 := cosiq*a9 + siniq*a10;
a5 := -siniq*a7 + cosiq*a8;
a6 := -siniq*a9 + cosiq*a10;
x1 := a1*cosomo + a2*sinomo;
x2 := a3*cosomo + a4*sinomo;
x3 := -a1*sinomo + a2*cosomo;
x4 := -a3*sinomo + a4*cosomo;
x5 := a5*sinomo;
x6 := a6*sinomo;
x7 := a5*cosomo;
x8 := a6*cosomo;
z31 := 12*x1*x1 - 3*x3*x3;
z32 := 24*x1*x2 - 6*x3*x4;
z33 := 12*x2*x2 - 3*x4*x4;
z1 := 3*(a1*a1 + a2*a2) + z31*eqsq;
z2 := 6*(a1*a3 + a2*a4) + z32*eqsq;
z3 := 3*(a3*a3 + a4*a4) + z33*eqsq;
z11 := -6*a1*a5 + eqsq*(-24*x1*x7 - 6*x3*x5);
z12 := -6*(a1*a6 + a3*a5)
+ eqsq*(-24*(x2*x7 + x1*x8) - 6*(x3*x6 + x4*x5));
z13 := -6*a3*a6 + eqsq*(-24*x2*x8 - 6*x4*x6);
z21 := 6*a2*a5 + eqsq*(24*x1*x5 - 6*x3*x7);
z22 := 6*(a4*a5 + a2*a6)
+ eqsq*(24*(x2*x5 + x1*x6) - 6*(x4*x7 + x3*x8));
z23 := 6*a4*a6 + eqsq*(24*x2*x6 - 6*x4*x8);
z1 := z1 + z1 + bsq*z31;
z2 := z2 + z2 + bsq*z32;
z3 := z3 + z3 + bsq*z33;
s3 := cc*xnoi;
s2 := -0.5*s3/rteqsq;
s4 := s3*rteqsq;
s1 := -15*eq*s4;
s5 := x1*x3 + x2*x4;
s6 := x2*x3 + x1*x4;
s7 := x2*x4 - x1*x3;
se := s1*zn*s5;
si := s2*zn*(z11 + z13);
sl := -zn*s3*(z1 + z3 - 14 - 6*eqsq);
sgh := s4*zn*(z31 + z33 - 6);
sh := -zn*s2*(z21 + z23);
if (xqncl < 5.2359877E-2) then
sh := 0;
ee2 := 2*s1*s6;
e3 := 2*s1*s7;
xi2 := 2*s2*z12;
xi3 := 2*s2*(z13 - z11);
xl2 := -2*s3*z2;
xl3 := -2*s3*(z3 - z1);
xl4 := -2*s3*(-21 - 9*eqsq)*ze;
xgh2 := 2*s4*z32;
xgh3 := 2*s4*(z33 - z31);
xgh4 := -18*s4*ze;
xh2 := -2*s2*z22;
xh3 := -2*s2*(z23 - z21);
case ls of
30 : Goto 30;
40 : Goto 40;
else
Halt;
end; {case}
{ Do lunar terms }
30: sse := se;
ssi := si;
ssl := sl;
ssh := sh/siniq;
ssg := sgh - cosiq*ssh;
se2 := ee2;
si2 := xi2;
sl2 := xl2;
sgh2 := xgh2;
sh2 := xh2;
se3 := e3;
si3 := xi3;
sl3 := xl3;
sgh3 := xgh3;
sh3 := xh3;
sl4 := xl4;
sgh4 := xgh4;
zcosg := zcosgl;
zsing := zsingl;
zcosi := zcosil;
zsini := zsinil;
zcosh := zcoshl*cosq + zsinhl*sinq;
zsinh := sinq*zcoshl - cosq*zsinhl;
zn := znl;
cc := c1l;
ze := zel;
zmo := zmol;
ls := 40; {assign 40 to ls}
Goto 20;
40: sse := sse + se;
ssi := ssi + si;
ssl := ssl + sl;
ssg := ssg + sgh - cosiq/siniq*sh;
ssh := ssh + sh/siniq;
{ Geopotential resonance initialization for 12 hour orbits }
iresfl := 0;
isynfl := 0;
if (xnq < 0.0052359877) and (xnq > 0.0034906585) then
Goto 70;
if (xnq < 8.26E-3) or (xnq > 9.24E-3) then
exit;
if (eq < 0.5) then
exit;
iresfl := 1;
eoc := eq*eqsq;
g201 := -0.306 - (eq - 0.64)*0.440;
if (eq > 0.65) then
Goto 45;
g211 := 3.616 - 13.247*eq + 16.290*eqsq;
g310 := -19.302 + 117.390*eq - 228.419*eqsq + 156.591*eoc;
g322 := -18.9068 + 109.7927*eq - 214.6334*eqsq + 146.5816*eoc;
g410 := -41.122 + 242.694*eq - 471.094*eqsq + 313.953*eoc;
g422 := -146.407 + 841.880*eq - 1629.014*eqsq + 1083.435*eoc;
g520 := -532.114 + 3017.977*eq - 5740*eqsq + 3708.276*eoc;
Goto 55;
45: g211 := -72.099 + 331.819*eq - 508.738*eqsq + 266.724*eoc;
g310 := -346.844 + 1582.851*eq - 2415.925*eqsq + 1246.113*eoc;
g322 := -342.585 + 1554.908*eq - 2366.899*eqsq + 1215.972*eoc;
g410 := -1052.797 + 4758.686*eq - 7193.992*eqsq + 3651.957*eoc;
g422 := -3581.69 + 16178.11*eq - 24462.77*eqsq + 12422.52*eoc;
if (eq > 0.715) then
Goto 50;
g520 := 1464.74 - 4664.75*eq + 3763.64*eqsq;
Goto 55;
50: g520 := -5149.66 + 29936.92*eq - 54087.36*eqsq + 31324.56*eoc;
55: if (eq >= (0.7)) then
Goto 60;
g533 := -919.2277 + 4988.61*eq - 9064.77*eqsq + 5542.21*eoc;
g521 := -822.71072 + 4568.6173*eq - 8491.4146*eqsq + 5337.524*eoc;
g532 := -853.666 + 4690.25*eq - 8624.77*eqsq + 5341.4*eoc;
Goto 65;
60: g533 := -37995.78 + 161616.52*eq - 229838.2*eqsq + 109377.94*eoc;
g521 := -51752.104 + 218913.95*eq - 309468.16*eqsq + 146349.42*eoc;
g532 := -40023.88 + 170470.89*eq - 242699.48*eqsq + 115605.82*eoc;
65: sini2 := siniq*siniq;
f220 := 0.75*(1 + 2*cosiq + cosq2);
f221 := 1.5*sini2;
f321 := 1.875*siniq*(1 - 2*cosiq - 3*cosq2);
f322 := -1.875*siniq*(1 + 2*cosiq - 3*cosq2);
f441 := 35*sini2*f220;
f442 := 39.3750*sini2*sini2;
f522 := 9.84375*siniq*(sini2*(1 - 2*cosiq - 5*cosq2)
+ 0.33333333*(-2 + 4*cosiq + 6*cosq2));
f523 := siniq*(4.92187512*sini2*(-2 - 4*cosiq + 10*cosq2)
+ 6.56250012*(1 + 2*cosiq - 3*cosq2));
f542 := 29.53125*siniq*(2 - 8*cosiq + cosq2*(-12 + 8*cosiq + 10*cosq2));
f543 := 29.53125*siniq*(-2 - 8*cosiq + cosq2*(12 + 8*cosiq - 10*cosq2));
xno2 := xnq*xnq;
ainv2 := aqnv*aqnv;
temp1 := 3*xno2*ainv2;
temp := temp1*root22;
d2201 := temp*f220*g201;
d2211 := temp*f221*g211;
temp1 := temp1*aqnv;
temp := temp1*root32;
d3210 := temp*f321*g310;
d3222 := temp*f322*g322;
temp1 := temp1*aqnv;
temp := 2*temp1*root44;
d4410 := temp*f441*g410;
d4422 := temp*f442*g422;
temp1 := temp1*aqnv;
temp := temp1*root52;
d5220 := temp*f522*g520;
d5232 := temp*f523*g532;
temp := 2*temp1*root54;
d5421 := temp*f542*g521;
d5433 := temp*f543*g533;
xlamo := xmao + xnodeo + xnodeo - thgr - thgr;
bfact := xlldot + xnodot + xnodot - thdt - thdt;
bfact := bfact + ssl + ssh + ssh;
Goto 80;
{ Synchronous resonance terms initialization }
70: iresfl := 1;
isynfl := 1;
g200 := 1 + eqsq*(-2.5 + 0.8125*eqsq);
g310 := 1 + 2*eqsq;
g300 := 1 + eqsq*(-6 + 6.60937*eqsq);
f220 := 0.75*(1 + cosiq)*(1 + cosiq);
f311 := 0.9375*siniq*siniq*(1 + 3*cosiq) - 0.75*(1 + cosiq);
f330 := 1 + cosiq;
f330 := 1.875*f330*f330*f330;
del1 := 3*xnq*xnq*aqnv*aqnv;
del2 := 2*del1*f220*g200*q22;
del3 := 3*del1*f330*g300*q33*aqnv;
del1 := del1*f311*g310*q31*aqnv;
fasx2 := 0.13130908;
fasx4 := 2.8843198;
fasx6 := 0.37448087;
xlamo := xmao + xnodeo + omegao - thgr;
bfact := xlldot + xpidot - thdt;
bfact := bfact + ssl + ssg + ssh;
80: xfact := bfact - xnq;
{ Initialize integrator }
xli := xlamo;
xni := xnq;
atime := 0;
stepp := 720;
stepn := -720;
step2 := 259200;
end; {dpinit}
dpsec : begin { Entrance for deep space secular effects }
xll := xll + ssl*t;
omgasm := omgasm + ssg*t;
xnodes := xnodes + ssh*t;
_em := eo + sse*t;
xinc := xincl + ssi*t;
if (xinc >= 0) then
Goto 90;
xinc := -xinc;
xnodes := xnodes + pi;
omgasm := omgasm - pi;
90: if (iresfl = 0) then
exit;
100: if (atime = 0) then
Goto 170;
if (t >= 0) and (atime < 0) then
Goto 170;
if (t < 0) and (atime >= 0) then
Goto 170;
105: if (Abs(t) >= Abs(atime)) then
Goto 120;
delt := stepp;
if (t >= 0) then
delt := stepn;
110: iret := 100; {assign 100 to iret}
Goto 160;
120: delt := stepn;
if (t > 0) then
delt := stepp;
125: if (Abs(t - atime) < stepp) then
Goto 130;
iret := 125; {assign 125 to iret}
Goto 160;
130: ft := t - atime;
iretn := 140; {assign 140 to iretn}
Goto 150;
140: xn := xni + xndot*ft + xnddt*ft*ft*0.5;
xl := xli + xldot*ft + xndot*ft*ft*0.5;
temp := -xnodes + thgr + t*thdt;
xll := xl - omgasm + temp;
if (isynfl = 0) then
xll := xl + temp + temp;
exit;
{ Dot terms calculated }
150: if (isynfl = 0) then
Goto 152;
xndot := del1*Sin(xli - fasx2) + del2*Sin(2*(xli - fasx4))
+ del3*Sin(3*(xli - fasx6));
xnddt := del1*Cos(xli - fasx2)
+ 2*del2*Cos(2*(xli - fasx4))
+ 3*del3*Cos(3*(xli - fasx6));
Goto 154;
152: xomi := omegaq + omgdt*atime;
x2omi := xomi + xomi;
x2li := xli + xli;
xndot := d2201*Sin(x2omi + xli - g22)
+ d2211*Sin(xli - g22)
+ d3210*Sin(xomi + xli - g32)
+ d3222*Sin(-xomi + xli - g32)
+ d4410*Sin(x2omi + x2li - g44)
+ d4422*Sin(x2li - g44)
+ d5220*Sin(xomi + xli - g52)
+ d5232*Sin(-xomi + xli - g52)
+ d5421*Sin(xomi + x2li - g54)
+ d5433*Sin(-xomi + x2li - g54);
xnddt := d2201*Cos(x2omi + xli - g22)
+ d2211*Cos(xli - g22)
+ d3210*Cos(xomi + xli - g32)
+ d3222*Cos(-xomi + xli - g32)
+ d5220*Cos(xomi + xli - g52)
+ d5232*Cos(-xomi + xli - g52)
+ 2*(d4410*Cos(x2omi + x2li - g44)
+ d4422*Cos(x2li - g44)
+ d5421*Cos(xomi + x2li - g54)
+ d5433*Cos(-xomi + x2li - g54));
154: xldot := xni + xfact;
xnddt := xnddt*xldot;
case iretn of
140 : Goto 140;
165 : Goto 165;
else
Halt;
end; {case}
{ Integrator }
160: iretn := 165; {assign 165 to iretn}
Goto 150;
165: xli := xli + xldot*delt + xndot*step2;
xni := xni + xndot*delt + xnddt*step2;
atime := atime + delt;
case iret of
100 : Goto 100;
125 : Goto 125;
else
Halt;
end; {case}
{ Epoch restart }
170: if (t >= 0) then
Goto 175;
delt := stepn;
Goto 180;
175: delt := stepp;
180: atime := 0;
xni := xnq;
xli := xlamo;
Goto 125;
end; {dpsec}
dpper : begin { Entrance for lunar-solar periodics }
sinis := Sin(xinc);
cosis := Cos(xinc);
if (Abs(savtsn - t) < 30) then
Goto 210;
savtsn := t;
zm := zmos + zns*t;
205: zf := zm + 2*zes*Sin(zm);
sinzf := Sin(zf);
f2 := 0.5*sinzf*sinzf - 0.25;
f3 := -0.5*sinzf*Cos(zf);
ses := se2*f2 + se3*f3;
sis := si2*f2 + si3*f3;
sls := sl2*f2 + sl3*f3 + sl4*sinzf;
sghs := sgh2*f2 + sgh3*f3 + sgh4*sinzf;
shs := sh2*f2 + sh3*f3;
zm := zmol + znl*t;
zf := zm + 2*zel*Sin(zm);
sinzf := Sin(zf);
f2 := 0.5*sinzf*sinzf - 0.25;
f3 := -0.5*sinzf*Cos(zf);
sel := ee2*f2 + e3*f3;
sil := xi2*f2 + xi3*f3;
sll := xl2*f2 + xl3*f3 + xl4*sinzf;
sghl := xgh2*f2 + xgh3*f3 + xgh4*sinzf;
sh1 := xh2*f2 + xh3*f3;
pe := ses + sel;
pinc := sis + sil;
pl := sls + sll;
210: pgh := sghs + sghl;
ph := shs + sh1;
xinc := xinc + pinc;
_em := _em + pe;
if (xqncl < 0.2) then
Goto 220;
Goto 218;
{ Apply periodics directly }
218: ph := ph/siniq;
pgh := pgh - cosiq*ph;
omgasm := omgasm + pgh;
xnodes := xnodes + ph;
xll := xll + pl;
Goto 230;
{ Apply periodics with Lyddane modification }
220: sinok := Sin(xnodes);
cosok := Cos(xnodes);
alfdp := sinis*sinok;
betdp := sinis*cosok;
dalf := ph*cosok + pinc*cosis*sinok;
dbet := -ph*sinok + pinc*cosis*cosok;
alfdp := alfdp + dalf;
betdp := betdp + dbet;
xls := xll + omgasm + cosis*xnodes;
dls := pl + pgh - pinc*xnodes*sinis;
xls := xls + dls;
xnodes := Actan(alfdp,betdp);
xll := xll + pl;
omgasm := xls - xll - Cos(xinc)*xnodes;
230: end; {dpper}
end; {case}
end; {Procedure Deep}
Procedure Call_dpinit(var eosq,sinio,cosio,betao,aodp,theta2,
sing,cosg,betao2,xmdot,omgdot,
xnodott,xnodpp : double);
begin
eqsq := eosq; siniq := sinio; cosiq := cosio; rteqsq := betao;
ao := aodp; cosq2 := theta2; sinomo := sing; cosomo := cosg;
bsq := betao2; xlldot := xmdot; omgdt := omgdot; xnodot := xnodott;
xnodp := xnodpp;
Deep(1);
eosq := eqsq; sinio := siniq; cosio := cosiq; betao := rteqsq;
aodp := ao; theta2 := cosq2; sing := sinomo; cosg := cosomo;
betao2 := bsq; xmdot := xlldot; omgdot := omgdt; xnodott := xnodot;
xnodpp := xnodp;
end; {Procedure Call_dpinit}
Procedure Call_dpsec(var xmdf,omgadf,xnode,emm,xincc,xnn,tsince : double);
begin
xll := xmdf; omgasm := omgadf; xnodes := xnode; {_em := emm;
xinc := xincc;} xn := xnn; t := tsince;
Deep(2);
xmdf := xll; omgadf := omgasm; xnode := xnodes; emm := _em;
xincc := xinc; xnn := xn; tsince := t;
end; {Procedure Call_dpsec}
Procedure Call_dpper(var e,xincc,omgadf,xnode,xmam : double);
begin
_em := e; xinc := xincc; omgasm := omgadf; xnodes := xnode;
xll := xmam;
Deep(3);
e := _em; xincc := xinc; omgadf := omgasm; xnode := xnodes;
xmam := xll;
end; {Procedure Call_dpper}
Procedure SDP8(tsince : double; var iflag : integer; var pos,vel : vector);
const rho: double=0.15696615;
var a1, cosi, theta2, tthmun, eosq, betao2, betao, del1, a0, delo, aodp, xnodp,b : double;
po, pom2, sini, sing, cosg, temp, sinio2, cosio2, theta4, unm5th, unmth2, a3cof : double;
pardt1, pardt2, pardt4, xmdt1, xgdt1, xhdt1, xlldot, omgdt,xnodot, tsi, eta, eta2,psim2, alpha2, eeta : double;
z1,z7,xmamdf,omgasm,xnodes, xn, em, xmam, zc2,sine,cose,zc5,cape : double;
cos2g, d1, d2, d3, d4, d5, b1, b2, b3, c0, c1, c4, c5, xndt, xndtn, edot : double;
am, beta2m,SINOS, cosos, AXNM, AYNM, PM, G1, G2,G3,G4, G5, BETA, SNF,CSF,FM : double;
dr, diwc, di, sni2du, xlamb, y4, y5, r, rdot, rvdot, snlamb, cslamb, ux, vx, uy, vy, uz, vz : double;
snfg, csfg, sn2f2g, cs2f2g, ecosf, g10, rm, aovr, g13, g14, sini2 : double;
i : integer;
label EndKepler;
begin;
A1:=power((XKE/XNO),TOTHRD);
COSI:=COS(XINCL);
THETA2:=COSI*COSI;
TTHMUN:=3*THETA2-1;
EOSQ:=EO*EO;
BETAO2:=1.-EOSQ;
BETAO:=SQRT(BETAO2);
DEL1:=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2);
AO:=A1*(1-DEL1*(.5*TOTHRD+DEL1*(1+134/81*DEL1)));
DELO:=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2);
AODP:=AO/(1.-DELO);
XNODP:=XNO/(1.+DELO);
B:=2.*BSTAR/RHO;
PO:=AODP*BETAO2;
POM2:=1./(PO*PO);
SINI:=SIN(XINCL);
SING:=SIN(OMEGAO);
COSG:=COS(OMEGAO);
TEMP:=0.5*XINCL;
SINIO2:=SIN(TEMP);
COSIO2:=COS(TEMP);
THETA4:=sqr(THETA2);
UNM5TH:=1-5*THETA2;
UNMTH2:=1-THETA2;
A3COF:=-XJ3/CK2*power(AE,3);
PARDT1:=3.*CK2*POM2*XNODP;
PARDT2:=PARDT1*CK2*POM2;
PARDT4:=1.25*CK4*POM2*POM2*XNODP;
XMDT1:=0.5*PARDT1*BETAO*TTHMUN;
XGDT1:=-0.5*PARDT1*UNM5TH;
XHDT1:=-PARDT1*COSI;
XLLDOT:=XNODP+XMDT1+0.0625*PARDT2*BETAO*(13-78*THETA2+137*THETA4);
OMGDT:=XGDT1+0.0625*PARDT2*(7-114*THETA2+395*THETA4)+PARDT4*(3-36*THETA2+49*THETA4);
XNODOT:=XHDT1+(0.5*PARDT2*(4-19*THETA2)+2*PARDT4*(3-7*THETA2))*COSI;
TSI:=1/(PO-S);
ETA:=EO*S*TSI;
ETA2:=sqr(ETA);
PSIM2:=ABS(1./(1.-ETA2));
ALPHA2:=1.+EOSQ;
EETA:=EO*ETA;
COS2G:=2*sqr(COSG)-1;
D5:=TSI*PSIM2;
D1:=D5/PO;
D2:=12.+ETA2*(36.+4.5*ETA2);
D3:=ETA2*(15.+2.5*ETA2);
D4:=ETA*(5.+3.75*ETA2);
B1:=CK2*TTHMUN;
B2:=-CK2*UNMTH2;
B3:=A3COF*SINI;
C0:=0.5*B*RHO*QOMS2T*XNODP*AODP*power(TSI,4)*power(PSIM2,3.5)/SQRT(ALPHA2);
C1:=1.5*XNODP*sqr(ALPHA2)*C0;
C4:=D1*D3*B2;
C5:=D5*D4*B3;
XNDT:=C1*((2+ETA2*(3+34*EOSQ)+5*EETA*(4+ETA2)+8.5*EOSQ)+D1*D2*B1+ C4*COS2G+C5*SING);
XNDTN:=XNDT/XNODP;
EDOT:=-TOTHRD*XNDTN*(1.-EO);
Call_DpINIT(EOSQ,SINI,COSI,BETAO,AODP,THETA2,SING,COSG,BETAO2,XLLDOT,OMGDT,XNODOT,XNODP);
Z1:=0.5*XNDT*TSINCE*TSINCE;
Z7:=3.5*TOTHRD*Z1/XNODP;
XMAMDF:=XMO+XLLDOT*TSINCE;
OMGASM:=OMEGAO+OMGDT*TSINCE+Z7*XGDT1;
XNODES:=XNODEO+XNODOT*TSINCE+Z7*XHDT1;
XN:=XNODP;
CALL_DPSEC(XMAMDF,OMGASM,XNODES,EM,XINC,XN,TSINCE);
XN:=XN+XNDT*TSINCE;
EM:=EM+EDOT*TSINCE;
XMAM:=XMAMDF+Z1+Z7*XMDT1;
CALL_DPPER(EM,XINC,OMGASM,XNODES,XMAM);
XMAM:=FMOD2P(XMAM);
ZC2:=XMAM+EM*SIN(XMAM)*(1+EM*COS(XMAM));
for I:=1 to 10 do
begin;
SINE:=SIN(ZC2);
COSE:=COS(ZC2);
ZC5:=1/(1-EM*COSE);
CAPE:=(XMAM+EM*SINE-ZC2)*ZC5+ZC2;
IF(ABS(CAPE-ZC2) < E6A) then GOTO Endkepler;
ZC2:=CAPE;
end;
EndKepler:
AM:=power((XKE/XN),TOTHRD);
BETA2M:=1.-EM*EM;
SINOS:=SIN(OMGASM);
COSOS:=COS(OMGASM);
AXNM:=EM*COSOS;
AYNM:=EM*SINOS;
PM:=AM*BETA2M;
G1:=1./PM;
G2:=0.5*CK2*G1;
G3:=G2*G1;
BETA:=SQRT(BETA2M);
G4:=0.25*A3COF*SINI;
G5:=0.25*A3COF*G1;
SNF:=BETA*SINE*ZC5;
CSF:=(COSE-EM)*ZC5;
FM:=ACTAN(SNF,CSF);
SNFG:=SNF*COSOS+CSF*SINOS;
CSFG:=CSF*COSOS-SNF*SINOS;
SN2F2G:=2.*SNFG*CSFG;
CS2F2G:=2.*sqr(CSFG)-1;
ECOSF:=EM*CSF;
G10:=FM-XMAM+EM*SNF;
RM:=PM/(1.+ECOSF);
AOVR:=AM/RM;
G13:=XN*AOVR;
G14:=-G13*AOVR;
DR:=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG;
DIWC:=3.*G3*SINI*CS2F2G-G5*AYNM;
DI:=DIWC*COSI;
SINI2:=SIN(0.5*XINC);
SNI2DU:=SINIO2*(G3*(0.5*(1-7*THETA2)*SN2F2G-3*UNM5TH*G10)-G5*SINI*CSFG*(2+ECOSF))-0.5*G5*THETA2*AXNM/COSIO2;
XLAMB:=FM+OMGASM+XNODES+G3*(0.5*(1+6*COSI-7*THETA2)*SN2F2G-3*(UNM5TH+2*COSI)*G10)+G5*SINI*(COSI*AXNM/(1+COSI)-
(2+ECOSF)*CSFG);
Y4:=SINI2*SNFG+CSFG*SNI2DU+0.5*SNFG*COSIO2*DI;
Y5:=SINI2*CSFG-SNFG*SNI2DU+0.5*CSFG*COSIO2*DI;
R:=RM+DR;
RDOT:=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG);
RVDOT:=XN*sqr(AM)*BETA/RM+G14*DR+AM*G13*SINI*DIWC;
SNLAMB:=SIN(XLAMB);
CSLAMB:=COS(XLAMB);
TEMP:=2.*(Y5*SNLAMB-Y4*CSLAMB);
UX:=Y4*TEMP+CSLAMB;
VX:=Y5*TEMP-SNLAMB;
TEMP:=2.*(Y5*CSLAMB+Y4*SNLAMB);
UY:=-Y4*TEMP+SNLAMB;
VY:=-Y5*TEMP+CSLAMB;
TEMP:=2.*SQRT(1.-Y4*Y4-Y5*Y5);
UZ:=Y4*TEMP;
VZ:=Y5*TEMP;
pos[1]:=R*UX;
pos[2]:=R*UY;
pos[3]:=R*UZ;
vel[1]:=RDOT*UX+RVDOT*VX;
vel[2]:=RDOT*UY+RVDOT*VY;
vel[3]:=RDOT*UZ+RVDOT*VZ;
end;
procedure SGP8call (time : double; var pos,vel : vector);
var
tsince : double;
begin
tsince := (time - julian_epoch) * xmnpda;
if ideep = 0 then
SGP8(tsince,iflag,pos,vel)
else
SDP8(tsince,iflag,pos,vel);
end; {Procedure SGP8call}
begin;
Define_Derived_Constants;
end.