New

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@848 4e78687f-474d-0410-85f9-8d5e500ac6b2

8 files changed, 1869 insertions, 0 deletions

diff --git a/src/cs2cs.c b/src/cs2cs.c
new file mode 100644
index 00000000..4c1005c8
--- /dev/null
+++ b/src/cs2cs.c
@@ -0,0 +1,408 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Mainline program sort of like ``proj'' for converting between
+ *           two coordinate systems.
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2000/07/06 23:32:27  warmerda
+ * New
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include <math.h>
+#include "projects.h"
+#include "emess.h"
+
+#define MAX_LINE 200
+#define MAX_PARGS 100
+
+static PJ   *fromProj, *toProj;
+
+static int
+reversein = 0,	/* != 0 reverse input arguments */
+reverseout = 0,	/* != 0 reverse output arguments */
+echoin = 0,	/* echo input data to output line */
+tag = '#';	/* beginning of line tag character */
+	static char
+*oform = (char *)0,	/* output format for x-y or decimal degrees */
+*oterr = "*\t*",	/* output line for unprojectable input */
+*usage =
+"%s\nusage: %s [ -beEfiIlorstvwW [args] ] [ +opts[=arg] ]\n"
+"                   [+to [+opts[=arg] [ files ]\n";
+
+static struct FACTORS facs;
+static double (*informat)(const char *, 
+                          char **); /* input data deformatter function */
+
+
+/************************************************************************/
+/*                              process()                               */
+/*                                                                      */
+/*      File processing function.                                       */
+/************************************************************************/
+static void process(FILE *fid) 
+
+{
+    char line[MAX_LINE+3], *s, pline[40];
+    projUV data;
+
+    for (;;) {
+        double z;
+
+        ++emess_dat.File_line;
+        if (!(s = fgets(line, MAX_LINE, fid)))
+            break;
+        if (!strchr(s, '\n')) { /* overlong line */
+            int c;
+            (void)strcat(s, "\n");
+				/* gobble up to newline */
+            while ((c = fgetc(fid)) != EOF && c != '\n') ;
+        }
+        if (*s == tag) {
+            fputs(line, stdout);
+            continue;
+        }
+
+        if (reversein) {
+            data.v = (*informat)(s, &s);
+            data.u = (*informat)(s, &s);
+        } else {
+            data.u = (*informat)(s, &s);
+            data.v = (*informat)(s, &s);
+        }
+
+        z = strtod( s, &s );
+
+        if (data.v == HUGE_VAL)
+            data.u = HUGE_VAL;
+
+        if (!*s && (s > line)) --s; /* assumed we gobbled \n */
+
+        if ( echoin) {
+            int t;
+            t = *s;
+            *s = '\0';
+            (void)fputs(line, stdout);
+            *s = t;
+            putchar('\t');
+        }
+
+        if (data.u != HUGE_VAL) {
+            if( pj_transform( fromProj, toProj, 1, 0, 
+                              &(data.u), &(data.v), &z ) != 0 )
+            {
+                data.u = HUGE_VAL;
+                data.v = HUGE_VAL;
+            }
+        }
+
+        if (data.u == HUGE_VAL) /* error output */
+            fputs(oterr, stdout);
+
+        else if (toProj->is_latlong && !oform) {	/*ascii DMS output */
+            if (reverseout) {
+                fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
+                putchar('\t');
+                fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
+            } else {
+                fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
+                putchar('\t');
+                fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
+            }
+
+        } else {	/* x-y or decimal degree ascii output */
+            if ( toProj->is_latlong ) {
+                data.v *= RAD_TO_DEG;
+                data.u *= RAD_TO_DEG;
+            }
+            if (reverseout) {
+                printf(oform,data.v); putchar('\t');
+                printf(oform,data.u);
+            } else {
+                printf(oform,data.u); putchar('\t');
+                printf(oform,data.v);
+            }
+        }
+
+        putchar(' ');
+        printf( "%.3f", z );
+        fputs("\n", stdout );
+    }
+}
+
+/************************************************************************/
+/*                                main()                                */
+/************************************************************************/
+
+int main(int argc, char **argv) 
+{
+    char *arg, **eargv = argv, *from_argv[MAX_PARGS], *to_argv[MAX_PARGS],
+        **iargv = argv;
+    FILE *fid;
+    int from_argc=0, to_argc=0, iargc = argc, eargc = 0, c, mon = 0;
+    int have_to_flag = 0, inverse = 0, i;
+
+    if (emess_dat.Prog_name = strrchr(*argv,DIR_CHAR))
+        ++emess_dat.Prog_name;
+    else emess_dat.Prog_name = *argv;
+    inverse = ! strncmp(emess_dat.Prog_name, "inv", 3);
+    if (argc <= 1 ) {
+        (void)fprintf(stderr, usage, pj_release, emess_dat.Prog_name);
+        exit (0);
+    }
+    /* process run line arguments */
+    while (--argc > 0) { /* collect run line arguments */
+        if(**++argv == '-') for(arg = *argv;;) {
+            switch(*++arg) {
+              case '\0': /* position of "stdin" */
+                if (arg[-1] == '-') eargv[eargc++] = "-";
+                break;
+              case 'v': /* monitor dump of initialization */
+                mon = 1;
+                continue;
+              case 'I': /* alt. method to spec inverse */
+                inverse = 1;
+                continue;
+              case 'E': /* echo ascii input to ascii output */
+                echoin = 1;
+                continue;
+              case 't': /* set col. one char */
+                if (arg[1]) tag = *++arg;
+                else emess(1,"missing -t col. 1 tag");
+                continue;
+              case 'l': /* list projections, ellipses or units */
+                if (!arg[1] || arg[1] == 'p' || arg[1] == 'P') {
+                    /* list projections */
+                    struct PJ_LIST *lp;
+                    int do_long = arg[1] == 'P', c;
+                    char *str;
+
+                    for (lp = pj_list ; lp->id ; ++lp) {
+                        (void)printf("%s : ", lp->id);
+                        if (do_long)  /* possibly multiline description */
+                            (void)puts(*lp->descr);
+                        else { /* first line, only */
+                            str = *lp->descr;
+                            while ((c = *str++) && c != '\n')
+                                putchar(c);
+                            putchar('\n');
+                        }
+                    }
+                } else if (arg[1] == '=') { /* list projection 'descr' */
+                    struct PJ_LIST *lp;
+
+                    arg += 2;
+                    for (lp = pj_list ; lp->id ; ++lp)
+                        if (!strcmp(lp->id, arg)) {
+                            (void)printf("%9s : %s\n", lp->id, *lp->descr);
+                            break;
+                        }
+                } else if (arg[1] == 'e') { /* list ellipses */
+                    struct PJ_ELLPS *le;
+
+                    for (le = pj_ellps; le->id ; ++le)
+                        (void)printf("%9s %-16s %-16s %s\n",
+                                     le->id, le->major, le->ell, le->name);
+                } else if (arg[1] == 'u') { /* list units */
+                    struct PJ_UNITS *lu;
+
+                    for (lu = pj_units; lu->id ; ++lu)
+                        (void)printf("%12s %-20s %s\n",
+                                     lu->id, lu->to_meter, lu->name);
+                } else if (arg[1] == 'd') { /* list datums */
+                    struct PJ_DATUMS *ld;
+
+                    printf("__datum_id__ __ellipse___ __definition/comments______________________________\n" );
+                    for (ld = pj_datums; ld->id ; ++ld)
+                    {
+                        printf("%12s %-12s %-30s\n",
+                               ld->id, ld->ellipse_id, ld->defn);
+                        if( ld->comments != NULL && strlen(ld->comments) > 0 )
+                            printf( "%25s %s\n", " ", ld->comments );
+                    }
+                } else
+                    emess(1,"invalid list option: l%c",arg[1]);
+                exit(0);
+                continue; /* artificial */
+              case 'e': /* error line alternative */
+                if (--argc <= 0)
+                    noargument:			   
+                emess(1,"missing argument for -%c",*arg);
+                oterr = *++argv;
+                continue;
+              case 'W': /* specify seconds precision */
+              case 'w': /* -W for constant field width */
+                if ((c = arg[1]) != 0 && isdigit(c)) {
+                    set_rtodms(c - '0', *arg == 'W');
+                    ++arg;
+                } else
+                    emess(1,"-W argument missing or non-digit");
+                continue;
+              case 'f': /* alternate output format degrees or xy */
+                if (--argc <= 0) goto noargument;
+                oform = *++argv;
+                continue;
+              case 'r': /* reverse input */
+                reversein = 1;
+                continue;
+              case 's': /* reverse output */
+                reverseout = 1;
+                continue;
+              default:
+                emess(1, "invalid option: -%c",*arg);
+                break;
+            }
+            break;
+
+        } else if (strcmp(*argv,"+to") == 0 ) {
+            have_to_flag = 1;
+
+        } else if (**argv == '+') { /* + argument */
+            if( have_to_flag )
+            {
+                if( to_argc < MAX_PARGS )
+                    to_argv[to_argc++] = *argv + 1;
+                else
+                    emess(1,"overflowed + argument table");
+            }
+            else 
+            {
+                if (from_argc < MAX_PARGS)
+                    from_argv[from_argc++] = *argv + 1;
+                else
+                    emess(1,"overflowed + argument table");
+            }
+        } else /* assumed to be input file name(s) */
+            eargv[eargc++] = *argv;
+    }
+    if (eargc == 0 ) /* if no specific files force sysin */
+        eargv[eargc++] = "-";
+
+    /* 
+     * If no to projection is provided, generate a latlong equivelent to
+     * the from projection.
+     */
+    if( to_argc == 0 )
+    {
+        to_argv[to_argc++] = "proj=latlong";
+        for( i = 0; i < from_argc; i++ )
+        {
+            if( strncmp(from_argv[i],"a=",2) == 0 
+                || strncmp(from_argv[i],"b=",2) == 0 
+                || strncmp(from_argv[i],"ra=",3) == 0 
+                || strncmp(from_argv[i],"ellps=",6) == 0 
+                || strncmp(from_argv[i],"datum=",6) == 0 )
+                to_argv[to_argc++] = from_argv[i];
+        }
+    }
+
+    /* 
+     * If the user has requested inverse, then just reverse the
+     * coordinate systems.
+     */
+    if( inverse )
+    {
+        int     argcount;
+        
+        for( i = 0; i < MAX_PARGS; i++ )
+        {
+            char *arg;
+
+            arg = from_argv[i];
+            from_argv[i] = to_argv[i];
+            to_argv[i] = arg;
+        }
+
+        argcount = from_argc;
+        from_argc = to_argc;
+        to_argc = argcount;
+    }
+
+    if (!(fromProj = pj_init(from_argc, from_argv)))
+    {
+        printf( "Using from definition: " );
+        for( i = 0; i < from_argc; i++ )
+            printf( "%s ", from_argv[i] );
+        printf( "\n" );
+
+        emess(3,"projection initialization failure\ncause: %s",
+              pj_strerrno(pj_errno));
+    }
+
+    if (!(toProj = pj_init(to_argc, to_argv)))
+    {
+        printf( "Using to definition: " );
+        for( i = 0; i < to_argc; i++ )
+            printf( "%s ", to_argv[i] );
+        printf( "\n" );
+
+        emess(3,"projection initialization failure\ncause: %s",
+              pj_strerrno(pj_errno));
+    }
+
+    if (mon) {
+        printf( "%c ---- From Coordinate System ----\n", tag );
+        pj_pr_list(fromProj);
+        printf( "%c ---- To Coordinate System ----\n", tag );
+        pj_pr_list(toProj);
+    }
+
+    /* set input formating control */
+    if( !fromProj->is_latlong )
+        informat = strtod;
+    else {
+        informat = dmstor;
+    }
+
+    if( !toProj->is_latlong && !oform )
+        oform = "%.2f";
+
+    /* process input file list */
+    for ( ; eargc-- ; ++eargv) {
+        if (**eargv == '-') {
+            fid = stdin;
+            emess_dat.File_name = "<stdin>";
+
+        } else {
+            if ((fid = fopen(*eargv, "rt")) == NULL) {
+                emess(-2, *eargv, "input file");
+                continue;
+            }
+            emess_dat.File_name = *eargv;
+        }
+        emess_dat.File_line = 0;
+        process(fid);
+        fclose(fid);
+        emess_dat.File_name = 0;
+    }
+    exit(0); /* normal completion */
+}
diff --git a/src/geocent.c b/src/geocent.c
new file mode 100644
index 00000000..e4fbad94
--- /dev/null
+++ b/src/geocent.c
@@ -0,0 +1,319 @@
+/***************************************************************************/
+/* RSC IDENTIFIER:  GEOCENTRIC
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between Geodetic coordinates (latitude,
+ *    longitude in radians and height in meters) and Geocentric coordinates
+ *    (X, Y, Z) in meters.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found, the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *      GEOCENT_NO_ERROR        : No errors occurred in function
+ *      GEOCENT_LAT_ERROR       : Latitude out of valid range
+ *                                 (-90 to 90 degrees)
+ *      GEOCENT_LON_ERROR       : Longitude out of valid range
+ *                                 (-180 to 360 degrees)
+ *      GEOCENT_A_ERROR         : Semi-major axis lessthan or equal to zero
+ *      GEOCENT_B_ERROR         : Semi-minor axis lessthan or equal to zero
+ *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
+ *
+ *
+ * REUSE NOTES
+ *
+ *    GEOCENTRIC is intended for reuse by any application that performs
+ *    coordinate conversions between geodetic coordinates and geocentric
+ *    coordinates.
+ *    
+ *
+ * REFERENCES
+ *    
+ *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
+ *    Ralph Toms, February 1996  UCRL-JC-123138.
+ *    
+ *    Further information on GEOCENTRIC can be found in the Reuse Manual.
+ *
+ *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
+ *                                 Geospatial Information Division
+ *                                 7701 Telegraph Road
+ *                                 Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    GEOCENTRIC has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    GEOCENTRIC was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC version 2.8.1
+ *    2. Windows 95 with MS Visual C++ version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    25-02-97          Original Code
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                               INCLUDES
+ */
+#include <math.h>
+#include "geocent.h"
+/*
+ *    math.h     - is needed for calls to sin, cos, tan and sqrt.
+ *    geocent.h  - is needed for Error codes and prototype error checking.
+ */
+
+
+/***************************************************************************/
+/*
+ *                               DEFINES
+ */
+#define PI         3.14159265358979323e0
+#define PI_OVER_2  (PI / 2.0e0)
+#define FALSE      0
+#define TRUE       1
+#define COS_67P5   0.38268343236508977  /* cosine of 67.5 degrees */
+#define AD_C       1.0026000            /* Toms region 1 constant */
+
+
+/***************************************************************************/
+/*
+ *                              GLOBAL DECLARATIONS
+ */
+/* Ellipsoid parameters, default to WGS 84 */
+double Geocent_a = 6378137.0;     /* Semi-major axis of ellipsoid in meters */
+double Geocent_b = 6356752.3142;  /* Semi-minor axis of ellipsoid           */
+
+double Geocent_a2 = 40680631590769.0;        /* Square of semi-major axis */
+double Geocent_b2 = 40408299984087.05;       /* Square of semi-minor axis */
+double Geocent_e2 = 0.0066943799901413800;   /* Eccentricity squared  */
+double Geocent_ep2 = 0.00673949675658690300; /* 2nd eccentricity squared */
+/*
+ * These state variables are for optimization purposes.  The only function
+ * that should modify them is Set_Geocentric_Parameters.
+ */
+
+
+/***************************************************************************/
+/*
+ *                              FUNCTIONS     
+ */
+
+
+long Set_Geocentric_Parameters (double a, 
+                                double b) 
+{ /* BEGIN Set_Geocentric_Parameters */
+/*
+ * The function Set_Geocentric_Parameters receives the ellipsoid parameters
+ * as inputs and sets the corresponding state variables.
+ *
+ *    a  : Semi-major axis, in meters.          (input)
+ *    b  : Semi-minor axis, in meters.          (input)
+ */
+  long Error_Code = GEOCENT_NO_ERROR;
+
+  if (a <= 0.0)
+    Error_Code |= GEOCENT_A_ERROR;
+  if (b <= 0.0)
+    Error_Code |= GEOCENT_B_ERROR;
+  if (a < b)
+    Error_Code |= GEOCENT_A_LESS_B_ERROR;
+  if (!Error_Code)
+  {
+    Geocent_a = a;
+    Geocent_b = b;
+    Geocent_a2 = a * a;
+    Geocent_b2 = b * b;
+    Geocent_e2 = (Geocent_a2 - Geocent_b2) / Geocent_a2;
+    Geocent_ep2 = (Geocent_a2 - Geocent_b2) / Geocent_b2;
+  }
+  return (Error_Code);
+} /* END OF Set_Geocentric_Parameters */
+
+
+void Get_Geocentric_Parameters (double *a, 
+                                double *b)
+{ /* BEGIN Get_Geocentric_Parameters */
+/*
+ * The function Get_Geocentric_Parameters returns the ellipsoid parameters
+ * to be used in geocentric coordinate conversions.
+ *
+ *    a  : Semi-major axis, in meters.          (output)
+ *    b  : Semi-minor axis, in meters.          (output)
+ */
+
+  *a = Geocent_a;
+  *b = Geocent_b;
+} /* END OF Get_Geocentric_Parameters */
+
+
+long Convert_Geodetic_To_Geocentric (double Latitude,
+                                     double Longitude,
+                                     double Height,
+                                     double *X,
+                                     double *Y,
+                                     double *Z) 
+{ /* BEGIN Convert_Geodetic_To_Geocentric */
+/*
+ * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
+ * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
+ * according to the current ellipsoid parameters.
+ *
+ *    Latitude  : Geodetic latitude in radians                     (input)
+ *    Longitude : Geodetic longitude in radians                    (input)
+ *    Height    : Geodetic height, in meters                       (input)
+ *    X         : Calculated Geocentric X coordinate, in meters    (output)
+ *    Y         : Calculated Geocentric Y coordinate, in meters    (output)
+ *    Z         : Calculated Geocentric Z coordinate, in meters    (output)
+ *
+ */
+  long Error_Code = GEOCENT_NO_ERROR;
+  double Rn;            /*  Earth radius at location  */
+  double Sin_Lat;       /*  sin(Latitude)  */
+  double Sin2_Lat;      /*  Square of sin(Latitude)  */
+  double Cos_Lat;       /*  cos(Latitude)  */
+
+  if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
+  { /* Latitude out of range */
+    Error_Code |= GEOCENT_LAT_ERROR;
+  }
+  if ((Longitude < -PI) || (Longitude > (2*PI)))
+  { /* Longitude out of range */
+    Error_Code |= GEOCENT_LON_ERROR;
+  }
+  if (!Error_Code)
+  { /* no errors */
+    if (Longitude > PI)
+      Longitude -= (2*PI);
+    Sin_Lat = sin(Latitude);
+    Cos_Lat = cos(Latitude);
+    Sin2_Lat = Sin_Lat * Sin_Lat;
+    Rn = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * Sin2_Lat));
+    *X = (Rn + Height) * Cos_Lat * cos(Longitude);
+    *Y = (Rn + Height) * Cos_Lat * sin(Longitude);
+    *Z = ((Rn * (1 - Geocent_e2)) + Height) * Sin_Lat;
+
+  }
+  return (Error_Code);
+} /* END OF Convert_Geodetic_To_Geocentric */
+
+void Convert_Geocentric_To_Geodetic (double X,
+                                     double Y, 
+                                     double Z,
+                                     double *Latitude,
+                                     double *Longitude,
+                                     double *Height)
+{ /* BEGIN Convert_Geocentric_To_Geodetic */
+/*
+ * The function Convert_Geocentric_To_Geodetic converts geocentric
+ * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, 
+ * and height), according to the current ellipsoid parameters.
+ *
+ *    X         : Geocentric X coordinate, in meters.         (input)
+ *    Y         : Geocentric Y coordinate, in meters.         (input)
+ *    Z         : Geocentric Z coordinate, in meters.         (input)
+ *    Latitude  : Calculated latitude value in radians.       (output)
+ *    Longitude : Calculated longitude value in radians.      (output)
+ *    Height    : Calculated height value, in meters.         (output)
+ *
+ * The method used here is derived from 'An Improved Algorithm for
+ * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
+ */
+
+/* Note: Variable names follow the notation used in Toms, Feb 1996 */
+
+  double W;        /* distance from Z axis */
+  double W2;       /* square of distance from Z axis */
+  double T0;       /* initial estimate of vertical component */
+  double T1;       /* corrected estimate of vertical component */
+  double S0;       /* initial estimate of horizontal component */
+  double S1;       /* corrected estimate of horizontal component */
+  double Sin_B0;   /* sin(B0), B0 is estimate of Bowring aux variable */
+  double Sin3_B0;  /* cube of sin(B0) */
+  double Cos_B0;   /* cos(B0) */
+  double Sin_p1;   /* sin(phi1), phi1 is estimated latitude */
+  double Cos_p1;   /* cos(phi1) */
+  double Rn;       /* Earth radius at location */
+  double Sum;      /* numerator of cos(phi1) */
+  int At_Pole;     /* indicates location is in polar region */
+
+  At_Pole = FALSE;
+  if (X != 0.0)
+  {
+    *Longitude = atan2(Y,X);
+  }
+  else
+  {
+    if (Y > 0)
+    {
+      *Longitude = PI_OVER_2;
+    }
+    else if (Y < 0)
+    {
+      *Longitude = -PI_OVER_2;
+    }
+    else
+    {
+      At_Pole = TRUE;
+      *Longitude = 0.0;
+      if (Z > 0.0)
+      {  /* north pole */
+        *Latitude = PI_OVER_2;
+      }
+      else if (Z < 0.0)
+      {  /* south pole */
+        *Latitude = -PI_OVER_2;
+      }
+      else
+      {  /* center of earth */
+        *Latitude = PI_OVER_2;
+        *Height = -Geocent_b;
+        return;
+      } 
+    }
+  }
+  W2 = X*X + Y*Y;
+  W = sqrt(W2);
+  T0 = Z * AD_C;
+  S0 = sqrt(T0 * T0 + W2);
+  Sin_B0 = T0 / S0;
+  Cos_B0 = W / S0;
+  Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
+  T1 = Z + Geocent_b * Geocent_ep2 * Sin3_B0;
+  Sum = W - Geocent_a * Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0;
+  S1 = sqrt(T1*T1 + Sum * Sum);
+  Sin_p1 = T1 / S1;
+  Cos_p1 = Sum / S1;
+  Rn = Geocent_a / sqrt(1.0 - Geocent_e2 * Sin_p1 * Sin_p1);
+  if (Cos_p1 >= COS_67P5)
+  {
+    *Height = W / Cos_p1 - Rn;
+  }
+  else if (Cos_p1 <= -COS_67P5)
+  {
+    *Height = W / -Cos_p1 - Rn;
+  }
+  else
+  {
+    *Height = Z / Sin_p1 + Rn * (Geocent_e2 - 1.0);
+  }
+  if (At_Pole == FALSE)
+  {
+    *Latitude = atan(Sin_p1 / Cos_p1);
+  }
+} /* END OF Convert_Geocentric_To_Geodetic */
diff --git a/src/geocent.h b/src/geocent.h
new file mode 100644
index 00000000..58fca48e
--- /dev/null
+++ b/src/geocent.h
@@ -0,0 +1,163 @@
+#ifndef GEOCENT_H
+#define GEOCENT_H
+
+/***************************************************************************/
+/* RSC IDENTIFIER:  GEOCENTRIC
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between Geodetic coordinates (latitude,
+ *    longitude in radians and height in meters) and Geocentric coordinates
+ *    (X, Y, Z) in meters.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found, the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *      GEOCENT_NO_ERROR        : No errors occurred in function
+ *      GEOCENT_LAT_ERROR       : Latitude out of valid range
+ *                                 (-90 to 90 degrees)
+ *      GEOCENT_LON_ERROR       : Longitude out of valid range
+ *                                 (-180 to 360 degrees)
+ *      GEOCENT_A_ERROR         : Semi-major axis less than or equal to zero
+ *      GEOCENT_B_ERROR         : Semi-minor axis less than or equal to zero
+ *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
+ *
+ *
+ * REUSE NOTES
+ *
+ *    GEOCENTRIC is intended for reuse by any application that performs
+ *    coordinate conversions between geodetic coordinates and geocentric
+ *    coordinates.
+ *    
+ *
+ * REFERENCES
+ *    
+ *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
+ *    Ralph Toms, February 1996  UCRL-JC-123138.
+ *    
+ *    Further information on GEOCENTRIC can be found in the Reuse Manual.
+ *
+ *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
+ *                                 Geospatial Information Division
+ *                                 7701 Telegraph Road
+ *                                 Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    GEOCENTRIC has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    GEOCENTRIC was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC version 2.8.1
+ *    2. Windows 95 with MS Visual C++ version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                              DEFINES
+ */
+#define GEOCENT_NO_ERROR        0x0000
+#define GEOCENT_LAT_ERROR       0x0001
+#define GEOCENT_LON_ERROR       0x0002
+#define GEOCENT_A_ERROR         0x0004
+#define GEOCENT_B_ERROR         0x0008
+#define GEOCENT_A_LESS_B_ERROR  0x0010
+
+
+/***************************************************************************/
+/*
+ *                              FUNCTION PROTOTYPES
+ */
+
+/* ensure proper linkage to c++ programs */
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+  long Set_Geocentric_Parameters (double a, 
+                                  double b);
+/*
+ * The function Set_Geocentric_Parameters receives the ellipsoid parameters
+ * as inputs and sets the corresponding state variables.
+ *
+ *    a  : Semi-major axis, in meters.          (input)
+ *    b  : Semi-minor axis, in meters.          (input)
+ */
+
+
+  void Get_Geocentric_Parameters (double *a, 
+                                  double *b);
+/*
+ * The function Get_Geocentric_Parameters returns the ellipsoid parameters
+ * to be used in geocentric coordinate conversions.
+ *
+ *    a  : Semi-major axis, in meters.          (output)
+ *    b  : Semi-minor axis, in meters.          (output)
+ */
+
+
+  long Convert_Geodetic_To_Geocentric (double Latitude,
+                                       double Longitude,
+                                       double Height,
+                                       double *X,
+                                       double *Y,
+                                       double *Z);
+/*
+ * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
+ * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
+ * according to the current ellipsoid parameters.
+ *
+ *    Latitude  : Geodetic latitude in radians                     (input)
+ *    Longitude : Geodetic longitude in radians                    (input)
+ *    Height    : Geodetic height, in meters                       (input)
+ *    X         : Calculated Geocentric X coordinate, in meters.   (output)
+ *    Y         : Calculated Geocentric Y coordinate, in meters.   (output)
+ *    Z         : Calculated Geocentric Z coordinate, in meters.   (output)
+ *
+ */
+
+
+  void Convert_Geocentric_To_Geodetic (double X,
+                                       double Y, 
+                                       double Z,
+                                       double *Latitude,
+                                       double *Longitude,
+                                       double *Height);
+/*
+ * The function Convert_Geocentric_To_Geodetic converts geocentric
+ * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, 
+ * and height), according to the current ellipsoid parameters.
+ *
+ *    X         : Geocentric X coordinate, in meters.         (input)
+ *    Y         : Geocentric Y coordinate, in meters.         (input)
+ *    Z         : Geocentric Z coordinate, in meters.         (input)
+ *    Latitude  : Calculated latitude value in radians.       (output)
+ *    Longitude : Calculated longitude value in radians.      (output)
+ *    Height    : Calculated height value, in meters.         (output)
+ */
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* GEOCENT_H */
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
new file mode 100644
index 00000000..0f7e39a6
--- /dev/null
+++ b/src/pj_apply_gridshift.c
@@ -0,0 +1,270 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Apply datum shifts based on grid shift files (normally NAD27 to
+ *           NAD83 or the reverse).  This module is responsible for keeping
+ *           a list of loaded grids, and calling with each one that is 
+ *           allowed for a given datum (expressed as the nadgrids= parameter).
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2000/07/06 23:32:27  warmerda
+ * New
+ *
+ */
+
+#define PJ_LIB__
+
+#include <projects.h>
+#include <string.h>
+#include <math.h>
+
+#define GRID_MAX 100
+
+/* used only by pj_get_grid() and pj_deallocate_grids() */
+static int grid_count = 0;
+static char **grid_names = NULL;
+static struct CTABLE **grid_list = NULL;
+
+/* used only by pj_load_nadgrids() and pj_deallocate_grids() */
+static struct CTABLE **last_nadgrids_list = NULL;
+static char          *last_nadgrids = NULL;
+
+/************************************************************************/
+/*                        pj_deallocate_grids()                         */
+/*                                                                      */
+/*      Deallocate all loaded grids.                                    */
+/************************************************************************/
+
+void pj_deallocate_grids()
+
+{
+    if( grid_count > 0 )
+    {
+        int        i;
+
+        for( i = 0; i < grid_count; i++ )
+        {
+            if( grid_list[i] != NULL )
+                nad_free( grid_list[i] );
+            pj_dalloc( grid_names[i] );
+        }
+         
+        pj_dalloc( grid_names );
+        pj_dalloc( grid_list );
+
+        grid_names = NULL;
+        grid_list = NULL;
+
+        grid_count = 0;
+    }
+
+    if( last_nadgrids != NULL )
+    {
+        pj_dalloc( last_nadgrids );
+        last_nadgrids = NULL;
+
+        pj_dalloc( last_nadgrids_list );
+        last_nadgrids_list = NULL;
+    }
+}
+
+/************************************************************************/
+/*                            pj_get_grid()                             */
+/*                                                                      */
+/*      Find the requested grid in the list, or if not present, try     */
+/*      and load it.  On failure returns NULL and sets pj_errno.        */
+/************************************************************************/
+
+static struct CTABLE *pj_get_grid( const char *name )
+
+{
+    int   i;
+
+/* -------------------------------------------------------------------- */
+/*      First look in the existing list.                                */
+/* -------------------------------------------------------------------- */
+    for( i = 0; i < grid_count; i++ )
+    {
+        if( strcmp( grid_names[i], name ) == 0 )
+        {
+            if( grid_list[i] == NULL )
+                pj_errno = -38;
+
+            return grid_list[i];
+        }
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Add entry for this file in the grid list.                       */
+/* -------------------------------------------------------------------- */
+    if( grid_count == 0 )
+    {
+        grid_names = pj_malloc(sizeof(char *) * GRID_MAX);
+        memset( grid_names, 0, sizeof(char *) * GRID_MAX );
+        grid_list = pj_malloc(sizeof(struct CTABLE *) * GRID_MAX );
+        memset( grid_list, 0, sizeof(struct CTABLE *) * GRID_MAX );
+    }
+    else if( grid_count >= GRID_MAX )
+    {
+        pj_errno = -38;
+        return NULL;
+    }
+
+    grid_count++;
+
+    grid_names[grid_count-1] = (char *) pj_malloc(strlen(name)+1);
+    strcpy( grid_names[grid_count-1], name );
+
+/* -------------------------------------------------------------------- */
+/*      Read the file.                                                  */
+/* -------------------------------------------------------------------- */
+    grid_list[grid_count-1] = nad_init( (char *) name );
+
+    return grid_list[grid_count-1];
+}
+
+/************************************************************************/
+/*                          pj_load_nadgrids()                          */
+/*                                                                      */
+/*      This functions loads the list of grids corresponding to a       */
+/*      particular nadgrids string into a list, and returns it.  The    */
+/*      list is kept around till a request is made with a different     */
+/*      string in order to cut down on the string parsing cost, and     */
+/*      the cost of building the list of tables each time.              */
+/************************************************************************/
+
+static struct CTABLE **pj_load_nadgrids( const char *nadgrids )
+
+{
+    int     nadgrids_count = 0;
+    const char *s;
+
+    pj_errno = 0;
+
+    if( last_nadgrids != NULL 
+        && strcmp(nadgrids,last_nadgrids) == 0 )
+        return last_nadgrids_list;
+
+/* -------------------------------------------------------------------- */
+/*      Free old one, if any, and make space for new list.              */
+/* -------------------------------------------------------------------- */
+    if( last_nadgrids != NULL )
+    {
+        pj_dalloc(last_nadgrids);
+    }
+    
+    last_nadgrids = (char *) pj_malloc(strlen(nadgrids)+1);
+    strcpy( last_nadgrids, nadgrids );
+
+    if( last_nadgrids_list == NULL )
+        last_nadgrids_list = (struct CTABLE **) 
+            pj_malloc(sizeof(struct CTABLE *) * GRID_MAX);
+
+/* -------------------------------------------------------------------- */
+/*      Loop processing names out of nadgrids one at a time.            */
+/* -------------------------------------------------------------------- */
+    for( s = nadgrids; *s != '\0'; s++ )
+    {
+        int   end_char;
+        char  name[128];
+
+        for( end_char = 0; 
+             s[end_char] != '\0' && s[end_char] != ','; 
+             end_char++ ) {}
+
+        if( end_char > sizeof(name) )
+        {
+            pj_errno = -38;
+            return NULL;
+        }
+        
+        strncpy( name, s, end_char );
+        name[end_char] = '\0';
+
+        s += end_char;
+        if( *s == ',' )
+            s++;
+
+        last_nadgrids_list[nadgrids_count] = pj_get_grid( name );
+        if( last_nadgrids_list[nadgrids_count] == NULL )
+            return NULL;
+
+        nadgrids_count++;
+    }
+
+    last_nadgrids_list[nadgrids_count] = NULL;
+
+    return last_nadgrids_list;
+}
+
+/************************************************************************/
+/*                         pj_apply_gridshift()                         */
+/************************************************************************/
+
+int pj_apply_gridshift( const char *nadgrids, int inverse, 
+                        long point_count, int point_offset,
+                        double *x, double *y, double *z )
+
+{
+    struct CTABLE   **tables = pj_load_nadgrids( nadgrids );
+    int  i;
+
+    if( tables == NULL )
+        return pj_errno;
+
+    for( i = 0; i < point_count; i++ )
+    {
+        long io = i * point_offset;
+        LP   input, output;
+        int  itable;
+
+        input.phi = y[io];
+        input.lam = x[io];
+
+        /* keep trying till we find a table that works */
+        for( itable = 0; tables[itable] != NULL; itable++ )
+        {
+            output = nad_cvt( input, inverse, tables[itable] );
+            if( output.lam != HUGE_VAL )
+                break;
+        }
+
+        if( output.lam == HUGE_VAL )
+        {
+            pj_errno = -38;
+            return pj_errno;
+        }
+        else
+        {
+            y[io] = output.phi;
+            x[io] = output.lam;
+        }
+    }
+
+    return 0;
+}
+
diff --git a/src/pj_datum_set.c b/src/pj_datum_set.c
new file mode 100644
index 00000000..417479a5
--- /dev/null
+++ b/src/pj_datum_set.c
@@ -0,0 +1,132 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Apply datum definition to PJ structure from initialization string.
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2000/07/06 23:32:27  warmerda
+ * New
+ *
+ */
+
+#include <projects.h>
+#include <string.h>
+
+/************************************************************************/
+/*                            pj_datum_set()                            */
+/************************************************************************/
+
+int pj_datum_set(paralist *pl, PJ *projdef)
+
+{
+    const char *name, *towgs84, *nadgrids;
+
+    projdef->datum_type = PJD_UNKNOWN;
+
+/* -------------------------------------------------------------------- */
+/*      Is there a datum definition in the parameters list?  If so,     */
+/*      add the defining values to the parameter list.  Note that       */
+/*      this will append the ellipse definition as well as the          */
+/*      towgs84= and related parameters.  It should also be pointed     */
+/*      out that the addition is permanent rather than temporary        */
+/*      like most other keyword expansion so that the ellipse           */
+/*      definition will last into the pj_ell_set() function called      */
+/*      after this one.                                                 */
+/* -------------------------------------------------------------------- */
+    if( (name = pj_param(pl,"sdatum").s) != NULL )
+    {
+        paralist *curr;
+        const char *s;
+        int i;
+
+        /* find the end of the list, so we can add to it */
+        for (curr = pl; curr && curr->next ; curr = curr->next) {}
+        
+        /* find the datum definition */
+        for (i = 0; (s = pj_datums[i].id) && strcmp(name, s) ; ++i) {}
+
+        if (!s) { pj_errno = -9; return 1; }
+
+        if( pj_datums[i].ellipse_id && strlen(pj_datums[i].ellipse_id) > 0 )
+        {
+            char	entry[100];
+            
+            strcpy( entry, "ellps=" );
+            strncat( entry, pj_datums[i].ellipse_id, 80 );
+            curr = curr->next = pj_mkparam(entry);
+        }
+        
+        if( pj_datums[i].defn && strlen(pj_datums[i].defn) > 0 )
+            curr = curr->next = pj_mkparam(pj_datums[i].defn);
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Check for nadgrids parameter.                                   */
+/* -------------------------------------------------------------------- */
+    if( (nadgrids = pj_param(pl,"snadgrids").s) != NULL )
+    {
+        /* We don't actually save the value separately.  It will continue
+           to exist int he param list for use in pj_apply_gridshift.c */
+
+        projdef->datum_type = PJD_GRIDSHIFT;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Check for towgs84 parameter.                                    */
+/* -------------------------------------------------------------------- */
+    else if( (towgs84 = pj_param(pl,"stowgs84").s) != NULL )
+    {
+        int    parm_count = 0;
+        const char *s;
+
+        memset( projdef->datum_params, 0, sizeof(double) * 7);
+
+        /* parse out the parameters */
+        s = towgs84;
+        for( s = towgs84; *s != '\0'; ) 
+        {
+            projdef->datum_params[parm_count++] = atof(s);
+            while( *s != '\0' && *s != ',' )
+                s++;
+            if( *s == ',' )
+                s++;
+        }
+
+        if( projdef->datum_params[3] != 0.0 
+            || projdef->datum_params[4] != 0.0 
+            || projdef->datum_params[5] != 0.0 
+            || projdef->datum_params[6] != 0.0 )
+            projdef->datum_type = PJD_7PARAM;
+        else 
+            projdef->datum_type = PJD_3PARAM;
+
+        /* Note that pj_init() will later switch datum_type to 
+           PJD_WGS84 if shifts are all zero, and ellipsoid is WGS84 or GRS80 */
+    }
+
+    return 0;
+}
diff --git a/src/pj_datums.c b/src/pj_datums.c
new file mode 100644
index 00000000..f8f8385e
--- /dev/null
+++ b/src/pj_datums.c
@@ -0,0 +1,57 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Built in datum list.
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2000/07/06 23:32:27  warmerda
+ * New
+ *
+ */
+
+#define PJ_DATUMS__
+
+#include <projects.h>
+
+/* 
+ * The ellipse code must match one from pj_ellps.c.  The datum id should
+ * be kept to 12 characters or less if possible.  Use the official OGC 
+ * datum name for the comments if available. 
+ */
+
+struct PJ_DATUMS pj_datums[] = {
+/* id       definition                               ellipse  comments */
+/* --       ----------                               -------  -------- */
+"WGS84",    "towgs84=0,0,0", 		             "WGS84", "",
+"GGRS87",   "towgs84=-199.87,74.79,246.62",          "GRS80", 
+				"Greek_Geodetic_Reference_System_1987",
+"NAD83",    "towgs84=0,0,0",                         "GRS80", 
+				"North_American_Datum_1983",
+"NAD27",    "nadgrids=conus",                        "clrk66", 
+				"North_American_Datum_1927",
+NULL,       NULL,                                    NULL,    NULL 
+};
diff --git a/src/pj_latlong.c b/src/pj_latlong.c
new file mode 100644
index 00000000..554691d0
--- /dev/null
+++ b/src/pj_latlong.c
@@ -0,0 +1,61 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Stub projection implementation for lat/long coordinates. We 
+ *           don't actually change the coordinates, but we want proj=latlong
+ *           to act sort of like a projection.
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2000/07/06 23:32:27  warmerda
+ * New
+ *
+ */
+
+/* very loosely based upon DMA code by Bradford W. Drew */
+#define PJ_LIB__
+#include	<projects.h>
+PROJ_HEAD(latlong, "Lat/long (Geodetic)")  "\n\t";
+
+FORWARD(forward);
+
+        xy.x = lp.lam / P->a;
+        xy.y = lp.phi / P->a;
+        return xy;
+}
+INVERSE(inverse);
+
+        lp.phi = xy.y * P->a;
+        lp.lam = xy.x * P->a;
+        return lp;
+}
+FREEUP; if (P) pj_dalloc(P); }
+ENTRY0(latlong)
+        P->is_latlong = 1;
+        P->x0 = 0.0;
+        P->y0 = 0.0;
+	P->inv = inverse; P->fwd = forward;
+ENDENTRY(P)
diff --git a/src/pj_transform.c b/src/pj_transform.c
new file mode 100644
index 00000000..9bdd4314
--- /dev/null
+++ b/src/pj_transform.c
@@ -0,0 +1,459 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Perform overall coordinate system to coordinate system 
+ *           transformations (pj_transform() function) including reprojection
+ *           and datum shifting.
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2000/07/06 23:32:27  warmerda
+ * New
+ *
+ */
+
+#include <projects.h>
+#include <string.h>
+#include <math.h>
+#include "geocent.h"
+
+#ifndef SRS_WGS84_SEMIMAJOR
+#define SRS_WGS84_SEMIMAJOR 6378137.0
+#endif
+
+#ifndef SRS_WGS84_ESQUARED
+#define SRS_WGS84_ESQUARED 0.006694379990
+#endif
+
+/************************************************************************/
+/*                            pj_transform()                            */
+/*                                                                      */
+/*      Currently this function doesn't recognise if two projections    */
+/*      are identical (to short circuit reprojection) because it is     */
+/*      difficult to compare PJ structures (since there are some        */
+/*      projection specific components).                                */
+/************************************************************************/
+
+int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
+                  double *x, double *y, double *z )
+
+{
+    long      i;
+    int       need_datum_shift;
+
+    pj_errno = 0;
+
+    if( point_offset == 0 )
+        point_offset = 1;
+
+/* -------------------------------------------------------------------- */
+/*      Transform source points to lat/long, if they aren't             */
+/*      already.                                                        */
+/* -------------------------------------------------------------------- */
+    if( !srcdefn->is_latlong )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            XY         projected_loc;
+            LP	       geodetic_loc;
+
+            projected_loc.u = x[point_offset*i];
+            projected_loc.v = y[point_offset*i];
+
+            geodetic_loc = pj_inv( projected_loc, srcdefn );
+            if( pj_errno != 0 )
+                return pj_errno;
+
+            x[point_offset*i] = geodetic_loc.u;
+            y[point_offset*i] = geodetic_loc.v;
+        }
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Convert datums if needed, and possible.                         */
+/* -------------------------------------------------------------------- */
+    if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset, 
+                            x, y, z ) != 0 )
+        return pj_errno;
+
+/* -------------------------------------------------------------------- */
+/*      Transform destination points to projection coordinates, if      */
+/*      desired.                                                        */
+/* -------------------------------------------------------------------- */
+    if( !dstdefn->is_latlong )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            XY         projected_loc;
+            LP	       geodetic_loc;
+
+            geodetic_loc.u = x[point_offset*i];
+            geodetic_loc.v = y[point_offset*i];
+
+            projected_loc = pj_fwd( geodetic_loc, dstdefn );
+            if( pj_errno != 0 )
+                return pj_errno;
+
+            x[point_offset*i] = projected_loc.u;
+            y[point_offset*i] = projected_loc.v;
+        }
+    }
+
+    return 0;
+}
+
+/************************************************************************/
+/*                     pj_geodetic_to_geocentric()                      */
+/************************************************************************/
+
+int pj_geodetic_to_geocentric( double a, double es, 
+                               long point_count, int point_offset,
+                               double *x, double *y, double *z )
+
+{
+    double b;
+    int    i;
+
+    if( es == 0.0 )
+        b = a;
+    else
+        b = a * sqrt(1-es);
+
+    if( Set_Geocentric_Parameters( a, b ) != 0 )
+    {
+        pj_errno = PJD_ERR_GEOCENTRIC;
+        return pj_errno;
+    }
+
+    for( i = 0; i < point_count; i++ )
+    {
+        long io = i * point_offset;
+
+        if( Convert_Geodetic_To_Geocentric( y[io], x[io], z[io], 
+                                            x+io, y+io, z+io ) != 0 )
+        {
+            pj_errno = PJD_ERR_GEOCENTRIC;
+            return PJD_ERR_GEOCENTRIC;
+        }
+    }
+
+    return 0;
+}
+
+/************************************************************************/
+/*                     pj_geodetic_to_geocentric()                      */
+/************************************************************************/
+
+int pj_geocentric_to_geodetic( double a, double es, 
+                               long point_count, int point_offset,
+                               double *x, double *y, double *z )
+
+{
+    double b;
+    int    i;
+
+    if( es == 0.0 )
+        b = a;
+    else
+        b = a * sqrt(1-es);
+
+    if( Set_Geocentric_Parameters( a, b ) != 0 )
+    {
+        pj_errno = PJD_ERR_GEOCENTRIC;
+        return pj_errno;
+    }
+
+    for( i = 0; i < point_count; i++ )
+    {
+        long io = i * point_offset;
+
+        Convert_Geocentric_To_Geodetic( x[io], y[io], z[io], 
+                                        y+io, x+io, z+io );
+    }
+
+    return 0;
+}
+
+/************************************************************************/
+/*                         pj_compare_datums()                          */
+/*                                                                      */
+/*      Returns TRUE if the two datums are identical, otherwise         */
+/*      FALSE.                                                          */
+/************************************************************************/
+
+int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
+
+{
+    if( srcdefn->datum_type != dstdefn->datum_type )
+    {
+        return 0;
+    }
+    else if( srcdefn->a != dstdefn->a 
+             || ABS(srcdefn->es - dstdefn->es) > 0.000000000050 )
+    {
+        /* the tolerence for es is to ensure that GRS80 and WGS84 are
+           considered identical */
+        return 0;
+    }
+    else if( srcdefn->datum_type == PJD_3PARAM )
+    {
+        return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
+                && srcdefn->datum_params[1] == dstdefn->datum_params[1]
+                && srcdefn->datum_params[2] == dstdefn->datum_params[2]);
+    }
+    else if( srcdefn->datum_type == PJD_7PARAM )
+    {
+        return (srcdefn->datum_params[0] == dstdefn->datum_params[0]
+                && srcdefn->datum_params[1] == dstdefn->datum_params[1]
+                && srcdefn->datum_params[2] == dstdefn->datum_params[2]
+                && srcdefn->datum_params[3] == dstdefn->datum_params[3]
+                && srcdefn->datum_params[4] == dstdefn->datum_params[4]
+                && srcdefn->datum_params[5] == dstdefn->datum_params[5]
+                && srcdefn->datum_params[6] == dstdefn->datum_params[6]);
+    }
+    else if( srcdefn->datum_type == PJD_GRIDSHIFT )
+    {
+        return strcmp( pj_param(srcdefn->params,"snadgrids").s,
+                       pj_param(dstdefn->params,"snadgrids").s ) == 0;
+    }
+    else
+        return 1;
+}
+
+/************************************************************************/
+/*                       pj_geocentic_to_wgs84()                        */
+/************************************************************************/
+
+int pj_geocentric_to_wgs84( PJ *defn, 
+                            long point_count, int point_offset,
+                            double *x, double *y, double *z )
+
+{
+    int       i;
+
+    pj_errno = 0;
+
+    if( defn->datum_type == PJD_3PARAM )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            long io = i * point_offset;
+            
+            x[io] = x[io] + defn->datum_params[0];
+            y[io] = y[io] + defn->datum_params[1];
+            z[io] = z[io] + defn->datum_params[2];
+        }
+    }
+    else if( defn->datum_type == PJD_7PARAM )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            long io = i * point_offset;
+            double x_out, y_out, z_out;
+            
+            x_out = x[io] + defn->datum_params[0] 
+                + y[io] * defn->datum_params[5]
+                + z[io] * defn->datum_params[4]
+                + x[io] * defn->datum_params[6];
+
+            y_out = y[io] + defn->datum_params[1]
+                + x[io] * defn->datum_params[5]
+                + z[io] * defn->datum_params[3]
+                + y[io] * defn->datum_params[6];
+
+            z_out = z[io] + defn->datum_params[2]
+                + x[io] * defn->datum_params[4]
+                + y[io] * defn->datum_params[3]
+                + z[io] * defn->datum_params[6];
+
+            x[io] = x_out;
+            y[io] = y_out;
+            z[io] = z_out;
+        }
+    }
+
+    return 0;
+}
+
+/************************************************************************/
+/*                      pj_geocentic_from_wgs84()                       */
+/************************************************************************/
+
+int pj_geocentric_from_wgs84( PJ *defn, 
+                              long point_count, int point_offset,
+                              double *x, double *y, double *z )
+
+{
+    int       i;
+
+    pj_errno = 0;
+
+    if( defn->datum_type == PJD_3PARAM )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            long io = i * point_offset;
+            
+            x[io] = x[io] - defn->datum_params[0];
+            y[io] = y[io] - defn->datum_params[1];
+            z[io] = z[io] - defn->datum_params[2];
+        }
+    }
+    else if( defn->datum_type == PJD_7PARAM )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            long io = i * point_offset;
+            double x_out, y_out, z_out;
+            
+            x_out = x[io] - defn->datum_params[0] 
+                - y[io] * defn->datum_params[5]
+                - z[io] * defn->datum_params[4]
+                - x[io] * defn->datum_params[6];
+
+            y_out = y[io] - defn->datum_params[1]
+                - x[io] * defn->datum_params[5]
+                - z[io] * defn->datum_params[3]
+                - y[io] * defn->datum_params[6];
+
+            z_out = z[io] - defn->datum_params[2]
+                - x[io] * defn->datum_params[4]
+                - y[io] * defn->datum_params[3]
+                - z[io] * defn->datum_params[6];
+
+            x[io] = x_out;
+            y[io] = y_out;
+            z[io] = z_out;
+        }
+    }
+
+    return 0;
+}
+
+/************************************************************************/
+/*                         pj_datum_transform()                         */
+/************************************************************************/
+
+int pj_datum_transform( PJ *srcdefn, PJ *dstdefn, 
+                        long point_count, int point_offset,
+                        double *x, double *y, double *z )
+
+{
+    double      src_a, src_es, dst_a, dst_es;
+
+    pj_errno = 0;
+
+/* -------------------------------------------------------------------- */
+/*      Short cut if the datums are identical.                          */
+/* -------------------------------------------------------------------- */
+    if( pj_compare_datums( srcdefn, dstdefn ) )
+        return 0;
+
+    src_a = srcdefn->a;
+    src_es = srcdefn->es;
+
+    dst_a = dstdefn->a;
+    dst_es = dstdefn->es;
+
+/* -------------------------------------------------------------------- */
+/*	If this datum requires grid shifts, then apply it to geodetic   */
+/*      coordinates.                                                    */
+/* -------------------------------------------------------------------- */
+    if( srcdefn->datum_type == PJD_GRIDSHIFT )
+    {
+        pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0, 
+                            point_count, point_offset, x, y, z );
+
+        if( pj_errno != 0 )
+            return pj_errno;
+
+        src_a = SRS_WGS84_SEMIMAJOR;
+        src_es = 0.006694379990;
+    }
+
+    if( dstdefn->datum_type == PJD_GRIDSHIFT )
+    {
+        dst_a = SRS_WGS84_SEMIMAJOR;
+        dst_es = 0.006694379990;
+    }
+        
+/* ==================================================================== */
+/*      Do we need to go through geocentric coordinates?                */
+/* ==================================================================== */
+    if( srcdefn->datum_type == PJD_3PARAM 
+        || srcdefn->datum_type == PJD_7PARAM
+        || dstdefn->datum_type == PJD_3PARAM 
+        || dstdefn->datum_type == PJD_7PARAM)
+    {
+/* -------------------------------------------------------------------- */
+/*      Convert to geocentric coordinates.                              */
+/* -------------------------------------------------------------------- */
+        pj_geodetic_to_geocentric( src_a, src_es,
+                                   point_count, point_offset, x, y, z );
+
+        if( pj_errno )
+            return pj_errno;
+
+/* -------------------------------------------------------------------- */
+/*      Convert between datums.                                         */
+/* -------------------------------------------------------------------- */
+        if( srcdefn->datum_type != PJD_UNKNOWN
+            && dstdefn->datum_type != PJD_UNKNOWN )
+        {
+            pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
+            if( pj_errno != 0 )
+                return pj_errno;
+            
+            pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
+            if( pj_errno != 0 )
+                return pj_errno;
+        }
+
+/* -------------------------------------------------------------------- */
+/*      Convert back to geodetic coordinates.                           */
+/* -------------------------------------------------------------------- */
+        pj_geocentric_to_geodetic( dst_a, dst_es,
+                                   point_count, point_offset, x, y, z );
+        
+        if( pj_errno )
+            return pj_errno;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Apply grid shift to destination if required.                    */
+/* -------------------------------------------------------------------- */
+    if( dstdefn->datum_type == PJD_GRIDSHIFT )
+    {
+        pj_apply_gridshift( pj_param(dstdefn->params,"snadgrids").s, 1,
+                            point_count, point_offset, x, y, z );
+
+        if( pj_errno != 0 )
+            return pj_errno;
+    }
+
+    return 0;
+}
+