This is sgp8.c in view mode; [Download] [Up]
/************************/
/* sgp8.c */
/* */
/* sdp8 orbital model */
/************************/
/***** description
*
* $Id: sgp8.c,v 1.3 1993/04/28 17:12:25 craig Exp $
*
*/
/***** modification history
*
* translated by f2c (version of 12 March 1993 7:07:21).
*
* $Log: sgp8.c,v $
* Revision 1.3 1993/04/28 17:12:25 craig
* changed atan2 call to matan2 call
*
* Revision 1.2 1993/04/21 19:53:51 craig
* changed ecnsts to pcnsts.
*
* Revision 1.1 1993/04/02 18:05:34 craig
* Initial revision
*
*
*/
/***** include files *****/
#include <math.h>
#include "satellite.h"
#include "satproto.h"
#include "aaproto.h"
/***** global variables *****/
extern struct ELEMENT element;
extern struct PCONSTANTS pcnsts;
extern struct MCONSTANTS mcnsts;
/******************/
/* SGP8 14 NOV 80 */
/******************/
int sgp8 (int *iflag, double tsince)
{
/* Initialized data */
static double rho = .15696615;
/* Local variables */
double alpha2, aodp, aovr, axnm, aynm, beta, beta2m, betao, betao2;
double cape, cs2f2g, csf, cslamb, del1, delo, diwc;
double ecosf, eeta, eosq, eta, eta2;
double omgasm, pardt1, pardt2, pardt4, pom2, psim2, rdot, rvdot;
double sn2f2g, snf, sni2du, snlamb, temp, theta4, tsi;
double xlamb, xmam, xndtn, xnodes;
double b, r;
double am, ao, di, dr, em, fm, pm, po, rm, xn;
double a1, b1, b2, b3, c0, c1, c4, c5, d1, d2, d3, d4, d5;
double g1, g2, g3, g4, g5, g10, g13, g14, ux, uy, uz;
double vx, vy, vz, y4, y5, z1, z7;
double zc2, zc5;
double csfg, snfg, cos2g;
double cose, sine, cosg, sing, sinos, cosos;
double aldtal, c0dtc0, c1dtc1, eddot, etddt, etdt, psdtps;
double temp1, tmnddt, tsddts, tsdtts, xnddt, xntrdt;
double c8, c9, d6, d7, d8, d9, d10, d11, d12, d13, d14, d15;
double d16, d17, d18, d19, d20, d23, d25;
double c4dt, c5dt, d1dt, d2dt, d3dt, d4dt, d5dt;
double sin2g;
static int isimp;
register int i;
static double a3cof, edot, omgdt, theta2, tthmun, unmth2, unm5th;
static double xgdt1, xhdt1, xlldot, xmdt1, xndt, xnodp, xnodot;
static double cosi, sini, cosio2, sinio2;
static double d1ddt, ed, gamma, ovgpp, pp, qq, xnd;
if (*iflag != 0)
{
/* RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP) */
/* FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT */
/* (B TERM) FROM INPUT B* DRAG TERM */
a1 = pow (pcnsts.xke / element.xno, mcnsts.tothrd);
cosi = cos (element.xincl);
theta2 = cosi * cosi;
tthmun = theta2 * 3. - 1.;
eosq = element.eo * element.eo;
betao2 = 1. - eosq;
betao = sqrt (betao2);
del1 = pcnsts.ck2 * 1.5 * tthmun / (a1 * a1 * betao * betao2);
ao = a1 * (1. - del1 * (mcnsts.tothrd * .5 +
del1 * (del1 * 134. / 81.)));
delo = pcnsts.ck2 * 1.5 * tthmun / (ao * ao * betao * betao2);
aodp = ao / (1. - delo);
xnodp = element.xno / (delo + 1.);
b = element.bstar * 2. / rho;
/* INITIALIZATION */
isimp = 0;
po = aodp * betao2;
pom2 = 1. / (po * po);
sini = sin (element.xincl);
sing = sin (element.omegao);
cosg = cos (element.omegao);
temp = element.xincl * .5;
sinio2 = sin (temp);
cosio2 = cos (temp);
theta4 = theta2 * theta2;
unm5th = 1. - theta2 * 5.;
unmth2 = 1. - theta2;
a3cof = -pcnsts.xj3 / pcnsts.ck2 * pow (pcnsts.ae, 3.);
pardt1 = pcnsts.ck2 * 3. * pom2 * xnodp;
pardt2 = pardt1 * pcnsts.ck2 * pom2;
pardt4 = pcnsts.ck4 * 1.25 * pom2 * pom2 * xnodp;
xmdt1 = pardt1 * .5 * betao * tthmun;
xgdt1 = pardt1 * -.5 * unm5th;
xhdt1 = -pardt1 * cosi;
xlldot = xnodp + xmdt1 + pardt2 * .0625 * betao *
(13. - theta2 * 78. + theta4 * 137.);
omgdt = xgdt1 +
pardt2 * .0625 * (7. - theta2 * 114. + theta4 * 395.) +
pardt4 * (3. - theta2 * 36. + theta4 * 49.);
xnodot = xhdt1 + (pardt2 * .5 * (4. - theta2 * 19.) +
pardt4 * 2. * (3. - theta2 * 7.)) * cosi;
tsi = 1. / (po - pcnsts.s);
eta = element.eo * pcnsts.s * tsi;
eta2 = eta * eta;
psim2 = fabs (1. / (1. - eta2));
alpha2 = eosq + 1.;
eeta = element.eo * eta;
cos2g = cosg * cosg * 2. - 1.;
d5 = tsi * psim2;
d1 = d5 / po;
d2 = eta2 * (eta2 * 4.5 + 36.) + 12.;
d3 = eta2 * (eta2 * 2.5 + 15.);
d4 = eta * (eta2 * 3.75 + 5.);
b1 = pcnsts.ck2 * tthmun;
b2 = -pcnsts.ck2 * unmth2;
b3 = a3cof * sini;
c0 = b * .5 * rho * pcnsts.qoms2t * xnodp * aodp *
pow (tsi, 4.) * pow (psim2, 3.5) / sqrt (alpha2);
c1 = xnodp * 1.5 * (alpha2 * alpha2) * c0;
c4 = d1 * d3 * b2;
c5 = d5 * d4 * b3;
xndt = c1 * (eta2 * (eosq * 34. + 3.) + 2.
+ eeta * 5. * (eta2 + 4.) +
eosq * 8.5 + d1 * d2 * b1 +
c4 * cos2g + c5 * sing);
xndtn = xndt / xnodp;
/* IF DRAG IS VERY SMALL, THE ISIMP FLAG IS SET AND THE */
/* EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN MEAN */
/* MOTION AND QUADRATIC VARIATION IN MEAN ANOMALY */
if (fabs (xndtn * pcnsts.xmnpda) < .00216)
{
isimp = 1;
edot = -mcnsts.tothrd * xndtn * (1. - element.eo);
}
else
{
d6 = eta * (eta2 * 22.5 + 30.);
d7 = eta * (eta2 * 12.5 + 5.);
d8 = eta2 * (eta2 + 6.75) + 1.;
c8 = d1 * d7 * b2;
c9 = d5 * d8 * b3;
edot = -c0 * (eta * (eta2 + 4. + eosq * (eta2 * 7. + 15.5))
+ element.eo * (eta2 * 15. + 5.) + d1 * d6 * b1
+ c8 * cos2g + c9 * sing);
d20 = mcnsts.tothrd * .5 * xndtn;
aldtal = element.eo * edot / alpha2;
tsdtts = aodp * 2. * tsi * (d20 * betao2 + element.eo * edot);
etdt = (edot + element.eo * tsdtts) * tsi * pcnsts.s;
psdtps = -eta * etdt * psim2;
sin2g = sing * 2. * cosg;
c0dtc0 = d20 + tsdtts * 4. - aldtal - psdtps * 7.;
c1dtc1 = xndtn + aldtal * 4. + c0dtc0;
d9 = eta * (eosq * 68. + 6.) + element.eo * (eta2 * 15. + 20.);
d10 = eta * 5. * (eta2 + 4.) + element.eo * (eta2 * 68. + 17.);
d11 = eta * (eta2 * 18. + 72.);
d12 = eta * (eta2 * 10. + 30.);
d13 = eta2 * 11.25 + 5.;
d14 = tsdtts - psdtps * 2.;
d15 = (d20 + element.eo * edot / betao2) * 2.;
d1dt = d1 * (d14 + d15);
d2dt = etdt * d11;
d3dt = etdt * d12;
d4dt = etdt * d13;
d5dt = d5 * d14;
c4dt = b2 * (d1dt * d3 + d1 * d3dt);
c5dt = b3 * (d5dt * d4 + d5 * d4dt);
d16 = d9 * etdt + d10 * edot + b1 * (d1dt * d2 + d1 * d2dt)
+ c4dt * cos2g + c5dt * sing + xgdt1
* (c5 * cosg - c4 * 2. * sin2g);
xnddt = c1dtc1 * xndt + c1 * d16;
eddot = c0dtc0 * edot - c0 * ((eta2 * 3. + 4. +
eeta * 30. + eosq * (eta2 * 21. + 15.5)) * etdt +
(eta2 * 15. + 5. + eeta * (eta2 * 14. + 31.)) * edot +
b1 * (d1dt * d6 + d1 * etdt * (eta2 * 67.5 + 30.)) +
b2 * (d1dt * d7 + d1 * etdt * (eta2 * 37.5 + 5.)) * cos2g +
b3 * (d5dt * d8 + d5 * etdt * eta * (eta2 * 4. + 13.5)) *
sing + xgdt1 * (c9 * cosg - c8 * 2. * sin2g));
d25 = edot * edot;
d17 = xnddt / xnodp - xndtn * xndtn;
tsddts = tsdtts * 2. * (tsdtts - d20) + aodp * tsi *
(mcnsts.tothrd * betao2 * d17 - d20 * 4. * element.eo *
edot + (d25 + element.eo * eddot) * 2.);
etddt = (eddot + edot * 2. * tsdtts) * tsi * pcnsts.s +
tsddts * eta;
d18 = tsddts - tsdtts * tsdtts;
d19 = -(psdtps * psdtps) / eta2 - eta * etddt * psim2 -
psdtps * psdtps;
d23 = etdt * etdt;
d1ddt = d1dt * (d14 + d15) + d1 * (d18 - d19 * 2. +
mcnsts.tothrd * d17 + (alpha2 * d25 / betao2 +
element.eo * eddot) * 2. / betao2);
xntrdt = xndt * (mcnsts.tothrd * 2. * d17 + (d25 +
element.eo * eddot) * 3. / alpha2 - aldtal * aldtal * 6.
+ d18 * 4. - d19 * 7.) + c1dtc1 * xnddt + c1 *
(c1dtc1 * d16 + d9 * etddt + d10 * eddot + d23 *
(eeta * 30. + 6. + eosq * 68.) + etdt * edot *
(eta2 * 30. + 40. + eeta * 272.) + d25 *
(eta2 * 68. + 17.) + b1 * (d1ddt * d2 + d1dt * 2. * d2dt +
d1 * (etddt * d11 + d23 * (eta2 * 54. + 72.))) +
b2 * (d1ddt * d3 + d1dt * 2. * d3dt + d1 * (etddt * d12 +
d23 * (eta2 * 30. + 30.))) * cos2g + b3 *
((d5dt * d14 + d5 * (d18 - d19 * 2.)) * d4 +
d4dt * 2. * d5dt + d5 * (etddt * d13 +
eta * 22.5 * d23)) * sing + xgdt1 *
((d20 * 7. + element.eo * 4. * edot / betao2) *
(c5 * cosg - c4 * 2. * sin2g) + (c5dt * 2. * cosg -
c4dt * 4. * sin2g - xgdt1 * (c5 * sing +
c4 * 4. * cos2g))));
tmnddt = xnddt * 1.0e9;
temp = tmnddt * tmnddt - xndt * 1.0e18 * xntrdt;
pp = (temp + tmnddt * tmnddt) / temp;
gamma = -xntrdt / (xnddt * (pp - 2.));
xnd = xndt / (pp * gamma);
qq = 1. - eddot / (edot * gamma);
ed = edot / (qq * gamma);
ovgpp = 1. / (gamma * (pp + 1.));
}
*iflag = 0;
}
/* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG */
xmam = fmod2p (element.xmo + xlldot * tsince);
omgasm = element.omegao + omgdt * tsince;
xnodes = element.xnodeo + xnodot * tsince;
if (isimp == 1)
{
xn = xnodp + xndt * tsince;
em = element.eo + edot * tsince;
z1 = xndt * .5 * tsince * tsince;
}
else
{
temp = 1. - gamma * tsince;
temp1 = pow (temp, pp);
xn = xnodp + xnd * (1. - temp1);
em = element.eo + ed * (1. - pow (temp, qq));
z1 = xnd * (tsince + ovgpp * (temp * temp1 - 1.));
}
z7 = mcnsts.tothrd * 3.5 * z1 / xnodp;
xmam = fmod2p (xmam + z1 + z7 * xmdt1);
omgasm += z7 * xgdt1;
xnodes += z7 * xhdt1;
/* SOLVE KEPLERS EQUATION */
zc2 = xmam + em * sin (xmam) * (em * cos (xmam) + 1.);
for (i = 1; i <= 10; ++i)
{
sine = sin (zc2);
cose = cos (zc2);
zc5 = 1. / (1. - em * cose);
cape = (xmam + em * sine - zc2) * zc5 + zc2;
if (fabs (cape - zc2) <= mcnsts.e6a)
{
break;
}
zc2 = cape;
}
/* SHORT PERIOD PRELIMINARY QUANTITIES */
am = pow (pcnsts.xke / xn, mcnsts.tothrd);
beta2m = 1. - em * em;
sinos = sin (omgasm);
cosos = cos (omgasm);
axnm = em * cosos;
aynm = em * sinos;
pm = am * beta2m;
g1 = 1. / pm;
g2 = pcnsts.ck2 * .5 * g1;
g3 = g2 * g1;
beta = sqrt (beta2m);
g4 = a3cof * .25 * sini;
g5 = a3cof * .25 * g1;
snf = beta * sine * zc5;
csf = (cose - em) * zc5;
fm = matan2 (snf, csf);
snfg = snf * cosos + csf * sinos;
csfg = csf * cosos - snf * sinos;
sn2f2g = snfg * 2. * csfg;
cs2f2g = csfg * csfg * 2. - 1.;
ecosf = em * csf;
g10 = fm - xmam + em * snf;
rm = pm / (ecosf + 1.);
aovr = am / rm;
g13 = xn * aovr;
g14 = -g13 * aovr;
dr = g2 * (unmth2 * cs2f2g - tthmun * 3.) - g4 * snfg;
diwc = g3 * 3. * sini * cs2f2g - g5 * aynm;
di = diwc * cosi;
/* UPDATE FOR SHORT PERIOD PERIODICS */
sni2du = sinio2 * (g3 * ((1. - theta2 * 7.) * .5 * sn2f2g - unm5th *
3. * g10) - g5 * sini * csfg * (ecosf + 2.)) - g5 * .5 *
theta2 * axnm / cosio2;
xlamb = fm + omgasm + xnodes + g3 * ((cosi * 6. + 1. - theta2 * 7.) *
.5 * sn2f2g - (unm5th + cosi * 2.) * 3. * g10) + g5 * sini *
(cosi * axnm / (cosi + 1.) - (ecosf + 2.) * csfg);
y4 = sinio2 * snfg + csfg * sni2du + snfg * .5 * cosio2 * di;
y5 = sinio2 * csfg - snfg * sni2du + csfg * .5 * cosio2 * di;
r = rm + dr;
rdot = xn * am * em * snf / beta + g14 *
(g2 * 2. * unmth2 * sn2f2g + g4 * csfg);
rvdot = xn * (am * am) * beta / rm + g14 * dr +
am * g13 * sini * diwc;
/* ORIENTATION VECTORS */
snlamb = sin (xlamb);
cslamb = cos (xlamb);
temp = (y5 * snlamb - y4 * cslamb) * 2.;
ux = y4 * temp + cslamb;
vx = y5 * temp - snlamb;
temp = (y5 * cslamb + y4 * snlamb) * 2.;
uy = -y4 * temp + snlamb;
vy = -y5 * temp + cslamb;
temp = sqrt (1. - y4 * y4 - y5 * y5) * 2.;
uz = y4 * temp;
vz = y5 * temp;
/* POSITION AND VELOCITY */
element.x = r * ux;
element.y = r * uy;
element.z = r * uz;
element.xdot = rdot * ux + rvdot * vx;
element.ydot = rdot * uy + rvdot * vy;
element.zdot = rdot * uz + rvdot * vz;
return (0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.