README-UTMupgrade

Readme for mapx library UTM update
Terry Haran 09 June 2004

This document describes a set of changes to the mapx library implemented
in the course of adding the Transverse Mercator (TM) and Universal
Transverse Mercator (UTM) projections and creating a set of unit
tests. Also included was the changing of all float variables and pointers
to double in all the mapx interface routines and data structures, and in
all applications distributed with mapx that use the mapx library. The
latest round of changes also included adding error checking, particularly
for the UTM projection. These changes have now been committed to the cvs
repository. If you don't have an existing directory structure for mapx
under your home directory, you'll have to first type:

cd
mkdir mapx

Then to retrieve mapx from the cvs repository and compile it, type:

cvs checkout maps
make allall

If you don't install mapx into your home directory, you'll have to modify
TOPDIR in Makefile. Also if you're not on a linux system and/or you want
to turn on debugging, you'll have to use one of the alternate definitions
of CONFIG_FLAGS in Makefile. Note that if you turn on debugging, you'll
get different results for the unit tests. The expected results were
produced with optimization set to -O and debugging turned off.

IMPORTANT NOTE: All code that calls mapx routines forward_grid() and
inverse_grid() will have to be modified. The (float) and (more importantly)
the (float *) parameters to these routines have been changed to (double)
and (double *), respectively.

Testing with macct indicates that these latest changes did improve the
accuracy (as reported by macct) significantly in all cases that could be
compared.

     For example, for Azimuthal Equal Area:
     macct N200correct.mpp:
	     Old: 
               102400 points,  0 bad points
               average error = 8.6325e-05 km
               std dev error = 1.3475e-04 km
               maximum error = 8.3848e-04 km
               max error was at 7.90N 37.81E
	     New:
               102400 points,  0 bad points
               average error = 1.8762e-05 km
               std dev error = 4.1020e-05 km
               maximum error = 1.3426e-04 km
               max error was at 14.95N 180.00W

      for Polar Stereographic Ellipsoid:
      macct Sps.mpp:
             Old:
               102400 points,  0 bad points
               average error = 1.2009e-04 km
               std dev error = 1.8439e-04 km
               maximum error = 1.5249e-03 km
               max error was at 29.22S 157.43W
	     New:
               102400 points,  320 bad points
               average error = 1.9059e-05 km
               std dev error = 4.0235e-05 km
               maximum error = 1.3426e-04 km
               max error was at 71.13S 180.00W

The above results are typical. Errors now appear to be less than 1 meter
for all projections over the entire portion of the globe for which the
projection is defined, except for Transverse Mercator Ellipsoid and
Universal Transverse Mercator, which is a special case of Transverse
Mercator Ellipsoid: the equations used for these projections appear to
have 1 meter accuracy only out to about +-9 degrees of longitude from the
central meridian of the projection; beyond this limit the accuracy
decreases rapidly. By default, error checking for the UTM projection is
implemented so that areas having an error of greater than 100 meters
effectively fall outside of the region for which the projection is
defined. This corresponds to about +-17.5 degrees of longitude at a
latitude of +-45 degrees.

      For example, for Universal Tranverse Mercator:
      macct linuxhp_universal_transverse_mercator_e00.gpd:
        Map Projection:                 Universal Transverse Mercator
        Map Equatorial Radius:          6378206.4
        Map Eccentricity Squared:       0.00676866
        Map UTM Zone:                   18
	  102400 points,  90760 bad points
          average error = 1.0886e-02 km
          std dev error = 2.0299e-02 km
          maximum error = 9.8196e-02 km
          max error was at 78.15S 55.86W

The following changes were made to all .c and .h files:

1) <module>_rcsid character strings were added where missing.

2) Any <module>_RCSID identifies were changed to <module>_rcsid.

3) The function id_<module>() was added to all .c files. This function
   returns a pointer to the rcsid string. This was done in order to
   suppress warnings as well as provide a means for applications to
   retrieve version information.

4) The following files were not modified since they were provided from
   third parties:
   a) From the MODIS Reprojection Tool 3.1, gctp directory:
      i)   cproj.h
      ii)  isin.h
      iii) isinfor.c (formerly isinusfor.c)
      iv)  isininv.c (formerly isinusinv.c)
      v)   proj.h
   b) From the U.S. Naval Oceanographic Office:
      i) wdbpltc.c
   The compilations of the above .c files still produce warnings.

In addition to the changes mentioned in 1-4, the following changes were
also made:

5) ppgc.html:
   a) Added statement that keywords are not case sensitive, and that
      a colon must delimit the end of a keyword without any intervening
      spaces, tabs, or newlines.
   b) Added the use of pound sign (#) to denote comments.
   c) Added statement that all text following a semi-colon or pound sign
      on a particular line is ignored.
   d) Added keywords Map Eccentricity Squared and Map Polar Radius.
   e) Provided expanded description of map eccentricity and radius.
   f) Added keywords Map Origin X and Y.
   g) Modified keywords Map False Easting and Map False Northing.
   k) Added keyword Map Center Scale.
   i) Fixed the spelling of Interrupted Homolosine Equal-Area. The
      mispelled form is still supported in the code.
   j) Added the following map projection names:
      i)   Polar Stereographic (ellipsoid)
      ii)  Transverse Mercator
      iii) Transverse Mercator (ellipsoid)
      iv)  Universal Transverse Mercator
   k) Added a separate table and explanation of map parameters for ISin:
      i)  Moved Map Isin Justify to separate ISin table.
      ii) Added Map Isin NZone parameter.
   l) Added a separate table and explanation of map parameters for UTM.
   m) Added keyword Map Maximum Error and an explanation of its use.

6) Makefile:
   a) Changed INSTALL from cp to cp -f so that old files are overwritten.
   b) Changed RANLIB from touch to ranlib to get around a problem on MAC
      OSX.
   c) Changed CONFIG_FLAGS to -O -DLSB1ST.
   d) Added transverse_mercator and universal_transverse_mercator to
      PROJECTION_SRCS and PROJECTION_OBJS, removed isinusfor and
      insinusiv from PROJECTION_SRCS and PROJECTION_OBJS, and removed
      isin.h from PROJECTION_HDRS. 
   e) Added targets:
      i)   allall
      ii)  appall
      iii) testall
      iv)  cleanall
      v)   cleanexes
   f) Under cdb_edit, changed "$(CP) cdb_edit.mpp $(MAPDIR)" to
      "$(INSTALL) cdb_edit.mpp $(MAPDIR)".
   g) Added build of xytest, crtest, and gacct.
   h) Added $(INSTALL) for each test program.
   i) Changed name of tar file to mapsnew.tar.gz temporarily during
      testing.
   j) Created GCTP_SRCS and GCTP_OBJS containing isinfor, isininv, report,
      and cproj, and created GCTP_HDRS containing isin.h, cproj.h, and
      proj.h. All of these files are from the gctp directory in the
      MODIS Reprojection Tool version 3.1.
   k) In the "tar" build, included new directory unit_test, the unit_test
      files utest.pl and utest_new_target.pl, and the unit_test
      subdirectories other, snyder, tilecalc, sgi, and linux.

7) The following changes were made to all projection .c files:
   a) For forward projection functions, the output parameters *u and *v were
      changed to *x and *y, respectively, and the following two lines of
      code:
	*u = current->T00*x + current->T01*y - current->u0;
	*v = current->T10*x + current->T11*y - current->v0;
      were replaced with:
	*x += current->false_easting;
	*y += current->false_northing;
      All input and output parameters were changed from (float) and (float
      *) to (double) and (double *), respectively.
   b) For inverse projection functions, the input parameters u and v were
      changed to x and y, respectively, and the following two lines of
      code:
        x =  current->T00*(u+current->u0) - current->T01*(v + current->v0);
        y = -current->T10*(u+current->u0) + current->T11*(v + current->v0);
      were replaced with
        x -= current->false_easting;
        y -= current->false_northing;
      All input and output parameters were changed from (float) and (float
      *) to (double) and (double *), respectively.

8) albers_conic_equal_area.c:
   a) Added include of proj.h for use of asinz() in place of asin().
   b) In inverse_albers_conic_equal_area() and 
      albers_conic_equal_area_ellipsoid(), in the calculation of theta,
      the arguments to atan2() are now negated when n is negative as
      specified in Snyder.
   c) In inverse_albers_conic_equal_area_ellipsoid(), the series equation
      3-18 from Snyder was replaced with the iterative solution based on
      equation 3-16 in Snyder. This, together with the float to double
      changes mentioned above, improved the overall accuracy by about
      a factor of about 2500:
      i)   macct sgi_albers_conic_equal_area_e00.gpd
             Old:
               102400 points,  336 bad points
               average error = 3.6677e-02 km
               std dev error = 9.4522e-02 km
               maximum error = 4.4933e+00 km
               max error was at 90.00S 124.70E
	     New:
               102400 points,  0 bad points
               average error = 1.3655e-05 km
               std dev error = 3.4695e-05 km
               maximum error = 1.3426e-04 km
               max error was at 73.07S 96.49E

