// sdutil.cpp -- Miscellaneous functions
//
// Written by Frederic Bouvier, started October 2001.
//
// Copyright (C) 2001  Frederic Bouvier - fredb@users.sourceforge.net
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//
// $Id: sdutil.cpp,v 1.12 2005/05/09 07:02:14 fredb Exp $

#include <sys/stat.h>
#ifdef _MSC_VER
#include <direct.h>
#define  _USE_MATH_DEFINES
#endif
#include <math.h>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <simgear/misc/sg_path.hxx>
#include "sdutil.hpp"
#include <FL/filename.H>

#ifdef _WIN32
#include <ctype.h>
#define CHAR_COMPARE(c)       tolower(c)
#else
#define CHAR_COMPARE(c)       c
#endif

std::string FGSD_Util::buildRelativePath( const char *ref, const char *abspath, bool &rel )
{
   size_t i = 0;
   size_t lastslash = (size_t)-1;
   while ( ref[ i ] != '\0' && abspath[ i ] != '\0' && CHAR_COMPARE( ref[ i ] ) == CHAR_COMPARE( abspath[ i ] ) )
      {
      if ( ref[ i ] == '/' )
         lastslash = i;

      i += 1;
      }

   std::string ret;
   if ( lastslash == (size_t)-1 )
      {
      ret = abspath;
      rel = false;
      }
   else
      {
      size_t nb = 0;
      i = lastslash + 1;
      while ( ref[ i ] != '\0' )
         {
         if ( ref[ i ] == '/' )
            nb += 1;
         i += 1;
         }
      for ( i = 0; i < nb; i++ )
         ret += "../";
      ret += &abspath[ lastslash + 1 ];
      rel = true;
      }

   return ret;
}

std::string FGSD_Util::buildAbsolutePath( const char *ref, const char *relpath )
{
   size_t nb = 0;
   size_t i = 0;
   while ( relpath[ i ] == '.' || relpath[ i ] == '/' )
      {
      if ( relpath[ i ] == '/' && i > 1 && relpath[ i - 1 ] == '.' && relpath[ i - 2 ] == '.' )
         nb += 1;

      i += 1;
      }
   size_t firstcar = i;

   size_t lastslash = (size_t)-1;
   nb += 1;
   i = strlen( ref );
   while ( nb > 0 && i != (size_t)-1 )
      {
      if ( ref[ i ] == '/' )
         {
         nb -= 1;
         lastslash = i;
         }

      i -= 1;
      }

   std::string ret;
   ret.append( ref, lastslash + 1 );
   ret.append( &relpath[ firstcar ] );
   return ret;
}

std::string FGSD_Util::getDirectoryFromFile( const char *fname )
{
   std::string fname_ = fname;
   size_t len = fname_.size();
   size_t pos = 0;
   size_t oldPos = (size_t)-1;
   while ( pos < len && ( pos = fname_.find( '/', pos ) ) != std::string::npos )
   {
      oldPos = pos;
      pos = pos + 1;
   }
   std::string dir;
   if ( oldPos > 0 )
      dir = fname_.substr( 0, oldPos );
   return dir;
}

double FGSD_Util::parsePosition( const char *_pos_, bool &_error_ )
{
   int sign = 1;
   size_t len = strlen( _pos_ );
   double hrs = 0.0; bool fhrs = false;
   double min = 0.0; bool fmin = false;
   double sec = 0.0; bool fsec = false;
   double dec = 0.0; bool fdec = false;
   double div = 1.0;
   double result = 0.0;
   bool error = false;

   if ( len > 0 )
      {
      if ( _pos_[ 0 ] == '+' || _pos_[ 0 ] == 'E' || _pos_[ 0 ] == 'N' )
         {
         len -= 1;
         _pos_ += 1;
         }
      else if ( _pos_[ 0 ] == '-' || _pos_[ 0 ] == 'W' || _pos_[ 0 ] == 'S' )
         {
         len -= 1;
         _pos_ += 1;
         sign = -1;
         }

      while ( _pos_[ 0 ] == ' ' )
         {
         len -= 1;
         _pos_ += 1;
         }

      while ( len > 0 && isdigit( _pos_[ 0 ] ) )
         {
         fhrs = true;
         hrs = hrs * 10 + _pos_[ 0 ] - '0';
         _pos_ += 1;
         len -= 1;
         }

      bool point = false;
      if ( len > 0 && _pos_[ 0 ] == ':' )
         {
         len -= 1;
         _pos_ += 1;

         while ( len > 0 && isdigit( _pos_[ 0 ] ) )
            {
            fmin = true;
            min = min * 10 + _pos_[ 0 ] - '0';
            _pos_ += 1;
            len -= 1;
            }
         if ( min < 60.0 )
            {
            if ( len > 0 && ( _pos_[ 0 ] == ':' || _pos_[ 0 ] == '\'' ) )
               {
               len -= 1;
               _pos_ += 1;

               while ( len > 0 && isdigit( _pos_[ 0 ] ) )
                  {
                  fsec = true;
                  sec = sec * 10 + _pos_[ 0 ] - '0';
                  _pos_ += 1;
                  len -= 1;
                  }
               if ( sec < 60.0 )
                  {
                  if ( len > 0 && _pos_[ 0 ] == '.' )
                     {
                     point = true;
                     }
                  else if ( len > 0 )
                     {
                     error = true;
                     }
                  }
               else
                  {
                  error = true;
                  }
               }
            else if ( len > 0 && _pos_[ 0 ] == '.' )
               {
               point = true;
               }
            else if ( len > 0 )
               {
               error = true;
               }
            }
         else
            {
            error = true;
            }
         }
      else if ( len > 0 && _pos_[ 0 ] == '.' )
         {
         point = true;
         }
      else if ( len > 0 )
         {
         error = true;
         }

      if ( !error && point )
         {
         len -= 1;
         _pos_ += 1;
         while ( len > 0 && isdigit( _pos_[ 0 ] ) )
            {
            fdec = true;
            dec = dec * 10 + _pos_[ 0 ] - '0';
            div = div * 10;
            len -= 1;
            _pos_ += 1;
            }

         if ( len > 0 )
            error = true;
         }

      if ( !error )
         {
         result = hrs;
         if ( fmin )
            {
            result = result + min / 60.0;

            if ( fsec )
               {
               result = result + sec / 3600.0;

               if ( fdec )
                  result = result + dec / ( 3600.0 * div );
               }
            else if ( fdec )
               {
               result = result + dec / ( 60 * div );
               }
            }
         else if ( fdec )
            {
            result = result + dec / div;
            }
         }
      }
   else
      error = true;

   _error_ = error;
   return sign * result;
}

