Created
April 30, 2018 09:50
-
-
Save schwehr/2a1389d8d5b711ff9796530887247d2d to your computer and use it in GitHub Desktop.
A first pass a trying to cleanup PJ_aeqd.c
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/****************************************************************************** | |
* Project: PROJ.4 | |
* Purpose: Implementation of the aeqd (Azimuthal Equidistant) projection. | |
* Author: Gerald Evenden | |
* | |
****************************************************************************** | |
* Copyright (c) 1995, Gerald Evenden | |
* | |
* Permission is hereby granted, free of charge, to any person obtaining a | |
* copy of this software and associated documentation files (the "Software"), | |
* to deal in the Software without restriction, including without limitation | |
* the rights to use, copy, modify, merge, publish, distribute, sublicense, | |
* and/or sell copies of the Software, and to permit persons to whom the | |
* Software is furnished to do so, subject to the following conditions: | |
* | |
* The above copyright notice and this permission notice shall be included | |
* in all copies or substantial portions of the Software. | |
* | |
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | |
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL | |
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | |
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER | |
* DEALINGS IN THE SOFTWARE. | |
*****************************************************************************/ | |
#define PJ_LIB__ | |
#include <errno.h> | |
#include <math.h> | |
#include "geodesic.h" | |
#include "proj.h" | |
#include "proj_api.h" | |
#include "proj_math.h" | |
#include "projects.h" | |
enum Mode { | |
N_POLE = 0, | |
S_POLE = 1, | |
EQUIT = 2, | |
OBLIQ = 3 | |
}; | |
struct pj_opaque { | |
double sinph0; | |
double cosph0; | |
double *en; | |
double M1; | |
double N1; | |
double Mp; | |
double He; | |
double G; | |
enum Mode mode; | |
struct geod_geodesic g; | |
}; | |
PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam"; | |
static const double EPS10 = 1.0e-10; | |
static const double TOL = 1.0e-14; | |
static void *destructor(PJ *P, int errlev) { /* Destructor */ | |
if (0==P) | |
return 0; | |
if (0==P->opaque) | |
return pj_default_destructor(P, errlev); | |
pj_dealloc(P->opaque->en); | |
return pj_default_destructor(P, errlev); | |
} | |
static XY e_guam_fwd(LP lp, PJ *P) { /* Guam elliptical */ | |
const struct pj_opaque *Q = P->opaque; | |
const double cosphi = cos(lp.phi); | |
const double sinphi = sin(lp.phi); | |
const double t = 1. / sqrt(1.0 - P->es * sinphi * sinphi); | |
const XY xy = { | |
lp.lam * cosphi * t, | |
pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->M1 + | |
0.5 * lp.lam * lp.lam * cosphi * sinphi * t}; | |
return xy; | |
} | |
static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ | |
XY xy = {0.0,0.0}; | |
const struct pj_opaque *Q = P->opaque; | |
double coslam = cos(lp.lam); | |
const double cosphi = cos(lp.phi); | |
const double sinphi = sin(lp.phi); | |
switch (Q->mode) { | |
case N_POLE: | |
coslam = -coslam; | |
/*-fallthrough*/ | |
case S_POLE: | |
{ | |
const double rho = | |
fabs(Q->Mp - pj_mlfn(lp.phi, sinphi, cosphi, Q->en)); | |
xy.x = rho * sin(lp.lam); | |
xy.y = rho * coslam; | |
} | |
break; | |
case EQUIT: | |
case OBLIQ: | |
if (fabs(lp.lam) < EPS10 && fabs(lp.phi - P->phi0) < EPS10) { | |
xy.x = 0.0; | |
xy.y = 0.0; | |
break; | |
} | |
{ | |
const double phi1 = P->phi0 / DEG_TO_RAD; | |
const double lam1 = P->lam0 / DEG_TO_RAD; | |
const double phi2 = lp.phi / DEG_TO_RAD; | |
const double lam2 = (lp.lam+P->lam0) / DEG_TO_RAD; | |
double azi1 = 0.0; | |
double azi2 = 0.0; | |
double s12 = 0.0; | |
geod_inverse(&Q->g, phi1, lam1, phi2, lam2, &s12, &azi1, &azi2); | |
/* axi2 and sl2 not used after set. */ | |
azi1 *= DEG_TO_RAD; | |
xy.x = s12 * sin(azi1) / P->a; | |
xy.y = s12 * cos(azi1) / P->a; | |
break; | |
} | |
} | |
return xy; | |
} | |
static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ | |
XY xy = {0.0,0.0}; | |
struct pj_opaque *Q = P->opaque; | |
double coslam = cos(lp.lam); | |
const double cosphi = cos(lp.phi); | |
const double sinphi = sin(lp.phi); | |
switch (Q->mode) { | |
case EQUIT: | |
case OBLIQ: | |
xy.y = | |
Q->mode == EQUIT | |
? cosphi * coslam | |
: Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam; | |
if (fabs(fabs(xy.y) - 1.0) < TOL) { | |
if (xy.y < 0.0) { | |
proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); | |
return xy; | |
} else { | |
xy.x = 0.0; | |
xy.y = 0.0; | |
} | |
} else { | |
xy.y = acos(xy.y); | |
xy.y /= sin(xy.y); | |
xy.x = xy.y * cosphi * sin(lp.lam); | |
xy.y *= Q->mode == EQUIT ? sinphi : | |
Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam; | |
} | |
break; | |
case N_POLE: | |
lp.phi = -lp.phi; | |
coslam = -coslam; | |
/*-fallthrough*/ | |
case S_POLE: | |
if (fabs(lp.phi - M_HALFPI) < EPS10) { | |
proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); | |
return xy; | |
} | |
xy.y = M_HALFPI + lp.phi; | |
xy.x = xy.y * sin(lp.lam); | |
xy.y *= coslam; | |
break; | |
} | |
return xy; | |
} | |
static LP e_guam_inv(XY xy, PJ *P) { /* Guam elliptical */ | |
const struct pj_opaque *Q = P->opaque; | |
const double x2 = 0.5 * xy.x * xy.x; | |
LP lp = {0.0 /* lam */, P->phi0}; | |
int i; /* Not used after for */ | |
double t = 0.0; | |
for (i = 0; i < 3; ++i) { | |
const double t_tmp = P->e * sin(lp.phi); | |
t = sqrt(1. - t_tmp * t_tmp); | |
lp.phi = pj_inv_mlfn(P->ctx, Q->M1 + xy.y - | |
x2 * tan(lp.phi) * t, P->es, Q->en); | |
} | |
lp.lam = xy.x * t / cos(lp.phi); | |
return lp; | |
} | |
static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ | |
LP lp = {0.0,0.0}; | |
struct pj_opaque *Q = P->opaque; | |
const double c = hypot(xy.x, xy.y); | |
if (c < EPS10) { | |
lp.phi = P->phi0; | |
lp.lam = 0.0; | |
return lp; | |
} | |
if (Q->mode == OBLIQ || Q->mode == EQUIT) { | |
const double x2 = xy.x * P->a; | |
const double y2 = xy.y * P->a; | |
const double lat1 = P->phi0 / DEG_TO_RAD; | |
const double lon1 = P->lam0 / DEG_TO_RAD; | |
const double azi1 = atan2(x2, y2) / DEG_TO_RAD; | |
const double s12 = sqrt(x2 * x2 + y2 * y2); | |
double lat2 = 0.0; | |
double lon2 = 0.0; | |
double azi2 = 0.0; | |
geod_direct(&Q->g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2); | |
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, Q->mode == N_POLE ? Q->Mp - c : Q->Mp + c, | |
P->es, Q->en); | |
lp.lam = atan2(xy.x, Q->mode == N_POLE ? -xy.y : xy.y); | |
} | |
return lp; | |
} | |
static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ | |
LP lp = {0.0,0.0}; | |
struct pj_opaque *Q = P->opaque; | |
double c_rh = hypot(xy.x, xy.y); | |
if (c_rh > M_PI) { | |
if (c_rh - EPS10 > M_PI) { | |
proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); | |
return lp; | |
} | |
c_rh = M_PI; | |
} else if (c_rh < EPS10) { | |
lp.phi = P->phi0; | |
lp.lam = 0.; | |
return lp; | |
} | |
if (Q->mode == OBLIQ || Q->mode == EQUIT) { | |
const double sinc = sin(c_rh); | |
const double cosc = cos(c_rh); | |
if (Q->mode == EQUIT) { | |
lp.phi = aasin(P->ctx, xy.y * sinc / c_rh); | |
xy.x *= sinc; | |
xy.y = cosc * c_rh; | |
} else { | |
lp.phi = aasin(P->ctx,cosc * Q->sinph0 + xy.y * sinc * Q->cosph0 / | |
c_rh); | |
xy.y = (cosc - Q->sinph0 * sin(lp.phi)) * c_rh; | |
xy.x *= sinc * Q->cosph0; | |
} | |
lp.lam = xy.y == 0.0 ? 0.0 : atan2(xy.x, xy.y); | |
} else if (Q->mode == N_POLE) { | |
lp.phi = M_HALFPI - c_rh; | |
lp.lam = atan2(xy.x, -xy.y); | |
} else { | |
lp.phi = c_rh - M_HALFPI; | |
lp.lam = atan2(xy.x, xy.y); | |
} | |
return lp; | |
} | |
PJ *PROJECTION(aeqd) { | |
struct pj_opaque *Q = pj_calloc(1, sizeof (struct pj_opaque)); | |
if (0==Q) | |
return pj_default_destructor(P, ENOMEM); | |
P->opaque = Q; | |
P->destructor = destructor; | |
geod_init(&Q->g, P->a, P->es / (1 + sqrt(P->one_es))); | |
if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) { | |
Q->mode = P->phi0 < 0.0 ? S_POLE : N_POLE; | |
Q->sinph0 = P->phi0 < 0.0 ? -1.0 : 1.0; | |
Q->cosph0 = 0.0; | |
} else if (fabs(P->phi0) < EPS10) { | |
Q->mode = EQUIT; | |
Q->sinph0 = 0.0; | |
Q->cosph0 = 1.0; | |
} else { | |
Q->mode = OBLIQ; | |
Q->sinph0 = sin(P->phi0); | |
Q->cosph0 = cos(P->phi0); | |
} | |
if (P->es == 0.0) { | |
P->inv = s_inverse; | |
P->fwd = s_forward; | |
} else { | |
if (!(Q->en = pj_enfn(P->es))) | |
return pj_default_destructor(P, 0); | |
if (pj_param(P->ctx, P->params, "bguam").i) { | |
Q->M1 = pj_mlfn(P->phi0, Q->sinph0, Q->cosph0, Q->en); | |
P->inv = e_guam_inv; | |
P->fwd = e_guam_fwd; | |
} else { | |
switch (Q->mode) { | |
case N_POLE: | |
Q->Mp = pj_mlfn(M_HALFPI, 1.0, 0.0, Q->en); | |
break; | |
case S_POLE: | |
Q->Mp = pj_mlfn(-M_HALFPI, -1.0, 0.0, Q->en); | |
break; | |
case EQUIT: | |
case OBLIQ: | |
P->inv = e_inverse; | |
P->fwd = e_forward; | |
Q->N1 = 1.0 / sqrt(1.0 - P->es * Q->sinph0 * Q->sinph0); | |
Q->He = P->e / sqrt(P->one_es); | |
Q->G = Q->sinph0 * Q->He; | |
Q->He *= Q->cosph0; | |
break; | |
} | |
P->inv = e_inverse; | |
P->fwd = e_forward; | |
} | |
} | |
return P; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment