the Bertin 1953 projection

5 files changed, 112 insertions, 2 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 84846c1c..e8bea068 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -55,7 +55,7 @@ libproj_la_SOURCES = \
 	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 \
-	PJ_chamb.c PJ_hammer.c PJ_lagrng.c PJ_larr.c \
+	PJ_bertin1953.c PJ_chamb.c PJ_hammer.c PJ_lagrng.c PJ_larr.c \
 	PJ_lask.c PJ_latlong.c PJ_nocol.c PJ_ob_tran.c PJ_oea.c \
 	PJ_tpeqd.c PJ_vandg.c PJ_vandg2.c PJ_vandg4.c \
 	PJ_wag7.c PJ_lcca.c PJ_geos.c proj_etmerc.c \
diff --git a/src/PJ_bertin1953.c b/src/PJ_bertin1953.c
new file mode 100644
index 00000000..19aa4d53
--- /dev/null
+++ b/src/PJ_bertin1953.c
@@ -0,0 +1,108 @@
+/*
+  Created by Jacques Bertin in 1953, this projection was the go-to choice
+  of the French cartographic school when they wished to represent phenomena
+  on a global scale.
+  
+  Formula designed by Philippe Rivière, 2017.
+  https://visionscarto.net/bertin-projection-1953
+
+  Port to PROJ by Philippe Rivière, 21 September 2018
+*/
+
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj.h"
+#include "projects.h"
+
+PROJ_HEAD(bertin1953, "Bertin 1953")
+    "\n\tMisc Sph no inv.";
+
+struct pj_opaque {
+    double w;
+    double m, rm;
+    double cosDeltaPhi, sinDeltaPhi, cosDeltaGamma, sinDeltaGamma, deltaLambda;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    // Projection constants
+    double fu = 1.4, k = 12., w = 1.68;
+
+    // Aliases
+    double lambda = lp.lam, phi = lp.phi;
+    
+    // Variable
+    double d;
+
+    // Apply rotation
+	lambda += Q->deltaLambda;
+
+    double cosphi = cos(phi),
+      x = cos(lambda) * cosphi,
+      y = sin(lambda) * cosphi,
+      z = sin(phi),
+      z0 = z * Q->cosDeltaPhi + x * Q->sinDeltaPhi;
+    lambda = atan2(y * Q->cosDeltaGamma - z0 * Q->sinDeltaGamma,
+       x * Q->cosDeltaPhi - z * Q->sinDeltaPhi);
+    z0 = z0 * Q->cosDeltaGamma + y * Q->sinDeltaGamma;
+    phi = asin(z0);
+
+    lambda = adjlon(lambda);
+
+    // Adjust pre-projection
+    if (lambda + phi < -fu) {
+      double u = (lambda - phi + 1.6) * (lambda + phi + fu) / 8.;
+      lambda += u;
+      phi -= 0.8 * u * sin(phi + M_PI / 2.);
+    }
+
+    // Project with Hammer (1.68,2)
+    cosphi = cos(phi);
+    d = sqrt(2./(1. + cosphi * cos(lambda / 2.)));
+    xy.x = w * d * cosphi * sin(lambda / 2.);
+    xy.y = d * sin(phi);
+
+    // Adjust post-projection
+    d = (1. - cos(lambda * phi)) / k;
+    if (xy.y < 0.) {
+      xy.x *= 1. + d;
+    }
+    if (xy.y > 0.) {
+      xy.x *= 1. + d / 1.5 * xy.x * xy.x;
+    }
+
+    return xy;
+}
+
+
+PJ *PROJECTION(bertin1953) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+
+    // force +lon_0 = -16.5
+    double deltaLambda = P->lam0 -16.5 / 180. * M_PI;
+
+    // force +lat_0=-42
+    double deltaPhi = -42. / 180. * M_PI,
+      deltaGamma = 0.  / 180. * M_PI;
+
+    Q->deltaLambda = deltaLambda;
+    Q->cosDeltaPhi = cos(deltaPhi);
+    Q->sinDeltaPhi = sin(deltaPhi);
+    Q->cosDeltaGamma = cos(deltaGamma);
+    Q->sinDeltaGamma = sin(deltaGamma);
+
+    P->es = 0.;
+    P->fwd = s_forward;
+    // P->inv = s_inverse;
+
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 0779568a..29f3e94f 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -45,6 +45,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_august.c
         PJ_axisswap.c
         PJ_bacon.c
+        PJ_bertin1953.c
         PJ_bipc.c
         PJ_boggs.c
         PJ_bonne.c
diff --git a/src/makefile.vc b/src/makefile.vc
index cf9d878c..fc6a6861 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -21,7 +21,7 @@ cylinder = \
 	PJ_gstmerc.obj proj_etmerc.obj PJ_comill.obj
 
 misc = \
-	PJ_airy.obj PJ_aitoff.obj PJ_august.obj PJ_bacon.obj \
+	PJ_airy.obj PJ_aitoff.obj PJ_august.obj PJ_bacon.obj PJ_bertin1953.obj \
 	PJ_chamb.obj PJ_hammer.obj PJ_lagrng.obj PJ_larr.obj \
 	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 \
diff --git a/src/pj_list.h b/src/pj_list.h
index 0a8d9d1b..b4fc2306 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -14,6 +14,7 @@ PROJ_HEAD(apian, "Apian Globular I")
 PROJ_HEAD(august, "August Epicycloidal")
 PROJ_HEAD(axisswap, "Axis ordering")
 PROJ_HEAD(bacon, "Bacon Globular")
+PROJ_HEAD(bertin1953, "Bertin 1953")
 PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
 PROJ_HEAD(boggs, "Boggs Eumorphic")
 PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")