std::string FGSD_Util::positionToString( double pos, bool lon )
{
   char sign;
   if ( pos < 0 )
      {
      sign = lon ? 'W' : 'S';
      pos = -pos;
      }
   else
      sign = lon ? 'E' : 'N';
   int deg = (int)pos;
   pos -= deg; pos *= 60;
   int min = (int)pos;
   pos -= min; pos *= 60;
   int sec = (int)pos;
   pos -= sec; pos *= 10;
   int milli = (int)pos;
   std::stringstream str;
   str << sign << " " << std::setw( lon ? 3 : 2 ) << std::setfill( '0' ) << deg << ":"
                      << std::setw( 2 ) << std::setfill( '0' ) << min << "'"
                      << std::setw( 2 ) << std::setfill( '0' ) << sec << "."
                      << std::setw( 1 ) << std::setfill( '0' ) << milli << std::ends;

   return str.str();
}

bool FGSD_Util::createPath( const char *path, bool lastIsFile )
{
   bool ret = false;
   std::string path_ = path;
   if ( path_.empty() )
      return ret;
   size_t pos = path_.rfind( '/' );
   if ( pos != std::string::npos && lastIsFile )
      {
      path_.erase( pos );
      pos = path_.rfind( '/' );
      }
   struct stat info;
   if ( stat( path_.c_str(), &info ) != 0 )
      {
      std::string parent = path_.substr( 0, pos );
      ret = createPath( parent.c_str(), false );
      if ( ret )
         {
#ifdef _MSC_VER
         ret = ( _mkdir( path_.c_str() ) == 0 );
#else
         ret = ( mkdir( path_.c_str(), 0755 ) == 0 );
#endif
         }
      }
   else if ( info.st_mode & S_IFDIR )
      ret = true;

    return ret;
}

bool FGSD_Util::copyToFolder( const char *file, const char *folder )
{
   bool ret = false;
#ifdef _MSC_VER
   std::ifstream in( file, std::ios_base::in | std::ios_base::binary );
#else
   std::ifstream in( file );
#endif
   if ( in.is_open() )
      {
      std::string file_ = file;
      size_t pos = file_.rfind( '/' );
      if ( pos != std::string::npos )
         {
         std::string basename = file_.substr( pos + 1, file_.size() - ( pos + 1 ) );
         SGPath folder_; folder_ = folder;
         folder_.append( basename );

#ifdef _MSC_VER
         std::ofstream out( folder_.str().c_str(), std::ios_base::out | std::ios_base::trunc | std::ios_base::binary );
#else
         std::ofstream out( folder_.str().c_str() );
#endif

         if ( out.is_open() )
            {
            ret = true;
            char buffer[ 1024 ];
            size_t nb = sizeof buffer;
            while ( nb == sizeof buffer && in.good() )
               {
               in.read( buffer, sizeof buffer );
               nb = in.gcount();
               out.write( buffer, nb );
               }
            }
         }
      }

   return ret;
}

void FGSD_Util::splitPath( const std::string &path, std::list<std::string> &pathList, char sep )
{
   std::string curdir = path,
               token;
   while ( !curdir.empty() )
      {
      size_t pos = curdir.find( sep );
      if ( pos == std::string::npos )
         {
         token = curdir;
         curdir = "";
         }
      else
         {
         token = curdir.substr( 0, pos );
         curdir.erase( 0, pos + 1 );
         }
      pathList.push_back( token );
      }
}

bool FGSD_Util::isDir( const char *name )
{
   return fl_filename_isdir( name ) != 0;
}

/*
 * Structure used to pass datum parameters between functions.
 */
struct datum  {
	double a;		// Semimajor ellipsoid radius (equatorial radius)
	double b;		// Semiminor ellipsoid radius (polar radius)
	double e_2;		// Eccentricity squared,  e^2 = 2*f - f*f = 1 - b*b / a*a
	double f_inv;		// Inverse flattening,  1/f = a / (a - b)
	double k0;		// scale factor along central meridian
	double a0;		// First coefficient in Redfearn integral expansion
	double a2;		// Second coefficient in Redfearn integral expansion
	double a4;		// Third coefficient in Redfearn integral expansion
	double a6;		// Fourth coefficient in Redfearn integral expansion
};

/*
 * These are the parameters for the Clarke 1866 ellipsoid, which is used with the
 * North American Datum of 1927 (NAD-27) datum.  The NAD-27 used a point on
 * Meade Ranch in Kansas as its reference origin.
 */
