simgear/math/sg_geodesy.hxx

   1 #ifndef _SG_GEODESY_HXX
   2 #define _SG_GEODESY_HXX
   3
   4 #include "SGMath.hxx"
   5
   6 // Compatibility header.
   7 // Please use the SGGeodesy and SGMath functions directly.
   8
   9 /**
  10  * Convert from geodetic coordinates to geocentric coordinates.
  11  * WARNING: this function is non-reversible.  Due to the fact that
  12  * "up" is a different direction for geocentric and geodetic frames,
  13  * you can not simply add your "alt" parameter to the "sl_radius"
  14  * result and get back (via sgGeodToGeoc()) to the coordinates you
  15  * started with.  The error under normal conditions will be of
  16  * centimeter order; whether that is important or not is application
  17  * dependent. Consider using sgGeodToCart() instead.
  18  *
  19  * @param lat_geod (in) Geodetic latitude, radians, + = North
  20  * @param alt (in) C.G. altitude above mean sea level (meters)
  21  * @param sl_radius (out) SEA LEVEL radius to earth center (meters)
  22  * @param lat_geoc (out) Geocentric latitude, radians, + = North
  23  */
  24 inline void sgGeodToGeoc(double lat_geod, double alt,
  25                          double *sl_radius, double *lat_geoc)
  26 {
  27   SGVec3<double> cart;
  28   SGGeod geod = SGGeod::fromRadM(0, lat_geod, alt);
  29   SGGeodesy::SGGeodToCart(geod, cart);
  30   SGGeoc geoc;
  31   SGGeodesy::SGCartToGeoc(cart, geoc);
  32   *lat_geoc = geoc.getLatitudeRad();
  33   *sl_radius = SGGeodesy::SGGeodToSeaLevelRadius(geod);
  34 }
  35
  36
  37 /**
  38  * Convert a cartesian point to a geodetic lat/lon/altitude.
  39  *
  40  * @param xyz (in) Pointer to cartesian point.
  41  * @param lat (out) Latitude, in radians
  42  * @param lon (out) Longitude, in radians
  43  * @param alt (out) Altitude, in meters above the WGS84 ellipsoid
  44  */
  45 inline void sgCartToGeod(const double* xyz, double* lat, double* lon, double* alt)
  46 {
  47   SGGeod geod;
  48   SGGeodesy::SGCartToGeod(SGVec3<double>(xyz), geod);
  49   *lat = geod.getLatitudeRad();
  50   *lon = geod.getLongitudeRad();
  51   *alt = geod.getElevationM();
  52 }
  53
  54 /**
  55  * Convert a geodetic lat/lon/altitude to a cartesian point.
  56  *
  57  * @param lat (in) Latitude, in radians
  58  * @param lon (in) Longitude, in radians
  59  * @param alt (in) Altitude, in meters above the WGS84 ellipsoid
  60  * @param xyz (out) Pointer to cartesian point.
  61  */
  62 inline void sgGeodToCart(double lat, double lon, double alt, double* xyz)
  63 {
  64   SGVec3<double> cart;
  65   SGGeodesy::SGGeodToCart(SGGeod::fromRadM(lon, lat, alt), cart);
  66   xyz[0] = cart(0);
  67   xyz[1] = cart(1);
  68   xyz[2] = cart(2);
  69 }
  70
  71 /**
  72  * Given a starting position and an offset radial and distance,
  73  * calculate an ending positon on a wgs84 ellipsoid.
  74  * @param alt (in) meters (unused)
  75  * @param lat1 (in) degrees
  76  * @param lon1 (in) degrees
  77  * @param az1 (in) degrees
  78  * @param s (in) distance in meters
  79  * @param lat2 (out) degrees
  80  * @param lon2 (out) degrees
  81  * @param az2 (out) return course in degrees
  82  */
  83 inline int geo_direct_wgs_84 ( double lat1, double lon1, double az1,
  84                                double s, double *lat2, double *lon2,
  85                                double *az2 )
  86 {
  87   SGGeod p2;
  88   if (!SGGeodesy::direct(SGGeod::fromDeg(lon1, lat1), az1, s, p2, *az2))
  89     return 1;
  90   *lat2 = p2.getLatitudeDeg();
  91   *lon2 = p2.getLongitudeDeg();
  92   return 0;
  93 }
  94 inline int geo_direct_wgs_84 ( double alt, double lat1,
  95                         double lon1, double az1,
  96                         double s, double *lat2, double *lon2,
  97                         double *az2 )
  98 { return geo_direct_wgs_84(lat1, lon1, az1, s, lat2, lon2, az2); }
  99
 100 /**
 101  * Given a starting position and an offset radial and distance,
 102  * calculate an ending positon on a wgs84 ellipsoid.
 103  * @param p1 (in) geodetic position
 104  * @param az1 (in) degrees
 105  * @param s (in) distance in meters
 106  * @param p2 (out) geodetic position
 107  * @param az2 (out) return course in degrees
 108  */
 109 inline int geo_direct_wgs_84(const SGGeod& p1, double az1,
 110                              double s, SGGeod& p2, double *az2 )
 111 {
 112   return !SGGeodesy::direct(p1, az1, s, p2, *az2);
 113 }
 114
 115 /**
 116  * Given an altitude and two sets of (lat, lon) calculate great circle
 117  * distance between them as well as the starting and ending azimuths.
 118  * @param alt (in) meters (unused)
 119  * @param lat1 (in) degrees
 120  * @param lon1 (in) degrees
 121  * @param lat2 (in) degrees
 122  * @param lon2 (in) degrees
 123  * @param az1 (out) start heading degrees
 124  * @param az2 (out) end heading degrees
 125  * @param s (out) distance meters
 126  */
 127 inline int geo_inverse_wgs_84( double lat1, double lon1, double lat2,
 128                                double lon2, double *az1, double *az2,
 129                                double *s )
 130 {
 131   return !SGGeodesy::inverse(SGGeod::fromDeg(lon1, lat1),
 132                              SGGeod::fromDeg(lon2, lat2), *az1, *az2, *s);
 133 }
 134 inline int geo_inverse_wgs_84( double alt, double lat1,
 135                                double lon1, double lat2,
 136                                double lon2, double *az1, double *az2,
 137                                double *s )
 138 { return geo_inverse_wgs_84(lat1, lon1, lat2, lon2, az1, az2, s); }
 139
 140
 141 /**
 142  * Given an altitude and two sets of (lat, lon) calculate great circle
 143  * distance between them as well as the starting and ending azimuths.
 144  * @param p1 (in) first position
 145  * @param p2 (in) fsecond position
 146  * @param az1 (out) start heading degrees
 147  * @param az2 (out) end heading degrees
 148  * @param s (out) distance meters
 149  */
 150 inline int geo_inverse_wgs_84(const SGGeod& p1, const SGGeod& p2,
 151                               double *az1, double *az2, double *s )
 152 {
 153   return !SGGeodesy::inverse(p1, p2, *az1, *az2, *s);
 154 }
 155
 156 #endif // _SG_GEODESY_HXX