This is sgp.c in view mode; [Download] [Up]
/************************/
/* sgp.c */
/* */
/* sgp orbital model */
/************************/
/***** description
*
* $Id: sgp.c,v 1.3 1993/04/28 17:14:02 craig Exp $
*
*/
/***** modification history
*
* translated by f2c (version of 12 March 1993 7:07:21).
*
* $Log: sgp.c,v $
* Revision 1.3 1993/04/28 17:14:02 craig
* changed atan2 call to matan2 call.
*
* Revision 1.2 1993/04/21 19:51:53 craig
* changed ecnsts to pcnsts.
*
* Revision 1.1 1993/04/02 18:05:57 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;
/*****************/
/* SGP 31 OCT 80 */
/*****************/
int sgp (int *iflag, double tsince)
{
/* Local variables */
double axnsl, aynsl, ecose, el2, eo1, esine;
double omgas, pl2, po2no, rdot, rvdot, temp, tem2, tem5;
double xinck, xls, xmx, xmy, xnodek, xnodes;
double cosik, sinik, cosio, sinio, cosu, sinu, cos2u, sin2u;
double coseo1, sineo1, sinuk, cosuk, cosnok, sinnok;
double a, e, p, r, u;
double a1, c1, c2, c3, c4, d1;
double pl, po, rk, su, uk, ux, uy, uz, vx, vy, vz, xl;
register int i;
static double ao, c5, c6, omgdt, qo, xlo, xnodot;
static double d1o, d2o, d3o, d4o;
if (*iflag != 0)
{
/* INITIALIZATION */
c1 = pcnsts.ck2 * 1.5;
c2 = pcnsts.ck2 / 4.;
c3 = pcnsts.ck2 / 2.;
c4 = pcnsts.xj3 * pow (pcnsts.ae, 3.) / (pcnsts.ck2 * 4.);
cosio = cos (element.xincl);
sinio = sin (element.xincl);
a1 = pow (pcnsts.xke / element.xno, mcnsts.tothrd);
d1 = c1 / a1 / a1 * (cosio * 3. * cosio - 1.) /
pow ((1. - element.eo * element.eo), 2.5);
ao = a1 * (1. - d1 * 1. / 3. - d1 * d1 - d1 *
134. / 81. * d1 * d1);
po = ao * (1. - element.eo * element.eo);
qo = ao * (1. - element.eo);
xlo = element.xmo + element.omegao + element.xnodeo;
d1o = c3 * sinio * sinio;
d2o = c2 * (cosio * 7. * cosio - 1.);
d3o = c1 * cosio;
d4o = d3o * sinio;
po2no = element.xno / (po * po);
omgdt = c1 * po2no * (cosio * 5. * cosio - 1.);
xnodot = d3o * -2. * po2no;
c5 = c4 * .5 * sinio * (cosio * 5. + 3.) / (cosio + 1.);
c6 = c4 * sinio;
*iflag = 0;
}
/* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG */
a = element.xno + (element.xndt2o * 2. + element.xndd6o * 3.
* tsince) * tsince;
a = ao * pow (element.xno / a, mcnsts.tothrd);
e = mcnsts.e6a;
if (a > qo)
{
e = 1. - qo / a;
}
p = a * (1. - e * e);
xnodes = element.xnodeo + xnodot * tsince;
omgas = element.omegao + omgdt * tsince;
xls = fmod2p (xlo + (element.xno + omgdt + xnodot +
(element.xndt2o + element.xndd6o * tsince) * tsince) * tsince);
/* LONG PERIOD PERIODICS */
axnsl = e * cos (omgas);
aynsl = e * sin (omgas) - c6 / p;
xl = fmod2p (xls - c5 / p * axnsl);
/* SOLVE KEPLERS EQUATION */
u = fmod2p (xl - xnodes);
eo1 = u;
tem5 = 1.;
for (i = 0; i < 11; i++)
{
sineo1 = sin (eo1);
coseo1 = cos (eo1);
if (fabs (tem5) < mcnsts.e6a)
{
break;
}
tem5 = 1. - coseo1 * axnsl - sineo1 * aynsl;
tem5 = (u - aynsl * coseo1 + axnsl * sineo1 - eo1) / tem5;
tem2 = fabs (tem5);
if (tem2 > 1.)
{
tem5 = tem2 / tem5;
}
eo1 += tem5;
}
sineo1 = sin (eo1);
coseo1 = cos (eo1);
/* SHORT PERIOD PRELIMINARY QUANTITIES */
ecose = axnsl * coseo1 + aynsl * sineo1;
esine = axnsl * sineo1 - aynsl * coseo1;
el2 = axnsl * axnsl + aynsl * aynsl;
pl = a * (1. - el2);
pl2 = pl * pl;
r = a * (1. - ecose);
rdot = pcnsts.xke * sqrt (a) / r * esine;
rvdot = pcnsts.xke * sqrt (pl) / r;
temp = esine / (sqrt (1. - el2) + 1.);
sinu = a / r * (sineo1 - aynsl - axnsl * temp);
cosu = a / r * (coseo1 - axnsl + aynsl * temp);
su = matan2 (sinu, cosu);
/* UPDATE FOR SHORT PERIODICS */
sin2u = (cosu + cosu) * sinu;
cos2u = 1. - sinu * 2. * sinu;
rk = r + d1o / pl * cos2u;
uk = su - d2o / pl2 * sin2u;
xnodek = xnodes + d3o * sin2u / pl2;
xinck = element.xincl + d4o / pl2 * cos2u;
/* ORIENTATION VECTORS */
sinuk = sin (uk);
cosuk = cos (uk);
sinnok = sin (xnodek);
cosnok = cos (xnodek);
sinik = sin (xinck);
cosik = cos (xinck);
xmx = -sinnok * cosik;
xmy = cosnok * cosik;
ux = xmx * sinuk + cosnok * cosuk;
uy = xmy * sinuk + sinnok * cosuk;
uz = sinik * sinuk;
vx = xmx * cosuk - cosnok * sinuk;
vy = xmy * cosuk - sinnok * sinuk;
vz = sinik * cosuk;
/* POSITION AND VELOCITY */
element.x = rk * ux;
element.y = rk * uy;
element.z = rk * uz;
element.xdot = rdot * ux;
element.ydot = rdot * uy;
element.zdot = rdot * uz;
element.xdot = rvdot * vx + element.xdot;
element.ydot = rvdot * vy + element.ydot;
element.zdot = rvdot * vz + element.zdot;
return (0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.