{
        error_code |= DATUM_7PARAM_OVERFLOW_ERROR;
      }
    }
    fclose(fp_7param);
    Datum_7Param_Count = index;
    index = 0;

    while ((!feof(fp_3param)) && (!error_code))
    {
      if (index < MAX_3PARAM)
      { /* build 3-parameter datum table entries */
        if (fscanf(fp_3param, "%s ", Datum_Table_3Param[index].Code) <= 0)
          error_code |= DATUM_3PARAM_FILE_PARSING_ERROR;
		else
		{
		  if (Datum_Table_3Param[index].Code[0] == '*')
		  {
			long i;
			Datum_Table_3Param[index].User_Defined = TRUE;
		    for (i = 0; i < DATUM_CODE_LENGTH; i++)
			  Datum_Table_3Param[index].Code[i] = Datum_Table_3Param[index].Code[i+1];
		  }
		  else
		    Datum_Table_3Param[index].User_Defined = FALSE;
		}
        if (fscanf(fp_3param, "\"%32[^\"]\"", Datum_Table_3Param[index].Name) <= 0)
          Datum_Table_3Param[index].Name[0] = '\0';
        if (fscanf(fp_3param, " %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf ",
                   Datum_Table_3Param[index].Ellipsoid_Code,
                   &(Datum_Table_3Param[index].Parameters[0]),
                   &(Datum_Table_3Param[index].Sigma_X),
                   &(Datum_Table_3Param[index].Parameters[1]),
                   &(Datum_Table_3Param[index].Sigma_Y),
                   &(Datum_Table_3Param[index].Parameters[2]),
                   &(Datum_Table_3Param[index].Sigma_Z),
                   &(Datum_Table_3Param[index].South_latitude),
                   &(Datum_Table_3Param[index].North_latitude),
                   &(Datum_Table_3Param[index].West_longitude),
                   &(Datum_Table_3Param[index].East_longitude)) <= 0)
        {
          error_code |= DATUM_3PARAM_FILE_PARSING_ERROR;
        }
        else
        {
          Datum_Table_3Param[index].Parameters[3] = 0.0;
          Datum_Table_3Param[index].Parameters[4] = 0.0;
          Datum_Table_3Param[index].Parameters[5] = 0.0;
          Datum_Table_3Param[index].Parameters[6] = 1.0;
          Datum_Table_3Param[index].South_latitude *= (PI / 180.0);
          Datum_Table_3Param[index].North_latitude *= (PI / 180.0);
          Datum_Table_3Param[index].West_longitude *= (PI / 180.0);
          Datum_Table_3Param[index].East_longitude *= (PI / 180.0);
          Datum_Table_3Param[index].Type = Three_Param_Datum;
        }
        index++;
      }
      else
      {
        error_code |= DATUM_3PARAM_OVERFLOW_ERROR;
      }
    }
    fclose(fp_3param);
    Datum_3Param_Count = index;

    /* Initialize array of pointers to datums */
    if (!error_code)
    {
      /* set total number of datums available */
      Number_of_Datums = Datum_3Param_Count + Datum_7Param_Count + 2;

      /* build WGS84 and WGS72 datum table entries */
      WGS84.Type = WGS84_Datum;
      strcpy(WGS84.Name,"World Geodetic System 1984");
      strcpy(WGS84.Code,"WGE");
      strcpy(WGS84.Ellipsoid_Code,"WE");
      WGS72.Type = WGS72_Datum;
      strcpy(WGS72.Name,"World Geodetic System 1972");
      strcpy(WGS72.Code,"WGC");
      strcpy(WGS72.Ellipsoid_Code,"WD");
      for (i=0; i<6; i++)
      {
        WGS84.Parameters[i] = 0.0;
        WGS72.Parameters[i] = 0.0;
      }
      WGS84.Parameters[6] = 1.0;
      WGS72.Parameters[6] = 1.0;

      WGS84.Sigma_X = 0.0;
      WGS84.Sigma_Y = 0.0;
      WGS84.Sigma_Z = 0.0;
      WGS84.South_latitude = -PI / 2.0;
      WGS84.North_latitude = +PI / 2.0;
      WGS84.West_longitude = -PI;
      WGS84.East_longitude = +PI;

      WGS72.Sigma_X = 0.0;
      WGS72.Sigma_Y = 0.0;
      WGS72.Sigma_Z = 0.0;
      WGS72.South_latitude = -PI / 2.0;
      WGS72.North_latitude = +PI / 2.0;
      WGS72.West_longitude = -PI;
      WGS72.East_longitude = +PI;

      Datum_WGS84_Index = 1;
      Datum_Table[Datum_WGS84_Index - 1] = &WGS84;
      Datum_WGS72_Index = 2;
      Datum_Table[Datum_WGS72_Index - 1] = &WGS72;
      index = 2;
      for (i = 0; i < Datum_7Param_Count; i++)
      {
        Datum_Table[index++] = &(Datum_Table_7Param[i]);
      }
      for (i = 0; i < Datum_3Param_Count; i++)
      {
        Datum_Table[index++] = &(Datum_Table_3Param[i]);
      }

      if (error_code)
      {
        error_code |= DATUM_ELLIPSE_ERROR;
        Datum_Initialized = 0;
        Number_of_Datums = 0;
      }
      else
        Datum_Initialized = 1;      
    }
  }
  return (error_code);
} /* End Initialize_Datums */