9) azimuthal_equal_area.c:
   a) Added include of proj.h for use of asinz() in place of asin().
   b) In init_azimuthal_equal_area_ellipsoid(), an erroneous
      recalculation of qp was removed as well as redundant calculations
      of q1, cos_phi1, and sin_phi1.
   c) In azimuthal_equal_area_ellipsoid(), for the south polar case,
      a test of fabs(qp - q) was changed to fabs(qp + q) prior to the
      calculation of rho.
   d) In inverse_azimuthal_equal_area_ellipsoid(), the series equation
      3-18 from Snyder was replaced with the iterative solution based on
      equation 3-16 in Snyder. This, together with the float to double
      changes mentioned above, improved the average accuracy for the
      north polar case by about a factor of 280 (the old oblique and south
      polar cases had bugs):
      i)   macct sgi_azimuthal_equal_area_e01.gpd
	     Old:
               102400 points,  150 bad points
	       average error = 4.6800e-03 km
	       std dev error = 9.9562e-02 km
	       maximum error = 3.5940e+00 km
	       max error was at 90.00S 51.35W
	     New:
               102400 points,  320 bad points
               average error = 1.6654e-05 km
               std dev error = 3.6599e-05 km
               maximum error = 1.3426e-04 km
               max error was at 73.07S 176.61W

10) cdb.c:
    a) Changed float to double throughout.
    b) In init_cb(),

      changed:
        this->index_order = this->header->index_order;
      to:        
        this->index_order = (cdb_index_sort)(this->header->index_order);

      in order to suppress a compiler warning.

11) cdb.h:
    a) Changed float to double throughout.

12) cdb_byteswap.h:
    a) added conditional compile testing definition of cdb_byteswap_h_.
    b) added definition of cdb_byteswap_h_.

13) cdb_edit.c:
    a) Changed float to double throughout.
    b) In main(), removed declaration of unused variables "i", "ios", and
       "extent".
    c) In finish_new_file(), removed declaration of unused variable
       "command".
    d) In join_map(), append_candidate(), and cdb_index_entry(), added a
       (byte4) typecast for NULL used in integer comparisons to suppress
       compiler warnings.

14) cdb_list.c:
    a) In main(), removed unused variable "i".

15) cproj.c:
    a) New file copied from the MODIS Reprojection Tool 3.1, gctp
       directory. Contains utilities needed for ISin routines.

16) cproj.h:
    a) New file copied from the MODIS Reprojection Tool 3.1, gctp
       directory. Needed by cproj.c.

17) cylindrical_equal_area.c:
    a) In cylindrical_equal_area() and cylindrical_equal_area_ellipsoid(),
       changed "dlon" from float to double.
    b) In init_cylindrical_equal_area_ellipsoid(), changed use of
       current->lat0 (Reference Latitude) to current->lat1 (Second
       Reference Latitude) to match init_cylindrical_equal_area() as well
       as the documentation. Thus Reference Latitude is effectively ignored
       for both the spherical and ellipsoidal cylindrical equal area
       projections.

18) cylindrical_equidistant.c:
    a) In cylindrical_equidistant(), changed "dlon" from float to double.

19) gridloc.c:
    a) Changed float variables to double as needed for call to
       inverse_grid().
    b) Added -D option to optionally produce output in double precision.
    c) Modifed behavior of -o option so that when -D is specified, the
       output file has the name output_name.WIDTHxHEIGHTxNBANDS.double.

20) grids.c:
    a) Changed float to double throughout.
    b) In decode_gpd(), changed call to new_mapx() to include new "quiet"
       parameter se to TRUE to suppress messages about an unknown projection.
    c) Removed static function next_line_from_buffer() so that modified
       version of next_line_from_buffer() in mapx is now used by
       old_fixed_format_decode_gpd().
    d) In gtest main(), changed "enter r s: " prompt to "enter col row: ".
    e) Added crtest, an interactive test similar to gtest, but produces
       output that is compatible with utest.pl.
    f) Added gacct, a grid accuracy test similar to macct, but for .gpd
       files rather than .mpp files. gacct tests the accuracy for each
       point in the grid.

21) grids.h:
    a) Changed float to double throughout.

22) integerized_sinusoidal.c:
    a) In init_integerized_sinusoidal(), replaced calculation of nrows
       with use of Map ISin NZone (isin_nzone) parameter value.
    b) In init_integerized_sinusoidal(), added passing of false easting
       and false northing values to Isin_inv_init().

23) interupted_homolosine_equal_area.c:
    a) In interupted_homolosine_equal_area(), changed max_it from 30 to
       35 to fix non-convergence problem near latitude = -90.0 on linux.

24) irregrid.c:
    a) Changed all float variables to double except for variables
       defining buffers for input and output data: these are still float.

