#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

/*
  BOR_CORD VERSION 2.2
  
  This program (developed for BOREAS) performs the geographic coordinate 
  conversions over the BOREAS study area of 51 N --> 60 N latitude and
  111 W --> 93 W longitude.
      
  The BOREAS grid is based on the ellipsoidal version of the Albers Equal-Area
  Conic (AEAC) projection as defined within the North American Datum 1983
  (NAD83).  The origin of the grid is at 111 W, 51 N and the standard
  parallels are set to 52.5 N and 58.5 N as prescribed in 'Map Projections -
  A Working Manual', USGS Professional Paper 1395, John P. Snyder, 1987.
  The projection equations used in this program were all taken from
  this manual.
  
  The software will take input coordinates as Universal Transverse Mercator
  (UTM) northing and easting; decimal degrees of latitude and longitude;
  or BOREAS (x,y) grid.
  
  Output coordinates for the BOREAS grid are always given under the NAD83
  datum.  Output coordinates for UTM and lat, lon coordinates are given under
  both NAD27 and NAD83 datums.

  The software handles conversion between NAD27 and NAD83 by using the
  function F to convert lat, lon coordinates between the two datums.
  The c_gridint function accesses datum shift information contained in the
  file datmshft.dat.  This file must be locatable for the software to work as
  designed.
  
  The BOREAS project thanks the Geodetic Survey of Canada for the FORTRAN-77
  source code and datum shift information files that made this possible.
===============================================================================
  The original development of this program was the result of collaboration
  between
              Jeffrey Newcomer, Scott Goetz, and Joe Lewthwaite
              Hughes STX Corporation    September 1992
-------------------------------------------------------------------------------
   The current version 2.0 is the result of efforts by:
        Jeffrey Newcomer retrofitted the version 1.3 code to procedurally
        handle the change over between datums.
        
        Fred Irani adapted the software in the function c_gridint that handles
        conversion between NAD27 and NAD83 datums from the National
        Transformation Program.
        
 -------------------------------------------------------------------------------
 Update: 10/94 F. Irani, HSTX: Added precision to longitude and latitude output.
 Update: 12/94 F. Irani, HSTX: X,Y Output now fits input format.
===============================================================================
*/

/*
  Make some definitions for machine portability.
*/
#ifndef TRUE
#define TRUE 1
#endif
#ifndef FALSE
#define FALSE 0
#endif

/*
  If you are running on a PC or VAX system, you will want to have
  the following definition activated so that reading of the file
  datmshft.dat by the c_gridint function works properly:

          #define  NRECBYTES 22


  If you are running on a Silicon Graphics or SUN work-station under 
  UNIX or on a Macintosh system, you will want to have the following 
  definition activated so that reading of the file datmshft.dat by 
  the c_gridint function works properly:

          #define  NRECBYTES 21
*/

#define  NRECBYTES 22 /* VAX, PC, etc */


/*
    North American Datum 27 (NAD27) (also Clarke 1866) parameters
*/
#define   A_27      6378206.4      /* NAD27 equatorial radius */
#define   E_27      0.0822719      /* NAD27 eccentricity */
#define   E2_27     (E_27*E_27)    /* NAD27 eccentricity squared */

/*
    North American Datum 83 (NAD83) parameters
*/
#define   A_83      6378137.0      /* NAD83 equatorial radius */
#define   E_83      0.0818192      /* NAD83 eccentricity */
#define   E2_83     (E_83*E_83)    /* NAD83 eccentricity squared */

/*
    Other needed map parameters
*/
#define  RAD2DEG  57.29577951308  /* Radian to degree conversion */
#define  DEG2RAD   0.01745329252  /* Degree to radian conversion */
#define  K0        0.9996         /* Scale factor on UTM Central Meridian */
#define  M0        0.0000         /* Distance along meridian */
                                  /* from equator to point */
/*
  Function to output software version banner
*/
void version_num (void);

/*
  Functions that control coordinate transformations
*/
void process_coords (FILE *, FILE *, FILE *, long int, char *, char *, char *, char *);
long int process_boreas (double, double, double, double,
                         FILE *, FILE *, FILE *);
long int process_utm27 (double, double, double, double,
                        FILE *, FILE *, FILE *);
long int process_utm83 (double, double, double, double,
                        FILE *, FILE *, FILE *);
long int process_latlon27 (double, double, double, double,
                           FILE *, FILE *, FILE *);
long int process_latlon83 (double, double, double, double,
                           FILE *, FILE *, FILE *);

/*
  Functions to read input coordinate information
*/
long int read_boreas (FILE *, FILE *, double *, double *);
long int read_ll (FILE *, FILE *, double *, double *);
long int read_utm (FILE *, FILE *, double *, double *, long int *);

/*
  Functions to perform basic coordinate transformations
*/
long int albers_2_boreas (double, double, double *, double *);
long int albers_2_ll (double, double, double, double, long int, double *,
                      double *, double *, double *);
long int boreas_2_albers (double, double, double *, double *);
void ll_2_albers (double, double, double, double, long int, double *,
                      double *, double *, double *);
long int ll_2_utm (double *, double *, double, double, double, long int,
                   double *, double *, long int *);
long int utm_2_ll (double *, double *, long int *, double, double, double,
                   long int, double *, double *);
              
/*
  Function to perform datum coordinate conversions
*/
long int c_gridint (FILE *, char *, long int, double *, double *,
                    double *, double *);

/*
  Function to output coordinate information
*/
long int write_coords (FILE *, double, double, double, double, double, double,
                       double, double, long int, double, double, long int);

/*
 Function to check user response for quit 
*/ 
long int check_response (char *);

/*
===============================================================================
   S T A R T   O F   M A I N   F U N C T I O N
===============================================================================
*/ 
main ()
{  

  char checkfile[100], filename[100], input_str[20];
  char title0[200], title1[200], title2[200], title3[200];
 
  FILE *datum_fp, *in_fp = NULL, *out_fp = NULL;

  long int choice, cnt, len;

  int  stat;
  
  /*
    Output banner and determine the coordinate conversions to be performed
  */
  version_num ();

  /*
    As long as the user wants to perform coordinate transformations
    keep prompting them for what transformations to perform.
  */
  for ( ; ; )
  {
    /*
      Keep presenting the menu and trying to get a good response from the user.
      Note: All files should be closed at this point.
    */
    do
    {
      puts ("\n\n  Input Coordinate Selection Menu.");
      puts ("      1) BOREAS X, Y.");
      puts ("      2) UTM  Easting, Northing");
      puts ("      3) Longitude, Latitude");
      puts ("      4) Exit.");
      printf (" Enter choice < 1 - 4 >: ");
      gets (input_str);
      len = strlen (input_str);
      if (len > 0)
      {
        cnt = sscanf (input_str, "%ld", &choice);
        if (cnt != 1)
        {
          puts ("*** Problem interpreting input ***");
          puts ("*** Please re-enter response ***");
          len = 0;
        }
      }
    } while (len == 0);

    if (choice == 4)
    {
      puts ("\n\n< Boreas_Coords Completed >");
      exit (0);
    }
	
    /*
      Do until a useable filename has been entered, <cr> or q(uit).
    */
    stat = TRUE;
    for (; ;)
    {
     /*
        Determine how the coordinate values are to be entered
        and if needed attempt to open the input file
      */
      puts ("\n\nEnter Input Filename or RETURN for keyboard input (q = quit)");
      printf (" --> ");
      gets (filename);
      if (!check_response(filename))
      {
	stat = FALSE;
	break;
     }
      else if (strlen (filename) == 0)
      {
	in_fp = (FILE *) NULL;
        break;
      }
      else if ((in_fp = fopen (filename, "r")) == NULL)
      {
	printf ("\n*** Unable to open %s for input ***\n", filename);
	continue;
      }
      break;
    }
    if (!stat) continue;

    /*
      Do forever until able to 
      determine where to send the converted values and
      attempt to open the output file if needed
    */
    stat = TRUE;
    for (; ;)
    {
      puts ("\n\nEnter Output filename or RETURN for screen output (q = quit)");
      printf (" --> ");
      gets (filename);
      if (!check_response(filename))
      {
	stat = FALSE;
	break;
      }
      else if (strlen (filename) == 0)
      {
       out_fp = (FILE *) NULL;
       break;
      }
      else if ((out_fp = fopen (filename, "w")) == NULL)
      {
	    printf ("\n*** Unable to  open %s for output ***\n", filename);
        continue;
      }
      /*
      Save valid user response for output file title.
      ----------------------------------------------- */
      sprintf(title0, "\nFILE: %s\n\n", filename);
      break;
    }
    if (!stat) continue;
    
    /*
      Open the file that contains the datum shift information
      needed by the c_gridint function for conversions between
      NAD27 and NAD83.
    */
    datum_fp = (FILE *) NULL;
    stat = TRUE;
    while (TRUE)
    {
      strcpy(filename, "datmshft.dat");
      printf ("\n\nEnter datum shift file name (Default = %s; q = quit)\n", filename);
      printf (" --> ");
      gets (filename);
      cnt = sscanf (filename, "%s", checkfile);
      if (cnt < 1) strcpy(filename, "datmshft.dat");
      else if ( !check_response(filename) )
      {
	stat = FALSE;
	break;
      }

      if ((datum_fp = fopen (filename, "r")) == NULL)
      {
	printf ("\n*** Unable to open %s for input ***\n", filename);
        continue;
      }
      else break;
    }
    if (!stat) continue;

    /*
      If an output file was specified, then set up the column titles to
      the top three records of the output file for later use 
   */
    if (out_fp != (FILE *) NULL)
    {
      title1[0] = title2[0] = title3[0] = '\0';
      sprintf(title1, 
      			"\n      * = INPUT\n\n   BOREAS GRID       NAD27 COORDINATES       ");
      strcat(title1,  "NAD83 COORDINATES      NAD27 COORDINATES");
      strcat(title1,  "           NAD83 COORDINATES     \n");

      sprintf(title2, " X (KM)    Y (KM)    LONGITUDE  LATITUDE    ");
      sprintf(title3, " --------  -------   ---------- --------    ");

      strcat(title2,  " LONGITUDE   LATITUDE   EASTING   NORTHING   ZONE");
      strcat(title3,  " ----------  --------   --------  ---------  ----");

      strcat(title2,  "   EASTING   NORTHING   ZONE\n");
      strcat(title3,  "   --------  ---------  ----\n");
    }
    
    /*
      Perform the specified coordinate conversion operations
    */
    process_coords (in_fp, out_fp, datum_fp, choice, title0, title1, title2, title3);

    /*
      Before continuing loop, close up the files
    */
    if (in_fp    != (FILE *) NULL) fclose (in_fp);
    if (out_fp   != (FILE *) NULL) fclose (out_fp);
    if (datum_fp != (FILE *) NULL) fclose (datum_fp);
  	 
  }  /* end of forever for loop */

} /* END OF MAIN */

