sgp4unit.cpp   sgp4unit.cpp 
skipping to change at line 218 skipping to change at line 218
double si3, double sl2, double sl3, double sl4, double t, double si3, double sl2, double sl3, double sl4, double t,
double xgh2, double xgh3, double xgh4, double xh2, double xh3, double xgh2, double xgh3, double xgh4, double xh2, double xh3,
double xi2, double xi3, double xl2, double xl3, double xl4, double xi2, double xi3, double xl2, double xl3, double xl4,
double zmol, double zmos, double inclo, double zmol, double zmos, double inclo,
char init, char init,
double& ep, double& inclp, double& nodep, double& argpp, double& mp, double& ep, double& inclp, double& nodep, double& argpp, double& mp,
char opsmode char opsmode
) )
{ {
/* --------------------- local variables ------------------------ */ /* --------------------- local variables ------------------------ */
const double twopi = 2.0 * pi; const double twopi = 2.0 * M_PI;
double alfdp, betdp, cosip, cosop, dalf, dbet, dls, double alfdp, betdp, cosip, cosop, dalf, dbet, dls,
f2, f3, pe, pgh, ph, pinc, pl , f2, f3, pe, pgh, ph, pinc, pl ,
sel, ses, sghl, sghs, shll, shs, sil, sel, ses, sghl, sghs, shll, shs, sil,
sinip, sinop, sinzf, sis, sll, sls, xls, sinip, sinop, sinzf, sis, sll, sls, xls,
xnoh, zf, zm, zel, zes, znl, zns; xnoh, zf, zm, zel, zes, znl, zns;
/* kill warning */ /* kill warning */
inclo = 0.; inclo = 0.;
/* ---------------------- constants ----------------------------- */ /* ---------------------- constants ----------------------------- */
skipping to change at line 320 skipping to change at line 320
nodep = nodep + twopi; nodep = nodep + twopi;
xls = mp + argpp + cosip * nodep; xls = mp + argpp + cosip * nodep;
dls = pl + pgh - pinc * nodep * sinip; dls = pl + pgh - pinc * nodep * sinip;
xls = xls + dls; xls = xls + dls;
xnoh = nodep; xnoh = nodep;
nodep = atan2(alfdp, betdp); nodep = atan2(alfdp, betdp);
// sgp4fix for afspc written intrinsic functions // sgp4fix for afspc written intrinsic functions
// nodep used without a trigonometric function ahead // nodep used without a trigonometric function ahead
if ((nodep < 0.0) && (opsmode == 'a')) if ((nodep < 0.0) && (opsmode == 'a'))
nodep = nodep + twopi; nodep = nodep + twopi;
if (fabs(xnoh - nodep) > pi){ if (fabs(xnoh - nodep) > M_PI){
if (nodep < xnoh) if (nodep < xnoh)
nodep = nodep + twopi; nodep = nodep + twopi;
else else
nodep = nodep - twopi; nodep = nodep - twopi;
} }
mp = mp + pl; mp = mp + pl;
argpp = xls - mp - cosip * nodep; argpp = xls - mp - cosip * nodep;
} }
} // if init == 'n' } // if init == 'n'
skipping to change at line 434 skipping to change at line 434
{ {
/* -------------------------- constants ------------------------- */ /* -------------------------- constants ------------------------- */
const double zes = 0.01675; const double zes = 0.01675;
const double zel = 0.05490; const double zel = 0.05490;
const double c1ss = 2.9864797e-6; const double c1ss = 2.9864797e-6;
const double c1l = 4.7968065e-7; const double c1l = 4.7968065e-7;
const double zsinis = 0.39785416; const double zsinis = 0.39785416;
const double zcosis = 0.91744867; const double zcosis = 0.91744867;
const double zcosgs = 0.1945905; const double zcosgs = 0.1945905;
const double zsings = -0.98088458; const double zsings = -0.98088458;
const double twopi = 2.0 * pi; const double twopi = 2.0 * M_PI;
/* --------------------- local variables ------------------------ */ /* --------------------- local variables ------------------------ */
int lsflg; int lsflg;
double a1 , a2 , a3 , a4 , a5 , a6 , a7 , double a1 , a2 , a3 , a4 , a5 , a6 , a7 ,
a8 , a9 , a10 , betasq, cc , ctem , stem , a8 , a9 , a10 , betasq, cc , ctem , stem ,
x1 , x2 , x3 , x4 , x5 , x6 , x7 , x1 , x2 , x3 , x4 , x5 , x6 , x7 ,
x8 , xnodce, xnoi , zcosg , zcosgl, zcosh , zcoshl, x8 , xnodce, xnoi , zcosg , zcosgl, zcosh , zcoshl,
zcosi , zcosil, zsing , zsingl, zsinh , zsinhl, zsini , zcosi , zcosil, zsing , zsingl, zsinh , zsinhl, zsini ,
zsinil, zx , zy; zsinil, zx , zy;
skipping to change at line 704 skipping to change at line 704
double& nm, double& nodem, double& nm, double& nodem,
int& irez, int& irez,
double& atime, double& d2201, double& d2211, double& d3210, double& d3222, double& atime, double& d2201, double& d2211, double& d3210, double& d3222,
double& d4410, double& d4422, double& d5220, double& d5232, double& d5421, double& d4410, double& d4422, double& d5220, double& d5232, double& d5421,
double& d5433, double& dedt, double& didt, double& dmdt, double& dndt, double& d5433, double& dedt, double& didt, double& dmdt, double& dndt,
double& dnodt, double& domdt, double& del1, double& del2, double& del3, double& dnodt, double& domdt, double& del1, double& del2, double& del3,
double& xfact, double& xlamo, double& xli, double& xni double& xfact, double& xlamo, double& xli, double& xni
) )
{ {
/* --------------------- local variables ------------------------ */ /* --------------------- local variables ------------------------ */
const double twopi = 2.0 * pi; const double twopi = 2.0 * M_PI;
double ainv2 , aonv=0.0, cosisq, eoc, f220 , f221 , f311 , double ainv2 , aonv=0.0, cosisq, eoc, f220 , f221 , f311 ,
f321 , f322 , f330 , f441 , f442 , f522 , f523 , f321 , f322 , f330 , f441 , f442 , f522 , f523 ,
f542 , f543 , g200 , g201 , g211 , g300 , g310 , f542 , f543 , g200 , g201 , g211 , g300 , g310 ,
g322 , g410 , g422 , g520 , g521 , g532 , g533 , g322 , g410 , g422 , g520 , g521 , g532 , g533 ,
ses , sgs , sghl , sghs , shs , shll , sis , ses , sgs , sghl , sghs , shs , shll , sis ,
sini2 , sls , temp , temp1 , theta , xno2 , q22 , sini2 , sls , temp , temp1 , theta , xno2 , q22 ,
q31 , q33 , root22, root44, root54, rptim , root32, q31 , q33 , root22, root44, root54, rptim , root32,
root52, x2o3 , xke , znl , emo , zns , emsqo, root52, x2o3 , xke , znl , emo , zns , emsqo,
tumin, mu, radiusearthkm, j2, j3, j4, j3oj2; tumin, mu, radiusearthkm, j2, j3, j4, j3oj2;
skipping to change at line 746 skipping to change at line 746
if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5)) if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5))
irez = 2; irez = 2;
/* ------------------------ do solar terms ------------------- */ /* ------------------------ do solar terms ------------------- */
ses = ss1 * zns * ss5; ses = ss1 * zns * ss5;
sis = ss2 * zns * (sz11 + sz13); sis = ss2 * zns * (sz11 + sz13);
sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq); sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
sghs = ss4 * zns * (sz31 + sz33 - 6.0); sghs = ss4 * zns * (sz31 + sz33 - 6.0);
shs = -zns * ss2 * (sz21 + sz23); shs = -zns * ss2 * (sz21 + sz23);
// sgp4fix for 180 deg incl // sgp4fix for 180 deg incl
if ((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2)) if ((inclm < 5.2359877e-2) || (inclm > M_PI - 5.2359877e-2))
shs = 0.0; shs = 0.0;
if (sinim != 0.0) if (sinim != 0.0)
shs = shs / sinim; shs = shs / sinim;
sgs = sghs - cosim * shs; sgs = sghs - cosim * shs;
/* ------------------------- do lunar terms ------------------ */ /* ------------------------- do lunar terms ------------------ */
dedt = ses + s1 * znl * s5; dedt = ses + s1 * znl * s5;
didt = sis + s2 * znl * (z11 + z13); didt = sis + s2 * znl * (z11 + z13);
dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq); dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
sghl = s4 * znl * (z31 + z33 - 6.0); sghl = s4 * znl * (z31 + z33 - 6.0);
shll = -znl * s2 * (z21 + z23); shll = -znl * s2 * (z21 + z23);
// sgp4fix for 180 deg incl // sgp4fix for 180 deg incl
if ((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2)) if ((inclm < 5.2359877e-2) || (inclm > M_PI - 5.2359877e-2))
shll = 0.0; shll = 0.0;
domdt = sgs + sghl; domdt = sgs + sghl;
dnodt = shs; dnodt = shs;
if (sinim != 0.0) if (sinim != 0.0)
{ {
domdt = domdt - cosim / sinim * shll; domdt = domdt - cosim / sinim * shll;
dnodt = dnodt + shll / sinim; dnodt = dnodt + shll / sinim;
} }
/* ----------- calculate deep space resonance effects -------- */ /* ----------- calculate deep space resonance effects -------- */
skipping to change at line 782 skipping to change at line 782
em = em + dedt * t; em = em + dedt * t;
inclm = inclm + didt * t; inclm = inclm + didt * t;
argpm = argpm + domdt * t; argpm = argpm + domdt * t;
nodem = nodem + dnodt * t; nodem = nodem + dnodt * t;
mm = mm + dmdt * t; mm = mm + dmdt * t;
// sgp4fix for negative inclinations // sgp4fix for negative inclinations
// the following if statement should be commented out // the following if statement should be commented out
//if (inclm < 0.0) //if (inclm < 0.0)
// { // {
// inclm = -inclm; // inclm = -inclm;
// argpm = argpm - pi; // argpm = argpm - M_PI;
// nodem = nodem + pi; // nodem = nodem + M_PI;
// } // }
/* -------------- initialize the resonance terms ------------- */ /* -------------- initialize the resonance terms ------------- */
if (irez != 0) if (irez != 0)
{ {
aonv = pow(nm / xke, x2o3); aonv = pow(nm / xke, x2o3);
/* ---------- geopotential resonance for 12 hour orbits ------ */ /* ---------- geopotential resonance for 12 hour orbits ------ */
if (irez == 2) if (irez == 2)
{ {
skipping to change at line 992 skipping to change at line 992
double d2201, double d2211, double d3210, double d3222, double d4410, double d2201, double d2211, double d3210, double d3222, double d4410,
double d4422, double d5220, double d5232, double d5421, double d5433, double d4422, double d5220, double d5232, double d5421, double d5433,
double dedt, double del1, double del2, double del3, double didt, double dedt, double del1, double del2, double del3, double didt,
double dmdt, double dnodt, double domdt, double argpo, double argpdot, double dmdt, double dnodt, double domdt, double argpo, double argpdot,
double t, double tc, double gsto, double xfact, double xlamo, double t, double tc, double gsto, double xfact, double xlamo,
double no, double no,
double& atime, double& em, double& argpm, double& inclm, double& xli, double& atime, double& em, double& argpm, double& inclm, double& xli,
double& mm, double& xni, double& nodem, double& dndt, double& nm double& mm, double& xni, double& nodem, double& dndt, double& nm
) )
{ {
const double twopi = 2.0 * pi; const double twopi = 2.0 * M_PI;
int iretn , iret; int iretn , iret;
double delt, ft, theta, x2li, x2omi, xl, xldot , xnddt, xndt, xomi, g2 2, g32, double delt, ft, theta, x2li, x2omi, xl, xldot , xnddt, xndt, xomi, g2 2, g32,
g44, g52, g54, fasx2, fasx4, fasx6, rptim , step2, stepn , stepp; g44, g52, g54, fasx2, fasx4, fasx6, rptim , step2, stepn , stepp;
fasx2 = 0.13130908; fasx2 = 0.13130908;
fasx4 = 2.8843198; fasx4 = 2.8843198;
fasx6 = 0.37448087; fasx6 = 0.37448087;
g22 = 5.7686396; g22 = 5.7686396;
g32 = 0.95240898; g32 = 0.95240898;
g44 = 1.8014998; g44 = 1.8014998;
skipping to change at line 1200 skipping to change at line 1200
/* --------------------- local variables ------------------------ */ /* --------------------- local variables ------------------------ */
double ak, d1, del, adel, po, x2o3, j2, xke, double ak, d1, del, adel, po, x2o3, j2, xke,
tumin, mu, radiusearthkm, j3, j4, j3oj2; tumin, mu, radiusearthkm, j3, j4, j3oj2;
/* kill warning */ /* kill warning */
satn = 0; satn = 0;
// sgp4fix use old way of finding gst // sgp4fix use old way of finding gst
double ds70; double ds70;
double ts70, tfrac, c1, thgr70, fk5r, c1p2p; double ts70, tfrac, c1, thgr70, fk5r, c1p2p;
const double twopi = 2.0 * pi; const double twopi = 2.0 * M_PI;
/* ----------------------- earth constants ---------------------- */ /* ----------------------- earth constants ---------------------- */
// sgp4fix identify constants and allow alternate values // sgp4fix identify constants and allow alternate values
getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j 3oj2 ); getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j 3oj2 );
x2o3 = 2.0 / 3.0; x2o3 = 2.0 / 3.0;
/* ------------- calculate auxillary epoch quantities ---------- */ /* ------------- calculate auxillary epoch quantities ---------- */
eccsq = ecco * ecco; eccsq = ecco * ecco;
omeosq = 1.0 - eccsq; omeosq = 1.0 - eccsq;
rteosq = sqrt(omeosq); rteosq = sqrt(omeosq);
skipping to change at line 1516 skipping to change at line 1516
if (fabs(cosio+1.0) > 1.5e-12) if (fabs(cosio+1.0) > 1.5e-12)
satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / ( 1.0 + cosio); satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / ( 1.0 + cosio);
else else
satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / t emp4; satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / t emp4;
satrec.aycof = -0.5 * j3oj2 * sinio; satrec.aycof = -0.5 * j3oj2 * sinio;
satrec.delmo = pow((1.0 + satrec.eta * cos(satrec.mo)), 3); satrec.delmo = pow((1.0 + satrec.eta * cos(satrec.mo)), 3);
satrec.sinmao = sin(satrec.mo); satrec.sinmao = sin(satrec.mo);
satrec.x7thm1 = 7.0 * cosio2 - 1.0; satrec.x7thm1 = 7.0 * cosio2 - 1.0;
/* --------------- deep space initialization ------------- */ /* --------------- deep space initialization ------------- */
if ((2*pi / satrec.no) >= 225.0) if ((2*M_PI / satrec.no) >= 225.0)
{ {
satrec.method = 'd'; satrec.method = 'd';
satrec.isimp = 1; satrec.isimp = 1;
tc = 0.0; tc = 0.0;
inclm = satrec.inclo; inclm = satrec.inclo;
dscom dscom
( (
epoch, satrec.ecco, satrec.argpo, tc, satrec.inclo, satr ec.nodeo, epoch, satrec.ecco, satrec.argpo, tc, satrec.inclo, satr ec.nodeo,
satrec.no, snodm, cnodm, sinim, cosim,sinomm, cosom m, satrec.no, snodm, cnodm, sinim, cosim,sinomm, cosom m,
skipping to change at line 1721 skipping to change at line 1721
xmdf , xmx , xmy , nodedf, xnode , nodep, tc , dndt, xmdf , xmx , xmy , nodedf, xnode , nodep, tc , dndt,
twopi, x2o3 , j2 , j3 , tumin, j4 , xke , j3oj2, radiusear thkm, twopi, x2o3 , j2 , j3 , tumin, j4 , xke , j3oj2, radiusear thkm,
mu, vkmpersec; mu, vkmpersec;
int ktr; int ktr;
/* ------------------ set mathematical constants --------------- */ /* ------------------ set mathematical constants --------------- */
// sgp4fix divisor for divide by zero check on inclination // sgp4fix divisor for divide by zero check on inclination
// the old check used 1.0 + cos(pi-1.0e-9), but then compared it to // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
// 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency
const double temp4 = 1.5e-12; const double temp4 = 1.5e-12;
twopi = 2.0 * pi; twopi = 2.0 * M_PI;
x2o3 = 2.0 / 3.0; x2o3 = 2.0 / 3.0;
// sgp4fix identify constants and allow alternate values // sgp4fix identify constants and allow alternate values
getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j 3oj2 ); getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j 3oj2 );
vkmpersec = radiusearthkm * xke/60.0; vkmpersec = radiusearthkm * xke/60.0;
/* --------------------- clear sgp4 error flag ----------------- */ /* --------------------- clear sgp4 error flag ----------------- */
satrec.t = tsince; satrec.t = tsince;
satrec.error = 0; satrec.error = 0;
/* ------- update for secular gravity and atmospheric drag ----- */ /* ------- update for secular gravity and atmospheric drag ----- */
skipping to change at line 1850 skipping to change at line 1850
satrec.sl4, satrec.t, satrec.xgh2, satrec.sl4, satrec.t, satrec.xgh2,
satrec.xgh3, satrec.xgh4, satrec.xh2, satrec.xgh3, satrec.xgh4, satrec.xh2,
satrec.xh3, satrec.xi2, satrec.xi3, satrec.xh3, satrec.xi2, satrec.xi3,
satrec.xl2, satrec.xl3, satrec.xl4, satrec.xl2, satrec.xl3, satrec.xl4,
satrec.zmol, satrec.zmos, satrec.inclo, satrec.zmol, satrec.zmos, satrec.inclo,
'n', ep, xincp, nodep, argpp, mp, satrec.operationmode 'n', ep, xincp, nodep, argpp, mp, satrec.operationmode
); );
if (xincp < 0.0) if (xincp < 0.0)
{ {
xincp = -xincp; xincp = -xincp;
nodep = nodep + pi; nodep = nodep + M_PI;
argpp = argpp - pi; argpp = argpp - M_PI;
} }
if ((ep < 0.0 ) || ( ep > 1.0)) if ((ep < 0.0 ) || ( ep > 1.0))
{ {
// printf("# error ep %f\n", ep); // printf("# error ep %f\n", ep);
satrec.error = 3; satrec.error = 3;
// sgp4fix add return // sgp4fix add return
return false; return false;
} }
} // if method = d } // if method = d
skipping to change at line 2011 skipping to change at line 2011
* *
* references : * references :
* vallado 2004, 191, eq 3-45 * vallado 2004, 191, eq 3-45
* ------------------------------------------------------------------------- -- */ * ------------------------------------------------------------------------- -- */
double gstime double gstime
( (
double jdut1 double jdut1
) )
{ {
const double twopi = 2.0 * pi; const double twopi = 2.0 * M_PI;
const double deg2rad = pi / 180.0; const double deg2rad = M_PI / 180.0;
double temp, tut1; double temp, tut1;
tut1 = (jdut1 - 2451545.0) / 36525.0; tut1 = (jdut1 - 2451545.0) / 36525.0;
temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 +
(876600.0*3600 + 8640184.812866) * tut1 + 67310.54841; // sec (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841; // sec
temp = fmod(temp * deg2rad / 240.0, twopi); //360/86400 = 1/240, to de g, to rad temp = fmod(temp * deg2rad / 240.0, twopi); //360/86400 = 1/240, to de g, to rad
// ------------------------ check quadrants --------------------- // ------------------------ check quadrants ---------------------
if (temp < 0.0) if (temp < 0.0)
temp += twopi; temp += twopi;
 End of changes. 13 change blocks. 
16 lines changed or deleted 16 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/