#define NAD27_SEMIMAJOR (6378206.4)		// 1866 Clarke ellipsoid, equatorial radius in meters
#define NAD27_SEMIMINOR (6356583.8)		// 1866 Clarke ellipsoid, polar radius in meters
#define NAD27_E_SQUARED (0.006768658)		// 1866 Clarke ellipsoid, eccentricity squared	(e^2 = 2*f - f*f = 1 - b^2 / a^2)
#define NAD27_F_INV (294.9786982)		// 1866 Clarke ellipsoid, inverse flattening	(1/f = a / (a - b))
#define NAD27_A0 (0.99830568187775514389)	// First constant in meridian integral expansion:  1 - (e^2 / 4) - (3 * e^4 / 64) - (5 * e^6 / 256)
#define NAD27_A2 (0.00254255550867060247)	// First constant in meridian integral expansion:  (3/8) * (e^2 + e^4 / 4 + 15 e^6 / 128)
#define NAD27_A4 (0.00000269808452963108)	// First constant in meridian integral expansion:  (15 / 256) * (e^4 + 3 * e^6 / 4)
#define NAD27_A6 (0.00000000353308874387)	// First constant in meridian integral expansion:  35 * e^6 / 3072

/*
 * These are the parameters for the Geodetic Reference System (GRS) 1980 ellipsoid,
 * which is used with the North American Datum of 1983 (NAD-83) datum.  This datum
 * is based on a reference point at the center of the earth, and is defined based on
 * satellite measurements.
 *
 * Other parameters:
 *  Polar radius of curvature (c) 6399593.6259 m
 *  angular velocity (w) 7292115e-11 radians/s
 *  Gravitational Constant (G) 986005e8 m^3/s^2
 *  Flattening (f) 0.00335281068118
 */
#define NAD83_SEMIMAJOR (6378137.0)		// GRS80 ellipsoid, equatorial radius
#define NAD83_SEMIMINOR (6356752.3141)		// GRS80 ellipsoid, polar radius
#define NAD83_E_SQUARED (0.00669438002290)	// GRS80 ellipsoid, e*e	(e^2 = 2*f - f*f = 1 - b^2 / a^2)
#define NAD83_F_INV (298.257222101)		// GRS80 ellipsoid, inverse flattening	(1/f = a / (a - b))
#define NAD83_A0 (0.99832429844458494622)	// First constant in meridian integral expansion:  1 - (e^2 / 4) - (3 * e^4 / 64) - (5 * e^6 / 256)
#define NAD83_A2 (0.00251460707284452333)	// First constant in meridian integral expansion:  (3/8) * (e^2 + e^4 / 4 + 15 e^6 / 128)
#define NAD83_A4 (0.00000263904662023027)	// First constant in meridian integral expansion:  (15 / 256) * (e^4 + 3 * e^6 / 4)
#define NAD83_A6 (0.00000000341804613677)	// First constant in meridian integral expansion:  35 * e^6 / 3072

/*
 * These are the parameters for the World Geodetic System (WGS) 1984
 * ellipsoid.
 * (The WGS-84 ellipsoid is virtually identical to GRS-80.)
 *
 * Other parameters:
 *  Polar radius of curvature (c) 6399593.6258 m
 *  angular velocity (w) 7292115e-11 radians/s
 *  Gravitational Constant (G) 986005e8 m^3/s^2
 *  Flattening (f) 0.00335281066474
 */
#define WGS84_SEMIMAJOR (6378137.0)		// WGS-84 ellipsoid, equatorial radius
#define WGS84_SEMIMINOR (6356752.3142)		// WGS-84 ellipsoid, polar radius
#define WGS84_E_SQUARED (0.00669437999013)	// WGS-84 ellipsoid, e*e	(e^2 = 2*f - f*f = 1 - b^2 / a^2)
#define WGS84_F_INV (298.257223563)		// WGS-84 ellipsoid, inverse flattening	(1/f = a / (a - b))
#define WGS84_A0 (0.99832429845279809866)	// First constant in meridian integral expansion:  1 - (e^2 / 4) - (3 * e^4 / 64) - (5 * e^6 / 256)
#define WGS84_A2 (0.00251460706051444693)	// First constant in meridian integral expansion:  (3/8) * (e^2 + e^4 / 4 + 15 e^6 / 128)
#define WGS84_A4 (0.00000263904659432867)	// First constant in meridian integral expansion:  (15 / 256) * (e^4 + 3 * e^6 / 4)
#define WGS84_A6 (0.00000000341804608657)	// First constant in meridian integral expansion:  35 * e^6 / 3072

/*
 * These are the parameters for Krassovsky
 * ellipsoid.
 * Other parameters:
 * Semimajor Axis: 6378245 
 * Semiminor Axis: 6356863.0187999997
 *
 * DELTA A ( nebo DA ): -108
 * DELTA F ( nebo DF ): +0.0048076
 * DELTA X ( nebo DX ): +23
 * DELTA Y ( nebo DY ): -124
 * DELTA Z ( nebo DZ ): -84
 *
 */
#define S42_SEMIMAJOR (6378245.0)		// S42 ellipsoid, equatorial radius
#define S42_SEMIMINOR (6356863.0187999997)		// S42 ellipsoid, polar radius
#define S42_E_SQUARED (0.006693421614543)	// S42 ellipsoid, e*e	(e^2 = 2*f - f*f = 1 - b^2 / a^2)
#define S42_F_INV (298.300000376)		// S42 ellipsoid, inverse flattening	(1/f = a / (a - b))
#define S42_A0 (0.9983245444684)	// First constant in meridian integral expansion:  1 - (e^2 / 4) - (3 * e^4 / 64) - (5 * e^6 / 256)
#define S42_A2 (0.002514233371121)	// First constant in meridian integral expansion:  (3/8) * (e^2 + e^4 / 4 + 15 e^6 / 128)
#define S42_A4 (2.625199120151e-06)	// First constant in meridian integral expansion:  (15 / 256) * (e^4 + 3 * e^6 / 4)
#define S42_A6 (2.286859905315e-11)	
// First constant in meridian integral expansion:  35 * e^6 / 3072



/*
 * To convert to/from UTM coordinates we need to know the scale factor on the central meridian.
 * For UTM, this is always 0.9996.
 */