25) isin.h:
    a) New version copied from the MODIS Reprojection Tool 3.1, gctp
       directory.

26) isinfor.c:
    a) New version (replacing isinusfor.c) copied from the MODIS
       Reprojection Tool 3.1, gctp directory.

27) isininv.c:
    a) New version (replacing isinusinv.c) copied from the MODIS
       Reprojection Tool 3.1, gctp directory.

28) keyval.c:
    a) In get_label_keyval(), added feature such that for each line,
       any characters following a semi-colon or a pound sign are replaced
       with blanks.
    b) In get_field_keyval(), added feature such that keyword and label
       being searched are both converted to uppercase before the search
       is performed, and that a valid keyword must be followed immediately
       by a colon.
    c) In get_value_keyval(), modified behavior such that for format
       strings "%lat" and "%lon", the value field returned by
       lat_lon_keyval() is now parsed as a double rather than a float.
    d) In get_value_keyval(), added feature such that if default_string is
       "KEYVAL_UNINITIALIZED" and format is "%f" or "%lf", and the keyword
       is not found, then the value KEYVAL_UNINITIALIZED is returned.
    e) In lat_lon_keyval(), changed the value parameter from (float *) to
       (double *).

29) keyval.h:
    a) Changed float to double throughout.
    b) Added definition of KEYVAL_UNINITIALIZED to be FLT_MAX, if FLT_MAX
       is defined, or 9e30 if it isn't.

30) lambert_conic_conformal.c:
    a) In init_lambert_conic_conformal_ellipsoid(), added use of Map
       Origin Latitude (center_lat) in addition to Map Reference Latitude
       (lat0) and Map Second Reference Latitude (lat1) as is done in
       init_albers_conic_equal_area().
    b) In inverse_lambert_conic_conformal_ellipsoid(), replaced lower
       accuracy series equation 3-5 from Snyder with higher accuracy
       iteration involving equation 7-9.

31) mapenum.c:
    a) Changed float to double throughout.

32) maps.c:
    a) Changed float to double throughout.

33) maps.h:
    a) Changed float to double throughout.

