From 7b30a05ab90cf122c5818c8eb04101687d4dd2f3 Mon Sep 17 00:00:00 2001
From: Charles Karney <ckarney@karney.com>
Date: Thu, 1 Oct 2015 17:21:45 -0400
Subject: Take 2 of the fix to dmstor.

I'm making this pull request soon after the release of 4.9.2.  It will
cause the results that people get out of proj.4 to change very slightly.
If there are problems, we'll get a chance to iron them out well before
the next release.

The important change is to use DEG_TO_RAD for degree to radian
conversions in dmstor.c instead of the slightly inaccurate number used
earlier.

This necessitates a change to geod_interface.c (which previously had to
work aroung the previous bad behavior).  PJ_aeqd.c now does conversions
in a compatible manner.

In src/proj_api.h, DEG_TO_RAD and RAD_TO_DEG are both given to 17
significant figures (this is just a cosmetic change).

I've "fixed" the testvarious tests so that they still pass (on my system
at least).  Everyone should be suitably skeptical of my fixes.

(1) Minor changes to "Test bug 244" and only ask for nanometer (instead
of picometer) accuracy on positions.

(2) 2 lon_wrap tests now return 0dE instead of 360dE (now it's tight?)

(3) The results for the forward healpix projection of (-180, +90) and
(-180, -90) are now different.  I have put the new values into
tv_out.dist.  I'm fairly confident that the new values are OK, since
this projection has various cuts which meet at the poles.  It would be
good if someone who knows about this projection can verify this.
---
 src/PJ_aeqd.c        | 17 ++++++++---------
 src/geod_interface.c | 22 ++++++++--------------
 src/proj_api.h       |  4 ++--
 3 files changed, 18 insertions(+), 25 deletions(-)

(limited to 'src')

diff --git a/src/PJ_aeqd.c b/src/PJ_aeqd.c
index 739dfdd6..22a75ac8 100644
--- a/src/PJ_aeqd.c
+++ b/src/PJ_aeqd.c
@@ -44,7 +44,6 @@ PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam";
 
 #define EPS10 1.e-10
 #define TOL 1.e-14
-#define RHO 57.295779513082320876798154814105
 
 #define N_POLE	0
 #define S_POLE	1
@@ -85,11 +84,11 @@ FORWARD(e_forward); /* elliptical */
 			break;
 		}
 
-		phi1 = P->phi0*RHO; lam1 = P->lam0*RHO;
-		phi2 = lp.phi*RHO;  lam2 = (lp.lam+P->lam0)*RHO;
+		phi1 = P->phi0 / DEG_TO_RAD; lam1 = P->lam0 / DEG_TO_RAD;
+		phi2 = lp.phi / DEG_TO_RAD;  lam2 = (lp.lam+P->lam0) / DEG_TO_RAD;
 
 		geod_inverse(&P->g, phi1, lam1, phi2, lam2, &s12, &azi1, &azi2);
-		azi1 /= RHO;
+		azi1 *= DEG_TO_RAD;
 		xy.x = s12 * sin(azi1) / P->a;
 		xy.y = s12 * cos(azi1) / P->a;
 		break;
@@ -160,13 +159,13 @@ INVERSE(e_inverse); /* elliptical */
 
 		x2 = xy.x * P->a;
 		y2 = xy.y * P->a;
-		lat1 = P->phi0 * RHO;
-		lon1 = P->lam0 * RHO;
-		azi1 = atan2(x2, y2) * RHO;
+		lat1 = P->phi0 / DEG_TO_RAD;
+		lon1 = P->lam0 / DEG_TO_RAD;
+		azi1 = atan2(x2, y2) / DEG_TO_RAD;
 		s12 = sqrt(x2 * x2 + y2 * y2);
 		geod_direct(&P->g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2);
-		lp.phi = lat2 / RHO;
-		lp.lam = lon2 / RHO;
+		lp.phi = lat2 * DEG_TO_RAD;
+		lp.lam = lon2 * DEG_TO_RAD;
 		lp.lam -= P->lam0;
 	} else { /* Polar */
 		lp.phi = pj_inv_mlfn(P->ctx, P->mode == N_POLE ? P->Mp - c : P->Mp + c,
diff --git a/src/geod_interface.c b/src/geod_interface.c
index 63b16b1e..a30377ac 100644
--- a/src/geod_interface.c
+++ b/src/geod_interface.c
@@ -1,20 +1,14 @@
 #include "projects.h"
 #include "geod_interface.h"
 
-/* DEG_IN is a crock to work around the problem that dmstor.c uses the wrong
- * value for pi/180 (namely .0174532925199433 which is an inaccurately
- * truncated version of DEG_TO_RAD).
- */
-#define DEG_IN .0174532925199433
-#define DEG_OUT DEG_TO_RAD;
-
 void geod_ini(void) {
   geod_init(&GlobalGeodesic, geod_a, geod_f);
 }
 
 void geod_pre(void) {
   double
-    lat1 = phi1 / DEG_IN, lon1 = lam1 / DEG_IN, azi1 = al12 / DEG_IN;
+    lat1 = phi1 / DEG_TO_RAD, lon1 = lam1 / DEG_TO_RAD,
+    azi1 = al12 / DEG_TO_RAD;
   geod_lineinit(&GlobalGeodesicLine, &GlobalGeodesic, lat1, lon1, azi1, 0U);
 }
 
@@ -23,17 +17,17 @@ void geod_for(void) {
     s12 = geod_S, lat2, lon2, azi2;
   geod_position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
   azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
-  phi2 = lat2 * DEG_OUT;
-  lam2 = lon2 * DEG_OUT;
-  al21 = azi2 * DEG_OUT;
+  phi2 = lat2 * DEG_TO_RAD;
+  lam2 = lon2 * DEG_TO_RAD;
+  al21 = azi2 * DEG_TO_RAD;
 }
 
 void geod_inv(void) {
   double
-    lat1 = phi1 / DEG_IN, lon1 = lam1 / DEG_IN,
-    lat2 = phi2 / DEG_IN, lon2 = lam2 / DEG_IN,
+    lat1 = phi1 / DEG_TO_RAD, lon1 = lam1 / DEG_TO_RAD,
+    lat2 = phi2 / DEG_TO_RAD, lon2 = lam2 / DEG_TO_RAD,
     azi1, azi2, s12;
   geod_inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
   azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
-  al12 = azi1 * DEG_OUT; al21 = azi2 * DEG_OUT; geod_S = s12;
+  al12 = azi1 * DEG_TO_RAD; al21 = azi2 * DEG_TO_RAD; geod_S = s12;
 }
diff --git a/src/proj_api.h b/src/proj_api.h
index 99faeafb..9e1e49a8 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -46,8 +46,8 @@ extern "C" {
 
 extern char const pj_release[]; /* global release id string */
 
-#define RAD_TO_DEG	57.29577951308232
-#define DEG_TO_RAD	.0174532925199432958
+#define RAD_TO_DEG	57.295779513082321
+#define DEG_TO_RAD	.017453292519943296
 
 
 extern int pj_errno;	/* global error return code */
-- 
cgit v1.2.3