From a0e79cb553264d1ef1555eaab782cef971676075 Mon Sep 17 00:00:00 2001
From: Matthew Bloch <masiyou@gmail.com>
Date: Wed, 5 May 2021 10:49:59 -0400
Subject: [PATCH] v0.0.30

---
 CHANGELOG.md              |  5 +++-
 package.json              |  2 +-
 src/projections/cupola.js | 57 ++++++++++++++++++++++-----------------
 3 files changed, 38 insertions(+), 26 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index cdc7296..2d34532 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
+v0.0.30
+* Updates to the Cupola projection.
+
 v0.0.29
-* Added forward spherical formula for the "Cupola" projection.
+* Added forward spherical formula for the Cupola projection.
 
 v0.0.27
 * Fixed bug related to reading UTM South .prj files
diff --git a/package.json b/package.json
index baf8ca7..1c5ef53 100644
--- a/package.json
+++ b/package.json
@@ -1,6 +1,6 @@
 {
 	"name": "mproj",
-	"version": "0.0.29",
+	"version": "0.0.30",
 	"description": "A JavaScript port of the Proj.4 cartographic projections library",
 	"keywords": [
 		"projections",
diff --git a/src/projections/cupola.js b/src/projections/cupola.js
index 51d2adf..6dcbf41 100644
--- a/src/projections/cupola.js
+++ b/src/projections/cupola.js
@@ -1,40 +1,49 @@
 pj_add(pj_cupola, 'cupola', 'Cupola', '\n\tMisc., Sph., NoInv.');
 
-// Source: https://github.com/OSGeo/PROJ/issues/2706
+// Source: https://www.tandfonline.com/eprint/EE7Y8RK4GXA4ITWUTQPY/full?target=10.1080/23729333.2020.1862962
 // See also: http://www.at-a-lanta.nl/weia/cupola.html
 
 function pj_cupola(P) {
-  // parameters for Cupola
-  var c1 = 0.5253;  // part of the equator on intermediate sphere, default = 1
-  var c2 = 0.7264;  // sin of angle of polarline; default = 1
-  var c3 = 0.4188;  // height of the equator, can be negative, default = 0
-  var c4 = 22.00;   // phi of centre projection in degrees, default = 0
-  var c5 = 0.9701;  // stretch in plane, default = 1
-  var c6 = 11.023;  // central meridian 11.023 degrees, also defines border of the map
-  var c7 = 180.00;  // degrees to the right of c6, default = 180
-  var r2 = 1.61885660611815;
-  var w123 = 0.258701109423297;
-  var c4r = 0.383972435438752;
-  var pc = 0.559562341761853;
+  var de = 0.5253;  // part of the equator on intermediate sphere, default = 1
+  var dp = 0.7264;  // sin of angle of polar line; default = 1
+  var ri = 1 / Math.sqrt(de * dp);
+  // height of the equator, can be negative, default = 0
+  var he = 0.4188;
+  // phi of projection center
+  var phi0 = 22 * DEG_TO_RAD;
+
+  // NOTE: Proj applies the +lon_0 (central meridian) parameter before passing
+  // coordinates to the projection function.
+  // TODO: Use 11.023 as default central meridian
+  // // var lam0 = 11.023 * DEG_TO_RAD;
+
+  var se = 0.9701; // stretch in plane, default = 1
+  // center of projection on intermediate sphere
+  var pc = calcP(phi0);
+  var qc = calcQ(0);
   var spc = sin(pc);
   var cpc = cos(pc);
-  var c6r = DEG_TO_RAD * (c6-c7+180);  // c6 in [rad], v of center of projection
-  var qc = c1*c6r;
 
   P.es = 0;
   P.fwd = s_fwd;
 
+  function calcP(phi) {
+    return asin(dp * sin(phi) + he * sqrt(de * dp));
+  }
+
+  function calcQ(lam) {
+    return de * lam;
+  }
+
   function s_fwd(lp, xy) {
-    // p, q: radians on intermediate sphere
-    var p = asin(c2 * sin(lp.phi) + w123);
+    var p = calcP(lp.phi);
+    var q = calcQ(lp.lam);
     var sp = sin(p);
     var cp = cos(p);
-    var q = c1 * lp.lam;
-    var sqqc = sin(q-qc);
-    var cqqc = cos(q-qc);
-    // k is Snyder's constant, taken as product with r2:
-    var r2k = r2 * sqrt(2/(1+spc*sp+cpc*cp*cqqc));
-    xy.x = r2k*cp*sqqc*c5;
-    xy.y = r2k*(cpc*sp-spc*cp*cqqc)/c5;
+    var sqqc = sin(q - qc);
+    var cqqc = cos(q - qc);
+    var K = sqrt(2 / (1 + sin(pc) * sp + cpc * cp * cqqc));
+    xy.x = ri * K * cp * sqqc * se;
+    xy.y = ri * K * (cpc * sp - spc * cp * cqqc) / se;
   }
 }