This is sgp4.c in view mode; [Download] [Up]
/************************/ /* sgp4.c */ /* */ /* sgp4 orbital model */ /************************/ /***** description * * $Id: sgp4.c,v 1.3 1993/04/28 17:13:26 craig Exp $ * */ /***** modification history * * translated by f2c (version of 12 March 1993 7:07:21). * * $Log: sgp4.c,v $ * Revision 1.3 1993/04/28 17:13:26 craig * changed atan2 call to matan2 call. * * Revision 1.2 1993/04/21 19:52:56 craig * changed ecnsts to pcnsts * * Revision 1.1 1993/04/02 18:05:11 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; /******************/ /* SGP4 3 NOV 80 */ /******************/ int sgp4 (int *iflag, double tsince) { /* Local variables */ double a3ovk2, axn, ayn, aynl, beta, betal, betao, betao2; double capu, coef, coef1, ddummy, del1, delo; double ecose, eeta, elsq, eosq, epw, esine, etasq; double omgadf, perige, pinvsq, psisq, qoms24, rdot; double rdotk, rfdot, rfdotk, theta2, theta4, tsi, tsq; double x1m5th, xhdot1, xinck, xll, xlt, xmdf, xmx, xmy; double xnoddf, xnode, xnodek; double temp, temp1, temp2, temp3, temp4, temp5, temp6; double tempa, tempe, templ; double cosu, sinu, cos2u, sin2u, cosik, sinik, cosuk, sinuk; double cosnok, sinnok, cosepw, sinepw; double a, e, r, u; double a1, ao, c2, pl, rk, s4, uk, ux, uy, uz; double vx, vy, vz, xl, xn; double c1sq, c3, delm, delomg, omega, tcube, tfour, xmp; static int isimp; register int i; static double aycof, omgdot, x1mth2, x3thm1, x7thm1, xlcof, xmdot; static double xnodcf, xnodot, xnodp; static double c1, c4, cosio, sinio, t2cof; static double aodp, eta; static double c5, d2, d3, d4, delmo, omgcof, sinmo, xmcof; static double t3cof, t4cof, t5cof; if (*iflag != 0) { /* RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP) */ /* FROM INPUT ELEMENTS */ ddummy = pcnsts.xke / element.xno; a1 = pow (ddummy, mcnsts.tothrd); cosio = cos (element.xincl); theta2 = cosio * cosio; x3thm1 = theta2 * 3. - 1.; eosq = element.eo * element.eo; betao2 = 1. - eosq; betao = sqrt (betao2); del1 = pcnsts.ck2 * 1.5 * x3thm1 / (a1 * a1 * betao * betao2); ao = a1 * (1. - del1 * (mcnsts.tothrd * .5 + del1 * (del1 * 134. / 81. + 1.))); delo = pcnsts.ck2 * 1.5 * x3thm1 / (ao * ao * betao * betao2); xnodp = element.xno / (delo + 1.); aodp = ao / (1. - delo); /* INITIALIZATION */ /* FOR PERIGEE LESS THAN 220 KILOMETERS, THE ISIMP FLAG IS SET AND */ /* THE EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN SQRT A AND */ /* QUADRATIC VARIATION IN MEAN ANOMALY. ALSO, THE C3 TERM, THE */ /* DELTA OMEGA TERM, AND THE DELTA M TERM ARE DROPPED. */ if (aodp * (1. - element.eo) / pcnsts.ae < 220. / pcnsts.xkmper + pcnsts.ae) { isimp = 1; } else { isimp = 0; } s4 = pcnsts.s; qoms24 = pcnsts.qoms2t; perige = (aodp * (1. - element.eo) - pcnsts.ae) * pcnsts.xkmper; /* FOR PERIGEE BELOW 156 KM, THE VALUES OF */ /* S AND QOMS2T ARE ALTERED */ if (perige < 156.) { if (perige > 98.) { s4 = perige - 78.; } else { s4 = 20.; } qoms24 = pow ((120. - s4) * pcnsts.ae / pcnsts.xkmper, 4.); s4 = s4 / pcnsts.xkmper + pcnsts.ae; } pinvsq = 1. / (aodp * aodp * betao2 * betao2); tsi = 1. / (aodp - s4); eta = aodp * element.eo * tsi; etasq = eta * eta; eeta = element.eo * eta; psisq = fabs (1. - etasq); coef = qoms24 * pow (tsi, 4.); coef1 = coef / pow (psisq, 3.5); c2 = coef1 * xnodp * (aodp * (etasq * 1.5 + 1. + eeta * (etasq + 4.)) + pcnsts.ck2 * .75 * tsi / psisq * x3thm1 * (etasq * 3. * (etasq + 8.) + 8.)); c1 = element.bstar * c2; sinio = sin (element.xincl); a3ovk2 = -pcnsts.xj3 / pcnsts.ck2 * pow (pcnsts.ae, 3.); c3 = coef * tsi * a3ovk2 * xnodp * pcnsts.ae * sinio / element.eo; x1mth2 = 1. - theta2; c4 = xnodp * 2. * coef1 * aodp * betao2 * (eta * (etasq * .5 + 2.) + element.eo * (etasq * 2. + .5) - pcnsts.ck2 * 2. * tsi / (aodp * psisq) * (x3thm1 * -3. * (1. - eeta * 2. + etasq * (1.5 - eeta * .5)) + x1mth2 * .75 * (etasq * 2. - eeta * (etasq + 1.)) * cos (element.omegao * 2.))); c5 = coef1 * 2. * aodp * betao2 * ((etasq + eeta) * 2.75 + 1. + eeta * etasq); theta4 = theta2 * theta2; temp1 = pcnsts.ck2 * 3. * pinvsq * xnodp; temp2 = temp1 * pcnsts.ck2 * pinvsq; temp3 = pcnsts.ck4 * 1.25 * pinvsq * pinvsq * xnodp; xmdot = xnodp + temp1 * .5 * betao * x3thm1 + temp2 * .0625 * betao * (13. - theta2 * 78. + theta4 * 137.); x1m5th = 1. - theta2 * 5.; omgdot = temp1 * -.5 * x1m5th + temp2 * .0625 * (7. - theta2 * 114. + theta4 * 395.) + temp3 * (3. - theta2 * 36. + theta4 * 49.); xhdot1 = -temp1 * cosio; xnodot = xhdot1 + (temp2 * .5 * (4. - theta2 * 19.) + temp3 * 2. * (3. - theta2 * 7.)) * cosio; omgcof = element.bstar * c3 * cos (element.omegao); xmcof = -mcnsts.tothrd * coef * element.bstar * pcnsts.ae / eeta; xnodcf = betao2 * 3.5 * xhdot1 * c1; t2cof = c1 * 1.5; xlcof = a3ovk2 * .125 * sinio * (cosio * 5. + 3.) / (cosio + 1.); aycof = a3ovk2 * .25 * sinio; delmo = pow (1. + eta * cos (element.xmo), 3.); sinmo = sin (element.xmo); x7thm1 = theta2 * 7. - 1.; if (isimp != 1) { c1sq = c1 * c1; d2 = aodp * 4. * tsi * c1sq; temp = d2 * tsi * c1 / 3.; d3 = (aodp * 17. + s4) * temp; d4 = temp * .5 * aodp * tsi * (aodp * 221. + s4 * 31.) * c1; t3cof = d2 + c1sq * 2.; t4cof = (d3 * 3. + c1 * (d2 * 12. + c1sq * 10.)) * .25; t5cof = (d4 * 3. + c1 * 12. * d3 + d2 * 6. * d2 + c1sq * 15. * (d2 * 2. + c1sq)) * .2; } *iflag = 0; } /* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG */ xmdf = element.xmo + xmdot * tsince; omgadf = element.omegao + omgdot * tsince; xnoddf = element.xnodeo + xnodot * tsince; omega = omgadf; xmp = xmdf; tsq = tsince * tsince; xnode = xnoddf + xnodcf * tsq; tempa = 1. - c1 * tsince; tempe = element.bstar * c4 * tsince; templ = t2cof * tsq; if (isimp != 1) { delomg = omgcof * tsince; delm = xmcof * (pow (1. + eta * cos (xmdf), 3.) - delmo); temp = delomg + delm; xmp = xmdf + temp; omega = omgadf - temp; tcube = tsq * tsince; tfour = tsince * tcube; tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour; tempe += element.bstar * c5 * (sin (xmp) - sinmo); templ = templ + t3cof * tcube + tfour * (t4cof + tsince * t5cof); } a = aodp * tempa * tempa; e = element.eo - tempe; xl = xmp + omega + xnode + xnodp * templ; beta = sqrt (1. - e * e); xn = pcnsts.xke / pow (a, 1.5); /* LONG PERIOD PERIODICS */ axn = e * cos (omega); temp = 1. / (a * beta * beta); xll = temp * xlcof * axn; aynl = temp * aycof; xlt = xl + xll; ayn = e * sin (omega) + aynl; /* SOLVE KEPLERS EQUATION */ capu = fmod2p (xlt - xnode); temp2 = capu; for (i = 1; i <= 10; ++i) { sinepw = sin (temp2); cosepw = cos (temp2); temp3 = axn * sinepw; temp4 = ayn * cosepw; temp5 = axn * cosepw; temp6 = ayn * sinepw; epw = (capu - temp4 + temp3 - temp2) / (1. - temp5 - temp6) + temp2; if (fabs (epw - temp2) <= mcnsts.e6a) { break; } temp2 = epw; } /* SHORT PERIOD PRELIMINARY QUANTITIES */ ecose = temp5 + temp6; esine = temp3 - temp4; elsq = axn * axn + ayn * ayn; temp = 1. - elsq; pl = a * temp; r = a * (1. - ecose); temp1 = 1. / r; rdot = pcnsts.xke * sqrt (a) * esine * temp1; rfdot = pcnsts.xke * sqrt (pl) * temp1; temp2 = a * temp1; betal = sqrt (temp); temp3 = 1. / (betal + 1.); cosu = temp2 * (cosepw - axn + ayn * esine * temp3); sinu = temp2 * (sinepw - ayn - axn * esine * temp3); u = matan2 (sinu, cosu); sin2u = sinu * 2. * cosu; cos2u = cosu * 2. * cosu - 1.; temp = 1. / pl; temp1 = pcnsts.ck2 * temp; temp2 = temp1 * temp; /* UPDATE FOR SHORT PERIODICS */ rk = r * (1. - temp2 * 1.5 * betal * x3thm1) + temp1 * .5 * x1mth2 * cos2u; uk = u - temp2 * .25 * x7thm1 * sin2u; xnodek = xnode + temp2 * 1.5 * cosio * sin2u; xinck = element.xincl + temp2 * 1.5 * cosio * sinio * cos2u; rdotk = rdot - xn * temp1 * x1mth2 * sin2u; rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + x3thm1 * 1.5); /* ORIENTATION VECTORS */ sinuk = sin (uk); cosuk = cos (uk); sinik = sin (xinck); cosik = cos (xinck); sinnok = sin (xnodek); cosnok = cos (xnodek); 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 = rdotk * ux + rfdotk * vx; element.ydot = rdotk * uy + rfdotk * vy; element.zdot = rdotk * uz + rfdotk * vz; return (0); }
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.