diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/PJ_aeqd.c | 17 | ||||
| -rw-r--r-- | src/geod_interface.c | 22 | ||||
| -rw-r--r-- | src/proj_api.h | 4 |
3 files changed, 18 insertions, 25 deletions
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 */ |