This is sdp8.c in view mode; [Download] [Up]
/************************/
/* sdp8.c */
/* */
/* sdp8 orbital model */
/************************/
/***** description
*
* $Id: sdp8.c,v 1.3 1993/04/28 17:09:14 craig Exp $
*
*/
/***** modification history
*
* translated by f2c (version of 12 March 1993 7:07:21).
*
* $Log: sdp8.c,v $
* Revision 1.3 1993/04/28 17:09:14 craig
* changed atan2 call to matan2 call.
*
* Revision 1.2 1993/04/21 19:50:33 craig
* changed approximation to psi term to real number.
* changed ecnsts to pcnsts
*
* Revision 1.1 1993/04/02 18:04:25 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;
/******************/
/* SDP8 14 NOV 80 */
/******************/
int sdp8 (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 xinc, xlamb, xmam, xmamdf, 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, sini2;
double cose, sine, cosg, sing, cosos, sinos;
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;
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. + 1.)));
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 */
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 = sqrt (1. - eta2);
psim2 = 1. / psim2 / psim2;
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;
edot = -mcnsts.tothrd * xndtn * (1. - element.eo);
*iflag = 0;
dpinit (eosq, sini, cosi, betao, aodp, theta2, sing, cosg,
betao2, xlldot, omgdt, xnodot, xnodp);
}
/* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG */
z1 = xndt * .5 * tsince * tsince;
z7 = mcnsts.tothrd * 3.5 * z1 / xnodp;
xmamdf = element.xmo + xlldot * tsince;
omgasm = element.omegao + omgdt * tsince + z7 * xgdt1;
xnodes = element.xnodeo + xnodot * tsince + z7 * xhdt1;
xn = xnodp;
dpsec (&xmamdf, &omgasm, &xnodes, &em, &xinc, &xn, tsince);
xn += xndt * tsince;
em += edot * tsince;
xmam = xmamdf + z1 + z7 * xmdt1;
dpper (&em, &xinc, &omgasm, &xnodes, &xmam);
xmam = fmod2p (xmam);
/* 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;
sini2 = sin (xinc * .5);
/* 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 = sini2 * snfg + csfg * sni2du + snfg * .5 * cosio2 * di;
y5 = sini2 * 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.