/*=============================================================================
  Function to control processing of coordinates

  Input parameters:
      in_fp == FILE pointer to file of input coordinate values
      out_fp == FILE pointer to file for output coordinate information
      datum_fp == FILE pointer to file containing datum shift information
      choice == The number that designates the user's
                selected input coordinate system

  Output parameters:
      Returned status value.
=============================================================================*/
void process_coords (in_fp, out_fp, datum_fp, choice, title0, title1, title2, title3)
FILE *in_fp, *out_fp, *datum_fp;
long int choice;
char *title0, *title1, *title2, *title3;
{
  double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
  
  long int cnt, in_datum, stat;
  
  char input_str[80];
  
/*
    Set values for the BOREAS Grid and Albers Equal Area Conic
    (AEAC) projection parameters
*/
  aeac_phi0 = 51.00;        /* Origin latitude for AEAC projection */
                            /*  (southern most) */
  aeac_lambda0 = -111.00;   /* Origin longitude for AEAC projection */
                            /*  (western most) */
  aeac_phi1 = 52.50;        /* Latitude of southern standard parallel */
  aeac_phi2 = 58.50;        /* Latitude of northern standard parallel */

  switch (choice)
  {
    case 1 :  /* BOREAS (x,y) */
    {
      puts ("\n\n\nNote that the BOREAS X and Y grid values are assumed");
      puts ("to have been calculated under the NAD83 datum");
      if (out_fp != (FILE *) NULL)
      {
        title1[20] = '*';
        stat = fprintf(out_fp, title0);
        stat = fprintf(out_fp, title1);
        stat = fprintf(out_fp, title2);
        stat = fprintf(out_fp, title3);
        if (stat < 0)
        {
          puts ("\n*** Error writing header records to output file ***");
          fclose (out_fp);
          if (in_fp != (FILE *) NULL) fclose (in_fp);
          fclose (datum_fp);
          datum_fp = in_fp = out_fp = (FILE *) NULL;
       	  return;
        }
      }
      process_boreas (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                      in_fp, datum_fp, out_fp);
      break;
    }
    case 2 :  /* UTM */
    {
      do
      {
        puts ("\n\nEnter the datum for the input UTM coordinates");
        puts ("as 27 for NAD27 or 83 for NAD83 (q = quit)");
        puts (" --> ");
        if (gets (input_str) == NULL)
        {
          puts ("*** Error reading from terminal ***");
          continue;
        }
	if ( !check_response (input_str) ) return;
    	cnt = sscanf (input_str, "%ld", &in_datum);
        if (cnt != 1)
        {
          puts ("*** Error extracting input datum from input record ***");
          continue;
        }
        if (in_datum != 27 && in_datum != 83)
            puts ("*** Datum values need to be given as 27 or 83 ***");
      } while (in_datum != 27 && in_datum !=83);
      if (in_datum == 27) 
      {
        if (out_fp != (FILE *) NULL)
        {
          title1[89] = '*';
          stat = fprintf(out_fp, title0);
          stat = fprintf(out_fp, title1);
          stat = fprintf(out_fp, title2);
          stat = fprintf(out_fp, title3);
          if (stat < 0)
          {
            puts ("\n*** Error writing header records to output file ***");
            fclose (out_fp);
            if (in_fp != (FILE *) NULL) fclose (in_fp);
            fclose (datum_fp);
            datum_fp = in_fp = out_fp = (FILE *) NULL;
       	    return;
          }
        }
      	process_utm27 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                       in_fp, datum_fp, out_fp);
      }
      else if (in_datum == 83)
      {
        if (out_fp != (FILE *) NULL)
        {
          title1[117] = '*';
          stat = fprintf(out_fp, title0);
          stat = fprintf(out_fp, title1);
          stat = fprintf(out_fp, title2);
          stat = fprintf(out_fp, title3);
          if (stat < 0)
          {
            puts ("\n*** Error writing header records to output file ***");
            fclose (out_fp);
            if (in_fp != (FILE *) NULL) fclose (in_fp);
            fclose (datum_fp);
            datum_fp = in_fp = out_fp = (FILE *) NULL;
       	    return;
          }
        }
      	process_utm83 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                       in_fp, datum_fp, out_fp);
      }
      break;
    }
    case 3 :  /* Latitude, longitude */
    {
      do
      {
        puts ("\n\nEnter the datum for the input latitude and longitude");
        puts ("coordinates as 27 for NAD27 or 83 for NAD83. (q = quit)");
        puts (" --> ");
        if (gets (input_str) == NULL)
        {
          puts ("*** Error reading from terminal ***");
          continue;
        }
        if ( !check_response (input_str) ) return;
        cnt = sscanf (input_str, "%ld", &in_datum);
        if (cnt != 1)
        {
          puts ("*** Error extracting input datum from input record ***");
          continue;
        } 
      } while (in_datum != 27 && in_datum !=83);
      if (in_datum == 27) 
      {
        if (out_fp != (FILE *) NULL)
        {
          title1[39] = '*';
          stat = fprintf(out_fp, title0);
          stat = fprintf(out_fp, title1);
          stat = fprintf(out_fp, title2);
          stat = fprintf(out_fp, title3);
          if (stat < 0)
          {
            puts ("\n*** Error writing header records to output file ***");
            fclose (out_fp);
            if (in_fp != (FILE *) NULL) fclose (in_fp);
            fclose (datum_fp);
            datum_fp = in_fp = out_fp = (FILE *) NULL;
       	    return;
          }
        }
     	process_latlon27 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                          in_fp, datum_fp, out_fp);
      }
      else if (in_datum == 83)
      {
        if (out_fp != (FILE *) NULL)
        {
          title1[63] = '*';
          stat = fprintf(out_fp, title0);
          stat = fprintf(out_fp, title1);
          stat = fprintf(out_fp, title2);
          stat = fprintf(out_fp, title3);
          if (stat < 0)
          {
            puts ("\n*** Error writing header records to output file ***");
            fclose (out_fp);
            if (in_fp != (FILE *) NULL) fclose (in_fp);
            fclose (datum_fp);
            datum_fp = in_fp = out_fp = (FILE *) NULL;
            return;
          }
        }
      	process_latlon83 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                          in_fp, datum_fp, out_fp);
      }        
      break;
    }
  }  /* end of switch */

}  /* END OF PROCESS_COORDS */

