added the Natural Earth projection

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

8 files changed, 142 insertions, 4 deletions

diff --git a/NEWS b/NEWS
index 4ad2f9f6..99ecdc5f 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,8 @@
 4.8.0 Release Notes
 -------------------
 
+ o Added the Natural Earth projection.
+
  o Added HEALPIX, rHEALPIX and Icosahedral Snyder Equal Area projections.
 
  o nad2bin now produces "CTable2" format grid shift files by default which
diff --git a/nad/testvarious b/nad/testvarious
index aadb8796..3886de56 100755
--- a/nad/testvarious
+++ b/nad/testvarious
@@ -371,6 +371,37 @@ $EXE +proj=latlong +ellps=sphere \
 3820138.08 2888664.15
 EOF
 #
+echo "##############################################################" >> ${OUT}
+echo "Test the natural earth projection" >> ${OUT}
+$EXE +proj=latlong +a=1 +lon_0=0 \
+ +to +proj=natearth +a=6371008.7714 +b=6371008.7714 -f '%.'7'f' \
+ -E >>${OUT} <<EOF
+0.0 0.0 0 0.0 0.0
+0.0 22.5 0 0.0 2525419.569383768
+0.0 45.0 0 0.0 5052537.389973222
+0.0 67.5 0 0.0 7400065.6562573705
+0.0 90.0 0 0.0 9062062.394736718
+45.0 0.0 0 4356790.016612169 0.0
+45.0 22.5 0 4253309.544984069 2525419.569383768
+45.0 45.0 0 3924521.5829515466 5052537.389973222
+45.0 67.5 0 3354937.47115583 7400065.6562573705
+45.0 90.0 0 2397978.2448443635 9062062.394736718
+90.0 0.0 0 8713580.033224339 0.0
+90.0 22.5 0 8506619.089968137 2525419.569383768
+90.0 45.0 0 7849043.165903093 5052537.389973222
+90.0 67.5 0 6709874.94231166 7400065.6562573705
+90.0 90.0 0 4795956.489688727 9062062.394736718
+135.0 0.0 0 1.3070370049836507E7 0.0
+135.0 22.5 0 1.2759928634952208E7 2525419.569383768
+135.0 45.0 0 1.177356474885464E7 5052537.389973222
+135.0 67.5 0 1.0064812413467491E7 7400065.6562573705
+135.0 90.0 0 7193934.734533091 9062062.394736718
+180.0 0.0 0 1.7427160066448677E7 0.0
+180.0 22.5 0 1.7013238179936275E7 2525419.569383768
+180.0 45.0 0 1.5698086331806187E7 5052537.389973222
+180.0 67.5 0 1.341974988462332E7 7400065.6562573705
+180.0 90.0 0 9591912.979377454 9062062.394736718
+EOF
 ##############################################################################
 # Done! 
 # do 'diff' with distribution results
diff --git a/nad/tv_out.dist b/nad/tv_out.dist
index 97c05ddf..912302e3 100644
--- a/nad/tv_out.dist
+++ b/nad/tv_out.dist
@@ -155,3 +155,30 @@ Test inv geos on a ellipsoid with alternate sweep
 -3779077.27 -2327750.87	43d11'47"W	22d54'30"S 0.000
 1605067.15 -3400461.47	18d25'26"E	33d55'31"S 0.000
 3820138.08 2888664.15	47d58'42"E	29d22'11"N 0.000
+##############################################################
+Test the natural earth projection
+0.0 0.0 0	0.0000000	0.0000000 0.0000000 0.0 0.0
+0.0 22.5 0	0.0000000	2525419.5693838 0.0000000 0.0 2525419.569383768
+0.0 45.0 0	0.0000000	5052537.3899732 0.0000000 0.0 5052537.389973222
+0.0 67.5 0	0.0000000	7400065.6562574 0.0000000 0.0 7400065.6562573705
+0.0 90.0 0	0.0000000	9062062.3947367 0.0000000 0.0 9062062.394736718
+45.0 0.0 0	4356790.0166122	0.0000000 0.0000000 4356790.016612169 0.0
+45.0 22.5 0	4253309.5449841	2525419.5693838 0.0000000 4253309.544984069 2525419.569383768
+45.0 45.0 0	3924521.5829515	5052537.3899732 0.0000000 3924521.5829515466 5052537.389973222
+45.0 67.5 0	3354937.4711558	7400065.6562574 0.0000000 3354937.47115583 7400065.6562573705
+45.0 90.0 0	2397978.2448444	9062062.3947367 0.0000000 2397978.2448443635 9062062.394736718
+90.0 0.0 0	8713580.0332243	0.0000000 0.0000000 8713580.033224339 0.0
+90.0 22.5 0	8506619.0899681	2525419.5693838 0.0000000 8506619.089968137 2525419.569383768
+90.0 45.0 0	7849043.1659031	5052537.3899732 0.0000000 7849043.165903093 5052537.389973222
+90.0 67.5 0	6709874.9423117	7400065.6562574 0.0000000 6709874.94231166 7400065.6562573705
+90.0 90.0 0	4795956.4896887	9062062.3947367 0.0000000 4795956.489688727 9062062.394736718
+135.0 0.0 0	13070370.0498365	0.0000000 0.0000000 1.3070370049836507E7 0.0
+135.0 22.5 0	12759928.6349522	2525419.5693838 0.0000000 1.2759928634952208E7 2525419.569383768
+135.0 45.0 0	11773564.7488546	5052537.3899732 0.0000000 1.177356474885464E7 5052537.389973222
+135.0 67.5 0	10064812.4134675	7400065.6562574 0.0000000 1.0064812413467491E7 7400065.6562573705
+135.0 90.0 0	7193934.7345331	9062062.3947367 0.0000000 7193934.734533091 9062062.394736718
+180.0 0.0 0	17427160.0664487	0.0000000 0.0000000 1.7427160066448677E7 0.0
+180.0 22.5 0	17013238.1799363	2525419.5693838 0.0000000 1.7013238179936275E7 2525419.569383768
+180.0 45.0 0	15698086.3318062	5052537.3899732 0.0000000 1.5698086331806187E7 5052537.389973222
+180.0 67.5 0	13419749.8846233	7400065.6562574 0.0000000 1.341974988462332E7 7400065.6562573705
+180.0 90.0 0	9591912.9793775	9062062.3947367 0.0000000 9591912.979377454 9062062.394736718
diff --git a/src/Makefile.am b/src/Makefile.am
index 97e9e833..621b0b4a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -57,7 +57,7 @@ libproj_la_SOURCES = \
 	pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \
 	pj_qsfn.c pj_strerrno.c pj_tsfn.c pj_units.c pj_ctx.c pj_log.c \
 	pj_zpoly1.c rtodms.c vector1.c pj_release.c pj_gauss.c \
-	PJ_healpix.c \
+	PJ_healpix.c PJ_natearth.c \
 	\
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 412bd0bd..61ef536d 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -102,7 +102,7 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
 	pj_open_lib.lo pj_param.lo pj_phi2.lo pj_pr_list.lo pj_qsfn.lo \
 	pj_strerrno.lo pj_tsfn.lo pj_units.lo pj_ctx.lo pj_log.lo \
 	pj_zpoly1.lo rtodms.lo vector1.lo pj_release.lo pj_gauss.lo \
-	PJ_healpix.lo nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \
+	PJ_healpix.lo PJ_natearth.lo nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \
 	pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \
 	pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \
 	pj_gridlist.lo jniproj.lo pj_mutex.lo pj_initcache.lo \
@@ -316,7 +316,7 @@ libproj_la_SOURCES = \
 	pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \
 	pj_qsfn.c pj_strerrno.c pj_tsfn.c pj_units.c pj_ctx.c pj_log.c \
 	pj_zpoly1.c rtodms.c vector1.c pj_release.c pj_gauss.c \
