sgp4ext.cpp   sgp4ext.cpp 
skipping to change at line 27 skipping to change at line 27
* fix sgn * fix sgn
* changes : * changes :
* 2 apr 07 david vallado * 2 apr 07 david vallado
* fix jday floor and str lengths * fix jday floor and str lengths
* updates for constants * updates for constants
* 14 aug 06 david vallado * 14 aug 06 david vallado
* original baseline * original baseline
* ---------------------------------------------------------------- */ * ---------------------------------------------------------------- */
#include "sgp4ext.h" #include "sgp4ext.h"
#ifdef _MSC_BUILD
#include "StelUtils.hpp" #include "StelUtils.hpp"
#endif
double sgn double sgn
( (
double x double x
) )
{ {
if (x < 0.0) if (x < 0.0)
{ {
return -1.0; return -1.0;
} }
skipping to change at line 168 skipping to change at line 170
* coupling : * coupling :
* dot dot product of two vectors * dot dot product of two vectors
* ------------------------------------------------------------------------- -- */ * ------------------------------------------------------------------------- -- */
double angle double angle
( (
double vec1[3], double vec1[3],
double vec2[3] double vec2[3]
) )
{ {
double small, undefined, magv1, magv2, temp; double sv, undefined, magv1, magv2, temp;
small = 0.00000001; sv = 0.00000001;
undefined = 999999.1; undefined = 999999.1;
magv1 = mag(vec1); magv1 = mag(vec1);
magv2 = mag(vec2); magv2 = mag(vec2);
if (magv1*magv2 > small*small) if (magv1*magv2 > sv*sv)
{ {
temp= dot(vec1,vec2) / (magv1*magv2); temp= dot(vec1,vec2) / (magv1*magv2);
if (fabs( temp ) > 1.0) if (fabs( temp ) > 1.0)
temp= sgn(temp) * 1.0; temp= sgn(temp) * 1.0;
return acos( temp ); return acos( temp );
} }
else else
return undefined; return undefined;
} // end angle } // end angle
skipping to change at line 227 skipping to change at line 229
* references : * references :
* vallado 2007, 85, alg 5 * vallado 2007, 85, alg 5
* ------------------------------------------------------------------------- -- */ * ------------------------------------------------------------------------- -- */
void newtonnu void newtonnu
( (
double ecc, double nu, double ecc, double nu,
double& e0, double& m double& e0, double& m
) )
{ {
double small, sine, cose; double sv, sine, cose;
// --------------------- implementation --------------------- // --------------------- implementation ---------------------
e0= 999999.9; e0= 999999.9;
m = 999999.9; m = 999999.9;
small = 0.00000001; sv = 0.00000001;
// --------------------------- circular ------------------------ // --------------------------- circular ------------------------
if ( fabs( ecc ) < small ) if ( fabs( ecc ) < sv )
{ {
m = nu; m = nu;
e0= nu; e0= nu;
} }
else else
// ---------------------- elliptical ----------------------- // ---------------------- elliptical -----------------------
if ( ecc < 1.0-small ) if ( ecc < 1.0-sv )
{ {
sine= ( sqrt( 1.0 -ecc*ecc ) * sin(nu) ) / ( 1.0 +ecc*cos(nu) ); sine= ( sqrt( 1.0 -ecc*ecc ) * sin(nu) ) / ( 1.0 +ecc*cos(nu) );
cose= ( ecc + cos(nu) ) / ( 1.0 + ecc*cos(nu) ); cose= ( ecc + cos(nu) ) / ( 1.0 + ecc*cos(nu) );
e0 = atan2( sine,cose ); e0 = atan2( sine,cose );
m = e0 - ecc*sin(e0); m = e0 - ecc*sin(e0);
} }
else else
// -------------------- hyperbolic -------------------- // -------------------- hyperbolic --------------------
if ( ecc > 1.0 + small ) if ( ecc > 1.0 + sv )
{ {
if ((ecc > 1.0 ) && (fabs(nu)+0.00001 < M_PI-acos(1.0 /ecc ))) if ((ecc > 1.0 ) && (fabs(nu)+0.00001 < M_PI-acos(1.0 /ecc )))
{ {
sine= ( sqrt( ecc*ecc-1.0 ) * sin(nu) ) / ( 1.0 + ec c*cos(nu) ); sine= ( sqrt( ecc*ecc-1.0 ) * sin(nu) ) / ( 1.0 + ec c*cos(nu) );
#ifdef _MSC_BUILD
e0 = StelUtils::asinh( sine ); e0 = StelUtils::asinh( sine );
#else
e0 = asinh( sine );
#endif
m = ecc*sinh(e0) - e0; m = ecc*sinh(e0) - e0;
} }
} }
else else
// ----------------- parabolic --------------------- // ----------------- parabolic ---------------------
if ( fabs(nu) < 168.0*M_PI/180.0 ) if ( fabs(nu) < 168.0*M_PI/180.0 )
{ {
e0= tan( nu*0.5 ); e0= tan( nu*0.5 );
m = e0 + (e0*e0*e0)/3.0; m = e0 + (e0*e0*e0)/3.0;
} }
skipping to change at line 341 skipping to change at line 347
* vallado 2007, 126, alg 9, ex 2-5 * vallado 2007, 126, alg 9, ex 2-5
* ------------------------------------------------------------------------- -- */ * ------------------------------------------------------------------------- -- */
void rv2coe void rv2coe
( (
double r[3], double v[3], double mu, double r[3], double v[3], double mu,
double& p, double& a, double& ecc, double& incl, double& omega, doub le& argp, double& p, double& a, double& ecc, double& incl, double& omega, doub le& argp,
double& nu, double& m, double& arglat, double& truelon, double& lonp er double& nu, double& m, double& arglat, double& truelon, double& lonp er
) )
{ {
double undefined, small, hbar[3], nbar[3], magr, magv, magn, ebar[3] , sme, double undefined, sv, hbar[3], nbar[3], magr, magv, magn, ebar[3], s me,
rdotv, infinite, temp, c1, hk, twopi, magh, halfpi, e; rdotv, infinite, temp, c1, hk, twopi, magh, halfpi, e;
int i; int i;
char typeorbit[3]; char typeorbit[3];
twopi = 2.0 * M_PI; twopi = 2.0 * M_PI;
halfpi = 0.5 * M_PI; halfpi = 0.5 * M_PI;
small = 0.00000001; sv = 0.00000001;
undefined = 999999.1; undefined = 999999.1;
infinite = 999999.9; infinite = 999999.9;
// ------------------------- implementation ----------------- // ------------------------- implementation -----------------
magr = mag( r ); magr = mag( r );
magv = mag( v ); magv = mag( v );
// ------------------ find h n and e vectors ---------------- // ------------------ find h n and e vectors ----------------
cross( r,v, hbar ); cross( r,v, hbar );
magh = mag( hbar ); magh = mag( hbar );
if ( magh > small ) if ( magh > sv )
{ {
nbar[0]= -hbar[1]; nbar[0]= -hbar[1];
nbar[1]= hbar[0]; nbar[1]= hbar[0];
nbar[2]= 0.0; nbar[2]= 0.0;
magn = mag( nbar ); magn = mag( nbar );
c1 = magv*magv - mu /magr; c1 = magv*magv - mu /magr;
rdotv = dot( r,v ); rdotv = dot( r,v );
for (i= 0; i <= 2; i++) for (i= 0; i <= 2; i++)
ebar[i]= (c1*r[i] - rdotv*v[i])/mu; ebar[i]= (c1*r[i] - rdotv*v[i])/mu;
ecc = mag( ebar ); ecc = mag( ebar );
// ------------ find a e and semi-latus rectum ---------- // ------------ find a e and semi-latus rectum ----------
sme= ( magv*magv*0.5 ) - ( mu /magr ); sme= ( magv*magv*0.5 ) - ( mu /magr );
if ( fabs( sme ) > small ) if ( fabs( sme ) > sv )
a= -mu / (2.0 *sme); a= -mu / (2.0 *sme);
else else
a= infinite; a= infinite;
p = magh*magh/mu; p = magh*magh/mu;
// ----------------- find inclination ------------------- // ----------------- find inclination -------------------
hk= hbar[2]/magh; hk= hbar[2]/magh;
incl= acos( hk ); incl= acos( hk );
// -------- determine type of orbit for later use -------- // -------- determine type of orbit for later use --------
// ------ elliptical, parabolic, hyperbolic inclined ------- // ------ elliptical, parabolic, hyperbolic inclined -------
strcpy(typeorbit,"ei"); strcpy(typeorbit,"ei");
if ( ecc < small ) if ( ecc < sv )
{ {
// ---------------- circular equatorial --------------- // ---------------- circular equatorial ---------------
if ((incl<small) | (fabs(incl-M_PI)<small)) if ((incl<sv) | (fabs(incl-M_PI)<sv))
strcpy(typeorbit,"ce"); strcpy(typeorbit,"ce");
else else
// -------------- circular inclined --------------- // -------------- circular inclined ---------------
strcpy(typeorbit,"ci"); strcpy(typeorbit,"ci");
} }
else else
{ {
// - elliptical, parabolic, hyperbolic equatorial -- // - elliptical, parabolic, hyperbolic equatorial --
if ((incl<small) | (fabs(incl-M_PI)<small)) if ((incl<sv) | (fabs(incl-M_PI)<sv))
strcpy(typeorbit,"ee"); strcpy(typeorbit,"ee");
} }
// ---------- find longitude of ascending node ------------ // ---------- find longitude of ascending node ------------
if ( magn > small ) if ( magn > sv )
{ {
temp= nbar[0] / magn; temp= nbar[0] / magn;
if ( fabs(temp) > 1.0 ) if ( fabs(temp) > 1.0 )
temp= sgn(temp); temp= sgn(temp);
omega= acos( temp ); omega= acos( temp );
if ( nbar[1] < 0.0 ) if ( nbar[1] < 0.0 )
omega= twopi - omega; omega= twopi - omega;
} }
else else
omega= undefined; omega= undefined;
skipping to change at line 448 skipping to change at line 454
{ {
arglat = angle( nbar,r ); arglat = angle( nbar,r );
if ( r[2] < 0.0 ) if ( r[2] < 0.0 )
arglat= twopi - arglat; arglat= twopi - arglat;
m = arglat; m = arglat;
} }
else else
arglat= undefined; arglat= undefined;
// -- find longitude of perigee - elliptical equatorial ---- // -- find longitude of perigee - elliptical equatorial ----
if (( ecc>small ) && (strcmp(typeorbit,"ee") == 0)) if (( ecc>sv ) && (strcmp(typeorbit,"ee") == 0))
{ {
temp= ebar[0]/ecc; temp= ebar[0]/ecc;
if ( fabs(temp) > 1.0 ) if ( fabs(temp) > 1.0 )
temp= sgn(temp); temp= sgn(temp);
lonper= acos( temp ); lonper= acos( temp );
if ( ebar[1] < 0.0 ) if ( ebar[1] < 0.0 )
lonper= twopi - lonper; lonper= twopi - lonper;
if ( incl > halfpi ) if ( incl > halfpi )
lonper= twopi - lonper; lonper= twopi - lonper;
} }
else else
lonper= undefined; lonper= undefined;
// -------- find true longitude - circular equatorial ------ // -------- find true longitude - circular equatorial ------
if (( magr>small ) && ( strcmp(typeorbit,"ce") == 0 )) if (( magr>sv ) && ( strcmp(typeorbit,"ce") == 0 ))
{ {
temp= r[0]/magr; temp= r[0]/magr;
if ( fabs(temp) > 1.0 ) if ( fabs(temp) > 1.0 )
temp= sgn(temp); temp= sgn(temp);
truelon= acos( temp ); truelon= acos( temp );
if ( r[1] < 0.0 ) if ( r[1] < 0.0 )
truelon= twopi - truelon; truelon= twopi - truelon;
if ( incl > halfpi ) if ( incl > halfpi )
truelon= twopi - truelon; truelon= twopi - truelon;
m = truelon; m = truelon;
 End of changes. 21 change blocks. 
18 lines changed or deleted 24 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/