/*=============================================================================
   Function to control conversion from BOREAS grid
   locations to latitude and longitude and UTM coordinates.
   
   Note that this function sets the in_datum parameter to a constant value
   of 83 since the BOREAS grid system was based on the NAD83 datum.
   
  Input parameters:
      aeac_phi0 == origin latitude for the projection (southern most)
      aeac_phi1 == latitude of southern standard parallel
      aeac_phi2 == latitude of northern standard parallel
      aeac_lambda0 == origin longitude for the projection (western most)
      in_fp == FILE pointer to file of input coordinate values
      datum_fp == FILE pointer to file containing datum shift information
      out_fp == FILE pointer to file for output coordinate information

  Output parameters:
      Returned status value.
=============================================================================*/
long int process_boreas (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                    in_fp, datum_fp, out_fp)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
FILE *in_fp, *datum_fp, *out_fp;
{

  double lat_27, lat_83, lon_27, lon_83;
  double easting_27, easting_83, northing_27, northing_83;
  double x, xgrid, y, ygrid;

  long int cnt, zone_27, zone_83;

  char datum_str[8], input_str[80];

  /*
    Set the input coordinate datum to NAD83 since the
    BOREAS grid is based on the NAD83 datum.
  */
  strcpy (datum_str, "NAD83");
  
  /*
    While there are coordinates to process, handle them accordingly
  */
  while ( read_boreas (in_fp, out_fp, &xgrid, &ygrid) )
  {
    /*
      Convert BOREAS x and y components into x and y AEAC meter
      values under NAD83
    */
    boreas_2_albers (xgrid, ygrid, &x, &y);

    /*
      Convert AEAC meters to latitude and longitude under NAD83
    */
    albers_2_ll (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, (long int) 1, &x, &y,
                 &lat_83, &lon_83);

    /*
      Convert NAD83 latitude and longitude to NAD83 UTM northing and easting
    */
    ll_2_utm (&lat_83, &lon_83, (double) A_83, (double) E_83, (double) E2_83,
              (long int) 1, &northing_83, &easting_83, &zone_83);
    
    /*
      Convert NAD83 latitude and longitude to NAD27 latitude and longitude    
    */
    if ( !c_gridint (datum_fp, datum_str, (long int) 1, &lat_83, &lon_83,
                     &lat_27, &lon_27) )  return (0);
    
    /*
      Convert NAD27 latitude and longitude to NAD27 UTM northing and easting
    */
    ll_2_utm (&lat_27, &lon_27, (double) A_27, (double) E_27, (double) E2_27,
              (long int) 1, &northing_27, &easting_27, &zone_27);
    
    /*
      Output the information as indicated
    */
    if (!write_coords (out_fp, xgrid, ygrid, lat_27, lon_27, lat_83, lon_83,
                       northing_27, easting_27, zone_27, northing_83,
                       easting_83, zone_83) )
    {
      puts ("*** Error writing coordinates to output ***");
      return (0);
    }
  } /* end of while loop */
  
  return (1);

}  /* END OF PROCESS_BOREAS */

/*=============================================================================
   Function to control conversion from Universal Transverse
   Mercator (UTM) under NAD27 to latitude, longitude and BOREAS grid

  Input parameters:
      aeac_phi0 == origin latitude for AEAC projection (southern most)
      aeac_phi1 == latitude of southern standard parallel for AEAC projection
      aeac_phi2 == latitude of northern standard parallel for AEAC projection
      aeac_lambda0 == origin longitude for the AEAC projection (western most)
      in_fp == FILE pointer to file of input coordinate values
      datum_fp == FILE pointer to file containing datum shift information
      out_fp == FILE pointer to file for output coordinate information

  Output parameters:
      Returned status value.
=============================================================================*/
long int process_utm27 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                   in_fp, datum_fp, out_fp)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
FILE *in_fp, *datum_fp, *out_fp;
{
  double lat_27, lat_83, lon_27, lon_83;
  double easting_27, easting_83, northing_27, northing_83;
  double temp_north, temp_east, x, xgrid, y, ygrid;

  long int cnt, zone_27, zone_83;

  char datum_str[8], input_str[80];

  /*
    Set the input coordinate datum to NAD27.
  */
  strcpy (datum_str, "NAD27");
  
  /*
    While there are coordinates to process, handle them accordingly
  */
  while ( read_utm (in_fp, out_fp, &northing_27, &easting_27, &zone_27) )
  {
    /*
      Convert UTM northing and easting to latitude and longitude
      under NAD27 datum
    */
    utm_2_ll ( &northing_27, &easting_27, &zone_27, (double) A_27,
               (double) E_27, (double) E2_27, (long int) 1, &lat_27, &lon_27 );

    /*
      Convert NAD27 latitude and longitude to NAD83 latitude and longitude
    */
    if ( !c_gridint (datum_fp, datum_str, (long int) 1, &lat_27, &lon_27,
                     &lat_83, &lon_83) ) return (0);
        
    /*
      Convert NAD83 latitude and longitude to NAD83 UTM northing and easting
    */
    ll_2_utm (&lat_83, &lon_83, (double) A_83, (double) E_83, (double) E2_83,
              (long int) 1, &northing_83, &easting_83, &zone_83);
    
    /*
      Convert NAD83 latitude and longitude values to AEAC meter values
    */
    ll_2_albers (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, (long int) 1,
                 &lat_83, &lon_83, &x, &y);

    /*
      Convert AEAC meters to BOREAS grid coordinates
    */
    albers_2_boreas (x, y, &xgrid, &ygrid);

    /*
      Output the information as indicated
    */
    if (!write_coords (out_fp, xgrid, ygrid, lat_27, lon_27, lat_83, lon_83,
                       northing_27, easting_27, zone_27, northing_83,
                       easting_83, zone_83) )
    {
      puts ("*** Error writing coordinates to output ***");
      return (0);
    }
  } /* end of while loop */
  
  return (1);

}  /* END OF PROCESS_UTM27 */

/*=============================================================================
   Function to control conversion from Universal Transverse
   Mercator (UTM) under NAD83 to latitude, longitude  and BOREAS grid

  Input parameters:
      aeac_phi0 == origin latitude for AEAC projection (southern most)
      aeac_phi1 == latitude of southern standard parallel for AEAC projection
      aeac_phi2 == latitude of northern standard parallel for AEAC projection
      aeac_lambda0 == origin longitude for the AEAC projection (western most)
      in_fp == FILE pointer to file of input coordinate values
      datum_fp == FILE pointer to file containing datum shift information
      out_fp == FILE pointer to file for output coordinate information

  Output parameters:
      Returned status value.
=============================================================================*/
long int process_utm83 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                   in_fp, datum_fp, out_fp)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
FILE *in_fp, *datum_fp, *out_fp;
{
  double lat_27, lat_83, lon_27, lon_83;
  double easting_27, easting_83, northing_27, northing_83;
  double temp_north, temp_east, x, xgrid, y, ygrid;

  long int cnt, zone_27, zone_83;

  char datum_str[8], input_str[80];

  /*
    Set the input coordinate datum to NAD83.
  */
  strcpy (datum_str, "NAD83");
  
  /*
    While there are coordinates to process, handle them accordingly
  */
  while ( read_utm (in_fp, out_fp, &northing_83, &easting_83, &zone_83) )
  {
    /*
      Convert UTM northing and easting to latitude and longitude
      under NAD83 datum
    */
    utm_2_ll ( &northing_83, &easting_83, &zone_83, (double) A_83,
               (double) E_83, (double) E2_83, (long int) 1, &lat_83, &lon_83 );

    /*
      Convert NAD83 latitude and longitude to NAD27 latitude and longitude
    */
    if ( !c_gridint (datum_fp, datum_str, (long int) 1, &lat_83, &lon_83,
                     &lat_27, &lon_27) )  return (0);
        
    /*
      Convert NAD27 latitude and longitude to NAD27 UTM northing and easting
    */
    ll_2_utm (&lat_27, &lon_27, (double) A_27, (double) E_27, (double) E2_27,
              (long int) 1, &northing_27, &easting_27, &zone_27);
    
    /*
      Convert NAD83 latitude and longitude values to AEAC meter values
    */
    ll_2_albers (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, (long int) 1,
                 &lat_83, &lon_83, &x, &y);

    /*
      Convert AEAC meters to BOREAS grid components
    */
    albers_2_boreas (x, y, &xgrid, &ygrid);

    /*
      Output the information as indicated
    */
    if (!write_coords (out_fp, xgrid, ygrid, lat_27, lon_27, lat_83, lon_83,
                       northing_27, easting_27, zone_27, northing_83,
                       easting_83, zone_83) )
    {
      puts ("*** Error writing coordinates to output ***");
      return (0);
    }
  } /* end of while loop */
  
  return (1);

}  /* END OF PROCESS_UTM83 */

/*=============================================================================
   Function to control conversion from latitude and
   longitude under NAD27 to BOREAS grid and UTM coordinates.

  Input parameters:
      aeac_phi0 == origin latitude for AEAC projection (southern most)
      aeac_phi1 == latitude of southern standard parallel for AEAC projection
      aeac_phi2 == latitude of northern standard parallel for AEAC projection
      aeac_lambda0 == origin longitude for the AEAC projection (western most)
      in_fp == FILE pointer to file of input coordinate values
      datum_fp == FILE pointer to file containing datum shift information
      out_fp == FILE pointer to file for output coordinate information

  Output parameters:
      Returned status value.
=============================================================================*/
long int process_latlon27 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                      in_fp, datum_fp, out_fp)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
FILE *in_fp, *datum_fp, *out_fp;
{
  double lat_27, lat_83, lon_27, lon_83;
  double easting_27, easting_83, northing_27, northing_83;
  double x, xgrid, y, ygrid;

  long int cnt, zone_27, zone_83;

  char datum_str[8], input_str[80];

  /*
    Set the input coordinate datum to NAD27.
  */
  strcpy (datum_str, "NAD27");
  
  /*
    While there are coordinates to process, handle them accordingly
  */
  while ( read_ll (in_fp, out_fp, &lat_27, &lon_27) )
  {
    /*
      Convert NAD27 latitude and longitude to NAD83 latitude and longitude
    */
    if ( !c_gridint (datum_fp, datum_str, (long int) 1, &lat_27, &lon_27,
                     &lat_83, &lon_83) )  return (0);
        
    /*
      Convert NAD83 latitude and longitude values to AEAC meter values
    */
    ll_2_albers (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, (long int) 1,
                 &lat_83, &lon_83, &x, &y);

    /*
      Convert AEAC meters to BOREAS grid components
    */
    albers_2_boreas (x, y, &xgrid, &ygrid);

    /*
      Convert NAD27 latitude and longitude to NAD27 UTM northing and easting
    */
    ll_2_utm (&lat_27, &lon_27, (double) A_27, (double) E_27, (double) E2_27,
              (long int) 1, &northing_27, &easting_27, &zone_27);
    
    /*
      Convert NAD83 latitude and longitude to NAD83 UTM northing and easting
    */
    ll_2_utm (&lat_83, &lon_83, (double) A_83, (double) E_83, (double) E2_83,
              (long int) 1, &northing_83, &easting_83, &zone_83);
    
    /*
      Output the information as indicated
    */
    if (!write_coords (out_fp, xgrid, ygrid, lat_27, lon_27, lat_83, lon_83,
                       northing_27, easting_27, zone_27, northing_83,
                       easting_83, zone_83) )
    {
      puts ("*** Error writing coordinates to output ***");
      return (0);
    }
  } /* end of while loop */
  
  return (1);

}  /* END OF PROCESS_LATLON27 */

