$a = 6378206.4;
      $b = 6356583.8;

##</code><code>##

// Reduced or Parametric Latitude
   Latitude Ellipsoid::
Eta(const Latitude & lat) const {
   // NOTE: sqrt(1-e^2) == b/a
   return Latitude(atan(b/a*tan(lat*deg2rad))*rad2deg);
}  // Eta


// Radius of Ellipsoid at Latitude.
   double Ellipsoid::
Radius(const Latitude & lat) const {
   Latitude eta = Eta(lat);
   double x = a * cos(eta*deg2rad);
   double y = b * sin(eta*deg2rad);
   return sqrt(x*x + y*y);
}  // Radius

##</code><code>##

// Given a Lat,Lon point, a distance (in meters), and an azimuth (in degrees),
// this routine returns a new Lat,Lon point.
   LatLon Ellipsoid::
NewLatLon(const LatLon & lat_lon, double c, double Az) const {
   // Convert to radians
   double lat = lat_lon.lat * deg2rad;
   // We do not need lon in radians.
   Az *= deg2rad;

   // Equations (5-5) and (5-6) in "Map Projections--A Working Manual", Page 31
   double cosAz = cos(Az);
   double sinAz = sin(Az);
   double sinlat = sin(lat);
   double coslat = cos(lat);

   c /= Radius(lat_lon.lat);
   double cosc = cos(c);
   double sinc = sin(c);

   double y = sinc*sinAz;
   double x = coslat*cosc - sinlat*sinc*cosAz;
   double at = (x == 0.0  &&  y == 0.0) ? 0.0 : atan2(y,x);
   lat = asin(sinlat*cosc + coslat*sinc*cosAz) * rad2deg;
   double lon = lat_lon.lon + at * rad2deg;

   return LatLon(lat, lon);
}  // NewLatLon