Add support for MISR SOM projection

5 files changed, 168 insertions, 2 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 73e3bdb1..871f8478 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -39,7 +39,7 @@ libproj_la_SOURCES = \
 	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_merc.c \
+	PJ_gall.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_tmerc.c \
 	PJ_airy.c PJ_aitoff.c PJ_august.c PJ_bacon.c \
diff --git a/src/PJ_misrsom.c b/src/PJ_misrsom.c
new file mode 100644
index 00000000..82a3275e
--- /dev/null
+++ b/src/PJ_misrsom.c
@@ -0,0 +1,164 @@
+/* based upon Snyder and Linck, USGS-NMD */
+#define PROJ_PARMS__ \
+    double a2, a4, b, c1, c3; \
+    double q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
+#define PJ_LIB__
+#include <projects.h>
+PROJ_HEAD(misrsom, "Space oblique for MISR")
+        "\n\tCyl, Sph&Ell\n\tpath=";
+#define TOL 1e-7
+#define PI_HALFPI 4.71238898038468985766
+#define TWOPI_HALFPI 7.85398163397448309610
+
+static void
+seraz0(double lam, double mult, PJ *P) {
+    double sdsq, h, s, fc, sd, sq, d__1;
+
+    lam *= DEG_TO_RAD;
+    sd = sin(lam);
+    sdsq = sd * sd;
+    s = P->p22 * P->sa * cos(lam) * sqrt((1. + P->t * sdsq) / ((
+            1. + P->w * sdsq) * (1. + P->q * sdsq)));
+    d__1 = 1. + P->q * sdsq;
+    h = sqrt((1. + P->q * sdsq) / (1. + P->w * sdsq)) * ((1. +
+            P->w * sdsq) / (d__1 * d__1) - P->p22 * P->ca);
+    sq = sqrt(P->xj * P->xj + s * s);
+    P->b += fc = mult * (h * P->xj - s * s) / sq;
+    P->a2 += fc * cos(lam + lam);
+    P->a4 += fc * cos(lam * 4.);
+    fc = mult * s * (h + P->xj) / sq;
+    P->c1 += fc * cos(lam);
+    P->c3 += fc * cos(lam * 3.);
+}
+FORWARD(e_forward); /* ellipsoid */
+    int l, nn;
+    double lamt, xlam, sdsq, c, d, s, lamdp, phidp, lampp, tanph,
+                lamtp, cl, sd, sp, fac, sav, tanphi;
+
+    if (lp.phi > HALFPI)
+        lp.phi = HALFPI;
+    else if (lp.phi < -HALFPI)
+        lp.phi = -HALFPI;
+    lampp = lp.phi >= 0. ? HALFPI : PI_HALFPI;
+    tanphi = tan(lp.phi);
+    for (nn = 0;;) {
+            sav = lampp;
+            lamtp = lp.lam + P->p22 * lampp;
+            cl = cos(lamtp);
+            if (fabs(cl) < TOL)
+                lamtp -= TOL;
+            fac = lampp - sin(lampp) * (cl < 0. ? -HALFPI : HALFPI);
+            for (l = 50; l; --l) {
+                    lamt = lp.lam + P->p22 * sav;
+                    if (fabs(c = cos(lamt)) < TOL)
+                        lamt -= TOL;
+                    xlam = (P->one_es * tanphi * P->sa + sin(lamt) * P->ca) / c;
+                    lamdp = atan(xlam) + fac;
+                    if (fabs(fabs(sav) - fabs(lamdp)) < TOL)
+                        break;
+                    sav = lamdp;
+            }
+            if (!l || ++nn >= 3 || (lamdp > P->rlm && lamdp < P->rlm2))
+                    break;
+            if (lamdp <= P->rlm)
+                lampp = TWOPI_HALFPI;
+            else if (lamdp >= P->rlm2)
+                lampp = HALFPI;
+    }
+    if (l) {
+            sp = sin(lp.phi);
+            phidp = aasin(P->ctx,(P->one_es * P->ca * sp - P->sa * cos(lp.phi) *
+                    sin(lamt)) / sqrt(1. - P->es * sp * sp));
+            tanph = log(tan(FORTPI + .5 * phidp));
+            sd = sin(lamdp);
+            sdsq = sd * sd;
+            s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)
+                     / ((1. + P->w * sdsq) * (1. + P->q * sdsq)));
+            d = sqrt(P->xj * P->xj + s * s);
+            xy.x = P->b * lamdp + P->a2 * sin(2. * lamdp) + P->a4 *
+                    sin(lamdp * 4.) - tanph * s / d;
+            xy.y = P->c1 * sd + P->c3 * sin(lamdp * 3.) + tanph * P->xj / d;
+    } else
+            xy.x = xy.y = HUGE_VAL;
+    return xy;
+}
+INVERSE(e_inverse); /* ellipsoid */
+    int nn;
+    double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
+
+    lamdp = xy.x / P->b;
+    nn = 50;
+    do {
+            sav = lamdp;
+            sd = sin(lamdp);
+            sdsq = sd * sd;
+            s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)
+                     / ((1. + P->w * sdsq) * (1. + P->q * sdsq)));
+            lamdp = xy.x + xy.y * s / P->xj - P->a2 * sin(
+                    2. * lamdp) - P->a4 * sin(lamdp * 4.) - s / P->xj * (
+                    P->c1 * sin(lamdp) + P->c3 * sin(lamdp * 3.));
+            lamdp /= P->b;
+    } while (fabs(lamdp - sav) >= TOL && --nn);
+    sl = sin(lamdp);
+    fac = exp(sqrt(1. + s * s / P->xj / P->xj) * (xy.y -
+            P->c1 * sl - P->c3 * sin(lamdp * 3.)));
+    phidp = 2. * (atan(fac) - FORTPI);
+    dd = sl * sl;
+    if (fabs(cos(lamdp)) < TOL)
+        lamdp -= TOL;
+    spp = sin(phidp);
+    sppsq = spp * spp;
+    lamt = atan(((1. - sppsq * P->rone_es) * tan(lamdp) *
+            P->ca - spp * P->sa * sqrt((1. + P->q * dd) * (
+            1. - sppsq) - sppsq * P->u) / cos(lamdp)) / (1. - sppsq
+            * (1. + P->u)));
+    sl = lamt >= 0. ? 1. : -1.;
+    scl = cos(lamdp) >= 0. ? 1. : -1;
+    lamt -= HALFPI * (1. - scl) * sl;
+    lp.lam = lamt - P->p22 * lamdp;
+    if (fabs(P->sa) < TOL)
+        lp.phi = aasin(P->ctx,spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
+    else
+            lp.phi = atan((tan(lamdp) * cos(lamt) - P->ca * sin(lamt)) /
+                    (P->one_es * P->sa));
+    return lp;
+}
+FREEUP; if (P) pj_dalloc(P); }
+ENTRY0(misrsom)
+    int path;
+    double lam, alf, esc, ess;
+
+    path = pj_param(P->ctx, P->params, "ipath").i;
+    if (path <= 0 || path > 233) E_ERROR(-29);
+    P->lam0 = DEG_TO_RAD * 129.3056 - TWOPI / 233. * path;
+    alf = 98.30382 * DEG_TO_RAD;
+    P->p22 = 98.88 / 1440.0;
+
+    P->sa = sin(alf);
+    P->ca = cos(alf);
+    if (fabs(P->ca) < 1e-9)
+        P->ca = 1e-9;
+    esc = P->es * P->ca * P->ca;
+    ess = P->es * P->sa * P->sa;
+    P->w = (1. - esc) * P->rone_es;
+    P->w = P->w * P->w - 1.;
+    P->q = ess * P->rone_es;
+    P->t = ess * (2. - P->es) * P->rone_es * P->rone_es;
+    P->u = esc * P->rone_es;
+    P->xj = P->one_es * P->one_es * P->one_es;
+    P->rlm = 0;
+    P->rlm2 = P->rlm + TWOPI;
+    P->a2 = P->a4 = P->b = P->c1 = P->c3 = 0.;
+    seraz0(0., 1., P);
+    for (lam = 9.; lam <= 81.0001; lam += 18.)
+        seraz0(lam, 4., P);
+    for (lam = 18; lam <= 72.0001; lam += 18.)
+        seraz0(lam, 2., P);
+    seraz0(90., 1., P);
+    P->a2 /= 30.;
+    P->a4 /= 60.;
+    P->b /= 30.;
+    P->c1 /= 15.;
+    P->c3 /= 45.;
+    P->inv = e_inverse; P->fwd = e_forward;
+ENDENTRY(P)
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 47cd8078..85207e87 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -88,6 +88,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_lcc.c
         PJ_loxim.c
         PJ_lsat.c
