/***************************************************************************/
/* RSC IDENTIFIER: Datum
 *
 * ABSTRACT
 *
 *    This component provides datum shifts for a large collection of local
 *    datums, WGS72, and WGS84.  A particular datum can be accessed by using its 
 *    standard 5-letter code to find its index in the datum table.  The index 
 *    can then be used to retrieve the name, type, ellipsoid code, and datum 
 *    shift parameters, and to perform shifts to or from that datum.
 *    
 *    By sequentially retrieving all of the datum codes and/or names, a menu
 *    of the available datums can be constructed.  The index values resulting
 *    from selections from this menu can then be used to access the parameters
 *    of the selected datum, or to perform datum shifts involving that datum.
 *
 *    This component supports both 3-parameter local datums, for which only X,
 *    Y, and Z translations relative to WGS 84 have been defined, and 
 *    7-parameter local datums, for which X, Y, and Z rotations, and a scale 
 *    factor, are also defined.  It also includes entries for WGS 84 (with an
 *    index of 0), and WGS 72 (with an index of 1), but no shift parameter 
 *    values are defined for these.
 *
 *    This component provides datum shift functions for both geocentric and
 *    geodetic coordinates.  WGS84 is used as an intermediate state when
 *    shifting from one local datum to another.  When geodetic coordinates are
 *    given Molodensky's method is used, except near the poles where the 3-step
 *    step method is used instead.  Specific algorithms are used for shifting 
 *    between WGS72 and WGS84.
 *
 *    This component depends on two data files, named 3_param.dat and 
 *    7_param.dat, which contain the datum parameter values.  Copies of these
 *    files must be located in the directory specified by the value of the 
 *    environment variable "DATUM_DATA", if defined, or else in the current 
 *    directory whenever a program containing this component is executed. 
 *
 *    Additional datums can be added to these files, either manually or using 
 *    the Create_Datum function.  However, if a large number of datums are 
 *    added, the datum table array sizes in this component will have to be 
 *    increased.
 *
 *    This component depends on two other components: the Ellipsoid component
 *    for access to ellipsoid parameters; and the Geocentric component for 
 *    conversions between geodetic and geocentric coordinates.
 *
 * ERROR HANDLING
 *
 *    This component checks for input file errors and input parameter errors.
 *    If an invalid value is found, the error code is combined with the current
 *    error code using the bitwise or.  This combining allows multiple error
 *    codes to be returned. The possible error codes are:
 *
 *  DATUM_NO_ERROR                  : No errors occurred in function
 *  DATUM_NOT_INITIALIZED_ERROR     : Datum module has not been initialized
 *  DATUM_7PARAM_FILE_OPEN_ERROR    : 7 parameter file opening error
 *  DATUM_7PARAM_FILE_PARSING_ERROR : 7 parameter file structure error
 *  DATUM_7PARAM_OVERFLOW_ERROR     : 7 parameter table overflow
 *  DATUM_3PARAM_FILE_OPEN_ERROR    : 3 parameter file opening error
 *  DATUM_3PARAM_FILE_PARSING_ERROR : 3 parameter file structure error
 *  DATUM_3PARAM_OVERFLOW_ERROR     : 3 parameter table overflow
 *  DATUM_INVALID_INDEX_ERROR       : Index out of valid range (less than one
 *                                      or more than Number_of_Datums)
 *  DATUM_INVALID_SRC_INDEX_ERROR   : Source datum index invalid
 *  DATUM_INVALID_DEST_INDEX_ERROR  : Destination datum index invalid
 *  DATUM_INVALID_CODE_ERROR        : Datum code not found in table
 *  DATUM_LAT_ERROR                 : Latitude out of valid range (-90 to 90)
 *  DATUM_LON_ERROR                 : Longitude out of valid range (-180 to
 *                                    360)
 *  DATUM_SIGMA_ERROR               : Standard error values must be positive
 *                                    (or -1 if unknown)
 *  DATUM_DOMAIN_ERROR              : Domain of validity not well defined
 *  DATUM_ELLIPSE_ERROR             : Error in ellipsoid module
 *  DATUM_NOT_USERDEF_ERROR         : Datum code is not user defined - cannot 
 *                                    be deleted
 *
 *
 * REUSE NOTES
 *
 *    Datum is intended for reuse by any application that needs access to 
 *    datum shift parameters relative to WGS 84.
 *
 *    
 * REFERENCES
 *
 *    Further information on Datum can be found in the Reuse Manual.
 *
 *    Datum originated from :  U.S. Army Topographic Engineering Center (USATEC)
 *                             Geospatial Information Division (GID)
 *                             7701 Telegraph Road
 *                             Alexandria, VA  22310-3864
 *
 * LICENSES
 *
 *    None apply to this component.
 *
 * RESTRICTIONS
 *
 *    Datum has no restrictions.
 *
 * ENVIRONMENT
 *
 *    Datum was tested and certified in the following environments:
 *
 *    1. Solaris 2.5 with GCC 2.8.1
 *    2. MS Windows 95 with MS Visual C++ 6
 *
 * MODIFICATIONS
 *
 *    Date              Description
 *    ----              -----------
 *    03/30/97          Original Code
 *    05/28/99          Added user-definable datums (for JMTK)
 *                      Added datum domain of validity checking (for JMTK)
 *                      Added datum shift accuracy calculation (for JMTK) 
 */