#define UTM_K0 (0.9996)			// UTM Scale factor on the central meridian

struct utm_zones  {
	long	zone;
	double	central_meridian;
	double	low_boundary;
	double	high_boundary;
} utm_zones[61] =  {
	 0,	   0.0,	   0.0,	   0.0,	// Dummy entry so that the zone indices are correct
	 1,	-177.0,	-180.0,	-174.0,
	 2,	-171.0,	-174.0,	-168.0,
	 3,	-165.0,	-168.0,	-162.0,
	 4,	-159.0,	-162.0,	-156.0,
	 5,	-153.0,	-156.0,	-150.0,
	 6,	-147.0,	-150.0,	-144.0,
	 7,	-141.0,	-144.0,	-138.0,
	 8,	-135.0,	-138.0,	-132.0,
	 9,	-129.0,	-132.0,	-126.0,
	10,	-123.0,	-126.0,	-120.0,
	11,	-117.0,	-120.0,	-114.0,
	12,	-111.0,	-114.0,	-108.0,
	13,	-105.0,	-108.0,	-102.0,
	14,	- 99.0,	-102.0,	- 96.0,
	15,	- 93.0,	-096.0,	- 90.0,
	16,	- 87.0,	-090.0,	- 84.0,
	17,	- 81.0,	-084.0,	- 78.0,
	18,	- 75.0,	-078.0,	- 72.0,
	19,	- 69.0,	-072.0,	- 66.0,
	20,	- 63.0,	-066.0,	- 60.0,
	21,	- 57.0,	-060.0,	- 54.0,
	22,	- 51.0,	-054.0,	- 48.0,
	23,	- 45.0,	-048.0,	- 42.0,
	24,	- 39.0,	-042.0,	- 36.0,
	25,	- 33.0,	-036.0,	- 30.0,
	26,	- 27.0,	-030.0,	- 24.0,
	27,	- 21.0,	-024.0,	- 18.0,
	28,	- 15.0,	-018.0,	- 12.0,
	29,	-  9.0,	-012.0,	-  6.0,
	30,	-  3.0,	-006.0,	   0.0,
	31,	   3.0,	 000.0,	   6.0,
	32,	   9.0,	 006.0,	  12.0,
	33,	  15.0,	 012.0,	  18.0,
	34,	  21.0,	 018.0,	  24.0,
	35,	  27.0,	 024.0,	  30.0,
	36,	  33.0,	 030.0,	  36.0,
	37,	  39.0,	 036.0,	  42.0,
	38,	  45.0,	 042.0,	  48.0,
	39,	  51.0,	 048.0,	  54.0,
	40,	  57.0,	 054.0,	  60.0,
	41,	  63.0,	 060.0,	  66.0,
	42,	  69.0,	 066.0,	  72.0,
	43,	  75.0,	 072.0,	  78.0,
	44,	  81.0,	 078.0,	  84.0,
	45,	  87.0,	 084.0,	  90.0,
	46,	  93.0,	 090.0,	  96.0,
	47,	  99.0,	 096.0,	 102.0,
	48,	 105.0,	 102.0,	 108.0,
	49,	 111.0,	 108.0,	 114.0,
	50,	 117.0,	 114.0,	 120.0,
	51,	 123.0,	 120.0,	 126.0,
	52,	 129.0,	 126.0,	 132.0,
	53,	 135.0,	 132.0,	 138.0,
	54,	 141.0,	 138.0,	 144.0,
	55,	 147.0,	 144.0,	 150.0,
	56,	 153.0,	 150.0,	 156.0,
	57,	 159.0,	 156.0,	 162.0,
	58,	 165.0,	 162.0,	 168.0,
	59,	 171.0,	 168.0,	 174.0,
	60,	 177.0,	 174.0,	 180.0,
};