34) mapx.c:
    a) Changed float to double throughout.
    b) added include isin.h
    c) Removed static declaration of next_line_from_buffer(). It is
       now declared externally in mapx.h so that the modified version
       may be used by old_fixed_format_decode_gpd() in grids.c.
    d) Added declarations of the following functions:
       i)   init_transverse_mercator()
       ii)  transverse_mercator()
       iii) inverse_transverse_mercator()
       iv)  init_transverse_mercator_ellipsoid()
       v)   transverse_mercator_ellipsoid()
       vi)  inverse_transverse_mercator_ellipsoid()
       vii) init_universal_transverse_mercator()
       viii)forward_xy_mapx_check()
       ix)  inverse_xy_mapx_check()
       x)   dist_latlon_map_units()
       xi)  dist_xy_map_units()
       Note that no forward or inverse functions for UTM are declared,
       since UTM uses transverse mercator forward and inverse functions.
    e) In init_mapx(), modified call to new_mapx() to include new "quiet"
       parameter set to FALSE so that a message about an unknown
       projection won't be suppressed.
    f) In new_mapx(), added new parameter "quiet" which, if set,
       suppresses displaying a message abount an unknown projection.
    g) In new_mapx(), added initialization for Transverse Mercator and
       Universal Transverse Mercator.
    h) In decode_mpp():
       i)   Added code to make Map Reference Latitude and Longitude
            optional keywords for UTM and ISin.
       ii)  Added code to process Map Origin X and Y keywords. These
            values are stored in x0 and y0 in the mapx structure.
       iii) Deferred assuming that Map Origin Latitude and Longitude
            take on Reference Latitude and Longitude values,
            respectively, when the former are not defined and the
            projection is UTM, until UTM initialization.
       iv)  Changed assignment of Map False Easting and Northing keyword
            values from u0, v0 to false_easting, false_northing,
            respectively, in the mapx structure.
       v)   Deferred assigning of default values to Map False Easting and
            Northing keywords when they are not defined and the
            projection is UTM until UTM initialization.
       vi)  Added code to process Map Center Scale, Map UTM Zone,
            Map ISin NZone, Map Eccentricity Squared, Map Polar Radius,
            and Map Maximum Error keywords.
       vii) Added code to assign default WGS-84 values for equatorial
            radius and eccentricity if map projection is UTM.
       viii)Added code to assign default value for equatorial
            radius if map projection is ISin.
       ix)  Added code to make eccentricity = 0 for spherical projections.
       x)   Added code to derive eccentricity from eccentricity squared
            as needed.
       xi)  Added code to derive eccentricity from equatorial radius and
            polar radius as needed.
       xii) Added code to derive equatorial radius from polar radius and
            eccentricity as needed.
    i) In old_fixed_format_decode_mpp():
        i)  Added code to set the default equatorial radius and
            eccentricity based on the map projection.
        ii) Added code to initialize the following new mpp parameters that
            are not supported by the old format:
              x0
              y0
	      false_easting
	      false_northing
	      center_scale
	      utm_zone
	      isin_nzone
	      isin_justify
              e2 (eccentricity squared)
              polar_radius
	      maximum_error
    j) In reinit_mapx():
       i)   Added code to derive polar radius from equatorial radius and
            eccentricity squared (after eccentricity squared has been
            computed from eccentricity).
       ii)  Added code to set the flattening (f).
       iii) If Map Origin X and Y not provided, then forward transform
            Map Origin Latitude and Longitude (center_lat and center_lon
            in the mapx structure) into Map Origin X and Y (x0 and y0
            in the mapx structure).
       iv)  Changed the initialization of u0 and v0 in the mapx structure.
            The old definitions of u0 and v0 were false easting and false
            northing, respectively. The new definition is the u and v values
            resulting from the forward rotation (if any) of (x0,y0). Thus
            (x0,y0) serves as the center of the rotation.
    k) In next_line_from_buffer(), added functionality such that lines
       beginning with # or ; are ignored. This allows comments to be
       inserted anywhere in old style gpd and mpp files which is needed
       for making unit tests using old style format.
    l) In forward_mapx():
       i)   Changed call to (*(this->geo_to_map))() so that pointers to
            x and y are passed rather than u and v. (x,y) is then rotated
            and translated to (u,v).
       ii)  Added call to forward_xy_mapx_check().
    m) In inverse_mapx():
       i)   Translate and rotate (u,v) to (x,y). Changed call to
            (*(this->map_to_geo))() so that x and y are passed rather than
            u and v.
       ii)  Added call to inverse_xy_mapx_check().
    n) Added forward_xy_mapx() which converts (lat,lon) to (x,y) and
       includes a call to forward_xy_mapx_check().
    o) Added inverse_xy_mapx() which converts (x,y) to (lat,lon) and
       includes a call to inverse_xy_mapx_check().
    p) Added the following functions:
       i)   forward_xy_mapx_check().
       ii)  inverse_xy_mapx_check().
       iii) dist_latlon_map_units().
       iv)  dist_xy_map_units().
    q) In standard_name():
       i)   Added support for correct spelling of "interrupted" for
            Goode Interrupted Homolosine Equal Area. Incorrect spelling
            is still supported.
       ii)  Added Transverse Mercator.
       iii) Added Transverse Mercator Ellipsoid.
       iv)  Added Universal Transverse Mercator.
            Note that Universal Transverse Mercator Ellipsoid is mapped
            to Universal Transverse Mercator.
    r) Added main program xytest, an interactive test for mapx routines
       that is similar to mtest, but works with original map coordinates
       (x,y) rather than rotated and translated coordinates (u,v).
    s) In mtest:
       i)   Changed declaration of readln[80] to readln[FILENAME_MAX].
       ii)  Added ability to specify the .mpp file name on the command
            line just as in gtest.

35) mapx.h:
    a) Changed float to double throughout.
    b) Added definitions of WGS-84 equatorial radius and eccentricity.
    c) Added definition of default sphere radius for ISin.
    d) Moved e2 from private to user specified constants in mapx structure.
    e) Added user specified constants to mapx structure:
       i)   x0
       ii)  y0
       iii) false_easting
       iv)  false_northing
       v)   center_scale
       vi)  utm_zone
       vii) isin_nzone
       viii)isin_justify
       ix)  polar_radius
       x)   maximum_error
    f) Added private projection constants to mapx structure:
       i)    t2
       ii)   e0
       iii)  e1p
       iv)   e2p
       v)    e3p
       vi)   ml0
       vii)  esp
       viii) f1
       ix)   f2
       x)    f3
       xi)   f4
       xii)  f
    g) Modified function prototype:
       i) new_mapx() now includes new "quiet" parameter
    h) Added function prototypes:
       i)   next_line_from_buffer() (formerly static)
       ii)  forward_xy_mapx()
       iii) inverse_xy_mapx()

36) mercator.c:
    a) In mercator(), changed "dlon" from float to double.