/***************************************************************************/
/*
 *                               INCLUDES
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include "ellipse.h"
#include "geocent.h"
#include "datum.h"
/* 
 *    stdio.h    - standard C input/output library
 *    stdlib.h   - standard C general utilities library 
 *    string.h   - standard C string handling library
 *    ctype.h    - standard C character handling library
 *    math.h     - standard C mathematics library
 *    ellipse.h  - used to get ellipsoid parameters
 *    geocent.h  - used to convert between geodetic and geocentric coordinates
 *    datum.h    - for prototype error ehecking and error codes
 */


/***************************************************************************/
/*
 *                               DEFINES
 */

#define SECONDS_PER_RADIAN     206264.8062471;   /* Seconds in a radian */
#define PI                          3.14159265358979323e0
#define MIN_LAT                  (-PI/2.0)
#define MAX_LAT                  (+PI/2.0)
#define MIN_LON                   -PI
#define MAX_LON                  (2.0 * PI)
#define DATUM_CODE_LENGTH           7
#define DATUM_NAME_LENGTH          33
#define ELLIPSOID_CODE_LENGTH       3
#define MAX_7PARAM                 25
#define MAX_3PARAM                250
#define MOLODENSKY_MAX            (89.75 * PI / 180.0) /* Polar limit */
#define FILENAME_LENGTH           128
#define FALSE 0
#define TRUE  1


/***************************************************************************/
/*
 *                          GLOBAL DECLARATIONS
 */

typedef struct Datum_Table_Row
{
  Datum_Type Type;
  char Code[DATUM_CODE_LENGTH];
  char Name[DATUM_NAME_LENGTH];
  char Ellipsoid_Code[ELLIPSOID_CODE_LENGTH];
  double Parameters[7];  /* interface for 3 and 7 Parameters */
  double Sigma_X;        /* standard error along X axis */
  double Sigma_Y;        /* standard error along Y axis */
  double Sigma_Z;        /* standard error along Z axis */
  double West_longitude; /* western boundary of validity rectangle */
  double East_longitude; /* eastern boundary of validity rectangle */
  double South_latitude; /* southern boundary of validity rectangle */
  double North_latitude; /* northern boundary of validity rectangle */
  long User_Defined;      /* Identifies a user defined datum */
} Datum_Row; /* defines a single entry in the datum table */

/* World Geodetic Systems */
static Datum_Row WGS84;
static Datum_Row WGS72;
const char *WGS84_Datum_Code = "WGE";
const char *WGS72_Datum_Code = "WGC";
static long Datum_WGS84_Index = 1; /* Index of WGS84 in datum table */
static long Datum_WGS72_Index = 2; /* Index of WGS72 in datum table */

static Datum_Row Datum_Table_3Param[MAX_3PARAM]; /* Array of 3Param datums */
static Datum_Row Datum_Table_7Param[MAX_7PARAM]; /* Array of 7Param datums */
static Datum_Row *Datum_Table[2 + MAX_3PARAM + MAX_7PARAM]; /* Datum pointer array, for sorting */
static long Datum_3Param_Count = 0;
static long Datum_7Param_Count = 0;
static long Number_of_Datums = 0; /* Count for datum table */

static long Datum_Initialized = 0; /* Indicates successful initialization */


/***************************************************************************/
/*
 *                              FUNCTIONS     
 */
/*Forward function declarations */
void Geodetic_Shift_WGS84_To_WGS72( const double WGS84_Lat,
                                    const double WGS84_Lon,
                                    const double WGS84_Hgt,
                                    double *WGS72_Lat,
                                    double *WGS72_Lon,
                                    double *WGS72_Hgt);

void Geodetic_Shift_WGS72_To_WGS84( const double WGS72_Lat,
                                    const double WGS72_Lon,
                                    const double WGS72_Hgt,
                                    double *WGS84_Lat,
                                    double *WGS84_Lon,
                                    double *WGS84_Hgt);

/* Index into 3 Param table function declaration */
long Datum_3Param_Index( const char *Code,
                         long *Index );