/*=============================================================================
   Function to control conversion from and
   longitude under NAD83 to BOREAS grid and UTM coordinates.

  Input parameters:
      aeac_phi0 == origin latitude for AEAC projection (southern most)
      aeac_phi1 == latitude of southern standard parallel for AEAC projection
      aeac_phi2 == latitude of northern standard parallel for AEAC projection
      aeac_lambda0 == origin longitude for the AEAC projection (western most)
      in_fp == FILE pointer to file of input coordinate values
      datum_fp == FILE pointer to file containing datum shift information
      out_fp == FILE pointer to file for output coordinate information

  Output parameters:
      Returned status value.
=============================================================================*/
long int process_latlon83 (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0,
                      in_fp, datum_fp, out_fp)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
FILE *in_fp, *datum_fp, *out_fp;
{
  double lat_27, lat_83, lon_27, lon_83;
  double easting_27, easting_83, northing_27, northing_83;
  double x, xgrid, y, ygrid;

  long int cnt, zone_27, zone_83;

  char datum_str[8], input_str[80];

  /*
    Set the input coordinate datum to NAD83.
  */
  strcpy (datum_str, "NAD83");
  
  /*
    While there are coordinates to process, handle them accordingly
  */
  while ( read_ll (in_fp, out_fp, &lat_83, &lon_83) )
  {
    /*
      Convert NAD83 latitude and longitude to NAD27 latitude and longitude
    */
    if ( !c_gridint (datum_fp, datum_str, (long int) 1, &lat_83, &lon_83,
                     &lat_27, &lon_27) )  return (0);
        
    /*
      Convert NAD83 latitude and longitude values to AEAC meter values
    */
    ll_2_albers (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, (long int) 1,
                 &lat_83, &lon_83, &x, &y);

    /*
      Convert AEAC meters to BOREAS grid components
    */
    albers_2_boreas (x, y, &xgrid, &ygrid);

    /*
      Convert NAD27 latitude and longitude to NAD27 UTM northing and easting
    */
    ll_2_utm (&lat_27, &lon_27, (double) A_27, (double) E_27, (double) E2_27,
              (long int) 1, &northing_27, &easting_27, &zone_27);
    
    /*
      Convert NAD83 latitude and longitude to NAD83 UTM northing and easting
    */
    ll_2_utm (&lat_83, &
lon_83, (double) A_83, (double) E_83, (double) E2_83,
              (long int) 1, &northing_83, &easting_83, &zone_83);
    
    /*
      Output the information as indicated
    */
    if (!write_coords (out_fp, xgrid, ygrid, lat_27, lon_27, lat_83, lon_83,
                       northing_27, easting_27, zone_27, northing_83,
                       easting_83, zone_83) )
    {
      puts ("*** Error writing coordinates to output ***");
      return (0);
    }
  } /* end of while loop */
  
  return (1);

}  /* END OF PROCESS_LATLON83 */

/*=============================================================================
  Function to convert Albers Equal-Area Conic (x,y) coordinates
  to latitude and longitude coordinates in NAD83 datum.
  
  All input latitude and longitide coordinates are assumed to be
  in decimal degree units.

  Input parameters:
      aeac_phi0 == origin latitude for the projection (southern most)
      aeac_phi1 == latitude of southern standard parallel
      aeac_phi2 == latitude of northern standard parallel
      aeac_lambda0 == origin longitude for the project (western most)
      ncoords == the number of latitude and longitude coordinate pairs
                 contained in the lat and lon arrays
      x == array of AEAC x values in meters
      y == array of AEAC y values in meters

  Output parameters:
      lat == array of latitude coordinates in degrees
      lon == array of longitude coordinates in degrees
=============================================================================*/
long int albers_2_ll (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, ncoords,
                 x, y, lat, lon)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
long int ncoords;
double *x, *y, *lat, *lon;
{

  double c, lambda, m, m_1, m_2, n, phi, q, q_0, q_1, q_2;
  double rho, rho_0, temp1, temp2, temp3, temp4, temp5, temp6, temp7, theta;

  long int i, j;

  /*
    Calculate the needed constants
  */
  temp1 = sin(aeac_phi0*DEG2RAD);
  temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
  q_0 = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));
  aeac_lambda0 = aeac_lambda0 * DEG2RAD;

  temp1 = sin(aeac_phi1*DEG2RAD);
  m_1 = cos(aeac_phi1*DEG2RAD)/sqrt(1.0-E2_83*temp1*temp1);
  temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
  q_1 = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));

  temp1 = sin(aeac_phi2*DEG2RAD);
  m_2 = cos(aeac_phi2*DEG2RAD)/sqrt(1.0-E2_83*temp1*temp1);
  temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
  q_2 = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));

  n = (m_1*m_1 - m_2*m_2)/(q_2 - q_1);
  c = m_1*m_1 + n*q_1;
  rho_0 = A_83*sqrt(c - n*q_0)/n;

  /*
    Calculate lat and lon coordinates from each Albers x, y pair
  */
  for (i = 0; i < ncoords; i++)
  {
    temp1 = x[i];
    temp2 = rho_0 - y[i];
    rho = sqrt( temp1*temp1 + temp2*temp2 );

    theta = atan(temp1/temp2);
  
    temp1 = rho*n/A_83;
    q = (c - temp1*temp1)/n;

    /*
      Set initial value of phi (i.e., lat) and derive iterative solution
    */
    phi = asin (q/2.0);
    lat[i] = phi * RAD2DEG;
    for (j = 0; j < 6; j++)
    {
      temp1 = sin(phi);
      temp2 = E_83*temp1;
      temp3 = temp2*temp2;
      temp4 = 1.0 - temp3;
      temp4 = (temp4*temp4)/(2.0*cos(phi));
      temp5 = q/(1.0 - E2_83);
      temp6 = temp1/(1.0 - temp3);
      temp7 = log( (1.0-temp2)/(1.0+temp2) );
      temp7 = temp7/(2.0*E_83);
      phi = phi + temp4*(temp5 - temp6 + temp7);
      lat[i] = phi * RAD2DEG;
    }

    lambda = aeac_lambda0 + theta/n;
    lon[i] = lambda * RAD2DEG;
  }
  return (0);

}  /* END OF ALBERS_2_LL */

/*=============================================================================
  Function to convert from latitude and longitude under NAD83
  to Albers Equal-Area Conic (x,y) coordinates.
  
  Note that all input latitude and longitide coordinates are assumed to be
  in decimal degree units.  Also the earth radius and eccentricity 
  parameters (A_83, E_83, and E2_83) are all hardcoded here since for the
  BOREAS project, the datum is assumed to have been NAD83.

  Input parameters:
      aeac_phi0 == origin latitude for the projection (southern most)
      aeac_phi1 == latitude of southern standard parallel
      aeac_phi2 == latitude of northern standard parallel
      aeac_lambda0 == origin longitude for the project (western most)
      ncoords == the number of latitude and longitude coordinate pairs
                 contained in the lat and lon arrays
      lat == array of latitude coordinates
      lon == array of longitude coordinates

  Output parameters:
      x == array of AEAC x values in meters
      y == array of AEAC y values in meters
=============================================================================*/
void ll_2_albers (aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0, ncoords,
                 lat, lon, x, y)