/*
 * The following two functions use Redfearn's formulas to calculate the
 * forward and inverse projections between UTM coordinates and geographic
 * (latitude/longitude) coordinates.
 *
 * Given some parameters for the selected ellipsoid, Redfearn's formulas
 * allow one to translate back and forth between UTM
 * coordinates and latitude/longitude coordinates.
 * Before we examine Redfearn's formulas, here are some preliminary notes.
 *
 * These formulas were apparently originally published in 1948:
 *
 *	Redfearn, J.C.B., "Transverse Mercator Formulae", Empire Survey Review, 69, 1948, 318-322.
 *
 * I was unable to find a copy of this reference for verification, but did
 * find several other documents that described the formulas.  From them, I
 * pieced together the formulas here.
 *
 * A good reference for projection calculations of all kinds is supposed to be:
 *      Snyder, John P., Map Projections -- A Working Manual:
 *      U.S. Geological Survey Professional Paper 1395,
 *      United States Government Printing Office, Washington D.C., 1987.
 * Although I haven't personnally had a chance to read it, it is frequently
 * recommended on the Internet.
 *
 * A software package for doing projections is the PROJ package.  It is
 * available on the Internet.  I chose not to use it because my needs are
 * limited, and it was easier to just write the software I need than integrate
 * PROJ.  It is, however, a fine package, as far as I can tell.
 *
 * As far as I know, most or all of the currently-available 7.5min USGS data assumes
 * the Clarke 1866 ellipsoid with the North American Datum of 1927 (nad-27).
 * According to the DEM standards document, they are nad-27 if they have the
 * old-format Type A record.  The new Type A record contains a field to specify the datum.
 * I gather that the data will ultimately be re-referenced to the new GRS80
 * ellipsoid and nad-83.
 *
 * The 1-degree DEMs may or may not be in the new WGS 84 datum.  The standards document
 * says that recomputed data have been made available to the USGS, but doesn't say
 * if these new data are what is available for download.  We probably don't care,
 * for the time being, because the 1-degree data come in latitude/longitude format,
 * and we don't convert them to any other form.
 *
 * Parameters for various ellipsoids are given in drawmap.h.
 *
 * Now, on to Redfearn's formulas.
 *
 *    Note:  In the equations that follow, y is the true northing, utm_y is the northing with
 *           the false northing (10,000,000m) added in, x is the true easting, utm_y is the
 *           easting with the false easting (500,000m) added in.
 *
 * They begin with the calculation of the length of arc of a meridian (a great
 * circle passing through the poles).  The formula (in Macsyma format) is:
 *
 *     m = a * (1 - e^2) * int(1 - (e^2 sin^2(lat))^(-3/2), lat, lat1, lat2)
 *
 * (where the caret sign, '^', represents exponentiation and 'int' represents
 * a definite integral of the given expression over the variable lat, from lat1 to lat2).
 *
 * We are interested in the case where lat1 = 0 (and the integral starts at the equator).
 *
 * This integral is normally calculated via a series expansion:
 *
 *     m = a * (A0 * lat - A2 * sin(2*lat) + A4 * sin(4*lat) - A6 * sin(6*lat))
 *
 *         A0 = 1 - (e^2 / 4) - (3 * e^4 / 64) - (5 * e^6 / 256)
 *         A2 = (3/8) * (e^2 + e^4 / 4 + 15 e^6 / 128)
 *         A4 = (15 / 256) * (e^4 + 3 * e^6 / 4)
 *         A6 = 35 * e^6 / 3072
 *
 * With the value of m in hand, we move on to calculate the foot-point latitude, lat', which is
 * the latitude for which m = y / k0
 * where y is the true northing (which is just the northing in the northern climes, but is the
 * nominal northing - 10,000,000 in the southern hemisphere) and k0 is the centeral-meridian
 * scale factor.
 *
 * The foot-point latitude is found as follows:
 *
 *      n = (a - b) / (a + b) = f / (2 - f)
 *      G = a * (1 - n) * (1 - n^2) * (1 + (9/4) * n^2 + (225/64) * n^4) * (pi / 180)
 *      sigma = (m * pi) / (180 * G)
 *      lat' = sigma + ((3*n/2) - (27*n^3/32)) * sin(2*sigma) + ((21*n^2/16) - (55*n^4/32)) * sin(4*sigma) + (151*n^3/96) * sin(6*sigma) + (1097*n^4/512) * sin(8*sigma)        in units of radians
 *
 * For the inverse projection, where the latitude is to-be-determined, there may
 * be some snazzy way to find lat', but I chose to do it iteratively using
 * Newton's nethod.  In my limited testing, this approach appears to converge
 * quite rapidly, in about 2 or 3 iterations.  (I have not, however, tested the convergence
 * rate rigorously.)
 *
 * Next, we need the radii of curvature, found from the formulas.
 *
 *      rho = a * (1 - e^2) / (1 - e^2 * sin^2(lat))^(3/2)
 *      nu  = a / (1 - e^2 * sin^2(lat))^(1/2)
 *      phi = nu / rho
 *
 * These are general formulas.  When we evaluate them specifically for the foot-point
 * latitude, then we prime each of them:
 *
 *      rho' = a * (1 - e^2) / (1 - e^2 * sin^2(lat'))^(3/2)
 *      nu'  = a / (1 - e^2 * sin^2(lat'))^(1/2)
 *      phi' = nu' / rho'
 *
 *
 * This completes the preliminary calculations.  Now, on to the actual conversions.
 *
 *
 * lat/long to UTM, performed by function redfearn():
 *    t = tan(lat)
 *    omega = longitude * pi / 180 - central_meridian
 *
 *    utm_x:
 *    x = k0 * nu * omega * cos(lat) * (1 +
 *                                       (omega^2 / 6) * cos^2(lat) * (phi - t^2) +
 *                                       (omega^4 / 120) * cos^4(lat) * (4 * phi^3 * (1 - 6 * t^2) + phi^2 * (1 + 8*t^2) - 2*phi*t^2 + t^4) +
 *                                       (omega^6 / 5040) * cos^6(lat) * (61 - 479*t^2 + 179*t^4 - t^6))
 *    utm_x = x + 500,000
 *
 *    utm_y:
 *    y = k0 * (m + (omega^2 / 2) * nu * sin(lat) * cos(lat) +
 *                   (omega^4 / 24) * nu * sin(lat) * cos^3(lat) * (4 * phi^2 + phi - t^2) +
 *                   (omega^6 / 720) * nu * sin(lat) * cos^5(lat) * (8 * phi^4 * (11 - 24*t^2) -
 *                    28 * phi^3 * (1 - 6*t^2) + phi^2 * (1 - 32*t^2) - 2*phi*t^2 + t^4) +
 *                   (omega^8 / 40320) * nu * sin(lat) * cos^7(lat) * (1385 - 3111*t^2 + 543*t^4 - t^6))
 *    utm_y = y + 10,000,000
 *
 *    grid convergence:
 *    gamma = -omega * sin(lat) -
 *            (omega^3 / 3) * sin(lat) * cos^2(lat) * (2*phi^2 - phi) -
 *            (omega^5 / 15) * sin(lat) * cos^4(lat) * (phi^4 * (11 - 24*t^2) - phi^3 * (11 - 36*t^2) +
 *             2*phi^2 * (1 - 7*t^2) + phi*t^2) -
 *            (omega^7 / 315) * sin(lat) * cos^6(lat) * (17 - 26*t^2 + 2*t^4)
 *
 *    point scale factor:
 *    k = k0 * (1 + (omega^2 / 2) * phi * cos^2(lat) +
 *                  (omega^4 / 24) * cos^4(lat) * (4*phi^3 * (1 - 6*t^2) + phi^2 * (1 + 24*t^2) - 4*phi*t^2) +
 *                  (omega^6 / 720) * cos^6(lat) * (61 - 148*t^2 + 16*t^4))
 *
 *
 * UTM to lat/long, performed by function redfearn_inverse():
 *    Note that the ugly construct, phi'^3, represents phi' taken to the third power.  It's
 *    ugly, but adding enough parentheses to clarify it would arguably be uglier still.
 *
 *    x = utm_x - 500,000
 *    d = x / (k0 * nu')
 *    t' = tan(lat')
 *
 *    latitude:
 *    lat = lat' - (nu' * t' / rho') * (d^2 / 2) +
 *                 (nu' * t' / rho') * (d^4 / 24) * (-4*phi'^2 + 9*phi' * (1 - t'^2) + 12*t'^2) -
 *                 (nu' * t' / rho') * (d^6 / 720) * (8*phi'^4 * (11 - 24*t'^2) - 12*phi'^3 * (21 - 71*t'^2) +
 *                  15*phi'^2 * (15 - 98*t'^2 + 15*t'^4) + 180*phi' * (5*t'^2 - 3*t'^4) + 360*t'^4) +
 *                 (nu' * t' / rho') * (d^8 / 40320) * (1385 + 3633*t'^2 + 4095*t'^4 + 1575*t'^6)
 *
 *    longitude:
 *    omega = d * sec(lat') -
 *            (d^3 / 6) * sec(lat') * (phi' + 2*t'^2) +
 *            (d^5 / 120) * sec(lat') * (-4*phi'^3 * (1 - 6*t'^2) + phi'^2 * (9 - 68*t'^2) + 72*phi'*t'^2 + 24*t'^4) -
 *            (d^7 / 5040) * sec(lat') * (61 + 662*t'^2 + 1320*t'^4 + 720*t'^6)
 *
 *    long = central_meridian + omega * 180 / pi
 *
 *    grid convergence:
 *    gamma = -t' * d +
 *            (t' * d^3 / 3) * (-2*phi'^2 + 3*phi' + t'^2) -
 *            (t' * d^5 / 15) * (phi'^4 * (11 - 24*t'^2) - 3*phi'^3 * (8 - 23*t'^2) +
 *             5*phi'^2 * (3 - 14*t'^2) + 30*phi'*t'^2 + 3*t'^4) +
 *            (t' * d^7 / 315) * (17 + 77*t'^2 + 105*t'^4 + 45*t'^6)
 *
 *    point scale factor:
 *    dd = x^2 / (k0^2 * rho' * nu')
 *    k = k0 * (1 + dd/2 + (dd^2 / 24) * (4*phi' * (1 - 6*t'^2) - 3 * (1 - 16*t'^2) - 24*t'^2 / phi') + dd^3 / 720)
 *
 * Now, on to the actual code.
 *
 * Note:  redfearn_inverse() returns 0 if the conversion appears to be successful, nonzero otherwise.
 *
 * Note further:  This function has not been tuned for efficiency.
 */