37) mollweide.c:
    a) In mollweide(), changed "dlon" from float to double.
    b) In mollweide(), changed the value of epsilon from 0.0025 to 1e-6.
       This, together with the float to double changes mentioned above,
       improved the average accuracy by about a factor of 209:
       i)   macct linuxhp_mollweide_s00.gpd
              Old:
	        102400 points,  0 bad points
                average error = 1.4900e-03 km
                std dev error = 3.0500e-03 km
                maximum error = 1.5307e-02 km
                max error was at 89.44S 89.72E
	      New:
                102400 points,  0 bad points
                average error = 7.1204e-06 km
                std dev error = 2.5006e-05 km
                maximum error = 9.4939e-05 km
                max error was at 45.99S 180.00W

38) orthographic.c:
    a) In orthographic(), changed the following temporary variables
       from float to double:
       i)   phi
       ii)  lam
       iii) sin_phi
       iv)  cos_phi
       v)   sin_lam
       vi)  cos_lam
       vii) cos_beta

39) polar_stereographic.c:
    a) In polar_stereographic(), changed the following temporary variables
       from float to double:
       i)   phi
       ii)  lam
       iii) rho
    b) In inverse_polar_stereographic(), fixed a bug by negating sin(phi)
       in the calculation of q for the south polar aspect.
    c) In init_polar_stereographic_ellipsoid(), changed the following
       temporary variables from float to double:
       i)   numerator
       ii)  denominator
    d) In polar_stereographic_ellipsoid(), changed the following
       temporary variables from float to double:
       i)   phi
       ii)  lam
       iii) rho
       iv)  t
       v)   numerator
       vi)  denominator
    e) Also in polar_stereographic_ellipsoid(), changed the use of
       equation 15-9 in Snyder to 15-9a.

40) proj.h:
    a) New file copied from the MODIS Reprojection Tool 3.1, gctp
       directory. Needed by cproj.h.

41) regrid.c:
    a) Changed all float variables to double except for variables
       defining buffers for input and output data: these are still float.

42) report.c:
    a) New file copied from the MODIS Reprojection Tool 3.1, gctp
       directory. Contains utilities needed by cproj.c.

43) resamp.c:
    a) Changed float to double throughout.

44) sinusoidal.c:
    a) In sinusoidal(), changed "dlon" from float to double.

45) Added transverse_mercator.c.

46) Added universal_transverse_mercator.c.

47) In wdbtocdb.c:
    a) Changed declaration of move_pu(), draw_pd(), and
       write_segment_data() from int to void to suppress warnings about
       not returning values.
    b) Changed declaration of thin from char to int to suppress a warning.

48) Added unit_test directory containing utest.pl, utest_new_target.pl,
    and several subdirectories containing gpd and mpp files used for unit
    tests.
    a) In the unit_test directory, typing utest.pl produces a rather
       lengthy usage message, which you can peruse if you're so inclined.

    b) On an sgi system, you can run all the sgi unit tests by typing:

         utest.pl sgi/sgi*

       on a linux system, try typing:

         utest.pl linux/linux*

       on a hp system running linux, type:

         utest.pl linuxhp/linuxhp*

       No screen output should be produced on any system. To produce
       screen output for each test regardless of the outcome, use the
       -v option. Note that you'll have to have compiled the mapx library
       with optimization set to -O and debugging turned off since this is
       how the expected results files were created (see CONFIG_FLAGS in
       the Makefile).

    c) You can also create a new set of unit tests for a new target (or
       regenerate the unit tests for an existing target) using
       the source unit tests in the "snyder", "tilecalc", and "other"
       subdirectories by means of the utest_new_target.pl script. This
       is how the sgi, linux, and linuxhp subdirectories were created on
       glacier, icewave, and megadune, respectively.

       Note that the -v option can also be used with utest_new_target.pl.

       Let's say you're on a sun, and you want to create a new set of unit
       tests. Then in the unit_test directory, you would type:

         utest_new_target.pl -v sun >&source2sun.log

       This will create a new set of unit tests in a new sun subdirectory.
       Once this is done, you can browse source2sun.log to verify that
       the only differences were very small. Note that all of the "snyder"
       and "tilecalc" xytests will have expected results; many of the
       "other" expected results in the source unit tests are specified as
       "dummy". In either case, expected results are replaced with actual
       results in creating the new expected results in the output files.

       Then to run the new unit tests, you would type:

         utest.pl sun/sun*

       No screen output should be produced.

       Alternatively, instead of using the source tests for expected
       results, you could create the sun tests using the sgi tests for
       the expected results. To do this, you would type:

         utest.pl -v -o sgi2sun -c sgi/sgi* >&sgi2sun.log

       The -o sgi2sun option says that a subdirectory called sgi2sun
       should be created containing the new expected results. The -c
       option indicates that a comment should be appended to each new
       xytest or crtest expected result indicating what the sgi expected
       result was.

       Again, you can browse sgi2sun.log to verify that the differences
       were small.

    d) Each of the unit test subdirectories mentioned above contain gpd
       and/or mpp files which contain "unit test" expected results to be
       used with utest.pl. Most of the files have names of the form:

         <tag>_<projection>_<s|e><nn>.gpd
	   where
             <tag> is the name of the subdirectory.
	     <projection> is the name of the map projection to be tested.
	     <s|e> is s for spherical projections,
                      e for ellipsoidal projections.
             <nn> is a two-digit serial number.

       Only the tag part of the filename has meaning to utest.pl, but even
       this behavior can be overridden by use of the -i option as
       described in the utest.pl usage.

    e) As an example, here is snyder/snyder_albers_conic_equal_area_e00.gpd:

# file: snyder_albers_conic_equal_area_e00.gpd
# Map Projections--A Working Manual, John P. Snyder
# Appendix A Numerical Examples
# Albers Conical Equal-Area Ellipsoid pp. 292-294
#
Map Projection:                 Albers Conic Equal-Area Ellipsoid
Map Equatorial Radius:          6378206.4
Map Eccentricity:               0.0822719
Map Reference Latitude:         29.5
Map Second Reference Latitude:  45.5
Map Origin Latitude:            23.0
Map Reference Longitude:        -96.0
Grid Width:                     300
Grid Height:                    300
Grid Map Units per Cell:        15000
Grid Map Origin Column:         149.5
Grid Map Origin Row:            149.5
#
# macct
#   dummy
#   dummy
#   dummy
#   dummy
#   dummy
#
# gacct
#   dummy
#   dummy
#   dummy
#   dummy
#   dummy
#
# xytest forward
#   lat,lon = 35.0 -75.0
#   x,y = 1885472.7       1535925.0
#
# xytest inverse
#   x,y = 1885472.7 1535925.0
#   lat,lon = 35.0000015 -75.0000003
#
# crtest forward
#   lat,lon = 35.0 -75.0
#   col,row = dummy dummy
#
# crtest inverse
#   col,row = 275.1981819 47.1050001
#   lat,lon = dummy dummy

    f) As mentioned above, to see produce output regardless of the
       results, use the -v option. For example on a linuxhp, typing:

         utest.pl -v snyder/snyder_albers_conic_equal_area_e00.gpd

       produces:

****************************************************************
file: snyder/snyder_albers_conic_equal_area_e00.gpd
  Map Projection:                 Albers Conic Equal-Area Ellipsoid
  Map Equatorial Radius:          6378206.4
  Map Eccentricity:               0.0822719
  Map Reference Latitude:         29.5
  Map Second Reference Latitude:  45.5
  Map Origin Latitude:            23.0
  Map Reference Longitude:        -96.0
  Grid Width:                     300
  Grid Height:                    300
  Grid Map Units per Cell:        15000
  Grid Map Origin Column:         149.5
  Grid Map Origin Row:            149.5
  **************************************************************
  macct expected results:
    dummy
    dummy
    dummy
    dummy
    dummy
  macct actual results:
    102400 points,  0 bad points
    average error = 4.8582e-06 km
    std dev error = 2.0920e-05 km
    maximum error = 9.4939e-05 km
    max error was at 40.34S 159.69W
  Expected and actual results differ
  **************************************************************
  gacct expected results:
    dummy
    dummy
    dummy
    dummy
    dummy
  gacct actual results:
    90000 points,  0 bad points
    average error = 4.6741e-12 pixels
    std dev error = 1.2855e-11 pixels
    maximum error = 7.5640e-11 pixels
    max error was at col: 257  row: 158   lat: 20.633791  lon: -80.890740
  Expected and actual results differ
  **************************************************************
  xytest forward results:
    lat,lon in:       35.0 -75.0
    x,y expected:     1885472.7       1535925.0

    x,y actual:       1885472.7282290 1535924.9987983
  Expected and actual results differ
  **************************************************************
  xytest inverse results:
    x,y in:       1885472.7 1535925.0
    lat,lon expected:     35.0000015 -75.0000003

    lat,lon actual:       35.0000001 -75.0000003
  Expected and actual results differ
  **************************************************************
  crtest forward results:
    lat,lon in:       35.0 -75.0
    col,row expected:     dummy dummy

    col,row actual:       275.1981819 47.1050001
  Expected and actual results differ
  **************************************************************
  crtest inverse results:
    col,row in:       275.1981819 47.1050001
    lat,lon expected:     dummy dummy

    lat,lon actual:       35.0000000 -75.0000000
  Expected and actual results differ
****************************************************************

    g) Again on a linuxhp, typing:

         utest.pl -v linuxhp/linuxhp_albers_conic_equal_area_e00.gpd

       produces:

