cass: add +hyperbolic switch for variant used by EPSG:3139 'Vanua Levu 1915 / Vanua Levu Grid'

3 files changed, 35 insertions, 1 deletions

diff --git a/docs/source/operations/projections/cass.rst b/docs/source/operations/projections/cass.rst
index 95684e51..5428260a 100644
--- a/docs/source/operations/projections/cass.rst
+++ b/docs/source/operations/projections/cass.rst
@@ -65,6 +65,10 @@ Options
 
 .. include:: ../options/R.rst
 
+.. option:: +hyperbolic
+
+    Use modified form of the standard Cassini-Soldner projection known as the Hyperbolic Cassini-Soldner.
+    This is used in particular for the "Vanua Levu Grid" of the the island of Vanua Levu, Fiji (EPSG:3139)
 
 
 Mathematical definition
diff --git a/src/projections/cass.cpp b/src/projections/cass.cpp
index aeddb27c..d45532ba 100644
--- a/src/projections/cass.cpp
+++ b/src/projections/cass.cpp
@@ -20,6 +20,7 @@ namespace { // anonymous namespace
 struct cass_data {
     double *en;
     double m0;
+    bool hyperbolic;
 };
 } // anonymous namespace
 
@@ -33,7 +34,8 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward
     const double cosphi = cos (lp.phi);
     const double M = pj_mlfn (lp.phi, sinphi, cosphi, Q->en);
 
-    const double nu  = 1./sqrt(1. - P->es * sinphi*sinphi);
+    const double nu_square = 1./(1. - P->es * sinphi*sinphi);
+    const double nu  = sqrt(nu_square);
     const double tanphi = tan(lp.phi);
     const double T = tanphi * tanphi;
     const double A = lp.lam * cosphi;
@@ -44,6 +46,11 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward
         (C1 - (8. - T + 8. * C) * A2 * C2));
     xy.y = M - Q->m0 + nu * tanphi * A2 *
         (.5 + (5. - T + 6. * C) * A2 * C3);
+    if( Q->hyperbolic )
+    {
+        const double rho = nu_square * (1. - P->es) * nu;
+        xy.y -= xy.y * xy.y * xy.y / (6 * rho * nu);
+    }
 
     return xy;
 }
@@ -74,6 +81,15 @@ static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
         (.5 - (1. + 3. * T1) * D2 * C3);
     lp.lam = D * (1. + T1 * D2 *
         (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1);
+
+    if( Q->hyperbolic )
+    {
+        // EPSG guidance note 7-2 suggests a custom approximation for the
+        // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
+        // generic inversion method
+        lp = pj_generic_inverse_2d(xy, P, lp);
+    }
+
     return lp;
 }
 
@@ -120,6 +136,8 @@ PJ *PROJECTION(cass) {
         return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
 
     Q->m0 = pj_mlfn (P->phi0,  sin (P->phi0),  cos (P->phi0), Q->en);
+    if (pj_param_exists(P->params, "hyperbolic"))
+        Q->hyperbolic = true;
     P->inv = cass_e_inverse;
     P->fwd = cass_e_forward;
 
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 3b4b098b..8c0cf6d2 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -850,6 +850,18 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+-------------------------------------------------------------------------------
+# Hyperbolic variant: test point from EPSG Guidance Note 7.2
+operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \
+          +lat_0=-16.25 +lon_0=179.33333333333333 +to_meter=20.1168 \
+          +x_0=251727.9155424 +y_0=334519.953768
+-------------------------------------------------------------------------------
+
+tolerance 0.1 mm
+accept  179.99433652777776 -16.841456527777776
+expect  16015.289017555102 13369.660053668682
+roundtrip 1
+
 ===============================================================================
 # Central Conic
 # 	Sph