1 files changed, 76 insertions, 0 deletions
diff --git a/src/projections/vandg2.cpp b/src/projections/vandg2.cpp
new file mode 100644
index 00000000..61d50044
--- /dev/null
+++ b/src/projections/vandg2.cpp
@@ -0,0 +1,76 @@
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "projects.h"
+
+namespace { // anonymous namespace
+struct pj_opaque {
+    int vdg3;
+};
+} // anonymous namespace
+
+PROJ_HEAD(vandg2, "van der Grinten II") "\n\tMisc Sph, no inv";
+PROJ_HEAD(vandg3, "van der Grinten III") "\n\tMisc Sph, no inv";
+
+#define TOL    1e-10
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    double x1, at, bt, ct;
+
+    bt = fabs(M_TWO_D_PI * lp.phi);
+    if ((ct = 1. - bt * bt) < 0.)
+        ct = 0.;
+    else
+        ct = sqrt(ct);
+    if (fabs(lp.lam) < TOL) {
+        xy.x = 0.;
+        xy.y = M_PI * (lp.phi < 0. ? -bt : bt) / (1. + ct);
+    } else {
+        at = 0.5 * fabs(M_PI / lp.lam - lp.lam / M_PI);
+        if (Q->vdg3) {
+            x1 = bt / (1. + ct);
+            xy.x = M_PI * (sqrt(at * at + 1. - x1 * x1) - at);
+            xy.y = M_PI * x1;
+        } else {
+            x1 = (ct * sqrt(1. + at * at) - at * ct * ct) /
+                (1. + at * at * bt * bt);
+            xy.x = M_PI * x1;
+            xy.y = M_PI * sqrt(1. - x1 * (x1 + 2. * at) + TOL);
+        }
+        if ( lp.lam < 0.) xy.x = -xy.x;
+        if ( lp.phi < 0.) xy.y = -xy.y;
+    }
+
+    return xy;
+}
+
+
+PJ *PROJECTION(vandg2) {
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+    if (nullptr==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+
+    Q->vdg3 = 0;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+PJ *PROJECTION(vandg3) {
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+    if (nullptr==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+
+    Q->vdg3 = 1;
+    P->es = 0.;
+    P->fwd = s_forward;
+
+    return P;
+}