long Create_Datum ( const char *Code,
                    const char *Name,
                    const char *Ellipsoid_Code,
                    double Delta_X,
                    double Delta_Y,
                    double Delta_Z,
                    double Sigma_X,
                    double Sigma_Y,
                    double Sigma_Z,
                    double South_latitude,
                    double North_latitude,
                    double West_longitude,
                    double East_longitude)
{ /* Begin Create_Datum */
/*
 *   Code           : 5-letter new datum code.                      (input)
 *   Name           : Name of the new datum                         (input)
 *   Ellipsoid_Code : 2-letter code for the associated ellipsoid    (input)
 *   Delta_X        : X translation to WGS84 in meters              (input)
 *   Delta_Y        : Y translation to WGS84 in meters              (input)
 *   Delta_Z        : Z translation to WGS84 in meters              (input)
 *   Sigma_X        : Standard error in X in meters                 (input)
 *   Sigma_Y        : Standard error in Y in meters                 (input)
 *   Sigma_Z        : Standard error in Z in meters                 (input)
 *   South_latitude : Southern edge of validity rectangle in radians(input)
 *   North_latitude : Northern edge of validity rectangle in radians(input)
 *   West_longitude : Western edge of validity rectangle in radians (input)
 *   East_longitude : Eastern edge of validity rectangle in radians (input)
 *
 * The function Create_Datum creates a new local (3-parameter) datum with the 
 * specified code, name, and axes.  If the datum table has not been initialized,
 * the specified code is already in use, or a new version of the 3-param.dat 
 * file cannot be created, an error code is returned, otherwise DATUM_NO_ERROR 
 * is returned.  Note that the indexes of all datums in the datum table may be 
 * changed by this function.
 */
  char datum_Code[DATUM_CODE_LENGTH];
  long error_code = DATUM_NO_ERROR;
  long index = 0;
  long ellipsoid_index = 0;
  long code_length = 0;
  char datum_name[DATUM_NAME_LENGTH+2];
  char *PathName;
  char FileName[FILENAME_LENGTH];
  FILE *fp_3param = NULL;

  if (!Datum_Initialized)
    error_code |= DATUM_NOT_INITIALIZED_ERROR;
  if (!(Datum_3Param_Count < MAX_3PARAM))
    error_code |= DATUM_3PARAM_OVERFLOW_ERROR;
  if (!(((Sigma_X > 0.0) || (Sigma_X == -1.0)) &&
        ((Sigma_Y > 0.0) || (Sigma_Y == -1.0)) &&
        ((Sigma_Z > 0.0) || (Sigma_Z == -1.0))))
    error_code |= DATUM_SIGMA_ERROR;
  if ((South_latitude < MIN_LAT) || (South_latitude > MAX_LAT))
    error_code |= DATUM_LAT_ERROR;
  if ((West_longitude < MIN_LON) || (West_longitude > MAX_LON))
    error_code |= DATUM_LON_ERROR;
  if ((North_latitude < MIN_LAT) || (North_latitude > MAX_LAT))
    error_code |= DATUM_LAT_ERROR;
  if ((East_longitude < MIN_LON) || (East_longitude > MAX_LON))
    error_code |= DATUM_LON_ERROR;
  if ((South_latitude >= North_latitude) || (West_longitude >= East_longitude))
    error_code |= DATUM_DOMAIN_ERROR;
  code_length = strlen(Code);
  if ((!Datum_Index(Code, &index)) || (code_length > (DATUM_CODE_LENGTH-1)))
    error_code |= DATUM_INVALID_CODE_ERROR;
  if (Ellipsoid_Index(Ellipsoid_Code, &ellipsoid_index))
    error_code |= DATUM_ELLIPSE_ERROR;
  if (!error_code)
  {
    long i;
    strcpy(datum_Code,Code);
    /* Convert code to upper case */
    for (i = 0; i < code_length; i++)
      datum_Code[i] = (char)toupper(datum_Code[i]);
    index = Datum_3Param_Count;
    strcpy(Datum_Table_3Param[index].Code, datum_Code);
    strcpy(Datum_Table_3Param[index].Name, Name);
    strcpy(Datum_Table_3Param[index].Ellipsoid_Code, Ellipsoid_Code);
    Datum_Table_3Param[index].Parameters[0] = Delta_X;
    Datum_Table_3Param[index].Parameters[1] = Delta_Y;
    Datum_Table_3Param[index].Parameters[2] = Delta_Z;
    Datum_Table_3Param[index].Parameters[3] = 0.0;
    Datum_Table_3Param[index].Parameters[4] = 0.0;
    Datum_Table_3Param[index].Parameters[5] = 0.0;
    Datum_Table_3Param[index].Parameters[6] = 1.0;
    Datum_Table_3Param[index].Sigma_X = Sigma_X;
    Datum_Table_3Param[index].Sigma_Y = Sigma_Y;
    Datum_Table_3Param[index].Sigma_Z = Sigma_Z;
    Datum_Table_3Param[index].South_latitude = South_latitude;
    Datum_Table_3Param[index].North_latitude = North_latitude;
    Datum_Table_3Param[index].West_longitude = West_longitude;
    Datum_Table_3Param[index].East_longitude = East_longitude;
    Datum_Table_3Param[index].Type = Three_Param_Datum;
    Datum_Table_3Param[index].User_Defined = TRUE;
    Datum_Table[Number_of_Datums] = &(Datum_Table_3Param[index]);
    Datum_3Param_Count++;
    Number_of_Datums++;

    /*output updated 3-parameter datum table*/
    PathName = getenv( "DATUM_DATA" );
    if (PathName != NULL)
    {
      strcpy( FileName, PathName );
      strcat( FileName, "/" );
    }
    else
    {
      strcpy( FileName, "./" );
    }
    strcat( FileName, "3_param.dat" );

    if ((fp_3param = fopen(FileName, "w")) == NULL)
    { /* fatal error */
      return DATUM_3PARAM_FILE_OPEN_ERROR;
    }

    /* write file */
    index = 0;
    while (index < Datum_3Param_Count)
    {
      strcpy( datum_name, "\"" );
      strcat( datum_name, Datum_Table_3Param[index].Name);
      strcat( datum_name, "\"" );
      if (Datum_Table_3Param[index].User_Defined)
        fprintf(fp_3param, "*%-6s %-33s%-2s %4.0f %4.0f %4.0f %4.0f %5.0f %4.0f %4.0f %4.0f %4.0f %4.0f \n",
                Datum_Table_3Param[index].Code,
                datum_name,
                Datum_Table_3Param[index].Ellipsoid_Code,
                Datum_Table_3Param[index].Parameters[0],
                Datum_Table_3Param[index].Sigma_X, 
                Datum_Table_3Param[index].Parameters[1],
                Datum_Table_3Param[index].Sigma_Y, 
                Datum_Table_3Param[index].Parameters[2],
                Datum_Table_3Param[index].Sigma_Z,
                (Datum_Table_3Param[index].South_latitude * 180.0 / PI),
                (Datum_Table_3Param[index].North_latitude * 180.0 / PI),
                (Datum_Table_3Param[index].West_longitude * 180.0 / PI),
                (Datum_Table_3Param[index].East_longitude * 180.0 / PI));
      else
        fprintf(fp_3param, "%-6s  %-33s%-2s %4.0f %4.0f %4.0f %4.0f %5.0f %4.0f %4.0f %4.0f %4.0f %4.0f \n",
                Datum_Table_3Param[index].Code,
                datum_name,
                Datum_Table_3Param[index].Ellipsoid_Code,
                Datum_Table_3Param[index].Parameters[0],
                Datum_Table_3Param[index].Sigma_X, 
                Datum_Table_3Param[index].Parameters[1],
                Datum_Table_3Param[index].Sigma_Y, 
                Datum_Table_3Param[index].Parameters[2],
                Datum_Table_3Param[index].Sigma_Z,
                (Datum_Table_3Param[index].South_latitude * 180.0 / PI),
                (Datum_Table_3Param[index].North_latitude * 180.0 / PI),
                (Datum_Table_3Param[index].West_longitude * 180.0 / PI),
                (Datum_Table_3Param[index].East_longitude * 180.0 / PI));

      index++;
    }
    fclose(fp_3param);
  }
  return (error_code);
} /* End Create_Datum */

