Introduce geodetic-geocentric conversions ... (#669)

* Introduce geodetic-geocentric conversions, as PJ_xxx style conversion step and as API entry points

* minor improvements and minor bug squashing

9 files changed, 119 insertions, 29 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 58514592..3fdc26d7 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -49,8 +49,8 @@ libproj_la_SOURCES = \
 	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \
 	PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
 	PJ_rpoly.c PJ_sconics.c proj_rouss.c \
-	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c \
-	PJ_gall.c PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
+	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c PJ_gall.c PJ_geoc.c \
+	PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
 	PJ_mill.c PJ_ocea.c PJ_omerc.c PJ_somerc.c \
 	PJ_tcc.c PJ_tcea.c PJ_times.c PJ_tmerc.c \
 	PJ_airy.c PJ_aitoff.c PJ_august.c PJ_bacon.c \
diff --git a/src/PJ_geoc.c b/src/PJ_geoc.c
new file mode 100644
index 00000000..865b1089
--- /dev/null
+++ b/src/PJ_geoc.c
@@ -0,0 +1,56 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Conversion from geographic to geocentric latitude and back.
+ * Author:   Thomas Knudsen (2017)
+ *
+ ******************************************************************************
+ * Copyright (c) 2017, SDFE, http://www.sdfe.dk
+ * Copyright (c) 2017, Thomas Knudsen
+ *
+ * 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.
+ *****************************************************************************/
+
+#define PJ_LIB__
+#include <proj.h>
+#include <errno.h>
+#include "projects.h"
+
+PROJ_HEAD(geoc, "Geocentric Latitude");
+
+/* Geographical to geocentric */
+static PJ_COORD forward(PJ_COORD coo, PJ *P) {
+    return proj_geoc_lat (P, PJ_FWD, coo);
+}
+
+/* Geocentric to geographical */
+static PJ_COORD inverse(PJ_COORD coo, PJ *P) {
+    return proj_geoc_lat (P, PJ_INV, coo);
+}
+
+
+static PJ *CONVERSION(geoc, 1) {
+    P->inv4d = inverse;
+    P->fwd4d = forward;
+
+    P->left   =  PJ_IO_UNITS_RADIANS;
+    P->right  =  PJ_IO_UNITS_RADIANS;
+
+    P->is_latlong = 1;
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 81de451e..52abc665 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -71,6 +71,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_fahey.c
         PJ_fouc_s.c
         PJ_gall.c
+        PJ_geoc.c
         PJ_geos.c
         PJ_gins8.c
         PJ_gnom.c
diff --git a/src/makefile.vc b/src/makefile.vc
index 6bc06bd1..f2acf982 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -26,7 +26,7 @@ misc = \
 	PJ_lask.obj PJ_nocol.obj PJ_ob_tran.obj PJ_oea.obj \
 	PJ_sch.obj PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj \
     PJ_vandg4.obj PJ_wag7.obj PJ_latlong.obj PJ_krovak.obj \
-    pj_geocent.obj PJ_healpix.obj PJ_qsc.obj
+    PJ_geoc.obj pj_geocent.obj PJ_healpix.obj PJ_qsc.obj
 
 pseudo = \
 	PJ_boggs.obj PJ_collg.obj PJ_crast.obj PJ_denoy.obj \
diff --git a/src/pj_init.c b/src/pj_init.c
index 534f0827..53fd1039 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -540,7 +540,6 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
             pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
             return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
         }
-
         PIN->a_orig = PIN->a;
         PIN->es_orig = PIN->es;
         if (pj_calc_ellps_params(PIN, PIN->a, PIN->es))
diff --git a/src/pj_list.h b/src/pj_list.h
index 0720b456..093355bd 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -42,6 +42,7 @@ PROJ_HEAD(fahey, "Fahey")
 PROJ_HEAD(fouc, "Foucaut")
 PROJ_HEAD(fouc_s, "Foucaut Sinusoidal")
 PROJ_HEAD(gall, "Gall (Gall Stereographic)")
+PROJ_HEAD(geoc, "Geocentric Latitude")
 PROJ_HEAD(geocent, "Geocentric")
 PROJ_HEAD(geos, "Geostationary Satellite View")
 PROJ_HEAD(gins8, "Ginsburg VIII (TsNIIGAiK)")
diff --git a/src/proj.def b/src/proj.def
index dd4a456c..c282525c 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -133,16 +133,17 @@ EXPORTS
 
     proj_torad                     @122
     proj_todeg                     @123
-    proj_rtodms                    @124
-    proj_dmstor                    @125
+    proj_geoc_lat                  @124
+    proj_rtodms                    @125
+    proj_dmstor                    @126
 
-    proj_derivatives               @126
-    proj_factors                   @127
+    proj_derivatives               @127
+    proj_factors                   @128
 
-    proj_list_operations           @128
-    proj_list_ellps                @129
-    proj_list_units                @130
-    proj_list_prime_meridians      @131
+    proj_list_operations           @129
+    proj_list_ellps                @130
+    proj_list_units                @131
+    proj_list_prime_meridians      @132
 
-    proj_angular_input             @132
-    proj_angular_output            @133
+    proj_angular_input             @133
+    proj_angular_output            @134
diff --git a/src/proj.h b/src/proj.h
index cb689621..e2ca6bea 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -370,8 +370,7 @@ int proj_angular_output (PJ *P, enum PJ_DIRECTION dir);
 
 
 PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coord);
