Use plib rad/degrees conversion constants.

This commit is contained in:
curt
2001-03-24 03:20:47 +00:00
parent 87e59ff876
commit a7459489ff
15 changed files with 61 additions and 87 deletions

View File

@@ -206,9 +206,9 @@ double SGBucket::get_width_m() const {
} else {
clat = (int)clat - 0.5;
}
double clat_rad = clat * DEG_TO_RAD;
double clat_rad = clat * SGD_DEGREES_TO_RADIANS;
double cos_lat = cos( clat_rad );
double local_radius = cos_lat * EQUATORIAL_RADIUS_M;
double local_radius = cos_lat * SG_EQUATORIAL_RADIUS_M;
double local_perimeter = 2.0 * local_radius * SGD_PI;
double degree_width = local_perimeter / 360.0;
@@ -218,7 +218,7 @@ double SGBucket::get_width_m() const {
// return height of the tile in meters
double SGBucket::get_height_m() const {
double perimeter = 2.0 * EQUATORIAL_RADIUS_M * SGD_PI;
double perimeter = 2.0 * SG_EQUATORIAL_RADIUS_M * SGD_PI;
double degree_height = perimeter / 360.0;
return SG_BUCKET_SPAN * degree_height;

View File

@@ -81,40 +81,14 @@
// Earth parameters for WGS 84, taken from LaRCsim/ls_constants.h
// Value of earth radius from [8]
#define EQUATORIAL_RADIUS_FT 20925650. // ft
#define EQUATORIAL_RADIUS_M 6378138.12 // meter
#define SG_EQUATORIAL_RADIUS_FT 20925650. // ft
#define SG_EQUATORIAL_RADIUS_M 6378138.12 // meter
// Radius squared
#define RESQ_FT 437882827922500. // ft
#define RESQ_M 40680645877797.1344 // meter
#if 0
// Value of earth flattening parameter from ref [8]
//
// Note: FP = f
// E = 1-f
// EPS = sqrt(1-(1-f)^2)
//
#define FP 0.003352813178
#define E 0.996647186
#define EPS 0.081819221
#define INVG 0.031080997
// Time Related Parameters
#define MJD0 2415020.0
#define J2000 (2451545.0 - MJD0)
#define SIDRATE .9972695677
#endif
#define SG_EQ_RAD_SQUARE_FT 437882827922500. // ft
#define SG_EQ_RAD_SQUARE_M 40680645877797.1344 // meter
// Conversions
// Degrees to Radians
#define DEG_TO_RAD 0.017453292 // deg*pi/180 = rad
// Radians to Degrees
#define RAD_TO_DEG 57.29577951 // rad*180/pi = deg
// Arc seconds to radians // (arcsec*pi)/(3600*180) = rad
#define ARCSEC_TO_RAD 4.84813681109535993589e-06

View File

@@ -62,7 +62,7 @@ void CelestialBody::updatePosition(double mjd, Star *ourSun)
actTime = sgCalcActTime(mjd);
// calcualate the angle bewteen ecliptic and equatorial coordinate system
ecl = DEG_TO_RAD * (23.4393 - 3.563E-7 *actTime);
ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
eccAnom = sgCalcEccAnom(M, e); //calculate the eccentric anomaly
xv = a * (cos(eccAnom) - e);
@@ -107,7 +107,7 @@ void CelestialBody::updatePosition(double mjd, Star *ourSun)
tmp = -1.0;
}
FV = RAD_TO_DEG * acos( tmp );
FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
}
/****************************************************************************
@@ -152,7 +152,7 @@ double CelestialBody::sgCalcEccAnom(double M, double e)
diff = fabs(E0 - E1);
E0 = E1;
}
while (diff > (DEG_TO_RAD * 0.001));
while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
return E0;
}
return eccAnom;

View File

@@ -150,10 +150,10 @@ inline CelestialBody::CelestialBody(double Nf, double Ns,
inline void CelestialBody::updateOrbElements(double mjd)
{
double actTime = sgCalcActTime(mjd);
M = DEG_TO_RAD * (MFirst + (MSec * actTime));
w = DEG_TO_RAD * (wFirst + (wSec * actTime));
N = DEG_TO_RAD * (NFirst + (NSec * actTime));
i = DEG_TO_RAD * (iFirst + (iSec * actTime));
M = SGD_DEGREES_TO_RADIANS * (MFirst + (MSec * actTime));
w = SGD_DEGREES_TO_RADIANS * (wFirst + (wSec * actTime));
N = SGD_DEGREES_TO_RADIANS * (NFirst + (NSec * actTime));
i = SGD_DEGREES_TO_RADIANS * (iFirst + (iSec * actTime));
e = eFirst + (eSec * actTime);
a = aFirst + (aSec * actTime);
}

View File

@@ -92,7 +92,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
// calculate the angle between ecliptic and equatorial coordinate system
// in Radians
ecl = ((DEG_TO_RAD * 23.4393) - (DEG_TO_RAD * 3.563E-7) * actTime);
ecl = ((SGD_DEGREES_TO_RADIANS * 23.4393) - (SGD_DEGREES_TO_RADIANS * 3.563E-7) * actTime);
eccAnom = sgCalcEccAnom(M, e); // Calculate the eccentric anomaly
xv = a * (cos(eccAnom) - e);
yv = a * (sqrt(1.0 - e*e) * sin(eccAnom));
@@ -116,7 +116,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
D = Lm - Ls;
F = Lm - N;
lonEcl += DEG_TO_RAD * (-1.274 * sin (M - 2*D)
lonEcl += SGD_DEGREES_TO_RADIANS * (-1.274 * sin (M - 2*D)
+0.658 * sin (2*D)
-0.186 * sin(ourSun->getM())
-0.059 * sin(2*M - 2*D)
@@ -129,7 +129,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
-0.015 * sin(2*F - 2*D)
+0.011 * sin(M - 4*D)
);
latEcl += DEG_TO_RAD * (-0.173 * sin(F-2*D)
latEcl += SGD_DEGREES_TO_RADIANS * (-0.173 * sin(F-2*D)
-0.055 * sin(M - F - 2*D)
-0.046 * sin(M + F - 2*D)
+0.033 * sin(F + 2*D)
@@ -151,8 +151,8 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
geoDec = atan2(ze, sqrt(xe*xe + ye*ye));
/* FG_LOG( FG_GENERAL, FG_INFO,
"(geocentric) geoRa = (" << (RAD_TO_DEG * geoRa)
<< "), geoDec= (" << (RAD_TO_DEG * geoDec) << ")" ); */
"(geocentric) geoRa = (" << (SGD_RADIANS_TO_DEGREES * geoRa)
<< "), geoDec= (" << (SGD_RADIANS_TO_DEGREES * geoDec) << ")" ); */
// Given the moon's geocentric ra and dec, calculate its
@@ -166,10 +166,10 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
// FG_LOG( FG_GENERAL, FG_INFO, "lat = " << f->get_Latitude() );
gclat = lat - 0.003358 *
sin (2 * DEG_TO_RAD * lat );
sin (2 * SGD_DEGREES_TO_RADIANS * lat );
// FG_LOG( FG_GENERAL, FG_INFO, "gclat = " << gclat );
rho = 0.99883 + 0.00167 * cos(2 * DEG_TO_RAD * lat);
rho = 0.99883 + 0.00167 * cos(2 * SGD_DEGREES_TO_RADIANS * lat);
// FG_LOG( FG_GENERAL, FG_INFO, "rho = " << rho );
if (geoRa < 0)
@@ -186,6 +186,6 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
declination = geoDec - mpar * rho * sin (gclat) * sin (g - geoDec) / sin(g);
/* FG_LOG( FG_GENERAL, FG_INFO,
"Ra = (" << (RAD_TO_DEG *rightAscension)
<< "), Dec= (" << (RAD_TO_DEG *declination) << ")" ); */
"Ra = (" << (SGD_RADIANS_TO_DEGREES *rightAscension)
<< "), Dec= (" << (SGD_RADIANS_TO_DEGREES *declination) << ")" ); */
}