static long
redfearn_inverse(struct datum *datum, double utm_x, double utm_y, long zone, double *latitude, double *longitude)
{
	double x, y;	// UTM coordinates with false easting and northing removed and scale factor applied.
	double lat_pm;	// foot-point latitude
	double d;
	double t_pm;
	double m, m_pm;
	double nu_pm;
	double rho_pm;
	double phi_pm;
	double slat, slat_2, clat, d_2, d_3, d_4, d_5, d_6, d_7, d_8;
	double t_pm_2, t_pm_4, t_pm_6, phi_pm_2, phi_pm_3, phi_pm_4;
	long i;

	x = (utm_x - 500000.0) / datum->k0;

	if ((zone > 60) || (zone == 0) || (zone < -60)) {
		return -1;
	}
	if (zone < 0)  {
		/* southern hemisphere */
		zone = -zone;
		y = (utm_y - 10000000.0) / datum->k0;
		lat_pm = -M_PI / 4.0;
	}
	else  {
		y = utm_y / datum->k0;
		lat_pm = M_PI / 4.0;
	}

	/*
	 * Find lat_pm, via iterative Newton's method.
	 * The goal is to find lat_pm, such that m == y, or equivalently
	 * to find a root of m-y.
	 */
	for (i = 0; i < 100; i++)  {
		m = datum->a * (datum->a0 * lat_pm - datum->a2 * sin(2.0 * lat_pm) + datum->a4 * sin(4.0 * lat_pm) - datum->a6 * sin(6.0 * lat_pm)) - y;
		m_pm = datum->a * (datum->a0 - datum->a2 * 2.0 * cos(2.0 * lat_pm) + datum->a4 * 4.0 * cos(4.0 * lat_pm) - datum->a6 * 6.0 * cos(6.0 * lat_pm));
		if (fabs(m / m_pm) < 1.0e-12)  {
			break;
		}
		lat_pm -= m / m_pm;
	}

	slat = sin(lat_pm);
	slat_2 = slat * slat;
	clat = sqrt(1.0 - slat_2);
	t_pm = slat / clat;

	nu_pm = datum->a / sqrt(1.0 - datum->e_2 * slat_2);
	rho_pm = datum->a * (1.0 - datum->e_2) / pow(1.0 - datum->e_2 * slat_2, 1.5);
	phi_pm = nu_pm / rho_pm;
	d = x / nu_pm;

	d_2 = d * d;
	d_3 = d_2 * d;
	d_4 = d_3 * d;
	d_5 = d_4 * d;
	d_6 = d_5 * d;
	d_7 = d_6 * d;
	d_8 = d_7 * d;
	t_pm_2 = t_pm * t_pm;
	t_pm_4 = t_pm_2 * t_pm_2;
	t_pm_6 = t_pm_2 * t_pm_4;
	phi_pm_2 = phi_pm * phi_pm;
	phi_pm_3 = phi_pm_2 * phi_pm;
	phi_pm_4 = phi_pm_3 * phi_pm;

	*latitude = lat_pm - (nu_pm * t_pm / rho_pm) * ((d_2 / 2.0) -
			 (d_4 / 24.0) * (-4.0 * phi_pm_2 + 9.0 * phi_pm * (1.0 - t_pm_2) + 12.0 * t_pm_2) +
			 (d_6 / 720.0) * (8.0 * phi_pm_4 * (11.0 - 24.0 * t_pm_2) - 12.0 * phi_pm_3 *
				(21.0 - 71.0 * t_pm_2) + 15.0 * phi_pm_2 * (15.0 - 98.0 * t_pm_2 + 15.0 * t_pm_4) +
				180.0 * phi_pm * (5.0 * t_pm_2 - 3.0 * t_pm_4) + 360.0 * t_pm_4) -
			 (d_8 / 40320.0) * (1385.0 + 3633.0 * t_pm_2 + 4095.0*t_pm_4 + 1575.0*t_pm_6));
	*longitude = d / clat - (d_3 / 6.0) * (phi_pm + 2.0 * t_pm_2) / clat +
		     (d_5 / 120.0) * (-4.0 * phi_pm_3 * (1.0 - 6.0 * t_pm_2) + phi_pm_2 * (9.0 - 68.0 * t_pm_2) +
			72.0 * phi_pm * t_pm_2 + 24.0 * t_pm_4) / clat -
		     (d_7 / 5040.0) * (61.0 + 662.0 * t_pm_2 + 1320.0 * t_pm_4 + 720.0 * t_pm_6) / clat;
	*latitude = *latitude * 180.0 / M_PI;
	*longitude = utm_zones[zone].central_meridian + *longitude * 180.0 / M_PI;

	return 0;
}

