added +sweep for +proj=geos (#146)

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

5 files changed, 172 insertions, 10 deletions

diff --git a/ChangeLog b/ChangeLog
index 05ff7f32..a80a9111 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,8 @@
 2012-02-26  Frank Warmerdam  <warmerdam@pobox.com>
 
+	* src/PJ_geos.c, nad/testvarious: Added GEOS +sweep and add GEOS 
+	to the test suite (#146)
+
 	* nad/CH: added swiss datum related definitions from strk (#145)
 
 	* src/Makefile.am, src/mutltistresstest.c: provide for building
diff --git a/nad/testvarious b/nad/testvarious
index 3505dd16..7a753404 100755
--- a/nad/testvarious
+++ b/nad/testvarious
@@ -295,6 +295,72 @@ $EXE +proj=latlong +a=1 +lon_0=0 +ellps=WGS84 \
 0 0.7853981633974483
 -1.5707963267948966 0
 EOF
+echo "##############################################################" >> ${OUT}
+echo "Test geos projection" >> ${OUT}
+echo "Test geos on a sphere" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere -E >>${OUT} <<EOF
+16d11'8" 58d35'31"
+-43d11'47" -22d54'30"
+18d25'26" -33d55'31"
+47d58'42" 29d22'11"
+EOF
+echo "Test geos on a ellipsoid" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 -E >>${OUT} <<EOF
+16d11'8" 58d35'31"
+-43d11'47" -22d54'30"
+18d25'26" -33d55'31"
+47d58'42" 29d22'11"
+EOF
+echo "Test inv geos on a sphere" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere -I -E >>${OUT} <<EOF
+849736.77 4960015.43
+-3780930.93 -2326595.36
+1608689.65 -3412115.56
+3825202.59 2885980.79
+EOF
+echo "Test inv geos on a ellipsoid" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 -I -E >>${OUT} <<EOF
+852862.53 4945122.70
+-3787026.57 -2314765.32
+1612331.00 -3397031.37
+3832522.65 2872185.29
+EOF
+echo "Test geos on a sphere with alternate sweep" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere +sweep=y -E >>${OUT} <<EOF
+16d11'8" 58d35'31"
+-43d11'47" -22d54'30"
+18d25'26" -33d55'31"
+47d58'42" 29d22'11"
+EOF
+echo "Test geos on a ellipsoid with alternate sweep" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 +sweep=y -E >>${OUT} <<EOF
+16d11'8" 58d35'31"
+-43d11'47" -22d54'30"
+18d25'26" -33d55'31"
+47d58'42" 29d22'11"
+EOF
+echo "Test inv geos on a sphere with alternate sweep" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere +sweep=y -I -E >>${OUT} <<EOF
+841586.28 4961396.21
+-3772913.22 -2339604.71
+1601377.77 -3415545.15
+3812722.89 2902474.62
+EOF
+echo "Test inv geos on a ellipsoid with alternate sweep" >> ${OUT}
+$EXE +proj=latlong +ellps=sphere \
+ +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 +sweep=y -I -E >>${OUT} <<EOF
+844731.03 4946509.59
+-3779077.27 -2327750.87
+1605067.15 -3400461.47
+3820138.08 2888664.15
+EOF
 #
 ##############################################################################
 # Done! 
diff --git a/nad/tv_out.dist b/nad/tv_out.dist
index b7365f25..e13367f6 100644
--- a/nad/tv_out.dist
+++ b/nad/tv_out.dist
@@ -110,3 +110,45 @@ Test forward projection on ellipsoid
 Test inverse projection on ellipsoid
 0 0.7853981633974483	0.00000	41.93785 0.00000
 -1.5707963267948966 0	-90.00000	0.00000 0.00000
+##############################################################
+Test geos projection
+Test geos on a sphere
+16d11'8" 58d35'31"	849736.77	4960015.43 0.00
+-43d11'47" -22d54'30"	-3780930.93	-2326595.36 0.00
+18d25'26" -33d55'31"	1608689.65	-3412115.56 0.00
+47d58'42" 29d22'11"	3825202.59	2885980.79 0.00
+Test geos on a ellipsoid
+16d11'8" 58d35'31"	852862.53	4945122.70 0.00
+-43d11'47" -22d54'30"	-3787026.57	-2314765.32 0.00
+18d25'26" -33d55'31"	1612331.00	-3397031.37 0.00
+47d58'42" 29d22'11"	3832522.65	2872185.29 0.00
+Test inv geos on a sphere
+849736.77 4960015.43	16d11'8"E	58d35'31"N 0.000
+-3780930.93 -2326595.36	43d11'47"W	22d54'30"S 0.000
+1608689.65 -3412115.56	18d25'26"E	33d55'31"S 0.000
+3825202.59 2885980.79	47d58'42"E	29d22'11"N 0.000
+Test inv geos on a ellipsoid
+852862.53 4945122.70	16d11'8"E	58d35'31"N 0.000
+-3787026.57 -2314765.32	43d11'47"W	22d54'30"S 0.000
+1612331.00 -3397031.37	18d25'26"E	33d55'31"S 0.000
+3832522.65 2872185.29	47d58'42"E	29d22'11"N 0.000
+Test geos on a sphere with alternate sweep
+16d11'8" 58d35'31"	841586.28	4961396.21 0.00
+-43d11'47" -22d54'30"	-3772913.22	-2339604.71 0.00
+18d25'26" -33d55'31"	1601377.77	-3415545.15 0.00
+47d58'42" 29d22'11"	3812722.89	2902474.62 0.00
+Test geos on a ellipsoid with alternate sweep
+16d11'8" 58d35'31"	844731.03	4946509.59 0.00
+-43d11'47" -22d54'30"	-3779077.27	-2327750.87 0.00
+18d25'26" -33d55'31"	1605067.15	-3400461.47 0.00
+47d58'42" 29d22'11"	3820138.08	2888664.15 0.00
+Test inv geos on a sphere with alternate sweep
+841586.28 4961396.21	16d11'8"E	58d35'31"N 0.000
+-3772913.22 -2339604.71	43d11'47"W	22d54'30"S 0.000
+1601377.77 -3415545.15	18d25'26"E	33d55'31"S 0.000
+3812722.89 2902474.62	47d58'42"E	29d22'11"N 0.000
+Test inv geos on a ellipsoid with alternate sweep
+844731.03 4946509.59	16d11'8"E	58d35'31"N 0.000
+-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
diff --git a/src/PJ_geos.c b/src/PJ_geos.c
index 9b9d36f4..2629d41d 100644
--- a/src/PJ_geos.c
+++ b/src/PJ_geos.c
@@ -2,6 +2,7 @@
 ** libproj -- library of cartographic projections
 **
 ** Copyright (c) 2004   Gerald I. Evenden
+** Copyright (c) 2012   Martin Raspaud
 */
 static const char
 LIBPROJ_ID[] = "$Id$";
@@ -35,7 +36,9 @@ LIBPROJ_ID[] = "$Id$";
 	double  radius_p_inv2; \
 	double  radius_g; \
 	double  radius_g_1; \
-	double  C;
+	double  C; \
+        char *  sweep_axis; \
+        int     flip_axis;
 #define PJ_LIB__
 #include	<projects.h>
 
@@ -54,8 +57,16 @@ FORWARD(s_forward); /* spheroid */
 	if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;
 /* Calculation based on view angles from satellite.*/
 	tmp = P->radius_g - Vx;
-	xy.x = P->radius_g_1 * atan(Vy / tmp);
-	xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
+        if(P->flip_axis)
+          {
+            xy.x = P->radius_g_1 * atan(Vy / hypot(Vz, tmp));
+            xy.y = P->radius_g_1 * atan(Vz / tmp);
+          }
+        else
+          {
+            xy.x = P->radius_g_1 * atan(Vy / tmp);
+            xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
+          }
 	return (xy);
 }
 FORWARD(e_forward); /* ellipsoid */
@@ -74,8 +85,16 @@ FORWARD(e_forward); /* ellipsoid */
 		F_ERROR;
 /* Calculation based on view angles from satellite. */
 	tmp = P->radius_g - Vx;
-	xy.x = P->radius_g_1 * atan (Vy / tmp);
-	xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
+        if(P->flip_axis)
+          {
+            xy.x = P->radius_g_1 * atan (Vy / hypot (Vz, tmp));
+            xy.y = P->radius_g_1 * atan (Vz / tmp);
+          }
+        else
+          {
+            xy.x = P->radius_g_1 * atan (Vy / tmp);
+            xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
+          }
 	return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
@@ -83,8 +102,16 @@ INVERSE(s_inverse); /* spheroid */
 
 /* Setting three components of vector from satellite to position.*/
 	Vx = -1.0;
-	Vy = tan (xy.x / (P->radius_g - 1.0));
-	Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
+	if(P->flip_axis)
+          {
+            Vz = tan (xy.y / (P->radius_g - 1.0));
+            Vy = tan (xy.x / (P->radius_g - 1.0)) * sqrt (1.0 + Vz * Vz);
+          }
+        else
+          {
+            Vy = tan (xy.x / (P->radius_g - 1.0));
+            Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
+          }
 /* Calculation of terms in cubic equation and determinant.*/
 	a   = Vy * Vy + Vz * Vz + Vx * Vx;
 	b   = 2 * P->radius_g * Vx;
@@ -104,8 +131,16 @@ INVERSE(e_inverse); /* ellipsoid */
 
 /* Setting three components of vector from satellite to position.*/
 	Vx = -1.0;
-	Vy = tan (xy.x / P->radius_g_1);
-	Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
+        if(P->flip_axis)
+          {
+            Vz = tan (xy.y / P->radius_g_1);
+            Vy = tan (xy.x / P->radius_g_1) * hypot(1.0, Vz);
+          }
+        else
+          {
+            Vy = tan (xy.x / P->radius_g_1);
+            Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
+          }
 /* Calculation of terms in cubic equation and determinant.*/
 	a = Vz / P->radius_p;
 	a   = Vy * Vy + a * a + Vx * Vx;
@@ -126,7 +161,22 @@ FREEUP; if (P) free(P); }
 ENTRY0(geos)
 	if ((P->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30);
 	if (P->phi0) E_ERROR(-46);
-	P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
+        P->sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
+        if (P->sweep_axis == NULL)
+          P->flip_axis = 0;
+        else
+          {
+            if (P->sweep_axis[1] != '\0' ||
+                (P->sweep_axis[0] != 'x' &&
+                 P->sweep_axis[0] != 'y'))
+              E_ERROR(-49);
+            if (P->sweep_axis[0] == 'y')
+              P->flip_axis = 1;
+            else
+              P->flip_axis = 0;
+          }
+        P->radius_g_1 = P->h / P->a;
+	P->radius_g = 1. + P->radius_g_1;
 	P->C  = P->radius_g * P->radius_g - 1.0;
 	if (P->es) {
 		P->radius_p      = sqrt (P->one_es);
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index 186a50e8..9d237039 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -54,6 +54,7 @@ pj_err_list[] = {
 	"unknown prime meridian conversion id",		/* -46 */
 	"illegal axis orientation combination",		/* -47 */
 	"point not within available datum shift grids", /* -48 */
+	"invalid sweep axis, choose x or y",            /* -49 */
 };
 	char *
 pj_strerrno(int err)