+        PJ_misrsom.c
         PJ_mbt_fps.c
         PJ_mbtfpp.c
         PJ_mbtfpq.c
diff --git a/src/makefile.vc b/src/makefile.vc
index ed7ab4c4..327d2c21 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -15,7 +15,7 @@ conic = \
 
 cylinder = \
 	PJ_cass.obj PJ_cc.obj PJ_cea.obj PJ_eqc.obj \
-	PJ_gall.obj PJ_labrd.obj PJ_lsat.obj PJ_merc.obj \
+	PJ_gall.obj PJ_labrd.obj PJ_lsat.obj PJ_misrsom.obj PJ_merc.obj \
 	PJ_mill.obj PJ_ocea.obj PJ_omerc.obj PJ_somerc.obj \
 	PJ_tcc.obj PJ_tcea.obj PJ_tmerc.obj PJ_geos.obj \
 	PJ_gstmerc.obj proj_etmerc.obj 
diff --git a/src/pj_list.h b/src/pj_list.h
index 20af2eb2..338eb2dd 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -79,6 +79,7 @@ PROJ_HEAD(mbtfps, "McBryde-Thomas Flat-Polar Sinusoidal")
 PROJ_HEAD(merc, "Mercator")
 PROJ_HEAD(mil_os, "Miller Oblated Stereographic")
 PROJ_HEAD(mill, "Miller Cylindrical")
+PROJ_HEAD(misrsom, "Space oblique for MISR")
 PROJ_HEAD(moll, "Mollweide")
 PROJ_HEAD(murd1, "Murdoch I")
 PROJ_HEAD(murd2, "Murdoch II")