static long
redfearn(struct datum *datum, double *utm_x, double *utm_y, int *zone,
		double latitude, double longitude, int east_west)
{
	double t;
	double m;
	double nu;
	double rho;
	double phi;
	double o;
	double t_2, t_4, t_6, o_2, o_3, o_4, o_5, o_6, o_7, o_8;
	double slat, slat_2, clat, clat_2, clat_3, clat_4, clat_5, clat_6, clat_7;
	double phi_2, phi_3, phi_4;
	long i;
	long south_flag = 0;	// If this flag is nonzero, the point is in the southern hemisphere.


	/*
	 * Note:  Originally the following check was
	 *
	 *    if ((latitude > 84.0) || (latitude < -80.0)) {
	 *
	 * because that is the valid range for UTM projections.
	 * However, in order to allow projections to be done for the
	 * GTOPO30 data, I changed the range to the full -90 to 90 span.
	 */
	if ((latitude > 90.0) || (latitude < -90.0)) {
		return -1;
	}
	if ((longitude > 180.0) || (longitude < -180.0)) {
		return -1;
	}

	/*
	 * Determine the zone.
	 *
	 * If the point falls on the boundary between zones,
	 * then use the parameter east_west to choose between zones.
	 * If east_west is nonzero, then choose the eastern zone.
	 * If east_west is 0, then choose the western zone.
	 * If we are on the boundary between zone 1 and zone 60, then ignore
	 * east_west and choose the zone based on the passed longitude.
	 */
	if (longitude == utm_zones[1].low_boundary)  {
		*zone = 1;
	}
	else if (longitude == utm_zones[60].high_boundary)  {
			*zone = 60;
	}
	else  {
		for (i = 1; i <= 60; i++)  {
			if (longitude == utm_zones[i].high_boundary)  {
				if (east_west == 0)  {
					*zone = i;
				}
				else  {
					*zone = i + 1;
				}
				break;
			}
			else if ((longitude > utm_zones[i].low_boundary) && (longitude < utm_zones[i].high_boundary))  {
				*zone = i;
				break;
			}
		}
	}

	o = (longitude - utm_zones[*zone].central_meridian) * M_PI / 180.0;

	latitude *= M_PI / 180.0;
	longitude *= M_PI / 180.0;
	slat = sin(latitude);
	slat_2 = slat * slat;
	clat = sqrt(1.0 - slat_2);	// cos(latitude)
	t = slat / clat;		// tan(latitude)

	t_2 = t * t;
	t_4 = t_2 * t_2;
	t_6 = t_2 * t_4;
	o_2 = o * o;
	o_3 = o_2 * o;
	o_4 = o_2 * o_2;
	o_5 = o_4 * o;
	o_6 = o_4 * o_2;
	o_7 = o_6 * o;
	o_8 = o_4 * o_4;
	clat_2 = clat * clat;
	clat_3 = clat_2 * clat;
	clat_4 = clat_2 * clat_2;
	clat_5 = clat_4 * clat;
	clat_6 = clat_4 * clat_2;
	clat_7 = clat_6 * clat;

	m = datum->a * (datum->a0 * latitude - datum->a2 * sin(2.0 * latitude) + datum->a4 * sin(4.0 * latitude) - datum->a6 * sin(6.0 * latitude));

	nu = datum->a / sqrt(1.0 - datum->e_2 * slat_2);
	rho = datum->a * (1.0 - datum->e_2) / pow(1.0 - datum->e_2 * slat_2, 1.5);
	phi = nu / rho;

	phi_2 = phi * phi;
	phi_3 = phi_2 * phi;
	phi_4 = phi_2 * phi_2;

	*utm_x = 500000.0 + datum->k0 * nu * clat * (o + (o_3 / 6.0) * clat_2 * (phi - t_2) +
			(o_5 / 120.0) * clat_4 * (4.0 * phi_3 * (1.0 - 6.0 * t_2) +
			phi_2 * (1.0 + 8.0 * t_2) - 2.0 * phi * t_2 + t_4) +
			(o_7 / 5040.0) * clat_6 * (61.0 - 479.0 * t_2 + 179.0 * t_4 - t_6));
	*utm_y = datum->k0 * (m + (o_2 / 2.0) * nu * slat * clat +
			(o_4 / 24.0) * nu * slat * clat_3 * (4.0 * phi_2 + phi - t_2) +
			(o_6 / 720.0) * nu * slat * clat_5 * (8.0 * phi_4 * (11.0 - 24.0 * t_2) -
			28.0 * phi_3 * (1.0 - 6.0 * t_2) + phi_2 * (1.0 - 32.0 * t_2) - 2.0 * phi * t_2 + t_4) +
			(o_8 / 40320.0) * nu * slat * clat_7 * (1385.0 - 3111.0*t_2 + 543.0*t_4 - t_6));

	if (latitude < 0)  {
		/* In the southern hemisphere, we return a negative zone number. */
		*zone = -*zone;
		*utm_y += 10000000.0;
	}

	return 0;
}