View File

@@ -85,7 +85,7 @@ void Star::updatePosition(double mjd)
updateOrbElements(mjd);
actTime = sgCalcActTime(mjd);
ecl = DEG_TO_RAD * (23.4393 - 3.563E-7 * actTime); // Angle in Radians
ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 * actTime); // Angle in Radians
eccAnom = sgCalcEccAnom(M, e); // Calculate the eccentric Anomaly (also known as solving Kepler's equation)
xv = cos(eccAnom) - e;

View File

@@ -42,15 +42,15 @@ if (argc == 8){
}
var = calc_magvar( DEG_TO_RAD * lat_deg, DEG_TO_RAD * lon_deg, h,
var = calc_magvar( SGD_DEGREES_TO_RADIANS * lat_deg, SGD_DEGREES_TO_RADIANS * lon_deg, h,
yymmdd_to_julian_days(yy,mm,dd), field );
fprintf(stdout,"%6.0lf %6.0lf %6.0lf\n", field[0], field[1], field[2] );
fprintf(stdout,"%6.0lf %6.0lf %6.0lf\n", field[3], field[4], field[5] );
fprintf(stdout,"%6.0lf %6.0lf %6.0lf %4.2lf %4.2lf \n",
field[3],field[4],field[5],
RAD_TO_DEG * (atan(field[5]/pow(field[3]*field[3]+field[4]*field[4],0.5))),
RAD_TO_DEG * var);
SGD_RADIANS_TO_DEGREES * (atan(field[5]/pow(field[3]*field[3]+field[4]*field[4],0.5))),
SGD_RADIANS_TO_DEGREES * var);
exit(0);
}