double aeac_phi0, aeac_phi1, aeac_phi2, aeac_lambda0;
long int ncoords;
double *lat, *lon, *x, *y;
{
  double c, lambda, m_1, m_2, n, phi, q, q_0, q_1, q_2;
  double rho, rho_0, temp1, temp2, theta;

  long int i;

  /*
    Calculate the needed constants
  */
  temp1 = sin(aeac_phi0*DEG2RAD);
  temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
  q_0 = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));
  aeac_lambda0 = aeac_lambda0 * DEG2RAD;

  temp1 = sin(aeac_phi1*DEG2RAD);
  m_1 = cos(aeac_phi1*DEG2RAD)/sqrt(1.0-E2_83*temp1*temp1);
  temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
  q_1 = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));

  temp1 = sin(aeac_phi2*DEG2RAD);
  m_2 = cos(aeac_phi2*DEG2RAD)/sqrt(1.0-E2_83*temp1*temp1);
  temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
  q_2 = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));

  n = (m_1*m_1 - m_2*m_2)/(q_2 - q_1);
  c = m_1*m_1 + n*q_1;

  rho_0 = A_83*sqrt(c - n*q_0)/n;

  /*
    Calculate the Albers x and y coordinates for each lat, lon pair
  */
  for (i = 0; i < ncoords; i++)
  {
    phi = lat[i] * DEG2RAD;
    lambda = lon[i] * DEG2RAD;
    theta = n*(lambda - aeac_lambda0);
    temp1 = sin(phi);
    temp2 = (1.0-E_83*temp1)/(1.0+E_83*temp1);
    q = (1.0-E2_83) * (temp1/(1.0-E2_83*temp1*temp1) - log(temp2)/(2.0*E_83));
    rho = A_83*sqrt(c - n*q)/n;
    x[i] = rho * sin(theta) + 0.5;
    y[i] = rho_0 - rho*cos(theta) + 0.5;
  }
}  /* END OF LL_2_ALBERS */

/*=============================================================================
  Function to convert from latitude and longitude
  to UTM northing and easting and zone
  See pp.48-65;269-270, USGS Map Projections Working Manual, J.P.Snyder
  
  Note that input latitude and longitide coordinates are assumed to be
  in decimal degree units.

  Input parameters:
      lat == array of latitude coordinates to be transformed
      lon == array of longitude coordinates to be transformed
      A == semimajor axis of the earth ellipsoid in meters
      E == eccentricity of the assumed ellipsoid
      E2 == eccentricity squared of the assumed ellipsoid
      ncoords == the number of latitude and longitude coordinate pairs
                 contained in the lat and lon arrays

  Output parameters:
      northing == array of northing coordinates
      easting == array of easting coordinates
      zone == array of UTM zones for the coordinates
=============================================================================*/
long int ll_2_utm (lat, lon, A, E, E2, ncoords, northing, easting, zone)
double *lat, *lon, A, E, E2;
long int ncoords;
double *northing, *easting;
long int *zone;
{
  /*
    The flow is from Snyder, solving in order equations:
    8.12 - 8.15, 3.22, 8.9, 8.10
  */

  double lat_rad;
  double e_prime, lambda0;
  double a, c, m, n, t;
  double coef1, coef2, coef3, coef4;

  long int i;

  /*
    Calculate some needed values outside the loop
  */
  e_prime = E2 / (1.0 - E2);  /* Equation 8.12 */
  
  /* The four following values are from Eq 3.21 */
  coef1 = 1.0 - pow(E,2.0)/4.0 - 3.0*pow(E,4.0)/64.0 - 5.0*pow(E,6.0)/256.0;
  coef1 = A * coef1 * DEG2RAD;
  coef2 = A * (3.0*pow(E,2.0)/8.0 + 3.0*pow(E,4.0)/32.0 + 45.0*pow(E,6.0)/1024.0);
  coef3 = A * (15.0*pow(E,4.0)/256.0 + 45.0*pow(E,6.0)/1024.0);
  coef4 = A * (35.0*pow(E,6.0)/3072.0);

  /*
    Calculate the Northing and Easting for each
    pair of input latitude and longitude values
  */
  for (i = 0; i < ncoords; i++)
  {
    /* 
      Assign appropriate central meridian to UTM zone based on input longitude 
    */
    zone[i] = (long int) ( (lon[i]+177.0)/6.0 + 1.5);
    lambda0 = (double) (zone[i]-1)*6 - 177;

    /*
      Equations 8.12 - 8.15
    */
    lat_rad = lat[i] * DEG2RAD;
    n = A / sqrt(1.0 - E2 * sin(lat_rad) * sin(lat_rad));
    t = tan(lat_rad) * tan(lat_rad);
    c = e_prime * cos(lat_rad) * cos(lat_rad);
    a = cos(lat_rad) * ((lon[i] - lambda0) * DEG2RAD);

    /*
      Equation 3.22, The distance along central meridian to given latitude
    m = (111132.0894 * lat[i]) - (16216.94 * sin(2.0 * lat_rad)) 
         + (17.21 * sin(4.0 * lat_rad)) - (0.02 * sin(6.0 * lat_rad));	
    */
    m = (coef1 * lat[i]) - (coef2 * sin(2.0 * lat_rad)) 
         + (coef3 * sin(4.0 * lat_rad)) - (coef4 * sin(6.0 * lat_rad));	

    /*
      Equations 8.9 and 8.10
    */
    easting[i] = K0 * n * (a + (1.0 - t + c) * pow(a,3.0) / 6.0
                    + (5.0 - 18.0 * t + pow(t,2.0) + 72.0 * c - 58.0 * e_prime)
                    * pow(a,5.0) / 120.0) + 500000;
	           
    northing[i] = K0 * (m + n * tan(lat_rad) 
                     * (pow(a,2.0) / 2.0 + (5.0 - t + 9.0 * c + 4.0 * pow(c,2.0))
                     * pow(a,4.0) / 24.0 + (61.0 - 58.0 * t + pow(t,2.0) 
                     + 600.0 * c - 330.0 * e_prime) * pow(a,6.0) / 720.0));
	            
  }  /* end of for */
  
  return(0);

}  /* END OF LL_2_UTM */

/*=============================================================================
  Function to calculate Latitude and Longitude coordinates from UTM coordinates
  See pp.48-65;270-271, USGS Map Projections Working Manual, J.P.Snyder

  Note that output latitude and longitide coordinates are in decimal
  degree units.

  Input parameters:
      northing == array of northing coordinates in decimal meters
      easting == array of easting coordinates in decimal meters
      zone == array of UTM zones for the coordinates
      A == semimajor axis of the earth ellipsoid in meters
      E == eccentricity of the assumed ellipsoid
      E2 == eccentricity squared of the assumed ellipsoid
      ncoords == the number of latitude and longitude coordinate pairs
                 contained in the lat and lon arrays

  Output parameters:
      lat == array of latitude coordinates to be transformed
      lon == array of longitude coordinates to be transformed
=============================================================================*/
long int utm_2_ll (northing, easting, zone, A, E, E2, ncoords, lat, lon)
double *northing, *easting;
long int *zone;
double A, E, E2;
long int ncoords;
double *lat, *lon;
{
  /*
    The flow is from Snyder, solving in order equations:
    8.12, 8.20, 3.24, 7.19, 3.26, 8.21-8.25, 8.17, 8.18 
  */

  double e_prime, east_temp, lambda0, m, mu;
  double phi1, phi1_deg, temp, temp2;
  double c1, d, e1, n1, t1, r1;
	
  long int i;
	
  /*
    Calculate some needed values outside the loop
  */
  e_prime = E2 / (1.0 - E2);
  e1 = (1.0 - sqrt(1.0 - E2)) / (1.0 + sqrt(1.0 - E2));
  temp = (A * (1.0 - E2/4.0 - 3.0*pow(E,4.0)/64.0 - 5.0*pow(E,6.0)/256.0));
  
  /*
    Calculate the latitude and longitude values for each pair of
    input Northing and Easting values
  */
  for (i = 0; i < ncoords; i++)
  {
    /*
      Calculate longitude at Central Meridian for given UTM zone
    */
    lambda0 = (double) (zone[i]-1)*6 - 177;
  
    /*
      Subtract "false" easting (Central Meridian value) from entered values
    */
    east_temp = easting[i] - 500000.0;

    /*
      Calculate values for Eqns 8.12, 8.20, 3.24, 7.19
    */
    m = M0 + northing[i] / K0;
    mu = m / temp;

    /*
      Calculate value for Eqn 3.26
    */
    phi1 = mu + (3.0*e1/2.0 - 27.0*pow(e1,3.0)/32.0) * sin(2.0*mu)
              + (21.0*e1*e1/16.0 - 55.0*pow(e1,4.0)/32.0) * sin(4.0*mu)
              + (151.0*pow(e1,3.0)/96.0) * sin(6.0*mu);
    phi1_deg = phi1*RAD2DEG;
  
    /*
      Calculate values for Eqns 8.21 - 8.25
    */
    c1 = e_prime * cos(phi1) * cos(phi1);
    t1 = tan(phi1) * tan(phi1);
    n1 = A / sqrt(1.0 - E2 * sin(phi1) * sin(phi1));
    temp2 = 1.0 - E2 * sin(phi1) * sin(phi1);
    r1 = A * (1.0 - E2) / pow(temp2,1.5);
    d = east_temp / (n1 * K0);

    /*
      Calculate values for Eqns 8.17 and 8.18
    */
    lat[i]=phi1_deg - (n1 * tan(phi1) / r1) * (RAD2DEG * (pow(d,2.0)/2.0 
                    - (5.0 + 3.0*t1 + 10.0*c1 - 4.0*pow(c1,2.0) 
                    - 9.0*pow(e_prime,2.0)) * pow(d,4.0)/24.0 
                    + (61.0 + 90.0*t1 + 298.0*c1 + 45.0*pow(t1,2.0) 
                    - 252.0*pow(e_prime,2.0) - 3.0*pow(c1,2.0)) * pow(d,6.0)/720.0));
			   
    lon[i] = lambda0 + (RAD2DEG * ( d - (1.0 + 2.0*t1 + c1) * pow(d,3.0)/6.0
                     + (5.0 - 2.0*c1 + 28.0*t1 - 3.0*pow(c1,2.0) + 8.0*e_prime
                     + 24.0*pow(t1,2.0) ) * pow(d,5.0) / 120.0) / cos(phi1));
  }

 return(0);

}  /* END OF UTM2LL */