-	PJ_healpix.c \
+	PJ_healpix.c PJ_natearth.c \
 	\
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
diff --git a/src/PJ_natearth.c b/src/PJ_natearth.c
new file mode 100644
index 00000000..f83a3ffd
--- /dev/null
+++ b/src/PJ_natearth.c
@@ -0,0 +1,77 @@
+/*
+The Natural Earth projection was designed by Tom Patterson, US National Park
+Service, in 2007, using Flex Projector. The shape of the original projection
+was defined at every 5 degrees and piece-wise cubic spline interpolation was
+used to compute the complete graticule.
+The code here uses polynomial functions instead of cubic splines and
+is therefore much simpler to program. The polynomial approximation was
+developed by Bojan Savric, in collaboration with Tom Patterson and Bernhard
+Jenny, Institute of Cartography, ETH Zurich. It slightly deviates from
+Patterson's original projection by adding additional curvature to meridians
+where they meet the horizontal pole line. This improvement is by intention
+and designed in collaboration with Tom Patterson.
+Port to PROJ.4 by Bernhard Jenny, 6 June 2011
+*/
+
+#define PJ_LIB__
+#include	<projects.h>
+PROJ_HEAD(natearth, "Natural Earth") "\n\tPCyl., Sph.";
+#define A0 0.8707
+#define A1 -0.131979
+#define A2 -0.013791
+#define A3 0.003971
+#define A4 -0.001529
+#define B0 1.007226
+#define B1 0.015085
+#define B2 -0.044475
+#define B3 0.028874
+#define B4 -0.005916
+#define C0 B0
+#define C1 (3 * B1)
+#define C2 (7 * B2)
+#define C3 (9 * B3)
+#define C4 (11 * B4)
+#define EPS 1e-11
+#define MAX_Y (0.8707 * 0.52 * PI)
+
+FORWARD(s_forward); /* spheroid */
+	double phi2, phi4;
+
+	phi2 = lp.phi * lp.phi;
+	phi4 = phi2 * phi2;
+	xy.x = lp.lam * (A0 + phi2 * (A1 + phi2 * (A2 + phi4 * phi2 * (A3 + phi2 * A4))));
+	xy.y = lp.phi * (B0 + phi2 * (B1 + phi4 * (B2 + B3 * phi2 + B4 * phi4)));
+	return (xy);
+}
+INVERSE(s_inverse); /* spheroid */
+	double yc, tol, y2, y4, f, fder;
+
+        /* make sure y is inside valid range */
+	if (xy.y > MAX_Y) {
+		xy.y = MAX_Y;
+	} else if (xy.y < -MAX_Y) {
+		xy.y = -MAX_Y;
+	}
+
+        /* latitude */
+	yc = xy.y;
+        for (;;) { /* Newton-Raphson */
+		y2 = yc * yc;
+		y4 = y2 * y2;
+		f = (yc * (B0 + y2 * (B1 + y4 * (B2 + B3 * y2 + B4 * y4)))) - xy.y;
+		fder = C0 + y2 * (C1 + y4 * (C2 + C3 * y2 + C4 * y4));
+		yc -= tol = f / fder;
+		if (fabs(tol) < EPS) {
+			break;
+		}
+	}
+	lp.phi = yc;
+
+        /* longitude */
+	y2 = yc * yc;
+	lp.lam = xy.x / (A0 + y2 * (A1 + y2 * (A2 + y2 * y2 * y2 * (A3 + y2 * A4))));
+
+	return (lp);
+}
+FREEUP; if (P) pj_dalloc(P); }
+ENTRY0(natearth) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
diff --git a/src/makefile.vc b/src/makefile.vc
index 220468ab..58c348f6 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -27,7 +27,7 @@ misc = \
 	PJ_lask.obj PJ_nocol.obj PJ_ob_tran.obj PJ_oea.obj \
 	PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj PJ_vandg4.obj \
 	PJ_wag7.obj pj_latlong.obj PJ_krovak.obj pj_geocent.obj \
-	PJ_healpix.obj
+	PJ_healpix.obj PJ_natearth.obj
 
 pseudo = \
 	PJ_boggs.obj PJ_collg.obj PJ_crast.obj PJ_denoy.obj \
diff --git a/src/pj_list.h b/src/pj_list.h
index ac88f9d8..8230f92f 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -82,6 +82,7 @@ PROJ_HEAD(moll, "Mollweide")
 PROJ_HEAD(murd1, "Murdoch I")
 PROJ_HEAD(murd2, "Murdoch II")
 PROJ_HEAD(murd3, "Murdoch III")
+PROJ_HEAD(natearth, "Natural Earth")
 PROJ_HEAD(nell, "Nell")
 PROJ_HEAD(nell_h, "Nell-Hammer")
 PROJ_HEAD(nicol, "Nicolosi Globular")