Generic coordinates (#431)

* Pipeline preliminaries

Introducing the PJ_OBSERVATION data type as the basis for generic geodetic
transformations.

Also introducing the elements of a new minimalistic API focused on
generic geodetic transformations. This API is documented in the new
proj.h header, and is orthogonal (non-intrusive) wrt. the existing API
from proj_api.h

Finally added a large amount of comments to the somewhat intractable
projects.h, and extended the PJ object with a number of additional
ellipsoidal parameters of general geodetic usefulness.

The PJ elements fwdobs and invobs extend fwd3d and inv3d in a homologous
way to how fwd3d and inv3d extend fwd and inv.

7 files changed, 669 insertions, 210 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 8acf2774..65944111 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,7 @@ EXTRA_PROGRAMS = multistresstest test228
 AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
 
-include_HEADERS = proj_api.h projects.h geodesic.h \
+include_HEADERS = proj.h proj_api.h projects.h geodesic.h \
 	org_proj4_Projections.h org_proj4_PJ.h
 
 EXTRA_DIST = makefile.vc proj.def bin_cs2cs.cmake \
@@ -72,7 +72,9 @@ libproj_la_SOURCES = \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
 	geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
 	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \
-	pj_strtod.c
+	pj_strtod.c \
+	\
+	pj_observation.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index a90cdac9..54e59ae7 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -36,7 +36,7 @@ endif(USE_THREAD AND NOT Threads_FOUND)
 
 
 ##############################################
-### librairie source list and include_list ###
+###  library source list and include_list  ###
 ##############################################
 SET(SRC_LIBPROJ_PJ
         nad_init.c
@@ -193,6 +193,7 @@ SET(SRC_LIBPROJ_CORE
         pj_mlfn.c
         pj_msfn.c
         pj_mutex.c
+        pj_observation.c
         pj_open_lib.c
         pj_param.c
         pj_phi2.c
@@ -217,6 +218,7 @@ SET(SRC_LIBPROJ_CORE
 set(HEADERS_LIBPROJ
         projects.h
         proj_api.h
+        proj.h
         geodesic.h
 )
 
@@ -334,7 +336,3 @@ boost_report_value(PROJ_CORE_TARGET)
 boost_report_value(PROJ_CORE_TARGET_OUTPUT_NAME)
 boost_report_value(PROJ_LIBRARY_TYPE)
 boost_report_value(PROJ_LIBRARIES)
-
-
-
-
diff --git a/src/makefile.vc b/src/makefile.vc
index 6f63bc61..f06a769a 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -58,9 +58,13 @@ support = \
 	pj_ctx.obj pj_fileapi.obj pj_log.obj pj_apply_vgridshift.obj \
 	pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj
 
+pipeline = \
+	pj_observation.obj
+
 geodesic = geodesic.obj
+
 LIBOBJ	=	$(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc) \
-	$(geodesic)
+	$(geodesic) $(pipeline)
 PROJEXE_OBJ	= proj.obj gen_cheb.obj p_series.obj emess.obj
 CS2CSEXE_OBJ	= cs2cs.obj gen_cheb.obj p_series.obj emess.obj
 GEODEXE_OBJ	= geod.obj geod_set.obj geod_interface.obj emess.obj
@@ -146,4 +150,3 @@ install: all
 	copy proj_api.h $(INSTDIR)\include
 	copy projects.h $(INSTDIR)\include
 	copy geodesic.h $(INSTDIR)\include
-
diff --git a/src/pj_observation.c b/src/pj_observation.c
new file mode 100644
index 00000000..0d1f7acf
--- /dev/null
+++ b/src/pj_observation.c
@@ -0,0 +1,114 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Implement some service routines for the PJ_OBSERVATION generic
+ *           geodetic data type
+ * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2016-06-09/2016-10-11
+ *
+ ******************************************************************************
+ * Copyright (c) 2016, Thomas Knudsen/SDFE
+ *
+ * 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.
+ *****************************************************************************/
+#include <proj.h>
+#include <projects.h>
+
+
+/* Constructor for the OBSERVATION object */
+PJ_OBSERVATION pj_observation (
+    double x, double y, double z, double t,
+    double o, double p, double k,
+    int id, unsigned int flags
+) {
+    PJ_OBSERVATION g;
+    g.coo.xyzt.x = x;
+    g.coo.xyzt.y = y;
+    g.coo.xyzt.z = z;
+    g.coo.xyzt.t = t;
+    g.anc.opk.o  = o;
+    g.anc.opk.p  = p;
+    g.anc.opk.k  = k;
+    g.id         = id;
+    g.flags      = flags;
+    return g;
+}
+
+
+PJ_OBSERVATION pj_apply (PJ *P, enum pj_direction direction, PJ_OBSERVATION obs) {
+
+    switch (direction) {
+        case 1:
+            obs.coo.xyz  =  pj_fwd3d (obs.coo.lpz, P);
+            return obs;
+        case -1:
+            obs.coo.lpz  =  pj_inv3d (obs.coo.xyz, P);
+            return obs;
+        case 0:
+            return obs;
+        default:
+            pj_ctx_set_errno (P->ctx, EINVAL);
+    }
+
+    return pj_observation (HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,0,0);
+}
+
+
+/* initial attempt, following a suggestion by Kristian Evers */
+double pj_roundtrip(PJ *P, enum pj_direction direction, int n, PJ_OBSERVATION obs) {
+    double  d;
+    PJ_OBSERVATION o;
+
+    /* multiple roundtrips not supported yet */
+    if (n > 1) {
+        pj_ctx_set_errno (P->ctx, EINVAL);
+        return HUGE_VAL;
+    }
+
+    switch (direction) {
+        case 1:
+            o.coo.xyz  =  pj_fwd3d (obs.coo.lpz, P);
+            break;
+        case -1:
+            o.coo.lpz  =  pj_inv3d (obs.coo.xyz, P);
+            break;
+        default:
+            pj_ctx_set_errno (P->ctx, EINVAL);
+            return HUGE_VAL;
+    }
+
+    o.coo.xyz = pj_fwd3d (obs.coo.lpz, P);
+    o.coo.lpz = pj_inv3d (o.coo.xyz, P);
+
+    /* Should add "both ways" here */
+    d = hypot (hypot (o.coo.v[0] - obs.coo.v[0], o.coo.v[1] - obs.coo.v[1]), o.coo.v[2] - obs.coo.v[2]);
+    return d;
+}
+
+
+int pj_show_triplet (FILE *stream, const char *banner, PJ_TRIPLET point) {
+    int i = 0;
+    if (banner)
+        i += fprintf (stream, "%s", banner);
+
+    i += fprintf(stream,  "%16.10f ", point.xyz.x);
+    i += fprintf(stream,  "%16.10f ", point.xyz.y);
+    i += fprintf(stream,  "%16.4f",   point.xyz.z);
+    if (banner)
+        i += fprintf(stream, "\n");
+    return i;
+}
diff --git a/src/proj.def b/src/proj.def
index 1f228a6c..dbf27f96 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -89,3 +89,7 @@ EXPORTS
     geod_polygon_testpoint          @87
     geod_polygon_clear              @88
     pj_run_selftests                @89
+    pj_observation                  @90
+    pj_show_triplet                 @91
+    pj_apply                        @92
+    
diff --git a/src/proj.h b/src/proj.h
new file mode 100644
index 00000000..3d8ba359
--- /dev/null
+++ b/src/proj.h
@@ -0,0 +1,208 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Revised, experimental, "bare essentials" API for PROJ.4.
+ *           Intended as the foundation for added geodetic functionality.
+ *
+ *           Introduces the OBSERVATION data type, for generic coordinate
+ *           and ancillary data handling.
+ *
+ *           Also introduces the PJ_SPATIOTEMPORAL and PJ_TRIPLET unions
+ *           making it possible to make explicit the previously used
+ *           "implicit type punning", where a XY is turned into a LP by
+ *           re#defining both as UV, behind the back of the user.
+ *
+ *           The bare essentials API presented here follows the PROJ.4
+ *           convention of sailing the coordinate to be reprojected, up on
+ *           the stack ("call by value"), and symmetrically returning the
+ *           result on the stack. Although the OBSERVATION object is 4 times
+ *           as large as the traditional XY and LP objects, timing results
+ *           have shown the overhead to be very reasonable.
+ *
+ * Author:   Thomas Knudsen, <thokn@sdfe.dk>
+ *
+ ******************************************************************************
+ * Copyright (c) 2016, Thomas Knudsen / SDFE
+ *
+ * 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.
+ *****************************************************************************/
+#ifdef _MSC_VER
+#define _USE_MATH_DEFINES
+#endif
+#include <math.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <errno.h>
+
+#ifndef PROJ_H
+#define PROJ_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+/******************************************************************************
+    proj.h is included by projects.h in order to determine the size of the
+    PJ_OBSERVATION object.
+
+    In order to stomp as little as possible on the traditional proj.4 name
+    space, proj.h is littered with inclusion guards, leaving only the minimum
+    possible implementation when included from projects.h
+
+    The PJ_OBSERVATION object is fully defined if proj.h is included alone or
+    in connection with *but before* projects.h (the latter may be needed in
+    some cases, where it is necessary to access this "bare essentials" API,
+    while still having direct access to PJ object internals)
+******************************************************************************/
+
+/* Data type for generic geodetic observations */
+struct PJ_OBSERVATION;
+typedef struct PJ_OBSERVATION PJ_OBSERVATION;
+
+/* Data type for generic geodetic 3D data */
+union PJ_TRIPLET;
+typedef union PJ_TRIPLET PJ_TRIPLET;
+
+/* Data type for generic geodetic 3D data plus epoch information */
+union PJ_SPATIOTEMPORAL;
+typedef union PJ_SPATIOTEMPORAL PJ_SPATIOTEMPORAL;
+
+#ifndef PROJECTS_H
+/* Data type for projection/transformation information */
+struct PJconsts;
+typedef struct PJconsts PJ;         /* the PJ object herself */
+
+/* Omega, Phi, Kappa: Rotations */
+typedef struct {double o, p, k;}  PJ_OPK;
+
+/* Easting, Northing, and some kind of height (orthometric or ellipsoidal) */
+typedef struct {double e, n, h;}  PJ_ENH;
+
+/* Geodetic spatiotemporal coordinate types */
+typedef struct { double   x,   y,  z, t; }  PJ_XYZT;
+typedef struct { double   e,   n,  h, t; }  PJ_ENHT;
+typedef struct { double   u,   v,  w, t; }  PJ_UVWT;
+typedef struct { double lam, phi,  z, t; }  PJ_LPZT;
+
+/* Classic proj.4 pair/triplet types */
+typedef struct { double   u,   v; }  UV;
+typedef struct { double   x,   y; }  XY;
+typedef struct { double lam, phi; }  LP;
+
+typedef struct { double   x,   y,  z; }  XYZ;
+typedef struct { double   u,   v,  w; }  UVW;
+typedef struct { double lam, phi,  z; }  LPZ;
+
+/* Ancillary pairs and triplets for geodetic computations */
+
+/* Degrees, minutes, and seconds */
+typedef struct { double d, m,  s; }  PJ_DMS;
+
+/* Geoid undulation (N) and deflections of the vertical (z, e) */
+typedef struct { double   N,   z,  e; }  PJ_NZE;
+
+/* Ellipsoidal parameters */
+typedef struct { double   a,   f; }  PJ_AF;
+#endif  /* ndef PROJECTS_H */
+
+union PJ_SPATIOTEMPORAL {
+#ifndef PROJECTS_H
+    PJ_XYZT xyzt;
+    PJ_UVWT uvwt;
+    PJ_ENHT enht;
+    PJ_LPZT lpzt;
+    PJ_ENH  enh;
+#endif
+    double v[4]; /* Who cares - it's just a vector! */
+    XYZ  xyz;
+    UVW  uvw;
+    LPZ  lpz;
+    XY   xy;
+    UV   uv;
+    LP   lp;
+};
+
+
+
+/* Avoid preprocessor renaming and implicit type-punning: Use a union to make it explicit */
+union PJ_TRIPLET {
+#ifndef PROJECTS_H
+    PJ_OPK  opk;
+    PJ_ENH  enh;
+    PJ_NZE  nze;
+    PJ_AF   af;
+#endif
+    double v[3]; /* Who cares - it's just a vector! */
+    XYZ    xyz;
+    LPZ    lpz;
+    UVW    uvw;
+    XY     xy;
+    LP     lp;
+    UV     uv;
+};
+
+struct PJ_OBSERVATION {
+    PJ_SPATIOTEMPORAL coo;  /* coordinate data */
+    PJ_TRIPLET anc;         /* ancillary data */
+    int id;                 /* integer ancillary data - e.g. observation number, EPSG code... */
+    unsigned int flags;     /* additional data, intended for flags */
+};
+
+
+
+#ifndef PROJECTS_H
+
+/* Direction: "+" forward, "-" reverse, 0: do nothing */
+enum pj_direction {
+    PJ_FWD   = -1,
+    PJ_IDENT =  0,
+    PJ_INV   =  1
+};
+
+/* Apply transformation to observation - in forward or inverse direction */
+PJ_OBSERVATION pj_apply (PJ *P, enum pj_direction direction, PJ_OBSERVATION obs);
+
+/* Measure internal consistency - in forward or inverse direction */
+double pj_roundtrip(PJ *P, enum pj_direction direction, int n, PJ_OBSERVATION obs);
+
+
+int pj_show_triplet (FILE *stream, const char *banner, PJ_TRIPLET point);
+
+/* Constructor for the OBSERVATION object */
+PJ_OBSERVATION pj_observation (
+    double x, double y, double z, double t,
+    double o, double p, double k,
+    int id, unsigned int flags
+);
+
+#ifndef TODEG
+#define TODEG(rad)  ((rad)*180.0/M_PI)
+#endif
+#ifndef TORAD
+#define TORAD(deg)  ((deg)*M_PI/180.0)
+#endif
+
+#endif /* ndef PROJECTS_H */
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* ndef PROJ_H */
diff --git a/src/projects.h b/src/projects.h
index 5b883133..f6ac6acf 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -37,7 +37,9 @@
 #    define _CRT_NONSTDC_NO_DEPRECATE
 #  endif
 /* enable predefined math constants M_* for MS Visual Studio workaround */
-#  define _USE_MATH_DEFINES
+#  ifndef _USE_MATH_DEFINES
+#     define _USE_MATH_DEFINES
+#  endif
 #endif
 
 /* standard inclusions */
@@ -135,199 +137,340 @@ extern double hypot(double, double);
 #define DIR_CHAR '/'
 #endif
 
-struct projFileAPI_t;
-
-/* proj thread context */
-typedef struct {
-    int     last_errno;
-    int     debug_level;
-    void    (*logger)(void *, int, const char *);
-    void    *app_data;
-    struct projFileAPI_t *fileapi;
-} projCtx_t;
-
-/* datum_type values */
-#define PJD_UNKNOWN   0
-#define PJD_3PARAM    1
-#define PJD_7PARAM    2
-#define PJD_GRIDSHIFT 3
-#define PJD_WGS84     4   /* WGS84 (or anything considered equivelent) */
-
-/* library errors */
-#define PJD_ERR_GEOCENTRIC          -45
-#define PJD_ERR_AXIS                -47
-#define PJD_ERR_GRID_AREA           -48
-#define PJD_ERR_CATALOG             -49
-
 #define USE_PROJUV
 
 typedef struct { double u, v; } projUV;
 typedef struct { double r, i; } COMPLEX;
 typedef struct { double u, v, w; } projUVW;
 
+/* If user explicitly includes proj.h, before projects.h, then avoid implicit type-punning */
+#ifndef PROJ_H
 #ifndef PJ_LIB__
 #define XY projUV
 #define LP projUV
 #define XYZ projUVW
 #define LPZ projUVW
+
+/* Yes, this is ridiculous, but a consequence of an old and bad decision about implicit type-punning through preprocessor abuse */
+typedef struct { double u, v; }        UV;
+typedef struct { double u, v, w; }     UVW;
+
 #else
-typedef struct { double x, y; }     XY;
-typedef struct { double lam, phi; } LP;
-typedef struct { double x, y, z; } XYZ;
+typedef struct { double x, y; }        XY;
+typedef struct { double x, y, z; }     XYZ;
+typedef struct { double lam, phi; }    LP;
 typedef struct { double lam, phi, z; } LPZ;
-#endif
+typedef struct { double u, v; }        UV;
+typedef struct { double u, v, w; }     UVW;
+#endif  /* ndef PJ_LIB__ */
+#endif  /* ndef PROJ_H   */
 
-typedef union { double  f; int  i; char *s; } PROJVALUE;
+
+/* Forward declarations and typedefs for stuff needed inside the PJ object */
 struct PJconsts;
+struct pj_opaque;
+struct ARG_list;
+struct FACTORS;
+struct PJ_REGION_S;
+typedef struct PJ_REGION_S  PJ_Region;
+typedef struct ARG_list paralist;   /* parameter list */
+enum pj_io_units {
+    PJ_IO_UNITS_CLASSIC = 0,   /* LEFT: Radians     RIGHT: Scaled meters */
+    PJ_IO_UNITS_METERS  = 1,   /* Meters  */
+    PJ_IO_UNITS_RADIANS = 2    /* Radians */
+};
+#ifndef PROJ_H
+typedef struct PJconsts PJ;         /* the PJ object herself */
+#endif
 
-struct PJ_LIST {
-    char             *id;                         /* projection keyword */
-    struct PJconsts  *(*proj)(struct PJconsts*);  /* projection entry point */
-    char * const     *descr;                      /* description text */
+struct PJ_REGION_S {
+    double ll_long;        /* lower left corner coordinates (radians) */
+    double ll_lat;
+    double ur_long;        /* upper right corner coordinates (radians) */
+    double ur_lat;
 };
 
-/* Merging this into the PJ_LIST infrastructure is tempting, but may imply ABI breakage. Perhaps at next major version? */
+
+#include "proj.h"  /* Need this for sizeof(PJ_OBSERVATION) */
+struct projCtx_t;
+typedef struct projCtx_t projCtx_t;
+
+
+/* base projection data structure */
+struct PJconsts {
+
+    /*************************************************************************************
+
+                            G E N E R A L   C O N T E X T
+
+    **************************************************************************************
+
+        TODO: Need some description here - especially about the thread context...
+
+    **************************************************************************************/
+
+    projCtx_t *ctx;
+    const char *descr;             /* From pj_list.h or individual PJ_*.c file */
+    paralist *params;              /* Parameter list */
+    struct pj_opaque *opaque;      /* Projection specific parameters, Defined in PJ_*.c */
+
+
+    /*************************************************************************************
+
+                          F U N C T I O N    P O I N T E R S
+
+    **************************************************************************************
+
+        For projection xxx, these are pointers to functions in the corresponding
+        PJ_xxx.c file.
+
+        pj_init() delegates the setup of these to pj_projection_specific_setup_xxx(),
+        a name which is currently hidden behind the magic curtain of the PROJECTION
+        macro.
+
+    **************************************************************************************/
+
+    XY  (*fwd)(LP,    PJ *);
+    LP  (*inv)(XY,    PJ *);
+    XYZ (*fwd3d)(LPZ, PJ *);
+    LPZ (*inv3d)(XYZ, PJ *);
+    PJ_OBSERVATION (*fwdobs)(PJ_OBSERVATION, PJ *);
+    PJ_OBSERVATION (*invobs)(PJ_OBSERVATION, PJ *);
+
+    void (*spc)(LP, PJ *, struct FACTORS *);
+
+    void (*pfree)(PJ *);
+
+
+    /*************************************************************************************
+
+                          E L L I P S O I D     P A R A M E T E R S
+
+    **************************************************************************************
+
+        Despite YAGNI, we add a large number of ellipsoidal shape parameters, which
+        are not yet set up in pj_init. They are, however, inexpensive to compute,
+        compared to the overall time taken for setting up the complex PJ object
+        (cf. e.g. https://en.wikipedia.org/wiki/Angular_eccentricity).
+
+        But during single point projections it will often be a useful thing to have
+        these readily available without having to recompute at every pj_fwd / pj_inv
+        call.
+
+        With this wide selection, we should be ready for quite a number of geodetic
+        algorithms, without having to incur further ABI breakage.
+
+    **************************************************************************************/
+
+    /* The linear parameters */
+
+    double  a;                         /* semimajor axis (radius if eccentricity==0) */
+    double  b;                         /* semiminor axis */
+    double  ra;                        /* 1/a */
+    double  rb;                        /* 1/b */
+
+    /* The eccentricities */
+
+    double  e;                         /* first  eccentricity */
+    double  es;                        /* first  eccentricity squared */
+    double  e2;                        /* second eccentricity */
+    double  e2s;                       /* second eccentricity squared */
+    double  e3;                        /* third  eccentricity */
+    double  e3s;                       /* third  eccentricity squared */
+    double  one_es;                    /* 1 - e^2 */
+    double  rone_es;                   /* 1/one_es */
+
+    /* The flattenings */
+
+    double  f;                         /* first  flattening */
+    double  f2;                        /* second flattening */
+    double  n;                         /* third  flattening */
+    double  rf;                        /* 1/f  */
+    double  rf2;                       /* 1/f2 */
+    double  rn;                        /* 1/n  */
+
+    /* This one's for GRS80 */
+    double  J;                         /* "Dynamic form factor" */
+
+    double  es_orig, a_orig;           /* es and a before any +proj related adjustment */
+
+
+    /*************************************************************************************
+
+                          C O O R D I N A T E   H A N D L I N G
+
+    **************************************************************************************/
+
+    int  over;                      /* Over-range flag */
+    int  geoc;                      /* Geocentric latitude flag */
+    int  is_latlong;                /* proj=latlong ... not really a projection at all */
+    int  is_geocent;                /* proj=geocent ... not really a projection at all */
+
+    enum pj_io_units left;          /* Flags for input/output coordinate types */
+    enum pj_io_units right;
+
+
+    /*************************************************************************************
+
+                       C A R T O G R A P H I C       O F F S E T S
+
+    **************************************************************************************/
+
+    double  lam0, phi0;                /* central longitude, latitude */
+    double  x0, y0;                    /* false easting and northing */
+
+
+    /*************************************************************************************
+
+                                    S C A L I N G
+
+    **************************************************************************************/
+
+    double  k0;                        /* General scaling factor - e.g. the 0.9996 of UTM */
+    double  to_meter, fr_meter;        /* Plane coordinate scaling. Internal unit [m] */
+    double  vto_meter, vfr_meter;      /* Vertical scaling. Internal unit [m] */
+
+
+    /*************************************************************************************
+
+                  D A T U M S   A N D   H E I G H T   S Y S T E M S
+
+    **************************************************************************************/
+
+    int     datum_type;                /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
+    double  datum_params[7];           /* Parameters for 3PARAM and 7PARAM */
+    struct _pj_gi **gridlist;          /* TODO: Description needed */
+    int     gridlist_count;
+
+    int     has_geoid_vgrids;          /* TODO: Description needed */
+    struct _pj_gi **vgridlist_geoid;   /* TODO: Description needed */
+    int     vgridlist_geoid_count;
+
+    double  from_greenwich;            /* prime meridian offset (in radians) */
+    double  long_wrap_center;          /* 0.0 for -180 to 180, actually in radians*/
+    int     is_long_wrap_set;
+    char    axis[4];                   /* TODO: Description needed */
+
+    /* New Datum Shift Grid Catalogs */
+    char   *catalog_name;
+    struct _PJ_GridCatalog *catalog;
+
+    double  datum_date;                 /* TODO: Description needed */
+
+    struct _pj_gi *last_before_grid;    /* TODO: Description needed */
+    PJ_Region     last_before_region;   /* TODO: Description needed */
+    double        last_before_date;     /* TODO: Description needed */
+
+    struct _pj_gi *last_after_grid;     /* TODO: Description needed */
+    PJ_Region     last_after_region;    /* TODO: Description needed */
+    double        last_after_date;      /* TODO: Description needed */
+};
+
+
+
+
+
+
+/* Parameter list (a copy of the +proj=... etc. parameters) */
+struct ARG_list {
+    paralist *next;
+    char used;
+    char param[1];
+};
+
+
+
+typedef union { double  f; int  i; char *s; } PROJVALUE;
+
+
 struct PJ_SELFTEST_LIST {
-    char    *id;                    /* projection keyword */
-    int     (* testfunc)(void);     /* projection entry point */
+    char    *id;                 /* projection keyword */
+    int     (* testfunc)(void);  /* projection entry point */
 };
 
 struct PJ_ELLPS {
-    char    *id;        /* ellipse keyword name */
-    char    *major;     /* a= value */
-    char    *ell;       /* elliptical parameter */
-    char    *name;      /* comments */
+    char    *id;           /* ellipse keyword name */
+    char    *major;        /* a= value */
+    char    *ell;          /* elliptical parameter */
+    char    *name;         /* comments */
 };
 struct PJ_UNITS {
-    char    *id;        /* units keyword */
-    char    *to_meter;  /* multiply by value to get meters */
-    char    *name;      /* comments */
+    char    *id;           /* units keyword */
+    char    *to_meter;     /* multiply by value to get meters */
+    char    *name;         /* comments */
 };
 
 struct PJ_DATUMS {
-    char    *id;        /* datum keyword */
-    char    *defn;      /* ie. "to_wgs84=..." */
-    char    *ellipse_id;/* ie from ellipse table */
-    char    *comments;  /* EPSG code, etc */
+    char    *id;           /* datum keyword */
+    char    *defn;         /* ie. "to_wgs84=..." */
+    char    *ellipse_id;   /* ie from ellipse table */
+    char    *comments;     /* EPSG code, etc */
 };
 
 struct PJ_PRIME_MERIDIANS {
-    char    *id;        /* prime meridian keyword */
-    char    *defn;      /* offset from greenwich in DMS format. */
+    char    *id;           /* prime meridian keyword */
+    char    *defn;         /* offset from greenwich in DMS format. */
 };
 
-typedef struct {
-    double ll_long;      /* lower left corner coordinates (radians) */
-    double ll_lat;
-    double ur_long;      /* upper right corner coordinates (radians) */
-    double ur_lat;
-} PJ_Region;
 
 struct DERIVS {
-    double x_l, x_p;    /* derivatives of x for lambda-phi */
-    double y_l, y_p;    /* derivatives of y for lambda-phi */
+    double x_l, x_p;       /* derivatives of x for lambda-phi */
+    double y_l, y_p;       /* derivatives of y for lambda-phi */
 };
 
 struct FACTORS {
     struct DERIVS der;
-    double h, k;        /* meridinal, parallel scales */
-    double omega, thetap;   /* angular distortion, theta prime */
-    double conv;        /* convergence */
-    double s;           /* areal scale factor */
-    double a, b;        /* max-min scale error */
-    int code;           /* info as to analytics, see following */
+    double h, k;           /* meridional, parallel scales */
+    double omega, thetap;  /* angular distortion, theta prime */
+    double conv;           /* convergence */
+    double s;              /* areal scale factor */
+    double a, b;           /* max-min scale error */
+    int code;              /* info as to analytics, see following */
 };
 
+
 #define IS_ANAL_XL_YL 01    /* derivatives of lon analytic */
 #define IS_ANAL_XP_YP 02    /* derivatives of lat analytic */
 #define IS_ANAL_HK    04    /* h and k analytic */
 #define IS_ANAL_CONV 010    /* convergence analytic */
 
-/* parameter list struct */
-typedef struct ARG_list {
-    struct ARG_list *next;
-    char used;
-    char param[1]; } paralist;
-
-/* base projection data structure */
-#ifdef PJ_LIB__
-    /* we need this forward declaration in order to be able to add a
-       pointer to struct opaque to the typedef struct PJconsts below */
-    struct pj_opaque;
-#endif
-
-typedef struct PJconsts {
-    projCtx_t *ctx;
-    XY  (*fwd)(LP, struct PJconsts *);
-    LP  (*inv)(XY, struct PJconsts *);
-    XYZ (*fwd3d)(LPZ, struct PJconsts *);
-    LPZ (*inv3d)(XYZ, struct PJconsts *);
-    void (*spc)(LP, struct PJconsts *, struct FACTORS *);
-    void (*pfree)(struct PJconsts *);
-
-    const char *descr;
-    paralist *params;           /* parameter list */
-    int over;                   /* over-range flag */
-    int geoc;                   /* geocentric latitude flag */
-    int is_latlong;             /* proj=latlong ... not really a projection at all */
-    int is_geocent;             /* proj=geocent ... not really a projection at all */
-    double a;                   /* major axis or radius if es==0 */
-    double a_orig;              /* major axis before any +proj related adjustment */
-    double es;                  /* e ^ 2 */
-    double es_orig;             /* es before any +proj related adjustment */
-    double e;                   /* eccentricity */
-    double ra;                  /* 1/A */
-    double one_es;              /* 1 - e^2 */
-    double rone_es;             /* 1/one_es */
-    double lam0, phi0;          /* central longitude, latitude */
-    double x0, y0;              /* easting and northing */
-    double k0;                  /* general scaling factor */
-    double to_meter, fr_meter;  /* cartesian scaling */
-
-    int     datum_type;         /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
-    double  datum_params[7];
-    struct _pj_gi **gridlist;
-    int     gridlist_count;
-
-    int     has_geoid_vgrids;
-    struct _pj_gi **vgridlist_geoid;
-    int     vgridlist_geoid_count;
-    double  vto_meter, vfr_meter;
-
-    double  from_greenwich;     /* prime meridian offset (in radians) */
-    double  long_wrap_center;   /* 0.0 for -180 to 180, actually in radians*/
-    int     is_long_wrap_set;
-    char    axis[4];
-
-    /* New Datum Shift Grid Catalogs */
-    char   *catalog_name;
-    struct _PJ_GridCatalog *catalog;
-
-    double   datum_date;
-
-    struct _pj_gi *last_before_grid;
-    PJ_Region     last_before_region;
-    double        last_before_date;
+/* datum_type values */
+#define PJD_UNKNOWN   0
+#define PJD_3PARAM    1
+#define PJD_7PARAM    2
+#define PJD_GRIDSHIFT 3
+#define PJD_WGS84     4   /* WGS84 (or anything considered equivelent) */
 
-    struct _pj_gi *last_after_grid;
-    PJ_Region     last_after_region;
-    double        last_after_date;
+/* library errors */
+#define PJD_ERR_GEOCENTRIC          -45
+#define PJD_ERR_AXIS                -47
+#define PJD_ERR_GRID_AREA           -48
+#define PJD_ERR_CATALOG             -49
 
-#ifdef PJ_LIB__
-        struct pj_opaque *opaque;
-#endif
+struct projFileAPI_t;
 
-#ifdef PROJ_PARMS__
-PROJ_PARMS__
-#endif /* end of optional extensions */
-} PJ;
+/* proj thread context */
+struct projCtx_t {
+    int     last_errno;
+    int     debug_level;
+    void    (*logger)(void *, int, const char *);
+    void    *app_data;
+    struct projFileAPI_t *fileapi;
+};
 
-/* public API */
+/* classic public API */
 #include "proj_api.h"
 
 
 /* Generate pj_list external or make list from include file */
 
+struct PJ_LIST {
+    char    *id;                 /* projection keyword */
+    PJ *(*proj)(PJ *);           /* projection entry point */
+    char    * const *descr;      /* description text */
+};
+
+
 #ifndef USE_PJ_LIST_H
 extern struct PJ_LIST pj_list[];
 extern struct PJ_SELFTEST_LIST pj_selftest_list[];
@@ -348,37 +491,20 @@ extern struct PJ_DATUMS pj_datums[];
 extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[];
 #endif
 
+
+
+
+
 #ifdef PJ_LIB__
-/* repetitive projection code */
 #define PROJ_HEAD(id, name) static const char des_##id [] = name
-#define ENTRYA(name) \
-    C_NAMESPACE_VAR const char * const pj_s_##name = des_##name; \
-    C_NAMESPACE PJ *pj_##name(PJ *P) { if (!P) { \
-    if( (P = (PJ*) pj_malloc(sizeof(PJ))) != NULL) { \
-        memset( P, 0, sizeof(PJ) ); \
-    P->pfree = freeup; P->fwd = 0; P->inv = 0; \
-        P->fwd3d = 0; P->inv3d = 0; \
-    P->spc = 0; P->descr = des_##name;
-#define ENTRYX } return P; } else {
-#define ENTRY0(name) ENTRYA(name) ENTRYX
-#define ENTRY1(name, a) ENTRYA(name) P->a = 0; ENTRYX
-#define ENTRY2(name, a, b) ENTRYA(name) P->a = 0; P->b = 0; ENTRYX
-#define ENDENTRY(p) } return (p); }
+
 #define E_ERROR(err) { pj_ctx_set_errno( P->ctx, err); freeup(P); return(0); }
 #define E_ERROR_0 { freeup(P); return(0); }
 #define F_ERROR { pj_ctx_set_errno( P->ctx, -20); return(xy); }
 #define F3_ERROR { pj_ctx_set_errno( P->ctx, -20); return(xyz); }
 #define I_ERROR { pj_ctx_set_errno( P->ctx, -20); return(lp); }
 #define I3_ERROR { pj_ctx_set_errno( P->ctx, -20); return(lpz); }
-#define FORWARD(name) static XY name(LP lp, PJ *P) { XY xy = {0.0,0.0}
-#define INVERSE(name) static LP name(XY xy, PJ *P) { LP lp = {0.0,0.0}
-#define FORWARD3D(name) static XYZ name(LPZ lpz, PJ *P) {XYZ xyz = {0.0, 0.0, 0.0}
-#define INVERSE3D(name) static LPZ name(XYZ xyz, PJ *P) {LPZ lpz = {0.0, 0.0, 0.0}
-#define FREEUP static void freeup(PJ *P) {
-#define SPECIAL(name) static void name(LP lp, PJ *P, struct FACTORS *fac)
-#define ELLIPSOIDAL(P) ((P->es==0)? (FALSE): (TRUE))
-
-/* cleaned up alternative to most of the "repetitive projection code" macros */
+
 #define PROJECTION(name)                                     \
 pj_projection_specific_setup_##name (PJ *P);                 \
 C_NAMESPACE_VAR const char * const pj_s_##name = des_##name; \
@@ -393,8 +519,7 @@ C_NAMESPACE PJ *pj_##name (PJ *P) {                          \
     return P;                                                \
 }                                                            \
 PJ *pj_projection_specific_setup_##name (PJ *P)
-
-#endif
+#endif /* def PJ_LIB__ */
 
 
 int pj_generic_selftest (
@@ -468,32 +593,36 @@ typedef struct _PJ_GridCatalog {
 /* procedure prototypes */
 double dmstor(const char *, char **);
 double dmstor_ctx(projCtx ctx, const char *, char **);
-void set_rtodms(int, int);
-char *rtodms(char *, double, int, int);
+void   set_rtodms(int, int);
+char  *rtodms(char *, double, int, int);
 double adjlon(double);
 double aacos(projCtx,double), aasin(projCtx,double), asqrt(double), aatan2(double, double);
+
 PROJVALUE pj_param(projCtx ctx, paralist *, const char *);
 paralist *pj_mkparam(char *);
+
 int pj_ell_set(projCtx ctx, paralist *, double *, double *);
 int pj_datum_set(projCtx,paralist *, PJ *);
 int pj_prime_meridian_set(paralist *, PJ *);
 int pj_angular_units_set(paralist *, PJ *);
+
 void pj_prepare (PJ *P, const char *description, void (*freeup)(struct PJconsts *), size_t sizeof_struct_opaque);
 
 paralist *pj_clone_paralist( const paralist* );
-paralist*pj_search_initcache( const char *filekey );
-void pj_insert_initcache( const char *filekey, const paralist *list);
+paralist *pj_search_initcache( const char *filekey );
+void      pj_insert_initcache( const char *filekey, const paralist *list);
 
 double *pj_enfn(double);
-double pj_mlfn(double, double, double, double *);
-double pj_inv_mlfn(projCtx, double, double, double *);
-double pj_qsfn(double, double, double);
-double pj_tsfn(double, double, double);
-double pj_msfn(double, double, double);
-double pj_phi2(projCtx, double, double);
-double pj_qsfn_(double, PJ *);
+double  pj_mlfn(double, double, double, double *);
+double  pj_inv_mlfn(projCtx, double, double, double *);
+double  pj_qsfn(double, double, double);
+double  pj_tsfn(double, double, double);
+double  pj_msfn(double, double, double);
+double  pj_phi2(projCtx, double, double);
+double  pj_qsfn_(double, PJ *);
 double *pj_authset(double);
-double pj_authlat(double, double *);
+double  pj_authlat(double, double *);
+
 COMPLEX pj_zpoly1(COMPLEX, COMPLEX *, int);
 COMPLEX pj_zpolyd1(COMPLEX, COMPLEX *, int, COMPLEX *);
 
@@ -513,25 +642,26 @@ typedef struct {    /* Chebyshev or Power series structure */
     int mu, mv;     /* maximum cu and cv index (+1 for count) */
     int power;      /* != 0 if power series, else Chebyshev */
 } Tseries;
+
 Tseries *mk_cheby(projUV, projUV, double, projUV *, projUV (*)(projUV), int, int, int);
-projUV bpseval(projUV, Tseries *);
-projUV bcheval(projUV, Tseries *);
-projUV biveval(projUV, Tseries *);
-void *vector1(int, int);
-void **vector2(int, int, int);
-void freev2(void **v, int nrows);
-int bchgen(projUV, projUV, int, int, projUV **, projUV(*)(projUV));
-int bch2bps(projUV, projUV, projUV **, int, int);
+projUV   bpseval(projUV, Tseries *);
+projUV   bcheval(projUV, Tseries *);
+projUV   biveval(projUV, Tseries *);
+void    *vector1(int, int);
+void   **vector2(int, int, int);
+void     freev2(void **v, int nrows);
+int      bchgen(projUV, projUV, int, int, projUV **, projUV(*)(projUV));
+int      bch2bps(projUV, projUV, projUV **, int, int);
 
 /* nadcon related protos */
-LP nad_intr(LP, struct CTABLE *);
-LP nad_cvt(LP, int, struct CTABLE *);
+LP             nad_intr(LP, struct CTABLE *);
+LP             nad_cvt(LP, int, struct CTABLE *);
 struct CTABLE *nad_init(projCtx ctx, char *);
 struct CTABLE *nad_ctable_init( projCtx ctx, PAFile fid );
-int nad_ctable_load( projCtx ctx, struct CTABLE *, PAFile fid );
+int            nad_ctable_load( projCtx ctx, struct CTABLE *, PAFile fid );
 struct CTABLE *nad_ctable2_init( projCtx ctx, PAFile fid );
-int nad_ctable2_load( projCtx ctx, struct CTABLE *, PAFile fid );
-void nad_free(struct CTABLE *);
+int            nad_ctable2_load( projCtx ctx, struct CTABLE *, PAFile fid );
+void           nad_free(struct CTABLE *);
 
 /* higher level handling of datum grid shift files */
 
@@ -553,8 +683,8 @@ PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx, const char *, int * );
 void pj_deallocate_grids();
 
 PJ_GRIDINFO *pj_gridinfo_init( projCtx, const char * );
-int pj_gridinfo_load( projCtx, PJ_GRIDINFO * );
-void pj_gridinfo_free( projCtx, PJ_GRIDINFO * );
+int          pj_gridinfo_load( projCtx, PJ_GRIDINFO * );
+void         pj_gridinfo_free( projCtx, PJ_GRIDINFO * );
 
 PJ_GridCatalog *pj_gc_findcatalog( projCtx, const char * );
 PJ_GridCatalog *pj_gc_readcatalog( projCtx, const char * );
@@ -574,20 +704,20 @@ PJ_GRIDINFO *pj_gc_findgrid( projCtx ctx,
 
 double pj_gc_parsedate( projCtx, const char * );
 
-void *proj_mdist_ini(double);
+void  *proj_mdist_ini(double);
 double proj_mdist(double, double, double, const void *);
 double proj_inv_mdist(projCtx ctx, double, const void *);
-void *pj_gauss_ini(double, double, double *,double *);
-LP pj_gauss(projCtx, LP, const void *);
-LP pj_inv_gauss(projCtx, LP, const void *);
+void  *pj_gauss_ini(double, double, double *,double *);
+LP     pj_gauss(projCtx, LP, const void *);
+LP     pj_inv_gauss(projCtx, LP, const void *);
 
 extern char const pj_release[];
 
-struct PJ_ELLPS *pj_get_ellps_ref( void );
-struct PJ_DATUMS *pj_get_datums_ref( void );
-struct PJ_UNITS *pj_get_units_ref( void );
-struct PJ_LIST  *pj_get_list_ref( void );
-struct PJ_SELFTEST_LIST  *pj_get_selftest_list_ref ( void );
+struct PJ_ELLPS            *pj_get_ellps_ref( void );
+struct PJ_DATUMS           *pj_get_datums_ref( void );
+struct PJ_UNITS            *pj_get_units_ref( void );
+struct PJ_LIST             *pj_get_list_ref( void );
+struct PJ_SELFTEST_LIST    *pj_get_selftest_list_ref ( void );
 struct PJ_PRIME_MERIDIANS  *pj_get_prime_meridians_ref( void );
 
 double pj_atof( const char* nptr );