-
-
+int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord);
 size_t proj_trans_generic (
     PJ *P,
     PJ_DIRECTION direction,
@@ -381,7 +380,7 @@ size_t proj_trans_generic (
     double *t, size_t st, size_t nt
 );
 
-int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord);
+
 
 /* Initializers */
 PJ_COORD proj_coord (double x, double y, double z, double t);
@@ -431,6 +430,9 @@ const PJ_PRIME_MERIDIANS  *proj_list_prime_meridians(void);
 double proj_torad (double angle_in_degrees);
 double proj_todeg (double angle_in_radians);
 
+/* Geographical to geocentric latitude - another of the "simple, but useful" */
+PJ_COORD proj_geoc_lat (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
+
 double proj_dmstor(const char *is, char **rs);
 char*  proj_rtodms(char *s, double r, int pos, int neg);
 
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 6e89fb60..7196b90c 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -334,6 +334,48 @@ size_t proj_trans_generic (
 
 
 /*************************************************************************************/
+PJ_COORD proj_geoc_lat (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+/**************************************************************************************
+    Convert geographical latitude to geocentric (or the other way round if
+    direction = PJ_INV)
+
+    The conversion involves a call to the tangent function, which goes through the
+    roof at the poles, so very close (the last few micrometers) to the poles no
+    conversion takes place and the input latitude is copied directly to the output.
+
+    Fortunately, the geocentric latitude converges to the geographical at the
+    poles, so the difference is negligible.
+
+    For the spherical case, the geographical latitude equals the geocentric, and
+    consequently, the input is copied directly to the output.
+**************************************************************************************/
+    const double limit = M_HALFPI - 1e-12;
+    PJ_COORD res = coo;
+    if ((coo.lp.phi > limit) || (coo.lp.phi < -limit) || (P->es==0))
+        return res;
+    if (direction==PJ_FWD)
+        res.lp.phi = atan (P->one_es * tan (coo.lp.phi) );
+    else
+        res.lp.phi = atan (P->rone_es * tan (coo.lp.phi) );
+
+    return res;
+}
+
+double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);}
+double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);}
+
+double proj_dmstor(const char *is, char **rs) {
+    return dmstor(is, rs);
+}
+
+char*  proj_rtodms(char *s, double r, int pos, int neg) {
+    return rtodms(s, r, pos, neg);
+}
+
+
+
+
+/*************************************************************************************/
 PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
 /**************************************************************************************
     Create a new PJ object in the context ctx, using the given definition. If ctx==0,
@@ -390,7 +432,6 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
     argv = (char **) calloc (argc, sizeof (char *));
     if (0==argv)
         return pj_dealloc (args);
-
     argv[0] = args;
     for (i = 0, j = 1;  i < n;  i++) {
         if (0==args[i])
@@ -841,14 +882,3 @@ const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void) {
     return pj_get_prime_meridians_ref();
 }
 
-double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);}
-double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);}
-
-
-double proj_dmstor(const char *is, char **rs) {
-    return dmstor(is, rs);
-}
-
-char*  proj_rtodms(char *s, double r, int pos, int neg) {
-    return rtodms(s, r, pos, neg);
-}