/*=============================================================================
  Function to convert Albers Equal-Area Conic
  coordinates into the BOREAS (x,y) grid components
  
  Input parameters:
      x == X coordinate in units of AEAC decimal meters
      y == Y coordinate in units of AEAC decimal meters

  Output parameters:
      xgrid == BOREAS X grid coordinate
      ygrid == BOREAS Y grid coordinate
=============================================================================*/
long int albers_2_boreas (x, y, xgrid, ygrid)
double x, y, *xgrid, *ygrid;
{

  *xgrid = x/1000.0;
  *ygrid = y/1000.0;

  return (1);

}  /* END OF ALBERS_2_BOREAS */

/*=============================================================================
  Function to convert BOREAS (x,y) grid component values
  into Albers Equal-Area Conic meters
  
  Input parameters:
      xgrid == BOREAS X grid coordinate in decimal kilometers
      ygrid == BOREAS Y grid coordinate in decimal kilometers

  Output parameters:
      x == X coordinate in units of AEAC decimal meters
      y == Y coordinate in units of AEAC decimal meters
=============================================================================*/
long int boreas_2_albers (xgrid, ygrid, x, y)
double xgrid, ygrid, *x, *y;
{

  *x = xgrid*1000.0;
  *y = ygrid*1000.0;

  return (1);

}  /* END OF BOREAS_2_ALBERS */

/*=============================================================================
  Function to get BOREAS (x,y) coordinate components

  Input parameters:
      infp == FILE pointer to input file of coordinates
      outfp == FILE pointer to output file of coordinates

  Output parameters:
      xgrid == BOREAS grid x coordinate in kilometers (east from origin)
      ygrid == BOREAS grid y coordinate in kilometers (north from origin)
=============================================================================*/
long int read_boreas (infp, outfp, xgrid, ygrid)
FILE *infp, *outfp;
double *xgrid, *ygrid;
{
  long int cnt;
  
  char input_str[80];
  
  /*
    The following do while loop allows reprompting of the
    user for erroneous values and handles returns to calling
    function if the input is a disk file
  */
  cnt = 0;
  do
  {
    /*
      Get an input record from the specified file or user
    */
    if (infp != (FILE *) NULL)
    {
      if (fgets (input_str, 80, infp) == NULL)
      {
        if ( ferror (infp) ) puts ("*** Error reading from input file ***");
        else if ( feof (infp) ) puts ("\n\nEnd of input file reached");
        return (0);
      }
    }
    else
    {
      puts ("\nEnter the X and Y BOREAS grid values as xxx.xx yyy.yy\n");
      puts ("(0 <= X <= 1000     0 <= Y <= 1000) ; q = quit");
      printf (" --> ");
      if (gets (input_str) == NULL)
      {
        puts ("*** Error reading from terminal ***");
        continue;
      }
      if ( !check_response (input_str) ) return (0);
    }

    /*
      See if anything was given in the input string
    */
    if ( strlen(input_str) == 0)
    {
      puts ("No coordinates found on input record");
      if (infp != (FILE *) NULL) return (0);
      else continue;
    }
    /*
      Since it seems like something is in the input string
      try to extract the expected coordinates from it
    */
    cnt = sscanf (input_str, "%lf %lf", xgrid, ygrid);
    if (cnt != 2)
    {
      puts ("*** Error extracting BOREAS coordinates from input record ***");
      cnt = 0;
      if (infp != (FILE *) NULL) return (0);
      else continue;
    }
  
    /*
      Check the ygrid value for validity
    */
    if (*ygrid > 1000.0 || *ygrid < 0.0)
    {
      if (outfp != (FILE *) NULL)
      {
        fputs ("*** Specified BOREAS Y coordinate is outside ***\n", outfp);
        fputs ("*** anticipated range.  Only BOREAS Y values ***\n", outfp);
        fputs ("*** between 0.0 and 1000.00 are accepted     ***\n", outfp);
        puts ("\n*** Specified BOREAS Y coordinate is outside ***");
        puts ("*** anticipated range.  Only BOREAS Y values ***");
        puts ("*** between 0.0 and 1000.00 are accepted     ***");
        return (0);
      }
      puts ("\n*** Specified BOREAS Y coordinate is outside ***");
      puts ("*** anticipated range.  Only BOREAS Y values ***");
      puts ("*** between 0.0 and 1000.00 are accepted     ***");
      cnt = 0;

    }
    /*
      Check the xgrid value for validity
    */
    if (*xgrid > 1000.00 || *xgrid < 0.0)
    {
      if (outfp != (FILE *) NULL)
      {
        fputs ("*** Specified BOREAS X coordinate is outside ***\n", outfp);
        fputs ("*** anticipated range.  Only BOREAS X values ***\n", outfp);
        fputs ("*** between 0.0 and 1000.00 are accepted     ***\n", outfp);
        puts ("\n*** Specified BOREAS X coordinate is outside ***");
        puts ("*** anticipated range.  Only BOREAS X values ***");
        puts ("*** between 0.0 and 1000.00 are accepted     ***");
        return (0);

      }
      puts ("\n*** Specified BOREAS X coordinate is outside ***");
      puts ("*** anticipated range.  Only BOREAS X values ***");
      puts ("*** between 0.0 and 1000.00 are accepted     ***");
      cnt = 0;
    }
    
  } while (cnt == 0);  /* end do while loop */
  
  return (1);  /* successful return */
 
}  /* END OF READ_BOREAS */

/*=============================================================================
  Function to get latitude and longitude values

  Input parameters:
      infp == FILE pointer to input file of coordinates
      outfp == FILE pointer to output file of coordinates

  Output parameters:
      lat == latitude coordinate in decimal degrees (+ N and - S)
      lon ==longitude coordinate in decimal degrees (+ E and - W)
=============================================================================*/
long int read_ll (infp, outfp, lat, lon)
FILE *infp, *outfp;
double *lat, *lon;
{
  long int cnt;
  
  char input_str[80];
  
  /*
    The following do while loop allows reprompting of the
    user for erroneous values and handles returns to calling
    function if the input is a disk file
  */
  cnt = 0;
  do
  {
    /*
      Get an input record from the specified file or user and
      extract the needed latitude and longitude values
    */
    if (infp != (FILE *) NULL)
    {
      if (fgets (input_str, 80, infp) == NULL)
      {
        if ( ferror (infp) ) puts ("*** Error reading from input file ***");
        else if ( feof (infp) ) puts ("\n\nEnd of input file reached");
        return (0);
      }
    }
    else
    {
      puts ("\n\nEnter longitude latitude in decimal degrees (xx.xx yy.yy)");
      puts ("(-111.00 <=lon<= -93.00 and 50.75 <=lat<= 60.00); q = quit");
      printf ("( --> ");
      if ( gets (input_str) == NULL) 
      {
        puts ("*** Error reading from terminal ***");
   		continue;
      }
      if ( !check_response (input_str) ) return (0);
    }

    /*
      See if anything was given in the input string
    */
    if ( strlen(input_str) == 0)
    {
      puts ("No coordinates found on input record");
      if (infp != (FILE *) NULL) return (0);
      else continue;
    }
    /*
      Since it seems like something is in the input string
      try to extract the expected coordinates from it
    */
    cnt = sscanf (input_str, "%lf %lf", lon, lat);
    if (cnt != 2)
    {
      puts ("*** Error extracting lon, lat from input record ***");
      if (infp != (FILE *) NULL) return (0);
      else continue;
    }
  
    /*
      Check the latitude value for validity
    */
    if (*lat > 60.00 || *lat < 50.75)
    {
      if (outfp != (FILE *) NULL)
      {
        fputs ("*** Specified latitude coordinate is outside ***\n", outfp);
        fputs ("*** anticipated range.  Only latitude values ***\n", outfp);
        fputs ("*** between 50.75 and 60.00 are accepted     ***\n", outfp);
        puts ("\n*** Specified latitude coordinate is outside     ***");
        puts ("*** the anticipated range.  Only latitude values ***");
        puts ("*** between 50.75 and 60.00 are accepted         ***");
        return (0);
      }
      else
      {
        puts ("\n*** Specified latitude coordinate is outside     ***");
        puts ("*** the anticipated range.  Only latitude values ***");
        puts ("*** between 50.75 and 60.00 are accepted         ***");
        cnt = 0;
      }
    }
    /*    
	  Check the longitude value for validity
    */
  if (*lon > -93.0 || *lon < -111.0)
    {
      if (outfp != (FILE *) NULL)
      {
        fputs ("*** Specified longitude coordinate is outside ***\n", outfp);
        fputs ("*** anticipated range.  Only longitude values ***\n", outfp);
        fputs ("*** between -93.00 and -111.00 are accepted   ***\n", outfp);
        puts ("\n*** Specified longitude coordinate is outside  ***");
        puts ("***  anticipated range.  Only longitude values ***");
        puts ("*** between -93.00 and -111.00 are accepted    ***");
        return (0);
      }
      else
      {
        puts ("\n*** Specified longitude coordinate is outside  ***");
        puts ("***  anticipated range.  Only longitude values ***");
        puts ("*** between -93.00 and -111.00 are accepted    ***");
        cnt = 0;
      }
    }

  } while (cnt == 0);  /* End of do while loop */

  return (1);  /* successful return */
  
}  /* END OF READ_LL */