long Delete_Datum (const char *Code)

{ /* Begin Delete_Datum */
/*
 *   Code           : 5-letter datum code.                      (input)
 *
 * The function Delete_Datum deletes a local (3-parameter) datum with the 
 * specified code.  If the datum table has not been initialized or a new 
 * version of the 3-param.dat file cannot be created, an error code is returned,  
 * otherwise DATUM_NO_ERROR is returned.  Note that the indexes of all datums 
 * in the datum table may be changed by this function.
 */

  char *PathName;
  char FileName[FILENAME_LENGTH];
  char datum_name[DATUM_NAME_LENGTH+2];
  FILE *fp_3param = NULL;
  long index = 0;
  long error_code = DATUM_NO_ERROR;

  if (!Datum_Initialized)
    error_code |= DATUM_NOT_INITIALIZED_ERROR;
  else
  {
    if (!Datum_3Param_Index(Code, &index))
    {
      if (!Datum_Table_3Param[index-1].User_Defined)
        error_code |= DATUM_NOT_USERDEF_ERROR;
    }
    else
      error_code |= DATUM_NOT_USERDEF_ERROR;
  }
  if (!error_code)
  {
    long i = 0;
    long count = 0;

    for (i = index-1; i < Datum_3Param_Count-1; i++)
      Datum_Table_3Param[i] = Datum_Table_3Param[i+1];

    if (Datum_3Param_Count != MAX_3PARAM)
      Datum_Table_3Param[i] = Datum_Table_3Param[i+1];
    else
    {
      Datum_Table_3Param[i].Type = ' ';
      strcpy(Datum_Table_3Param[i].Code, "");
      strcpy(Datum_Table_3Param[i].Name, "");
      strcpy(Datum_Table_3Param[i].Ellipsoid_Code, "");
      Datum_Table_3Param[i].Parameters[0] = 0;
      Datum_Table_3Param[i].Parameters[1] = 0;
      Datum_Table_3Param[i].Parameters[2] = 0;
      Datum_Table_3Param[i].Parameters[3] = 0;
      Datum_Table_3Param[i].Parameters[4] = 0;
      Datum_Table_3Param[i].Parameters[5] = 0;
      Datum_Table_3Param[i].Parameters[6] = 0;
      Datum_Table_3Param[i].Sigma_X = 0;
      Datum_Table_3Param[i].Sigma_Y = 0;
      Datum_Table_3Param[i].Sigma_Z = 0;
      Datum_Table_3Param[i].West_longitude = 0;
      Datum_Table_3Param[i].East_longitude = 0;
      Datum_Table_3Param[i].South_latitude = 0;
      Datum_Table_3Param[i].North_latitude = 0;
      Datum_Table_3Param[i].User_Defined = ' ';