simgear/math/sg_geodesy.hxx

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