/*=============================================================================
  Function to get UTM northing, easting, and zone values

  Input parameters:
      infp == FILE pointer to input file of coordinates
      outfp == FILE pointer to output file of coordinates

  Output parameters:
      northing == UTM northing in decimal meters
      easting == UTM easting in decimal meters
      zone == UTM zone
=============================================================================*/
long int read_utm (infp, outfp, northing, easting, zone)
FILE *infp, *outfp;
double *northing, *easting;
long int *zone;
{
  long int cnt;
  
  char input_str[80];
  
  /*
    The following do while loop allows reprompting of the
    user for erroneous values and handles returns to calling
    function if the input is a disk file
  */
  cnt = 0;
  do
  {
    /*
      Get an input record from the specified file or user and
      extract the needed northing, easting, and zone values
    */
    if (infp != (FILE *) NULL)
    {
      if (fgets (input_str, 80, infp) == NULL)
      {
        if ( ferror (infp) ) puts ("*** Error reading from input file ***");
        else if ( feof (infp) ) puts ("\n\nEnd of input file reached");
        return (0);
      }
    }
    else
    {
      puts ("\n\nEnter: Easting  Northing  Zone  (q = quit)");
      printf ("\r --> ");
      if (gets (input_str) == NULL)
      {
        puts ("*** Error reading from terminal ***");
        continue;
      }
      if ( !check_response (input_str) ) return (0);
    }

    /*
      See if anything was given in the input string
    */
    if ( strlen(input_str) == 0)
    {
      puts ("No coordinates found on input record");
      if (infp != (FILE *) NULL) return (0);
      else continue;
    }
    /*
      Since it seems like something is in the input string
      try to extract the expected coordinates from it
    */
    cnt = sscanf (input_str, "%lf %lf %ld", easting, northing, zone);
    if (cnt != 3)
    {
      puts("*** Error extracting  easting, northing, zone from input ***");
      if (infp != (FILE *) NULL) return (0);
      else continue;
    }

  
    /*
      Check the easting value for validity
    */
    if (*easting > 710496 || *easting < 289504)
    {
      if (outfp != (FILE *) NULL)
      {
        fputs ("*** Specified easting coordinate is outside ***\n", outfp);
        fputs ("*** anticipated range.  Only easting values ***\n", outfp);
        fputs ("*** between 289504 and 710496 are accepted  ***\n", outfp);
        puts ("\n*** Specified easting coordinate is outside ***");
        puts ("*** anticipated range.  Only easting values ***");
        puts ("*** between 289504 and 710496 are accepted  ***");
      }
      else
      {
        puts ("\n*** Specified easting coordinate is outside ***");
        puts ("*** anticipated range.  Only easting values ***");
        puts ("*** between 289504 and 710496 are accepted  ***");
        cnt = 0;
      }
    }
    /*
      Check the northing value for validity
    */
    if (*northing > 6655000 || *northing < 5649600)
    {
      if (outfp != (FILE *) NULL)
      {
        fputs ("*** Specified northing coordinate is outside ***\n", outfp);
        fputs ("*** anticipated range.  Only northing values ***\n", outfp);
        fputs ("*** between 5649600 and 6655000 are accepted ***\n", outfp);
        puts ("\n*** Specified northing coordinate is outside of ***");
        puts ("*** anticipated range.  Only northing values    ***");
        puts ("*** between 5649600 and 6655000 are accepted    ***");
      }
      else
      {
        puts ("\n*** Specified northing coordinate is outside of ***");
        puts ("*** anticipated range.  Only northing values    ***");
        puts ("*** between 5649600 and 6655000 are accepted    ***");
        cnt = 0;
      }
    }

} while (cnt == 0);  /* End do while loop */
  
  return (1);  /* Successful return */
  
}  /* END OF READ_UTM */