View File

@@ -45,11 +45,11 @@ double fgGeodAltFromCart(const Point3D& cp)
if( ( (SGD_PI_2 - lat_geoc) < SG_ONE_SECOND ) // near North pole
|| ( (SGD_PI_2 + lat_geoc) < SG_ONE_SECOND ) ) // near South pole
{
result = radius - EQUATORIAL_RADIUS_M*E;
result = radius - SG_EQUATORIAL_RADIUS_M*E;
} else {
t_lat = tan(lat_geoc);
x_alpha = E*EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
mu_alpha = atan2(sqrt(RESQ_M - x_alpha*x_alpha),E*x_alpha);
x_alpha = E*SG_EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
mu_alpha = atan2(sqrt(SG_EQ_RAD_SQUARE_M - x_alpha*x_alpha),E*x_alpha);
if (lat_geoc < 0) {
mu_alpha = - mu_alpha;
}

View File

@@ -61,13 +61,13 @@ void sgGeocToGeod( double lat_geoc, double radius, double
|| ( (SGD_PI_2 + lat_geoc) < SG_ONE_SECOND ) ) // near South pole
{
*lat_geod = lat_geoc;
*sea_level_r = EQUATORIAL_RADIUS_M*E;
*sea_level_r = SG_EQUATORIAL_RADIUS_M*E;
*alt = radius - *sea_level_r;
} else {
// cout << " lat_geoc = " << lat_geoc << endl;
t_lat = tan(lat_geoc);
// cout << " tan(t_lat) = " << t_lat << endl;
x_alpha = E*EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
x_alpha = E*SG_EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
#ifdef DOMAIN_ERR_DEBUG
if ( errno ) {
perror("fgGeocToGeod()");
@@ -75,12 +75,12 @@ void sgGeocToGeod( double lat_geoc, double radius, double
}
#endif
// cout << " x_alpha = " << x_alpha << endl;
double tmp = sqrt(RESQ_M - x_alpha * x_alpha);
double tmp = sqrt(SG_EQ_RAD_SQUARE_M - x_alpha * x_alpha);
if ( tmp < 0.0 ) { tmp = 0.0; }
#ifdef DOMAIN_ERR_DEBUG
if ( errno ) {
perror("fgGeocToGeod()");
FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" << RESQ_M - x_alpha * x_alpha
FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" << SG_EQ_RAD_SQUARE_M - x_alpha * x_alpha
<< ")" );
}
#endif
@@ -100,19 +100,19 @@ void sgGeocToGeod( double lat_geoc, double radius, double
1-EPS*EPS*sin_mu_a*sin_mu_a << ")" );
}
#endif
rho_alpha = EQUATORIAL_RADIUS_M*(1-EPS)/
rho_alpha = SG_EQUATORIAL_RADIUS_M*(1-EPS)/
(denom*denom*denom);
delt_mu = atan2(l_point*sin(delt_lambda),rho_alpha + *alt);
*lat_geod = mu_alpha - delt_mu;
lambda_sl = atan( E*E * tan(*lat_geod) ); // SL geoc. latitude
sin_lambda_sl = sin( lambda_sl );
*sea_level_r =
sqrt(RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
sqrt(SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
#ifdef DOMAIN_ERR_DEBUG
if ( errno ) {
perror("fgGeocToGeod()");
FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" <<
RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
<< ")" );
}
#endif
@@ -149,12 +149,12 @@ void sgGeodToGeoc( double lat_geod, double alt, double *sl_radius,
sin_mu = sin(lat_geod); // Geodetic (map makers') latitude
cos_mu = cos(lat_geod);
*sl_radius =
sqrt(RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
sqrt(SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
#ifdef DOMAIN_ERR_DEBUG
if ( errno ) {
perror("fgGeodToGeoc()");
FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" <<
RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
<< ")" );
}
#endif

View File

@@ -48,7 +48,7 @@ double CMetarStation::decodeDMS( char *b )
// Direction (E W N S)
if ( *b == 'W' || *b == 'S' ) r = -r;
}
return r * DEG_TO_RAD;
return r * SGD_DEGREES_TO_RADIANS;
}
// Constructor
@@ -86,10 +86,10 @@ CMetarStation::CMetarStation(
m_altitude = altitude;
s = t; t = strchr( s, ';' ); *t = 0; t++;
double ualtitude = atoi( s ) * FEET_TO_METER;
Point3D p( longitude, latitude, altitude+EQUATORIAL_RADIUS_M );
Point3D p( longitude, latitude, altitude+SG_EQUATORIAL_RADIUS_M );
m_locationPolar = p;
m_locationCart = sgPolarToCart3d( p );
Point3D up( ulongitude, ulatitude, ualtitude+EQUATORIAL_RADIUS_M );
Point3D up( ulongitude, ulatitude, ualtitude+SG_EQUATORIAL_RADIUS_M );
m_upperLocationPolar = up;
m_upperLocationCart = sgPolarToCart3d( up );
s = t;

View File

@@ -178,9 +178,9 @@ point_list calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
clat = (int)clat - 0.5;
}
double clat_rad = clat * DEG_TO_RAD;
double clat_rad = clat * SGD_DEGREES_TO_RADIANS;
double cos_lat = cos( clat_rad );
double local_radius = cos_lat * EQUATORIAL_RADIUS_M;
double local_radius = cos_lat * SG_EQUATORIAL_RADIUS_M;
double local_perimeter = 2.0 * local_radius * SGD_PI;
double degree_width = local_perimeter / 360.0;
@@ -191,7 +191,7 @@ point_list calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
// cout << "local_perimeter = " << local_perimeter << endl;
// cout << "degree_width = " << degree_width << endl;
double perimeter = 2.0 * EQUATORIAL_RADIUS_M * SGD_PI;
double perimeter = 2.0 * SG_EQUATORIAL_RADIUS_M * SGD_PI;
double degree_height = perimeter / 360.0;
// cout << "degree_height = " << degree_height << endl;

View File

@@ -61,14 +61,14 @@ void SGWayPoint::CourseAndDistance( const double cur_lon,
geo_inverse_wgs_84( cur_alt, cur_lat, cur_lon, target_lat, target_lon,
course, &reverse, distance );
} else if ( mode == SPHERICAL ) {
Point3D current( cur_lon * DEG_TO_RAD, cur_lat * DEG_TO_RAD, 0.0 );
Point3D target( target_lon * DEG_TO_RAD, target_lat * DEG_TO_RAD, 0.0 );
Point3D current( cur_lon * SGD_DEGREES_TO_RADIANS, cur_lat * SGD_DEGREES_TO_RADIANS, 0.0 );
Point3D target( target_lon * SGD_DEGREES_TO_RADIANS, target_lat * SGD_DEGREES_TO_RADIANS, 0.0 );
calc_gc_course_dist( current, target, course, distance );
*course = 360.0 - *course * RAD_TO_DEG;
*course = 360.0 - *course * SGD_RADIANS_TO_DEGREES;
} else if ( mode == CARTESIAN ) {
double dx = target_lon - cur_lon;
double dy = target_lat - cur_lat;
*course = -atan2( dy, dx ) * RAD_TO_DEG - 90;
*course = -atan2( dy, dx ) * SGD_RADIANS_TO_DEGREES - 90;
while ( *course < 0 ) {
*course += 360.0;
}

View File

@@ -64,7 +64,7 @@ GeoCoord::GeoCoord(const GeoCoord& other)
// GeoCoordVectorConstIterator i, nearest;
// for (i = data.begin(); i != data.end(); i++)
// {
// angle = RAD_TO_DEG * (*i)->getAngle(ref);
// angle = SGD_RADIANS_TO_DEGREES * (*i)->getAngle(ref);
// if (angle < maxAngle)
// {
// maxAngle = angle;

View File

@@ -62,9 +62,9 @@ public:
void set(float la, float lo) { lat = la; lon = lo; };
float getLat() const { return lat; };
float getLon() const { return lon; };
float getX() const { return cos(DEG_TO_RAD*lat) * cos(DEG_TO_RAD*lon); };
float getY() const { return cos(DEG_TO_RAD*lat) * sin(DEG_TO_RAD*lon); };
float getZ() const { return sin(DEG_TO_RAD*lat); };
float getX() const { return cos(SGD_DEGREES_TO_RADIANS*lat) * cos(SGD_DEGREES_TO_RADIANS*lon); };
float getY() const { return cos(SGD_DEGREES_TO_RADIANS*lat) * sin(SGD_DEGREES_TO_RADIANS*lon); };
float getZ() const { return sin(SGD_DEGREES_TO_RADIANS*lat); };
//double getAngle(const GeoCoord& other) const;

View File

@@ -62,7 +62,7 @@
#define DEGHR(x) ((x)/15.)
#define RADHR(x) DEGHR(x*RAD_TO_DEG)
#define RADHR(x) DEGHR(x*SGD_RADIANS_TO_DEGREES)
static const double MJD0 = 2415020.0;
@@ -89,7 +89,7 @@ SGTime::SGTime( double lon, double lat, const string& root )
<< zone.str() );
tzContainer = new TimezoneContainer( zone.c_str() );
GeoCoord location( RAD_TO_DEG * lat, RAD_TO_DEG * lon );
GeoCoord location( SGD_RADIANS_TO_DEGREES * lat, SGD_RADIANS_TO_DEGREES * lon );
GeoCoord* nearestTz = tzContainer->getNearest(location);
FGPath name( root );
@@ -135,7 +135,7 @@ static double sidereal_precise( double mjd, double lng )
mjd + MJD0, lng); */
// convert to required internal units
lng *= DEG_TO_RAD;
lng *= SGD_DEGREES_TO_RADIANS;
// compute LST and print
double gst = sgTimeCalcGST( mjd );
@@ -219,7 +219,7 @@ void SGTime::update( double lon, double lat, long int warp ) {
jd = mjd + MJD0;
FG_LOG( FG_EVENT, FG_DEBUG, " Current Julian Date = " << jd );
// printf(" Current Longitude = %.3f\n", FG_Longitude * RAD_TO_DEG);
// printf(" Current Longitude = %.3f\n", FG_Longitude * SGD_RADIANS_TO_DEGREES);
// Calculate local side real time
if ( gst_diff < -100.0 ) {
@@ -231,18 +231,18 @@ void SGTime::update( double lon, double lat, long int warp ) {
gst_diff = gst_precise - gst_course;
lst = sidereal_course( cur_time, gmt, -(lon * RAD_TO_DEG) ) + gst_diff;
lst = sidereal_course( cur_time, gmt, -(lon * SGD_RADIANS_TO_DEGREES) ) + gst_diff;
} else {
// course + difference should drift off very slowly
gst = sidereal_course( cur_time, gmt, 0.00 ) + gst_diff;
lst = sidereal_course( cur_time, gmt, -(lon * RAD_TO_DEG) ) + gst_diff;
lst = sidereal_course( cur_time, gmt, -(lon * SGD_RADIANS_TO_DEGREES) ) + gst_diff;
}
FG_LOG( FG_EVENT, FG_DEBUG,
" Current lon=0.00 Sidereal Time = " << gst );
FG_LOG( FG_EVENT, FG_DEBUG,
" Current LOCAL Sidereal Time = " << lst << " ("
<< sidereal_precise( mjd, -(lon * RAD_TO_DEG) )
<< sidereal_precise( mjd, -(lon * SGD_RADIANS_TO_DEGREES) )
<< ") (diff = " << gst_diff << ")" );
}
@@ -252,7 +252,7 @@ void SGTime::updateLocal( double lon, double lat, const string& root )
{
time_t currGMT;
time_t aircraftLocalTime;
GeoCoord location( RAD_TO_DEG * lat, RAD_TO_DEG * lon );
GeoCoord location( SGD_RADIANS_TO_DEGREES * lat, SGD_RADIANS_TO_DEGREES * lon );
GeoCoord* nearestTz = tzContainer->getNearest(location);
FGPath zone( root );
zone.append ( nearestTz->getDescription() );