bool FGSD_Util::getWGS84PositionFromUTM( Datum _datum, double utm_x, double utm_y, long zone, double &latitude, double &longitude )
{
   bool ret = true;
   struct datum datum;
   if (_datum == D_NAD27) {
      /* Fill in the datum parameters for NAD-27. */
      datum.a = NAD27_SEMIMAJOR;
      datum.b = NAD27_SEMIMINOR;
      datum.e_2 = NAD27_E_SQUARED;
      datum.f_inv = NAD27_F_INV;
      datum.k0 = UTM_K0;
      datum.a0 = NAD27_A0;
      datum.a2 = NAD27_A2;
      datum.a4 = NAD27_A4;
      datum.a6 = NAD27_A6;
   } else if (_datum == D_NAD83) {
      /* Fill in the datum parameters for NAD-83. */
      datum.a = NAD83_SEMIMAJOR;
      datum.b = NAD83_SEMIMINOR;
      datum.e_2 = NAD83_E_SQUARED;
      datum.f_inv = NAD83_F_INV;
      datum.k0 = UTM_K0;
      datum.a0 = NAD83_A0;
      datum.a2 = NAD83_A2;
      datum.a4 = NAD83_A4;
      datum.a6 = NAD83_A6;
   } else if (_datum == D_WGS84) {
      /* Fill in the datum parameters for WGS-84. */
      datum.a = WGS84_SEMIMAJOR;
      datum.b = WGS84_SEMIMINOR;
      datum.e_2 = WGS84_E_SQUARED;
      datum.f_inv = WGS84_F_INV;
      datum.k0 = UTM_K0;
      datum.a0 = WGS84_A0;
      datum.a2 = WGS84_A2;
      datum.a4 = WGS84_A4;
      datum.a6 = WGS84_A6;
   } else if (_datum == D_S42) {
      /* Fill in the datum parameters for S-42 Krassovsky elipsoid. */
      datum.a = S42_SEMIMAJOR;
      datum.b = S42_SEMIMINOR;
      datum.e_2 = S42_E_SQUARED;
      datum.f_inv = S42_F_INV;
      datum.k0 = 1.0;
      datum.a0 = S42_A0;
      datum.a2 = S42_A2;
      datum.a4 = S42_A4;
      datum.a6 = S42_A6;
   } else {
      ret = false;
   }

   if ( ret ) {
      redfearn_inverse( &datum, utm_x, utm_y, zone, &latitude, &longitude );
   }
   return ret;
}

bool FGSD_Util::getUTMPosition( Datum _datum, double &utm_x, double &utm_y, int &zone, double latitude, double longitude )
{
   bool ret = true;
   struct datum datum;
   if ( _datum == D_NAD27) {
      /* Fill in the datum parameters for NAD-27. */
      datum.a = NAD27_SEMIMAJOR;
      datum.b = NAD27_SEMIMINOR;
      datum.e_2 = NAD27_E_SQUARED;
      datum.f_inv = NAD27_F_INV;
      datum.k0 = UTM_K0;
      datum.a0 = NAD27_A0;
      datum.a2 = NAD27_A2;
      datum.a4 = NAD27_A4;
      datum.a6 = NAD27_A6;
   } else if (_datum == D_NAD83) {
      /* Fill in the datum parameters for NAD-83. */
      datum.a = NAD83_SEMIMAJOR;
      datum.b = NAD83_SEMIMINOR;
      datum.e_2 = NAD83_E_SQUARED;
      datum.f_inv = NAD83_F_INV;
      datum.k0 = UTM_K0;
      datum.a0 = NAD83_A0;
      datum.a2 = NAD83_A2;
      datum.a4 = NAD83_A4;
      datum.a6 = NAD83_A6;
   } else if (_datum == D_WGS84) {
      /* Fill in the datum parameters for WGS-84. */
      datum.a = WGS84_SEMIMAJOR;
      datum.b = WGS84_SEMIMINOR;
      datum.e_2 = WGS84_E_SQUARED;
      datum.f_inv = WGS84_F_INV;
      datum.k0 = UTM_K0;
      datum.a0 = WGS84_A0;
      datum.a2 = WGS84_A2;
      datum.a4 = WGS84_A4;
      datum.a6 = WGS84_A6;
   } else if (_datum == D_S42) {
      /* Fill in the datum parameters for S-42 Krassovsky elipsoid. */
      datum.a = S42_SEMIMAJOR;
      datum.b = S42_SEMIMINOR;
      datum.e_2 = S42_E_SQUARED;
      datum.f_inv = S42_F_INV;
      datum.k0 = 1.0;
      datum.a0 = S42_A0;
      datum.a2 = S42_A2;
      datum.a4 = S42_A4;
      datum.a6 = S42_A6;
   } else {
      ret = false;
   }

   if ( ret ) {
      redfearn( &datum, &utm_x, &utm_y, &zone, latitude, longitude, 0);
   }

   return ret;
}