/*=============================================================================
Function to compute the shifts in latitude and longitude
between the NAD27 and NAD83 datums.

  Input parameters:
      fp == FILE pointer to datum shift information file
      indatum == string giving input coordinate datum
      ncoors == number of coordinate pairs for which to calculate shifts
      inlat == array of input latitudes to be transformed
      inlong == array of input longitudes to be transformed 

  Output parameters:
      olat == array of shifted output latitudes
      olong == array of shifted output longitudes

Original program written  by Marc Veronneau          January 5, 1989
                          Modified by Mario Berube   May 9, 1989
                          Modified by Steve Farley   January 18 1990
                          Converted to C and modified
                          for BOREAS by Fred Irani, Hughes STX,   Feb 1994     

  This function is based on the FORTRAN-77 subroutine developed for the:
  National Transformation (for converting between NAD27 and NAD83 in Canada)
  by:
  
  		Geodetic Survey of Canada
  		615 Booth Street
  		Ottawa, Canada
  		K1A0E9
  		Phone: (613) 995-44120

=============================================================================*/
long int c_gridint(fp, indatum, ncoors, inlat, inlong, olat, olong)
FILE *fp;
char *indatum;
long int	ncoors;
double *inlat, *inlong, *olat, *olong;
{
  
  char astring[81];	/* Character buffer to read from transform shift file */
  
                         /* limits of current transformation grid: */
  double elong =  89.0,	 /* East longitute */
         llat  =  48.0,  /* Lower Latitude */
         ulat  =  60.0,  /* Upper Latitude */
         wlong = 112.0;  /* West longitude (updated from 110 on 3/23 */

                  /* Transformation grid shift values read from file */
  double alat,    /* First nearest Latitude shift */
         along,   /* First nearest longitude shift */
         blat,    /* Next latitude shift on same line of latitude */
         blong,   /* Next longitude shift on same line of latitude */
         clat,    /* First Nearest Latitude shift on next latitude line */
         clong,   /* First Nearest Longitude shift on next latitude line */
         dlat,    /* Next Latitude shift on upper latitude line */
         dlong,   /* Next Longitude shift on upper latitude line */
         diflat,  /* Calculated Latitude shift for specified lat/long */
         diflong, /* Calculated Longitude shift for specified lat/long */
         dir;     /* Direction of calculation: NAD27->83 or reverse (1, -1) */

  double n,      /* Number of grid shift values to desired latitude  */
         ngrid,  /* Number of grid shift values from south to north  */
         w,      /* Number of grid shift values to desired longitude */
         wgrid;  /* Number of grid shift values from east to west */

  double itemp, jtemp, lon_temp;  /* Temporary calculation values */

  long int i,        /* Target latitude index */
           iread,    /* target record index */
           irow,     /* Target row index */
           iscan,    /* sscanf status */
           j,        /* Target longitude index*/
           jread,    /*  Second read record index */
           num_long, /* Number of values to next longitude */
           pair,     /* Current user-specified lat long pair */
           ibytes,   /* Number of bytes to seek to desired record */
           jbytes;   /* Number of bytes to seek to desired record */
  	
  /*
    Initalize grid increment variables.
  */ 
  ngrid  =  (double) 5.0/60.0;
  wgrid  =  (double) 5.0/60.0;
    
  /*
    Calculate the dimensions of the grid point file
  */
  num_long = (long int) ((wlong - elong) / wgrid) + 1;	/* grid width  */
  
  /*
    Set dir as per datum shift direction.
  */
  if ( (strcmp(indatum, "NAD27") == 0) )
  {
    dir = 1.0;
  }
  else if ( (strcmp(indatum, "NAD83") == 0) )
  {
    dir = -1.0;
  }
  else
  {
    printf ("\n*** Invalid input datum (%s) ***\n", indatum);
    return (0);
  }

  /*
    Process the number of coordinates passed in
    
    Since the software requires a positive longitude value in its
    interpolation, assign the incoming negative longitude value
    to temp place holder variable
  */ 
  for ( pair = 0; pair < ncoors; ++pair)
  { 
 
    lon_temp = -inlong[pair];
    
    /*
    Check that the input lat and long are useable.  Report and return if not.
    ------------------------------------------------------------------------- */
    if (lon_temp < elong)
    {
    	printf("Input or calculated longitude of %lf is east of datum shift boundary.\n",
    	        lon_temp);
    	puts("This condition should not have occurred under the original design of\n");
    	puts("BOR_CORD for BOREAS purposes.\n");
    	return(0);
    }
    if (lon_temp > wlong)
    {
    	printf("Input or calculated longitude of %lf is west of datum shift boundary.\n",
    	        lon_temp);
    	puts("This condition should not have occurred under the original design of\n");
    	puts("BOR_CORD for BOREAS purposes.\n");
       	return(0);
    }
    if (inlat[pair] < llat)
    {
    	printf("Input or calculated latitude of %lf is south of datum shift boundary.\n",
    	        lon_temp);
    	puts("This condition should not have occurred under the original design of\n");
    	puts("BOR_CORD for BOREAS purposes.\n");
       	return(0);
    }
    /*
    For BOREAS, calculated latitude should be only slightly above 60 degrees north.
    If this is the case, treat the input latitude as if it were at 60 degrees 
    and report this to users on final output.
    ------------------------------------------------------------------------------- */
    if (inlat[pair] > ulat)
    {
     printf("Input or calculated latitude of %lf is north of 60.\n", inlat[pair]);
     puts("Latitude shifts will be handled as if at 60 degrees\n");
     puts("Note that datum shift file edge coordinates are interpolated");
     puts("differently than interior coordintes\n");
     it
emp = (double)  60.0 - llat;
    }
    else if (inlat[pair] == ulat)
    {
     printf(
     "Input or calculated latitude of %lf is at 60 degrees north.\n", inlat[pair]);
     puts("Note that datum shift file edge coordinates are interpolated");
     puts("differently than interior coordintes.\n");
     itemp = (double)  inlat[pair] - llat;

    }
    else
    { 
        itemp = (double)  inlat[pair] - llat;
    }
    
    itemp /= ngrid;
    i = (long int) itemp; /* Truncate result */

    jtemp = (double) lon_temp - elong;
    jtemp /= wgrid;
    j = (long int) jtemp; /* Like i, truncate result */
 
    /*
      Calculate where to seek to in the input file of datum shift values
      in order to start searching for the closest 4 grid points to
      compute needed latitude and longitude adjustments
    */
    irow = (long int) i * num_long;
    
   
    iread = irow  + j;  /* iread is: current row + # longs needed */
    jread = iread + 1;  /* jread is: get the next lat/long pair   */

    /*
      Read upper nearest lat long coords for interpolation.
    */ 
    ibytes =  iread * NRECBYTES;
    jbytes =  jread * NRECBYTES;

    fseek(fp, ibytes, 0);
    if (fgets(astring, 80, fp) == NULL)
    {
      if (ferror (fp)) puts ("*** Error reading from datum shift file ***");
      else if (feof (fp)) puts ("\n\nEnd of datum shift file reached");
      return (0);
    }
    iscan = sscanf(astring, "%lf %lf",&alat, &along);
    if ( (iscan != 2) && (inlat[pair] != ulat) )
    {
      puts ("*** Error extracting datum shift values from input 1 ***");
      return (0);
    }
    fseek(fp, jbytes, 0);
    if (fgets(astring, 80, fp) == NULL)
    {
      if (ferror (fp))
      {
       puts ("*** Error reading from datum shift file ***");
      }
      else if ( feof (fp) )
      {  
      	puts ("\n\nEnd of datum shift file reached");
      }
      return (0);
    }
    iscan = sscanf(astring, "%10lf %10lf", &blat, &blong);
    if (iscan != 2)
    {
      puts ("*** Error extracting datum shift values from input 2 ***");
      return (0);
    }
      
    /*
    Handle edge coordinates differently than interior coordinates.
    -------------------------------------------------------------- */
    if (inlat[pair] < 60.0)
    {
      /*
        Read upper nearest lat long coords for interpolation.
      */
      iread += num_long;
      jread += num_long;
      ibytes = iread * NRECBYTES;
      jbytes = jread * NRECBYTES;
      fseek(fp, ibytes, 0);
      if (fgets(astring, 80, fp) == NULL)
      { 
        if (ferror (fp))
        {
        	puts ("*** Error reading from datum shift file ***");
        }
        else if ( (feof (fp) ) && (inlat[pair] != ulat) )
        {  
        	puts ("\n\nEnd of datum shift file reached");
        }
        return (0);
      }
      iscan = sscanf(astring, "%lf %lf", &clat, &clong);
      if (iscan != 2)
      {
        puts ("*** Error extracting datum shift values from input 3 ***");
        return (0);
      }
      fseek(fp, jbytes, 0);
      if (fgets(astring, 80, fp) == NULL)
      {
        if (ferror (fp))
        {
        	puts ("*** Error reading from datum shift file ***");
        }
        else if ( (feof (fp) ) && (inlat[pair] != ulat) )
        {  
        	puts ("\n\nEnd of datum shift file reached");
        }
        return (0);
      }
      iscan = sscanf(astring, "%10lf %10lf", &dlat, &dlong);
      if (iscan != 2)
      {
        puts ("*** Error extracting datum shift values from input 4 ***");
        return (0);
      }
       n = (double) (( inlat[pair] - (llat  + ((i-1) * ngrid))) / ngrid);

    }  /* inlat[pair] < 60 */
    else 
    {
    	/*
    	Interpolate northern edge coordinates without northern values.
    	-------------------------------------------------------------- */
    	clat  = alat;
    	clong = along;
    	dlat  = blat;
    	dlong = blong;
    	n = (double) (( 60.0 - (llat  + ((i-1) * ngrid))) / ngrid);
    }
    
    w = (double) (( lon_temp - (elong + ((j-1) * wgrid))) / wgrid);

    diflat  = (double) (alat  + (clat  - alat ) * n + (blat - alat ) * w + 
                       (alat - blat  - clat  + dlat ) * n * w);

    diflong = (double) (along + (clong - along) * n + (blong - along) * w + 
                       (along - blong - clong + dlong) * n * w);

    /*
      Apply corrections and convert seconds to decimal degrees.
    */
    olat[pair] = (double) (inlat[pair] + (dir * ( diflat  / 3600.0 )));
    olong[pair] = (double) (lon_temp + (dir * ( diflong / 3600.0 )) );
    olong[pair] = -olong[pair];
  }  /* for each pair */

  return (1);

}  /* END OF C_GRIDINT */

/*=============================================================================
  Function to output the coordinate information

  Input parameters:
      outfp == FILE pointer to output file of coordinates
      xgrid == BOREAS X grid coordinate in decimal kilometers
      ygrid == BOREAS Y grid coordinate in decimal kilometers
      lat_27 == latitude coordinate in decimal degrees for NAD27
      lon_27 ==longitude coordinate in decimal degrees for NAD27
      lat_83 == latitude coordinate in decimal degrees for NAD83
      lon_83 ==longitude coordinate in decimal degrees for NAD83
      northing_27 == UTM northing in decimal meters for NAD27
      easting_27 == UTM easting in decimal meters for NAD27
      zone_27 == UTM zone for NAD27
      northing_83 == UTM northing in decimal meters for NAD83
      easting_83 == UTM easting in decimal meters for NAD83
      zone_83 == UTM zone for NAD83

  Output parameters:
      Returned status code.
=============================================================================*/
long int write_coords (outfp, xgrid, ygrid, lat_27, lon_27, lat_83, lon_83,
                  northing_27, easting_27, zone_27, northing_83, easting_83,
                  zone_83)
FILE *outfp;
double xgrid, ygrid, lat_27, lon_27, lat_83, lon_83;
double northing_27, easting_27;
long int zone_27;
double northing_83, easting_83;
long int zone_83;
{
  long int stat;
  
  /*
    Output the coordinate information to the user specified location
  */
  if (outfp == (FILE *) NULL)
  {
    stat=printf("\n\nNAD83: BOREAS X  =  %07.3f    BOREAS Y = %07.3f\n", xgrid, ygrid);
    stat=printf("NAD83  Longitude = %10.5lf  Latitude = %9.5lf\n", lon_83, lat_83);
    stat=printf("NAD83  Easting   = %9.1lf   Northing = %9.1lf  Zone = %2ld\n",
			    easting_83, northing_83, zone_83);
    stat=printf("NAD27  Longitude = %10.5lf  Latitude = %9.5lf\n", lon_27, lat_27);
    stat=printf("NAD27  Easting   = %9.1lf   Northing = %9.1lf  Zone = %2ld\n",
			    easting_27, northing_27, zone_27);
    if (stat < 0) return (0);
  }
  else
  {
    stat = fprintf(outfp, " %8.3lf %8.3lf   %10.5lf %8.5lf",
                            xgrid, ygrid, lon_27, lat_27);
    stat = fprintf(outfp, "     %10.5lf  %7.5lf %10.1lf  %9.1lf   %2.2ld   ",
                           lon_83, lat_83, easting_27, northing_27,  zone_27);
    stat = fprintf(outfp, "%9.1lf  %9.1lf   %2.2ld  \n",
                           easting_83, northing_83,  zone_83);
    if (stat < 0) return (0);
  }
  
  return (1);
  
}  /* END OF WRITE_COORDS */

/*=============================================================================
  Function to Check a user response for a quit and return FALSE
  if user does not wish to continue; TRUE if they do.
=============================================================================*/ 
long int check_response (response)
char	*response;
{
  if ( (strcmp(response, "q") == 0) || (strcmp(response, "Q") == 0) )
  {
    return(FALSE);
  }
  else
  {
    return(TRUE);
  }
     
} /* END OF CHECK_RESPONSE */

/*=============================================================================
  Function to output the software version number banner
=============================================================================*/ 
void version_num()
{
  puts ("\n *============================================*");
  puts   (" *  Version 2.2      BOREAS Coordinate        *");
  puts   (" *                Transformation Software     *");
  puts   (" *============================================*");
} /* END OF VERSION_NUM */