added goodes interrupted homolosine (#106)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1997 4e78687f-474d-0410-85f9-8d5e500ac6b2

5 files changed, 196 insertions, 5 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 1f5d0aa8..7a6b356d 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -42,7 +42,7 @@ libproj_la_SOURCES = \
 	PJ_boggs.c PJ_collg.c PJ_crast.c PJ_denoy.c \
 	PJ_eck1.c PJ_eck2.c PJ_eck3.c PJ_eck4.c \
 	PJ_eck5.c PJ_fahey.c PJ_fouc_s.c PJ_gins8.c PJ_gstmerc.c \
-	PJ_gn_sinu.c PJ_goode.c PJ_hatano.c PJ_loxim.c \
+	PJ_gn_sinu.c PJ_goode.c PJ_igh.c PJ_hatano.c PJ_loxim.c \
 	PJ_mbt_fps.c PJ_mbtfpp.c PJ_mbtfpq.c PJ_moll.c \
 	PJ_nell.c PJ_nell_h.c PJ_putp2.c PJ_putp3.c \
 	PJ_putp4p.c PJ_putp5.c PJ_putp6.c PJ_robin.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 6affea45..50369db2 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -88,7 +88,7 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
 	PJ_geos.lo proj_etmerc.lo PJ_boggs.lo PJ_collg.lo PJ_crast.lo \
 	PJ_denoy.lo PJ_eck1.lo PJ_eck2.lo PJ_eck3.lo PJ_eck4.lo \
 	PJ_eck5.lo PJ_fahey.lo PJ_fouc_s.lo PJ_gins8.lo PJ_gstmerc.lo \
-	PJ_gn_sinu.lo PJ_goode.lo PJ_hatano.lo PJ_loxim.lo \
+	PJ_gn_sinu.lo PJ_goode.lo PJ_igh.lo PJ_hatano.lo PJ_loxim.lo \
 	PJ_mbt_fps.lo PJ_mbtfpp.lo PJ_mbtfpq.lo PJ_moll.lo PJ_nell.lo \
 	PJ_nell_h.lo PJ_putp2.lo PJ_putp3.lo PJ_putp4p.lo PJ_putp5.lo \
 	PJ_putp6.lo PJ_robin.lo PJ_sts.lo PJ_urm5.lo PJ_urmfps.lo \
@@ -300,7 +300,7 @@ libproj_la_SOURCES = \
 	PJ_boggs.c PJ_collg.c PJ_crast.c PJ_denoy.c \
 	PJ_eck1.c PJ_eck2.c PJ_eck3.c PJ_eck4.c \
 	PJ_eck5.c PJ_fahey.c PJ_fouc_s.c PJ_gins8.c PJ_gstmerc.c \
-	PJ_gn_sinu.c PJ_goode.c PJ_hatano.c PJ_loxim.c \
+	PJ_gn_sinu.c PJ_goode.c PJ_igh.c PJ_hatano.c PJ_loxim.c \
 	PJ_mbt_fps.c PJ_mbtfpp.c PJ_mbtfpq.c PJ_moll.c \
 	PJ_nell.c PJ_nell_h.c PJ_putp2.c PJ_putp3.c \
 	PJ_putp4p.c PJ_putp5.c PJ_putp6.c PJ_robin.c \
@@ -504,6 +504,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_gstmerc.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_hammer.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_hatano.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_igh.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_imw_p.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_krovak.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_labrd.Plo@am__quote@
diff --git a/src/PJ_igh.c b/src/PJ_igh.c
new file mode 100644
index 00000000..40a6f89c
--- /dev/null
+++ b/src/PJ_igh.c
@@ -0,0 +1,185 @@
+#define PROJ_PARMS__ \
+        struct PJconsts* pj[12]; \
+        double dy0;
+#define PJ_LIB__
+#include	<projects.h>
+PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph.";
+	C_NAMESPACE PJ
+*pj_sinu(PJ *), *pj_moll(PJ *);
+
+static const double d4044118 = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD; // 40d 44' 11.8" [degrees]
+
+static const double d10  =  10 * DEG_TO_RAD;
+static const double d20  =  20 * DEG_TO_RAD;
+static const double d30  =  30 * DEG_TO_RAD;
+static const double d40  =  40 * DEG_TO_RAD;
+static const double d50  =  50 * DEG_TO_RAD;
+static const double d60  =  60 * DEG_TO_RAD;
+static const double d80  =  80 * DEG_TO_RAD;
+static const double d90  =  90 * DEG_TO_RAD;
+static const double d100 = 100 * DEG_TO_RAD;
+static const double d140 = 140 * DEG_TO_RAD;
+static const double d160 = 160 * DEG_TO_RAD;
+static const double d180 = 180 * DEG_TO_RAD;
+
+static const double EPSLN = 1.e-10; // allow a little 'slack' on zone edge positions
+
+FORWARD(s_forward); /* spheroid */
+        int z;
+        if (lp.phi >=  d4044118) {          // 1|2
+          z = (lp.lam <= -d40 ? 1: 2);
+        }
+        else if (lp.phi >=  0) {            // 3|4
+          z = (lp.lam <= -d40 ? 3: 4);
+        }
+        else if (lp.phi >= -d4044118) {     // 5|6|7|8
+               if (lp.lam <= -d100) z =  5; // 5
+          else if (lp.lam <=  -d20) z =  6; // 6
+          else if (lp.lam <=   d80) z =  7; // 7
+          else z = 8;                       // 8
+        }
+        else {                              // 9|10|11|12
+               if (lp.lam <= -d100) z =  9; // 9
+          else if (lp.lam <=  -d20) z = 10; // 10
+          else if (lp.lam <=   d80) z = 11; // 11
+          else z = 12;                      // 12
+        }
+
+        lp.lam -= P->pj[z-1]->lam0;
+        xy = P->pj[z-1]->fwd(lp, P->pj[z-1]);
+        xy.x += P->pj[z-1]->x0;
+        xy.y += P->pj[z-1]->y0;
+
+        return (xy);
+}
+INVERSE(s_inverse); /* spheroid */
+        const double y90 = P->dy0 + sqrt(2); // lt=90 corresponds to y=y0+sqrt(2)
+
+        int z = 0;
+        if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) // 0
+          z = 0;
+        else if (xy.y >=  d4044118)       // 1|2
+          z = (xy.x <= -d40? 1: 2);
+        else if (xy.y >=  0)              // 3|4
+          z = (xy.x <= -d40? 3: 4);
+        else if (xy.y >= -d4044118) {     // 5|6|7|8
+               if (xy.x <= -d100) z =  5; // 5
+          else if (xy.x <=  -d20) z =  6; // 6
+          else if (xy.x <=   d80) z =  7; // 7
+          else z = 8;                     // 8
+        }
+        else {                            // 9|10|11|12
+               if (xy.x <= -d100) z =  9; // 9
+          else if (xy.x <=  -d20) z = 10; // 10
+          else if (xy.x <=   d80) z = 11; // 11
+          else z = 12;                    // 12
+        }
+
+        if (z)
+        {
+          xy.x -= P->pj[z-1]->x0;
+          xy.y -= P->pj[z-1]->y0;
+          lp = P->pj[z-1]->inv(xy, P->pj[z-1]);
+          lp.lam += P->pj[z-1]->lam0;
+
+          int ok = 0;
+          switch (z) {
+            case  1: ok = (lp.lam >= -d180-EPSLN && lp.lam <=  -d40+EPSLN) ||
+                         ((lp.lam >=  -d40-EPSLN && lp.lam <=  -d10+EPSLN) &&
+                          (lp.phi >=   d60-EPSLN && lp.phi <=   d90+EPSLN)); break;
+            case  2: ok = (lp.lam >=  -d40-EPSLN && lp.lam <=  d180+EPSLN) ||
+                         ((lp.lam >= -d180-EPSLN && lp.lam <= -d160+EPSLN) &&
+                          (lp.phi >=   d50-EPSLN && lp.phi <=   d90+EPSLN)) ||
+                         ((lp.lam >=  -d50-EPSLN && lp.lam <=  -d40+EPSLN) &&
+                          (lp.phi >=   d60-EPSLN && lp.phi <=   d90+EPSLN)); break;
+            case  3: ok = (lp.lam >= -d180-EPSLN && lp.lam <=  -d40+EPSLN); break;
+            case  4: ok = (lp.lam >=  -d40-EPSLN && lp.lam <=  d180+EPSLN); break;
+            case  5: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break;
+            case  6: ok = (lp.lam >= -d100-EPSLN && lp.lam <=  -d20+EPSLN); break;
+            case  7: ok = (lp.lam >=  -d20-EPSLN && lp.lam <=   d80+EPSLN); break;
+            case  8: ok = (lp.lam >=   d80-EPSLN && lp.lam <=  d180+EPSLN); break;
+            case  9: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break;
+            case 10: ok = (lp.lam >= -d100-EPSLN && lp.lam <=  -d20+EPSLN); break;
+            case 11: ok = (lp.lam >=  -d20-EPSLN && lp.lam <=   d80+EPSLN); break;
+            case 12: ok = (lp.lam >=   d80-EPSLN && lp.lam <=  d180+EPSLN); break;
+          }
+
+          z = (!ok? 0: z); // projectable?
+        }
+     // if (!z) pj_errno = -15; // invalid x or y
+        if (!z) lp.lam = HUGE_VAL;
+        if (!z) lp.phi = HUGE_VAL;
+        return (lp);
+}
+FREEUP; int i;
+        if (P) {
+                for (i = 0; i < 12; ++i)
+                        if (P->pj[i]) free(P->pj[i]);
+                free(P);
+        }
+}
+ENTRY0(igh)
+/*
+  Zones:
+
+    -180            -40                       180
+      +--------------+-------------------------+    Zones 1,2,9,10,11 & 12:
+      |1             |2                        |      Mollweide projection
+      |              |                         |
+      +--------------+-------------------------+    Zones 3,4,5,6,7 & 8:
+      |3             |4                        |      Sinusoidal projection
+      |              |                         |
+    0 +-------+------+-+-----------+-----------+
+      |5      |6       |7          |8          |
+      |       |        |           |           |
+      +-------+--------+-----------+-----------+
+      |9      |10      |11         |12         |
+      |       |        |           |           |
+      +-------+--------+-----------+-----------+
+    -180    -100      -20         80          180
+*/
+
+#define SETUP(n, proj, x_0, y_0, lon_0) \
+    if (!(P->pj[n-1] = pj_##proj(0))) E_ERROR_0; \
+    if (!(P->pj[n-1] = pj_##proj(P->pj[n-1]))) E_ERROR_0; \
+    P->pj[n-1]->x0 = x_0; \
+    P->pj[n-1]->y0 = y_0; \
+    P->pj[n-1]->lam0 = lon_0; 
+
+        LP lp = { 0, d4044118 };
+        XY xy1;
+        XY xy3;
+
+        // sinusoidal zones
+        SETUP(3, sinu, -d100, 0, -d100);
+        SETUP(4, sinu,   d30, 0,   d30);
+        SETUP(5, sinu, -d160, 0, -d160);
+        SETUP(6, sinu,  -d60, 0,  -d60);
+        SETUP(7, sinu,   d20, 0,   d20);
+        SETUP(8, sinu,  d140, 0,  d140);
+
+        // mollweide zones
+        SETUP(1, moll, -d100, 0, -d100);
+
+        // y0 ?
+        xy1 = P->pj[0]->fwd(lp, P->pj[0]); // zone 1
+        xy3 = P->pj[2]->fwd(lp, P->pj[2]); // zone 3
+        // y0 + xy1.y = xy3.y for lt = 40d44'11.8"
+        P->dy0 = xy3.y - xy1.y;
+
+        P->pj[0]->y0 = P->dy0;
+
+        // mollweide zones (cont'd)
+        SETUP( 2, moll,   d30,  P->dy0,   d30);
+        SETUP( 9, moll, -d160, -P->dy0, -d160);
+        SETUP(10, moll,  -d60, -P->dy0,  -d60);
+        SETUP(11, moll,   d20, -P->dy0,   d20);
+        SETUP(12, moll,  d140, -P->dy0,  d140);
+
+        P->inv = s_inverse;
+        P->fwd = s_forward;
+        P->es = 0.;
+ENDENTRY(P)
+
+
+
diff --git a/src/pj_list.h b/src/pj_list.h
index 1e8df686..90b51abd 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -47,6 +47,7 @@ PROJ_HEAD(gs48, "Mod. Stererographics of 48 U.S.")
 PROJ_HEAD(gs50, "Mod. Stererographics of 50 U.S.")
 PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff")
 PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area")
+PROJ_HEAD(igh,  "Interrupted Goode Homolosine")
 PROJ_HEAD(imw_p, "Internation Map of the World Polyconic")
 PROJ_HEAD(kav5, "Kavraisky V")
 PROJ_HEAD(kav7, "Kavraisky VII")
diff --git a/src/projects.h b/src/projects.h
index 6cf90f4e..eb998d3a 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -454,8 +454,12 @@ struct PJ_LIST  *pj_get_list_ref( void );
 struct PJ_PRIME_MERIDIANS  *pj_get_prime_meridians_ref( void );
  
 #ifndef DISABLE_CVSID
-#  define PJ_CVSID(string)     static char pj_cvsid[] = string; \
-static char *cvsid_aw() { return( cvsid_aw() ? ((char *) NULL) : pj_cvsid ); }
+#  if defined(__GNUC__) && __GNUC__ >= 4
+#    define PJ_CVSID(string)     static char pj_cvsid[] __attribute__((used)) = string;
+#  else
+#    define PJ_CVSID(string)     static char pj_cvsid[] = string; \
+static char *cvsid_aw() { return( cvsid_aw() ? ((char *) NULL) : cpl_cvsid ); }
+#  endif
 #else
 #  define PJ_CVSID(string)
 #endif