****************************************************************
file: linuxhp/linuxhp_albers_conic_equal_area_e00.gpd
  Map Projection:                 Albers Conic Equal-Area Ellipsoid
  Map Equatorial Radius:          6378206.4
  Map Eccentricity:               0.0822719
  Map Reference Latitude:         29.5
  Map Second Reference Latitude:  45.5
  Map Origin Latitude:            23.0
  Map Reference Longitude:        -96.0
  Grid Width:                     300
  Grid Height:                    300
  Grid Map Units per Cell:        15000
  Grid Map Origin Column:         149.5
  Grid Map Origin Row:            149.5
  **************************************************************
  macct expected results:
    102400 points,  0 bad points
    average error = 4.8582e-06 km
    std dev error = 2.0920e-05 km
    maximum error = 9.4939e-05 km
    max error was at 40.34S 159.69W
  macct actual results:
    102400 points,  0 bad points
    average error = 4.8582e-06 km
    std dev error = 2.0920e-05 km
    maximum error = 9.4939e-05 km
    max error was at 40.34S 159.69W
  Expected results match actual results
  **************************************************************
  gacct expected results:
    90000 points,  0 bad points
    average error = 4.6741e-12 pixels
    std dev error = 1.2855e-11 pixels
    maximum error = 7.5640e-11 pixels
    max error was at col: 257  row: 158   lat: 20.633791  lon: -80.890740
  gacct actual results:
    90000 points,  0 bad points
    average error = 4.6741e-12 pixels
    std dev error = 1.2855e-11 pixels
    maximum error = 7.5640e-11 pixels
    max error was at col: 257  row: 158   lat: 20.633791  lon: -80.890740
  Expected results match actual results
  **************************************************************
  xytest forward results:
    lat,lon in:       35.0 -75.0
    x,y expected:     1885472.7282290 1535924.9987983
                  .vs 1885472.7       1535925.0 in snyder
    x,y actual:       1885472.7282290 1535924.9987983
  Expected results match actual results
  **************************************************************
  xytest inverse results:
    x,y in:       1885472.7 1535925.0
    lat,lon expected:     35.0000001 -75.0000003
                      .vs 35.0000015 -75.0000003 in snyder
    lat,lon actual:       35.0000001 -75.0000003
  Expected results match actual results
  **************************************************************
  crtest forward results:
    lat,lon in:       35.0 -75.0
    col,row expected:     275.1981819 47.1050001

    col,row actual:       275.1981819 47.1050001
  Expected results match actual results
  **************************************************************
  crtest inverse results:
    col,row in:       275.1981819 47.1050001
    lat,lon expected:     35.0000000 -75.0000000

    lat,lon actual:       35.0000000 -75.0000000
  Expected results match actual results
****************************************************************

    h) Some of the files in the "other" subdirectory contain gpd and
       mpp files using the old format. These have filenames of the form:

         other_<standard_name>.<gpd|mpp>

       For example, here is other/other_Na25.gpd:

# file: other_Na25.gpd
# Unit Test
# EASE-Grid Northern Hemisphere AVHRR 25 km
#
other/other_N200correct.mpp map projection parameters   # AVHRR 25km EASE-Grid
361 361         columns rows                            # Northern Hemisphere
8               grid cells per map unit                 # 25 km
180.0 180.0     map origin column,row                   #
#
# gacct
#   dummy
#   dummy
#   dummy
#   dummy
#   dummy
#
# crtest forward
#   lat,lon = 75.0 -150.0
#   col,row = dummy dummy
#
# crtest inverse
#   col,row = 146.8251201 122.5394225
#   lat,lon = dummy dummy

       Note that old format gpd files can only contain gacct and crtest
       tests; the macct and xytest tests are contained in the
       corresponding mpp file, in this case other/other_N200correct.mpp:

# file: other_N200corrent.mpp
# Unit Test
# EASE-Grid Northern Hemisphere
#
Azimuthal Equal-Area
90.0    0.0     lat0 lon0
0.0             rotation
200.5402        scale (km/pixel)
90.00   00.00   center lat lon
0.00   90.00    lat min max
-180.00  180.00 lon min max
15.00 30.00     grid
0.00    00.00   label lat lon
1 0 0           cil bdy riv
#
# macct
#   dummy
#   dummy
#   dummy
#   dummy
#   dummy
#
# xytest forward
#   lat,lon = 75.0 -150.0
#   x,y = dummy dummy
#
# xytest inverse
#   x,y = -4.1468600 7.1825722
#   lat,lon = dummy dummy

    i) The unit_test gpd and mpp files produced to date do not necessarily
       exercise all the paths through the code for each projection. More
       files could be added at some future date to do this.