void Assign_Datum_Row(Datum_Row *destination, Datum_Row *source)
{ /* Begin Assign_Datum_Row */
  /*
   * The function Assign_Datum_Row copies the data from source into destination.
   *
   * destination      : Pointer to the destination datum container.   (output)
   * source           : Pointer to the source datum container.        (input)
   */
  long i = 0;

  destination->Type = source->Type;
  strcpy(destination->Code, source->Code);
  strcpy(destination->Name, source->Name);
  strcpy(destination->Ellipsoid_Code, source->Ellipsoid_Code);

  for (i = 0; i < 7; i++)
  {
    destination->Parameters[i] = source->Parameters[i];
  }

  destination->Sigma_X = source->Sigma_X;
  destination->Sigma_Y = source->Sigma_Y;
  destination->Sigma_Z = source->Sigma_Z;

  destination->User_Defined = source->User_Defined;

} /* End Assign_Datum_Row */



long Initialize_Datums(void)
{ /* Begin Initialize_Datums */
/*
 * The function Initialize_Datums creates the datum table from two external
 * files.  If an error occurs, the initialization stops and an error code is
 * returned.  This function must be called before any of the other functions
 * in this component.
 */
  long index = 0, i = 0;
  char *PathName;
  char FileName7[FILENAME_LENGTH];
  FILE *fp_7param = NULL;
  FILE *fp_3param = NULL;
  char FileName3[FILENAME_LENGTH];
  long error_code = DATUM_NO_ERROR;

  if (Datum_Initialized)
  {
    return (error_code);
  }

  /*  Check the environment for a user provided path, else current directory;   */
  /*  Build a File Name, including specified or default path:                   */

  PathName = getenv( "DATUM_DATA" );
  if (PathName != NULL)
  {
    strcpy( FileName7, PathName );
    strcat( FileName7, "/" );
  }
  else
  {
    strcpy( FileName7, "./" );
  }
  strcat( FileName7, "7_param.dat" );

  /*  Open the File READONLY, or Return Error Condition:                        */

  if (( fp_7param = fopen( FileName7, "r" ) ) == NULL)
  {
    return ( DATUM_7PARAM_FILE_OPEN_ERROR);
  }

  if (PathName != NULL)
  {
    strcpy( FileName3, PathName );
    strcat( FileName3, "/" );
  }
  else
  {
    strcpy( FileName3, "./" );
  }
  strcat( FileName3, "3_param.dat" );

  /*  Open the File READONLY, or Return Error Condition:                        */

  if (( fp_3param = fopen( FileName3, "r" ) ) == NULL)
  {
    return ( DATUM_3PARAM_FILE_OPEN_ERROR);
  }

  if (!error_code)
  {
    while ((!feof(fp_7param)) && (!error_code))
    {
      if (index < MAX_7PARAM)
      { /* build 7-parameter datum table entries */
        if (fscanf(fp_7param, "%s ", Datum_Table_7Param[index].Code) <= 0)
          error_code |= DATUM_7PARAM_FILE_PARSING_ERROR;
        if (fscanf(fp_7param, "\"%32[^\"]\"", Datum_Table_7Param[index].Name) <= 0)
          Datum_Table_7Param[index].Name[0] = '\0';
        if (fscanf(fp_7param, " %s %lf %lf %lf %lf %lf %lf %lf ", 
                   Datum_Table_7Param[index].Ellipsoid_Code,
                   &(Datum_Table_7Param[index].Parameters[0]),
                   &(Datum_Table_7Param[index].Parameters[1]),
                   &(Datum_Table_7Param[index].Parameters[2]),
                   &(Datum_Table_7Param[index].Parameters[3]),
                   &(Datum_Table_7Param[index].Parameters[4]),
                   &(Datum_Table_7Param[index].Parameters[5]),
                   &(Datum_Table_7Param[index].Parameters[6])) <= 0)
        {
          error_code |= DATUM_7PARAM_FILE_PARSING_ERROR;
        }
        else
        { /* convert from degrees to radians */
          Datum_Table_7Param[index].Type = Seven_Param_Datum;
          Datum_Table_7Param[index].Parameters[3] /= SECONDS_PER_RADIAN;
          Datum_Table_7Param[index].Parameters[4] /= SECONDS_PER_RADIAN;
          Datum_Table_7Param[index].Parameters[5] /= SECONDS_PER_RADIAN;
          Datum_Table_7Param[index].Sigma_X = 0.0;
          Datum_Table_7Param[index].Sigma_Y = 0.0;
          Datum_Table_7Param[index].Sigma_Z = 0.0;
          Datum_Table_7Param[index].South_latitude = -PI / 2.0;
          Datum_Table_7Param[index].North_latitude = +PI / 2.0;
          Datum_Table_7Param[index].West_longitude = -PI;
          Datum_Table_7Param[index].East_longitude = +PI;

        }
        index++;
      }
      else