diff options
56 files changed, 1181 insertions, 810 deletions
diff --git a/src/PJ_aea.c b/src/PJ_aea.c index 524657b2..303a319e 100644 --- a/src/PJ_aea.c +++ b/src/PJ_aea.c @@ -28,7 +28,8 @@ *****************************************************************************/ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" # define EPS10 1.e-10 # define TOL7 1.e-7 @@ -78,12 +79,33 @@ struct pj_opaque { }; +static void *freeup_new (PJ *P) { /* Destructor */ + if (0==P) + return 0; + + if (0==P->opaque) + return pj_dealloc (P); + + pj_dealloc (P->opaque->en); + pj_dealloc (P->opaque); + return pj_dealloc(P); +} + + +static void freeup (PJ *P) { + freeup_new (P); + return; +} + static XY e_forward (LP lp, PJ *P) { /* Ellipsoid/spheroid, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - if ((Q->rho = Q->c - (Q->ellips ? Q->n * pj_qsfn(sin(lp.phi), - P->e, P->one_es) : Q->n2 * sin(lp.phi))) < 0.) F_ERROR + Q->rho = Q->c - (Q->ellips ? Q->n * pj_qsfn(sin(lp.phi), P->e, P->one_es) : Q->n2 * sin(lp.phi));; + if (Q->rho < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } Q->rho = Q->dd * sqrt(Q->rho); xy.x = Q->rho * sin( lp.lam *= Q->n ); xy.y = Q->rho0 - Q->rho * cos(lp.lam); @@ -104,8 +126,10 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoid/spheroid, inverse */ if (Q->ellips) { lp.phi = (Q->c - lp.phi * lp.phi) / Q->n; if (fabs(Q->ec - fabs(lp.phi)) > TOL7) { - if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL) - I_ERROR + if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } } else lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI; } else if (fabs(lp.phi = (Q->c - lp.phi * lp.phi) / Q->n2) <= 1.) @@ -121,24 +145,6 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoid/spheroid, inverse */ } -static void *freeup_new (PJ *P) { /* Destructor */ - if (0==P) - return 0; - - if (0==P->opaque) - return pj_dealloc (P); - - pj_dealloc (P->opaque->en); - pj_dealloc (P->opaque); - return pj_dealloc(P); -} - - -static void freeup (PJ *P) { - freeup_new (P); - return; -} - static PJ *setup(PJ *P) { double cosphi, sinphi; @@ -148,14 +154,17 @@ static PJ *setup(PJ *P) { P->inv = e_inverse; P->fwd = e_forward; - if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21); + if (fabs(Q->phi1 + Q->phi2) < EPS10) { + proj_errno_set(P, PJD_ERR_CONIC_LAT_EQUAL); + return freeup_new(P); + } Q->n = sinphi = sin(Q->phi1); cosphi = cos(Q->phi1); secant = fabs(Q->phi1 - Q->phi2) >= EPS10; if( (Q->ellips = (P->es > 0.))) { double ml1, m1; - if (!(Q->en = pj_enfn(P->es))) E_ERROR_0; + if (!(Q->en = pj_enfn(P->es))) return freeup_new(P); m1 = pj_msfn(sinphi, cosphi, P->es); ml1 = pj_qsfn(sinphi, P->e, P->one_es); if (secant) { /* secant cone */ diff --git a/src/PJ_aeqd.c b/src/PJ_aeqd.c index 225acfdd..4a46cf1c 100644 --- a/src/PJ_aeqd.c +++ b/src/PJ_aeqd.c @@ -27,7 +27,8 @@ #define PJ_LIB__ #include "geodesic.h" -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double sinph0; @@ -53,6 +54,24 @@ PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam"; #define OBLIQ 3 +static void *freeup_new (PJ *P) { /* Destructor */ + if (0==P) + return 0; + if (0==P->opaque) + return pj_dealloc (P); + + if (P->opaque->en) + pj_dealloc(P->opaque->en); + pj_dealloc (P->opaque); + return pj_dealloc(P); +} + + +static void freeup (PJ *P) { + freeup_new (P); + return; +} + static XY e_guam_fwd(LP lp, PJ *P) { /* Guam elliptical */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; @@ -124,8 +143,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam; oblcon: if (fabs(fabs(xy.y) - 1.) < TOL) - if (xy.y < 0.) - F_ERROR + if (xy.y < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } else xy.x = xy.y = 0.; else { @@ -141,7 +162,10 @@ oblcon: coslam = -coslam; /*-fallthrough*/ case S_POLE: - if (fabs(lp.phi - M_HALFPI) < EPS10) F_ERROR; + if (fabs(lp.phi - M_HALFPI) < EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = (xy.y = (M_HALFPI + lp.phi)) * sin(lp.lam); xy.y *= coslam; break; @@ -206,7 +230,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ double cosc, c_rh, sinc; if ((c_rh = hypot(xy.x, xy.y)) > M_PI) { - if (c_rh - EPS10 > M_PI) I_ERROR; + 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; @@ -238,25 +265,6 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ } -static void *freeup_new (PJ *P) { /* Destructor */ - if (0==P) - return 0; - if (0==P->opaque) - return pj_dealloc (P); - - if (P->opaque->en) - pj_dealloc(P->opaque->en); - pj_dealloc (P->opaque); - return pj_dealloc(P); -} - - -static void freeup (PJ *P) { - freeup_new (P); - return; -} - - PJ *PROJECTION(aeqd) { struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque)); if (0==Q) @@ -282,7 +290,7 @@ PJ *PROJECTION(aeqd) { P->inv = s_inverse; P->fwd = s_forward; } else { - if (!(Q->en = pj_enfn(P->es))) E_ERROR_0; + if (!(Q->en = pj_enfn(P->es))) return freeup_new(P); 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; diff --git a/src/PJ_airy.c b/src/PJ_airy.c index 3a91986e..2c58bb79 100644 --- a/src/PJ_airy.c +++ b/src/PJ_airy.c @@ -27,74 +27,79 @@ *****************************************************************************/ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(airy, "Airy") "\n\tMisc Sph, no inv.\n\tno_cut lat_b="; struct pj_opaque { - double p_halfpi; - double sinph0; - double cosph0; - double Cb; - int mode; - int no_cut; /* do not cut at hemisphere limit */ + double p_halfpi; + double sinph0; + double cosph0; + double Cb; + int mode; + int no_cut; /* do not cut at hemisphere limit */ }; # define EPS 1.e-10 -# define N_POLE 0 +# define N_POLE 0 # define S_POLE 1 -# define EQUIT 2 -# define OBLIQ 3 +# define EQUIT 2 +# define OBLIQ 3 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double sinlam, coslam, cosphi, sinphi, t, s, Krho, cosz; - - sinlam = sin(lp.lam); - coslam = cos(lp.lam); - switch (Q->mode) { - case EQUIT: - case OBLIQ: - sinphi = sin(lp.phi); - cosphi = cos(lp.phi); - cosz = cosphi * coslam; - if (Q->mode == OBLIQ) - cosz = Q->sinph0 * sinphi + Q->cosph0 * cosz; - if (!Q->no_cut && cosz < -EPS) - F_ERROR; - if (fabs(s = 1. - cosz) > EPS) { - t = 0.5 * (1. + cosz); - Krho = -log(t)/s - Q->Cb / t; - } else - Krho = 0.5 - Q->Cb; - xy.x = Krho * cosphi * sinlam; - if (Q->mode == OBLIQ) - xy.y = Krho * (Q->cosph0 * sinphi - - Q->sinph0 * cosphi * coslam); - else - xy.y = Krho * sinphi; - break; - case S_POLE: - case N_POLE: - lp.phi = fabs(Q->p_halfpi - lp.phi); - if (!Q->no_cut && (lp.phi - EPS) > M_HALFPI) - F_ERROR; - if ((lp.phi *= 0.5) > EPS) { - t = tan(lp.phi); - Krho = -2.*(log(cos(lp.phi)) / t + t * Q->Cb); - xy.x = Krho * sinlam; - xy.y = Krho * coslam; - if (Q->mode == N_POLE) - xy.y = -xy.y; - } else - xy.x = xy.y = 0.; - } - return xy; + double sinlam, coslam, cosphi, sinphi, t, s, Krho, cosz; + + sinlam = sin(lp.lam); + coslam = cos(lp.lam); + switch (Q->mode) { + case EQUIT: + case OBLIQ: + sinphi = sin(lp.phi); + cosphi = cos(lp.phi); + cosz = cosphi * coslam; + if (Q->mode == OBLIQ) + cosz = Q->sinph0 * sinphi + Q->cosph0 * cosz; + if (!Q->no_cut && cosz < -EPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + if (fabs(s = 1. - cosz) > EPS) { + t = 0.5 * (1. + cosz); + Krho = -log(t)/s - Q->Cb / t; + } else + Krho = 0.5 - Q->Cb; + xy.x = Krho * cosphi * sinlam; + if (Q->mode == OBLIQ) + xy.y = Krho * (Q->cosph0 * sinphi - + Q->sinph0 * cosphi * coslam); + else + xy.y = Krho * sinphi; + break; + case S_POLE: + case N_POLE: + lp.phi = fabs(Q->p_halfpi - lp.phi); + if (!Q->no_cut && (lp.phi - EPS) > M_HALFPI) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + if ((lp.phi *= 0.5) > EPS) { + t = tan(lp.phi); + Krho = -2.*(log(cos(lp.phi)) / t + t * Q->Cb); + xy.x = Krho * sinlam; + xy.y = Krho * coslam; + if (Q->mode == N_POLE) + xy.y = -xy.y; + } else + xy.x = xy.y = 0.; + } + return xy; } @@ -116,7 +121,7 @@ static void freeup (PJ *P) { PJ *PROJECTION(airy) { - double beta; + double beta; struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque)); if (0==Q) @@ -124,34 +129,34 @@ PJ *PROJECTION(airy) { P->opaque = Q; - Q->no_cut = pj_param(P->ctx, P->params, "bno_cut").i; - beta = 0.5 * (M_HALFPI - pj_param(P->ctx, P->params, "rlat_b").f); - if (fabs(beta) < EPS) - Q->Cb = -0.5; - else { - Q->Cb = 1./tan(beta); - Q->Cb *= Q->Cb * log(cos(beta)); - } - - if (fabs(fabs(P->phi0) - M_HALFPI) < EPS) - if (P->phi0 < 0.) { - Q->p_halfpi = -M_HALFPI; - Q->mode = S_POLE; - } else { - Q->p_halfpi = M_HALFPI; - Q->mode = N_POLE; - } - else { - if (fabs(P->phi0) < EPS) - Q->mode = EQUIT; - else { - Q->mode = OBLIQ; - Q->sinph0 = sin(P->phi0); - Q->cosph0 = cos(P->phi0); - } - } - P->fwd = s_forward; - P->es = 0.; + Q->no_cut = pj_param(P->ctx, P->params, "bno_cut").i; + beta = 0.5 * (M_HALFPI - pj_param(P->ctx, P->params, "rlat_b").f); + if (fabs(beta) < EPS) + Q->Cb = -0.5; + else { + Q->Cb = 1./tan(beta); + Q->Cb *= Q->Cb * log(cos(beta)); + } + + if (fabs(fabs(P->phi0) - M_HALFPI) < EPS) + if (P->phi0 < 0.) { + Q->p_halfpi = -M_HALFPI; + Q->mode = S_POLE; + } else { + Q->p_halfpi = M_HALFPI; + Q->mode = N_POLE; + } + else { + if (fabs(P->phi0) < EPS) + Q->mode = EQUIT; + else { + Q->mode = OBLIQ; + Q->sinph0 = sin(P->phi0); + Q->cosph0 = cos(P->phi0); + } + } + P->fwd = s_forward; + P->es = 0.; return P; } diff --git a/src/PJ_aitoff.c b/src/PJ_aitoff.c index 77d3f1c6..8b1d7f94 100644 --- a/src/PJ_aitoff.c +++ b/src/PJ_aitoff.c @@ -29,12 +29,13 @@ *****************************************************************************/ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { - double cosphi1; - int mode; + double cosphi1; + int mode; }; @@ -51,18 +52,18 @@ FORWARD(s_forward); /* spheroid */ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double c, d; - - if((d = acos(cos(lp.phi) * cos(c = 0.5 * lp.lam))) != 0.0) {/* basic Aitoff */ - xy.x = 2. * d * cos(lp.phi) * sin(c) * (xy.y = 1. / sin(d)); - xy.y *= d * sin(lp.phi); - } else - xy.x = xy.y = 0.; - if (Q->mode) { /* Winkel Tripel */ - xy.x = (xy.x + lp.lam * Q->cosphi1) * 0.5; - xy.y = (xy.y + lp.phi) * 0.5; - } - return (xy); + double c, d; + + if((d = acos(cos(lp.phi) * cos(c = 0.5 * lp.lam))) != 0.0) {/* basic Aitoff */ + xy.x = 2. * d * cos(lp.phi) * sin(c) * (xy.y = 1. / sin(d)); + xy.y *= d * sin(lp.phi); + } else + xy.x = xy.y = 0.; + if (Q->mode) { /* Winkel Tripel */ + xy.x = (xy.x + lp.lam * Q->cosphi1) * 0.5; + xy.y = (xy.y + lp.phi) * 0.5; + } + return (xy); } /*********************************************************************************** @@ -90,64 +91,64 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; int iter, MAXITER = 10, round = 0, MAXROUND = 20; - double EPSILON = 1e-12, D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y; - - if ((fabs(xy.x) < EPSILON) && (fabs(xy.y) < EPSILON )) { lp.phi = 0.; lp.lam = 0.; return lp; } - - /* initial values for Newton-Raphson method */ - lp.phi = xy.y; lp.lam = xy.x; - do { - iter = 0; - do { - sl = sin(lp.lam * 0.5); cl = cos(lp.lam * 0.5); - sp = sin(lp.phi); cp = cos(lp.phi); - D = cp * cl; - C = 1. - D * D; - D = acos(D) / pow(C, 1.5); - f1 = 2. * D * C * cp * sl; - f2 = D * C * sp; - f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl); - f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp; - f2p = sp * sp * cl / C + D * sl * sl * cp; - f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl); - if (Q->mode) { /* Winkel Tripel */ - f1 = 0.5 * (f1 + lp.lam * Q->cosphi1); - f2 = 0.5 * (f2 + lp.phi); - f1p *= 0.5; - f1l = 0.5 * (f1l + Q->cosphi1); - f2p = 0.5 * (f2p + 1.); - f2l *= 0.5; - } - f1 -= xy.x; f2 -= xy.y; - dl = (f2 * f1p - f1 * f2p) / (dp = f1p * f2l - f2p * f1l); - dp = (f1 * f2l - f2 * f1l) / dp; - dl = fmod(dl, M_PI); /* set to interval [-M_PI, M_PI] */ - lp.phi -= dp; lp.lam -= dl; - } while ((fabs(dp) > EPSILON || fabs(dl) > EPSILON) && (iter++ < MAXITER)); - if (lp.phi > M_PI_2) lp.phi -= 2.*(lp.phi-M_PI_2); /* correct if symmetrical solution for Aitoff */ - if (lp.phi < -M_PI_2) lp.phi -= 2.*(lp.phi+M_PI_2); /* correct if symmetrical solution for Aitoff */ - if ((fabs(fabs(lp.phi) - M_PI_2) < EPSILON) && (!Q->mode)) lp.lam = 0.; /* if pole in Aitoff, return longitude of 0 */ - - /* calculate x,y coordinates with solution obtained */ - if((D = acos(cos(lp.phi) * cos(C = 0.5 * lp.lam))) != 0.0) {/* Aitoff */ - x = 2. * D * cos(lp.phi) * sin(C) * (y = 1. / sin(D)); - y *= D * sin(lp.phi); - } else - x = y = 0.; - if (Q->mode) { /* Winkel Tripel */ - x = (x + lp.lam * Q->cosphi1) * 0.5; - y = (y + lp.phi) * 0.5; - } - /* if too far from given values of x,y, repeat with better approximation of phi,lam */ - } while (((fabs(xy.x-x) > EPSILON) || (fabs(xy.y-y) > EPSILON)) && (round++ < MAXROUND)); - - if (iter == MAXITER && round == MAXROUND) + double EPSILON = 1e-12, D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y; + + if ((fabs(xy.x) < EPSILON) && (fabs(xy.y) < EPSILON )) { lp.phi = 0.; lp.lam = 0.; return lp; } + + /* initial values for Newton-Raphson method */ + lp.phi = xy.y; lp.lam = xy.x; + do { + iter = 0; + do { + sl = sin(lp.lam * 0.5); cl = cos(lp.lam * 0.5); + sp = sin(lp.phi); cp = cos(lp.phi); + D = cp * cl; + C = 1. - D * D; + D = acos(D) / pow(C, 1.5); + f1 = 2. * D * C * cp * sl; + f2 = D * C * sp; + f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl); + f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp; + f2p = sp * sp * cl / C + D * sl * sl * cp; + f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl); + if (Q->mode) { /* Winkel Tripel */ + f1 = 0.5 * (f1 + lp.lam * Q->cosphi1); + f2 = 0.5 * (f2 + lp.phi); + f1p *= 0.5; + f1l = 0.5 * (f1l + Q->cosphi1); + f2p = 0.5 * (f2p + 1.); + f2l *= 0.5; + } + f1 -= xy.x; f2 -= xy.y; + dl = (f2 * f1p - f1 * f2p) / (dp = f1p * f2l - f2p * f1l); + dp = (f1 * f2l - f2 * f1l) / dp; + dl = fmod(dl, M_PI); /* set to interval [-M_PI, M_PI] */ + lp.phi -= dp; lp.lam -= dl; + } while ((fabs(dp) > EPSILON || fabs(dl) > EPSILON) && (iter++ < MAXITER)); + if (lp.phi > M_PI_2) lp.phi -= 2.*(lp.phi-M_PI_2); /* correct if symmetrical solution for Aitoff */ + if (lp.phi < -M_PI_2) lp.phi -= 2.*(lp.phi+M_PI_2); /* correct if symmetrical solution for Aitoff */ + if ((fabs(fabs(lp.phi) - M_PI_2) < EPSILON) && (!Q->mode)) lp.lam = 0.; /* if pole in Aitoff, return longitude of 0 */ + + /* calculate x,y coordinates with solution obtained */ + if((D = acos(cos(lp.phi) * cos(C = 0.5 * lp.lam))) != 0.0) {/* Aitoff */ + x = 2. * D * cos(lp.phi) * sin(C) * (y = 1. / sin(D)); + y *= D * sin(lp.phi); + } else + x = y = 0.; + if (Q->mode) { /* Winkel Tripel */ + x = (x + lp.lam * Q->cosphi1) * 0.5; + y = (y + lp.phi) * 0.5; + } + /* if too far from given values of x,y, repeat with better approximation of phi,lam */ + } while (((fabs(xy.x-x) > EPSILON) || (fabs(xy.y-y) > EPSILON)) && (round++ < MAXROUND)); + + if (iter == MAXITER && round == MAXROUND) { pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT ); /* fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl); */ } - return lp; + return lp; } @@ -168,10 +169,10 @@ static void freeup (PJ *P) { } static PJ *setup(PJ *P) { - P->inv = s_inverse; - P->fwd = s_forward; - P->es = 0.; - return P; + P->inv = s_inverse; + P->fwd = s_forward; + P->es = 0.; + return P; } @@ -181,7 +182,7 @@ PJ *PROJECTION(aitoff) { return freeup_new (P); P->opaque = Q; - Q->mode = 0; + Q->mode = 0; return setup(P); } @@ -192,13 +193,15 @@ PJ *PROJECTION(wintri) { return freeup_new (P); P->opaque = Q; - Q->mode = 1; - if (pj_param(P->ctx, P->params, "tlat_1").i) { - if ((Q->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f)) == 0.) - E_ERROR(-22) + Q->mode = 1; + if (pj_param(P->ctx, P->params, "tlat_1").i) { + if ((Q->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f)) == 0.) { + proj_errno_set(P, PJD_ERR_LAT_LARGER_THAN_90); + return freeup_new(P); + } } - else /* 50d28' or acos(2/pi) */ - Q->cosphi1 = 0.636619772367581343; + else /* 50d28' or acos(2/pi) */ + Q->cosphi1 = 0.636619772367581343; return setup(P); } diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c index 2d015de1..97284fcc 100644 --- a/src/PJ_bipc.c +++ b/src/PJ_bipc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(bipc, "Bipolar conic of western hemisphere") "\n\tConic Sph."; @@ -51,7 +52,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ sdlam = sin(sdlam); z = S20 * sphi + C20 * cphi * cdlam; if (fabs(z) > 1.) { - if (fabs(z) > ONEEPS) F_ERROR + if (fabs(z) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } else z = z < 0. ? -1. : 1.; } else z = acos(z); @@ -62,19 +66,31 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ } else { z = S45 * (sphi + cphi * cdlam); if (fabs(z) > 1.) { - if (fabs(z) > ONEEPS) F_ERROR + if (fabs(z) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } else z = z < 0. ? -1. : 1.; } else z = acos(z); Av = Azba; xy.y = -rhoc; } - if (z < 0.) F_ERROR; + if (z < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } r = F * (t = pow(tan(.5 * z), n)); - if ((al = .5 * (R104 - z)) < 0.) F_ERROR; + if ((al = .5 * (R104 - z)) < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } al = (t + pow(al, n)) / T; if (fabs(al) > 1.) { - if (fabs(al) > ONEEPS) F_ERROR + if (fabs(al) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } else al = al < 0. ? -1. : 1.; } else al = acos(al); @@ -125,7 +141,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ break; rl = r; } - if (! i) I_ERROR; + if (! i) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } Az = Av - Az / n; lp.phi = asin(s * cos(z) + c * sin(z) * cos(Az)); lp.lam = atan2(sin(Az), c / tan(z) - s * cos(Az)); diff --git a/src/PJ_bonne.c b/src/PJ_bonne.c index ca076f56..2a576c60 100644 --- a/src/PJ_bonne.c +++ b/src/PJ_bonne.c @@ -1,78 +1,85 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)") - "\n\tConic Sph&Ell\n\tlat_1="; -#define EPS10 1e-10 + "\n\tConic Sph&Ell\n\tlat_1="; +#define EPS10 1e-10 struct pj_opaque { - double phi1; - double cphi1; - double am1; - double m1; - double *en; + double phi1; + double cphi1; + double am1; + double m1; + double *en; }; static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double rh, E, c; + double rh, E, c; - rh = Q->am1 + Q->m1 - pj_mlfn(lp.phi, E = sin(lp.phi), c = cos(lp.phi), Q->en); - E = c * lp.lam / (rh * sqrt(1. - P->es * E * E)); - xy.x = rh * sin(E); - xy.y = Q->am1 - rh * cos(E); - return xy; + rh = Q->am1 + Q->m1 - pj_mlfn(lp.phi, E = sin(lp.phi), c = cos(lp.phi), Q->en); + E = c * lp.lam / (rh * sqrt(1. - P->es * E * E)); + xy.x = rh * sin(E); + xy.y = Q->am1 - rh * cos(E); + 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 E, rh; - - rh = Q->cphi1 + Q->phi1 - lp.phi; - if (fabs(rh) > EPS10) { - xy.x = rh * sin(E = lp.lam * cos(lp.phi) / rh); - xy.y = Q->cphi1 - rh * cos(E); - } else - xy.x = xy.y = 0.; - return xy; + double E, rh; + + rh = Q->cphi1 + Q->phi1 - lp.phi; + if (fabs(rh) > EPS10) { + xy.x = rh * sin(E = lp.lam * cos(lp.phi) / rh); + xy.y = Q->cphi1 - rh * cos(E); + } else + xy.x = xy.y = 0.; + return xy; } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double rh; - - rh = hypot(xy.x, xy.y = Q->cphi1 - xy.y); - lp.phi = Q->cphi1 + Q->phi1 - rh; - if (fabs(lp.phi) > M_HALFPI) I_ERROR; - if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) - lp.lam = 0.; - else - lp.lam = rh * atan2(xy.x, xy.y) / cos(lp.phi); - return lp; + double rh; + + rh = hypot(xy.x, xy.y = Q->cphi1 - xy.y); + lp.phi = Q->cphi1 + Q->phi1 - rh; + if (fabs(lp.phi) > M_HALFPI) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) + lp.lam = 0.; + else + lp.lam = rh * atan2(xy.x, xy.y) / 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; - double s, rh; - - rh = hypot(xy.x, xy.y = Q->am1 - xy.y); - lp.phi = pj_inv_mlfn(P->ctx, Q->am1 + Q->m1 - rh, P->es, Q->en); - if ((s = fabs(lp.phi)) < M_HALFPI) { - s = sin(lp.phi); - lp.lam = rh * atan2(xy.x, xy.y) * - sqrt(1. - P->es * s * s) / cos(lp.phi); - } else if (fabs(s - M_HALFPI) <= EPS10) - lp.lam = 0.; - else I_ERROR; - return lp; + double s, rh; + + rh = hypot(xy.x, xy.y = Q->am1 - xy.y); + lp.phi = pj_inv_mlfn(P->ctx, Q->am1 + Q->m1 - rh, P->es, Q->en); + if ((s = fabs(lp.phi)) < M_HALFPI) { + s = sin(lp.phi); + lp.lam = rh * atan2(xy.x, xy.y) * + sqrt(1. - P->es * s * s) / cos(lp.phi); + } else if (fabs(s - M_HALFPI) <= EPS10) + lp.lam = 0.; + else { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + return lp; } @@ -94,29 +101,32 @@ static void freeup (PJ *P) { PJ *PROJECTION(bonne) { - double c; + double c; struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque)); if (0==Q) return freeup_new (P); P->opaque = Q; - Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f; - if (fabs(Q->phi1) < EPS10) E_ERROR(-23); - if (P->es != 0.0) { - Q->en = pj_enfn(P->es); - Q->m1 = pj_mlfn(Q->phi1, Q->am1 = sin(Q->phi1), - c = cos(Q->phi1), Q->en); - Q->am1 = c / (sqrt(1. - P->es * Q->am1 * Q->am1) * Q->am1); - P->inv = e_inverse; - P->fwd = e_forward; - } else { - if (fabs(Q->phi1) + EPS10 >= M_HALFPI) - Q->cphi1 = 0.; - else - Q->cphi1 = 1. / tan(Q->phi1); - P->inv = s_inverse; - P->fwd = s_forward; - } + Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f; + if (fabs(Q->phi1) < EPS10) { + proj_errno_set(P, PJD_ERR_LAT1_IS_ZERO); + return freeup_new(P); + } + if (P->es != 0.0) { + Q->en = pj_enfn(P->es); + Q->m1 = pj_mlfn(Q->phi1, Q->am1 = sin(Q->phi1), + c = cos(Q->phi1), Q->en); + Q->am1 = c / (sqrt(1. - P->es * Q->am1 * Q->am1) * Q->am1); + P->inv = e_inverse; + P->fwd = e_forward; + } else { + if (fabs(Q->phi1) + EPS10 >= M_HALFPI) + Q->cphi1 = 0.; + else + Q->cphi1 = 1. / tan(Q->phi1); + P->inv = s_inverse; + P->fwd = s_forward; + } return P; } diff --git a/src/PJ_calcofi.c b/src/PJ_calcofi.c index ebffd4f5..25521eed 100644 --- a/src/PJ_calcofi.c +++ b/src/PJ_calcofi.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(calcofi, "Cal Coop Ocean Fish Invest Lines/Stations") "\n\tCyl, Sph&Ell"; @@ -46,7 +47,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ /* if the user has specified +lon_0 or +k0 for some reason, we're going to ignore it so that xy is consistent with point O */ lp.lam = lp.lam + P->lam0; - if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) F_ERROR; + if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = lp.lam; xy.y = -log(pj_tsfn(lp.phi, sin(lp.phi), P->e)); /* Mercator transform xy*/ oy = -log(pj_tsfn(PT_O_PHI, sin(PT_O_PHI), P->e)); @@ -74,7 +78,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ double l2; double ry; lp.lam = lp.lam + P->lam0; - if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) F_ERROR; + if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = lp.lam; xy.y = log(tan(M_FORTPI + .5 * lp.phi)); oy = log(tan(M_FORTPI + .5 * PT_O_PHI)); diff --git a/src/PJ_cc.c b/src/PJ_cc.c index 1b5d3152..d43f5a88 100644 --- a/src/PJ_cc.c +++ b/src/PJ_cc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(cc, "Central Cylindrical") "\n\tCyl, Sph"; #define EPS10 1.e-10 @@ -7,19 +8,22 @@ PROJ_HEAD(cc, "Central Cylindrical") "\n\tCyl, Sph"; static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; - if (fabs (fabs(lp.phi) - M_HALFPI) <= EPS10) F_ERROR; - xy.x = lp.lam; - xy.y = tan(lp.phi); - return xy; + if (fabs (fabs(lp.phi) - M_HALFPI) <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + xy.x = lp.lam; + xy.y = tan(lp.phi); + return xy; } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; - (void) P; - lp.phi = atan(xy.y); - lp.lam = xy.x; - return lp; + (void) P; + lp.phi = atan(xy.y); + lp.lam = xy.x; + return lp; } diff --git a/src/PJ_cea.c b/src/PJ_cea.c index d6dcd471..20f03547 100644 --- a/src/PJ_cea.c +++ b/src/PJ_cea.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double qp; @@ -44,7 +45,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ else lp.phi = asin(xy.y); lp.lam = xy.x / P->k0; - } else I_ERROR; + } else { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } return (lp); } @@ -76,14 +80,18 @@ PJ *PROJECTION(cea) { if (pj_param(P->ctx, P->params, "tlat_ts").i) { P->k0 = cos(t = pj_param(P->ctx, P->params, "rlat_ts").f); if (P->k0 < 0.) { - E_ERROR(-24); + proj_errno_set(P, PJD_ERR_LAT_TS_LARGER_THAN_90); + freeup_new(P); + return 0; } } if (P->es != 0.0) { t = sin(t); P->k0 /= sqrt(1. - P->es * t * t); P->e = sqrt(P->es); - if (!(Q->apa = pj_authset(P->es))) E_ERROR_0; + if (!(Q->apa = pj_authset(P->es))) { + return freeup_new(P); + } Q->qp = pj_qsfn(1., P->e, P->one_es); P->inv = e_inverse; P->fwd = e_forward; diff --git a/src/PJ_chamb.c b/src/PJ_chamb.c index e9e9848b..c6028a7a 100644 --- a/src/PJ_chamb.c +++ b/src/PJ_chamb.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" typedef struct { double r, Az; } VECT; struct pj_opaque { @@ -130,7 +131,10 @@ PJ *PROJECTION(chamb) { j = i == 2 ? 0 : i + 1; Q->c[i].v = vect(P->ctx,Q->c[j].phi - Q->c[i].phi, Q->c[i].cosphi, Q->c[i].sinphi, Q->c[j].cosphi, Q->c[j].sinphi, Q->c[j].lam - Q->c[i].lam); - if (Q->c[i].v.r == 0.0) E_ERROR(-25); + if (Q->c[i].v.r == 0.0) { + proj_errno_set(P, PJD_ERR_CONTROL_POINT_NO_DIST); + return freeup_new(P); + } /* co-linearity problem ignored for now */ } Q->beta_0 = lc(P->ctx,Q->c[0].v.r, Q->c[2].v.r, Q->c[1].v.r); diff --git a/src/PJ_collg.c b/src/PJ_collg.c index baac1af9..c646d99a 100644 --- a/src/PJ_collg.c +++ b/src/PJ_collg.c @@ -1,37 +1,43 @@ #define PJ_LIB__ -# include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(collg, "Collignon") "\n\tPCyl, Sph."; -#define FXC 1.12837916709551257390 -#define FYC 1.77245385090551602729 -#define ONEEPS 1.0000001 +#define FXC 1.12837916709551257390 +#define FYC 1.77245385090551602729 +#define ONEEPS 1.0000001 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; - (void) P; - if ((xy.y = 1. - sin(lp.phi)) <= 0.) - xy.y = 0.; - else - xy.y = sqrt(xy.y); - xy.x = FXC * lp.lam * xy.y; - xy.y = FYC * (1. - xy.y); - return (xy); + (void) P; + if ((xy.y = 1. - sin(lp.phi)) <= 0.) + xy.y = 0.; + else + xy.y = sqrt(xy.y); + xy.x = FXC * lp.lam * xy.y; + xy.y = FYC * (1. - xy.y); + return (xy); } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; - lp.phi = xy.y / FYC - 1.; - if (fabs(lp.phi = 1. - lp.phi * lp.phi) < 1.) - lp.phi = asin(lp.phi); - else if (fabs(lp.phi) > ONEEPS) I_ERROR - else lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI; - if ((lp.lam = 1. - sin(lp.phi)) <= 0.) - lp.lam = 0.; - else - lp.lam = xy.x / (FXC * sqrt(lp.lam)); - return (lp); + lp.phi = xy.y / FYC - 1.; + if (fabs(lp.phi = 1. - lp.phi * lp.phi) < 1.) + lp.phi = asin(lp.phi); + else if (fabs(lp.phi) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else { + lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI; + } + + if ((lp.lam = 1. - sin(lp.phi)) <= 0.) + lp.lam = 0.; + else + lp.lam = xy.x / (FXC * sqrt(lp.lam)); + return (lp); } diff --git a/src/PJ_eck2.c b/src/PJ_eck2.c index 763039a9..6ef2a96c 100644 --- a/src/PJ_eck2.c +++ b/src/PJ_eck2.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -# include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(eck2, "Eckert II") "\n\tPCyl. Sph."; @@ -28,9 +29,12 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ lp.lam = xy.x / (FXC * ( lp.phi = 2. - fabs(xy.y) / FYC) ); lp.phi = (4. - lp.phi * lp.phi) * C13; if (fabs(lp.phi) >= 1.) { - if (fabs(lp.phi) > ONEEPS) I_ERROR - else + if (fabs(lp.phi) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else { lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI; + } } else lp.phi = asin(lp.phi); if (xy.y < 0) diff --git a/src/PJ_eqc.c b/src/PJ_eqc.c index 6f0cfcf1..576e0e3a 100644 --- a/src/PJ_eqc.c +++ b/src/PJ_eqc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double rc; @@ -54,7 +55,10 @@ PJ *PROJECTION(eqc) { return freeup_new (P); P->opaque = Q; - if ((Q->rc = cos(pj_param(P->ctx, P->params, "rlat_ts").f)) <= 0.) E_ERROR(-24); + if ((Q->rc = cos(pj_param(P->ctx, P->params, "rlat_ts").f)) <= 0.) { + proj_errno_set(P, PJD_ERR_LAT_TS_LARGER_THAN_90); + return freeup_new(P); + } P->inv = s_inverse; P->fwd = s_forward; P->es = 0.; diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c index a83abede..eaf4db0b 100644 --- a/src/PJ_eqdc.c +++ b/src/PJ_eqdc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double phi1; @@ -96,9 +97,16 @@ PJ *PROJECTION(eqdc) { Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f; Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f; - if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21); + + if (fabs(Q->phi1 + Q->phi2) < EPS10) { + proj_errno_set(P, PJD_ERR_CONIC_LAT_EQUAL); + freeup_new(P); + return 0; + } + if (!(Q->en = pj_enfn(P->es))) - E_ERROR_0; + return freeup_new(P); + Q->n = sinphi = sin(Q->phi1); cosphi = cos(Q->phi1); secant = fabs(Q->phi1 - Q->phi2) >= EPS10; @@ -122,6 +130,7 @@ PJ *PROJECTION(eqdc) { Q->c = Q->phi1 + cos(Q->phi1) / Q->n; Q->rho0 = Q->c - P->phi0; } + P->inv = e_inverse; P->fwd = e_forward; P->spc = special; diff --git a/src/PJ_fouc_s.c b/src/PJ_fouc_s.c index f5df5d84..32eeeb4f 100644 --- a/src/PJ_fouc_s.c +++ b/src/PJ_fouc_s.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(fouc_s, "Foucaut Sinusoidal") "\n\tPCyl., Sph."; @@ -71,8 +72,10 @@ PJ *PROJECTION(fouc_s) { P->opaque = Q; Q->n = pj_param(P->ctx, P->params, "dn").f; - if (Q->n < 0. || Q->n > 1.) - E_ERROR(-99) + if (Q->n < 0. || Q->n > 1.) { + proj_errno_set(P, PJD_ERR_N_OUT_OF_RANGE); + return freeup_new(P); + } Q->n1 = 1. - Q->n; P->es = 0; P->inv = s_inverse; diff --git a/src/PJ_geos.c b/src/PJ_geos.c index 56357c2b..b929c06b 100644 --- a/src/PJ_geos.c +++ b/src/PJ_geos.c @@ -28,7 +28,8 @@ */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double h; @@ -91,8 +92,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ Vz = r * sin (lp.phi); /* Check visibility. */ - if (((Q->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * Q->radius_p_inv2) < 0.) - F_ERROR; + if (((Q->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * Q->radius_p_inv2) < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } /* Calculation based on view angles from satellite. */ tmp = Q->radius_g - Vx; @@ -127,7 +130,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ /* Calculation of terms in cubic equation and determinant.*/ a = Vy * Vy + Vz * Vz + Vx * Vx; b = 2 * Q->radius_g * Vx; - if ((det = (b * b) - 4 * a * Q->C) < 0.) I_ERROR; + if ((det = (b * b) - 4 * a * Q->C) < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } /* Calculation of three components of vector from satellite to position.*/ k = (-b - sqrt(det)) / (2 * a); @@ -163,7 +169,10 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ a = Vz / Q->radius_p; a = Vy * Vy + a * a + Vx * Vx; b = 2 * Q->radius_g * Vx; - if ((det = (b * b) - 4 * a * Q->C) < 0.) I_ERROR; + if ((det = (b * b) - 4 * a * Q->C) < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } /* Calculation of three components of vector from satellite to position.*/ k = (-b - sqrt(det)) / (2. * a); @@ -202,18 +211,24 @@ PJ *PROJECTION(geos) { return freeup_new (P); P->opaque = Q; - if ((Q->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30); + if ((Q->h = pj_param(P->ctx, P->params, "dh").f) <= 0.){ + proj_errno_set(P, PJD_ERR_H_LESS_THAN_ZERO); + return freeup_new(P); + } - if (P->phi0 != 0.0) E_ERROR(-46); + if (P->phi0 != 0.0) { + proj_errno_set(P, PJD_ERR_UNKNOWN_PRIME_MERIDIAN); + return freeup_new(P); + } Q->sweep_axis = pj_param(P->ctx, P->params, "ssweep").s; if (Q->sweep_axis == NULL) Q->flip_axis = 0; else { - if (Q->sweep_axis[1] != '\0' || - (Q->sweep_axis[0] != 'x' && - Q->sweep_axis[0] != 'y')) - E_ERROR(-49); + if (Q->sweep_axis[1] != '\0' || (Q->sweep_axis[0] != 'x' && Q->sweep_axis[0] != 'y')) { + proj_errno_set(P, PJD_ERR_INVALID_SWEEP_AXIS); + return freeup_new(P); + } if (Q->sweep_axis[0] == 'x') Q->flip_axis = 1; else diff --git a/src/PJ_gn_sinu.c b/src/PJ_gn_sinu.c index 68a4c936..57bbecc8 100644 --- a/src/PJ_gn_sinu.c +++ b/src/PJ_gn_sinu.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(gn_sinu, "General Sinusoidal Series") "\n\tPCyl, Sph.\n\tm= n="; PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)") "\n\tPCyl, Sph&Ell"; @@ -36,7 +37,7 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ } else if ((s - EPS10) < M_HALFPI) { lp.lam = 0.; } else { - I_ERROR; + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); } return lp; @@ -60,8 +61,11 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ if (fabs(V) < LOOP_TOL) break; } - if (!i) - F_ERROR + if (!i) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + } xy.x = Q->C_x * lp.lam * (Q->m + cos(lp.phi)); xy.y = Q->C_y * lp.phi; diff --git a/src/PJ_gnom.c b/src/PJ_gnom.c index a925d43e..099f32fb 100644 --- a/src/PJ_gnom.c +++ b/src/PJ_gnom.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph."; @@ -40,7 +41,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ break; } - if (xy.y <= EPS10) F_ERROR; + if (xy.y <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam); switch (Q->mode) { diff --git a/src/PJ_goode.c b/src/PJ_goode.c index fd63fa83..fff12a78 100644 --- a/src/PJ_goode.c +++ b/src/PJ_goode.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(goode, "Goode Homolosine") "\n\tPCyl, Sph."; @@ -71,12 +72,12 @@ PJ *PROJECTION(goode) { P->es = 0.; if (!(Q->sinu = pj_sinu(0)) || !(Q->moll = pj_moll(0))) - E_ERROR_0; + return freeup_new(P); Q->sinu->es = 0.; Q->sinu->ctx = P->ctx; Q->moll->ctx = P->ctx; if (!(Q->sinu = pj_sinu(Q->sinu)) || !(Q->moll = pj_moll(Q->moll))) - E_ERROR_0; + return freeup_new(P); P->fwd = s_forward; P->inv = s_inverse; diff --git a/src/PJ_hammer.c b/src/PJ_hammer.c index 891baa00..58182398 100644 --- a/src/PJ_hammer.c +++ b/src/PJ_hammer.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff") "\n\tMisc Sph, \n\tW= M="; @@ -33,7 +34,7 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ if (fabs(2.*z*z-1.) < EPS) { lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; - pj_errno = -14; + proj_errno_set(P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT); } else { lp.lam = aatan2(Q->w * xy.x * z,2. * z * z - 1)/Q->w; lp.phi = aasin(P->ctx,z * xy.y); @@ -65,11 +66,17 @@ PJ *PROJECTION(hammer) { P->opaque = Q; if (pj_param(P->ctx, P->params, "tW").i) { - if ((Q->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) E_ERROR(-27); + if ((Q->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) { + proj_errno_set(P, PJD_ERR_W_OR_M_ZERO_OR_LESS); + return freeup_new(P); + } } else Q->w = .5; if (pj_param(P->ctx, P->params, "tM").i) { - if ((Q->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) E_ERROR(-27); + if ((Q->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) { + proj_errno_set(P, PJD_ERR_W_OR_M_ZERO_OR_LESS); + return freeup_new(P); + } } else Q->m = 1.; diff --git a/src/PJ_hatano.c b/src/PJ_hatano.c index 9e681a3d..d75a96a3 100644 --- a/src/PJ_hatano.c +++ b/src/PJ_hatano.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph."; @@ -43,7 +44,8 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ th = xy.y * ( xy.y < 0. ? RYCS : RYCN); if (fabs(th) > 1.) { if (fabs(th) > ONETOL) { - I_ERROR; + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; } else { th = th > 0. ? M_HALFPI : - M_HALFPI; } @@ -56,7 +58,8 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN); if (fabs(lp.phi) > 1.) { if (fabs(lp.phi) > ONETOL) { - I_ERROR; + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; } else { lp.phi = lp.phi > 0. ? M_HALFPI : - M_HALFPI; } diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c index 7e1245a6..ef53e58f 100644 --- a/src/PJ_healpix.c +++ b/src/PJ_healpix.c @@ -29,7 +29,8 @@ * SOFTWARE. *****************************************************************************/ # define PJ_LIB__ -# include <projects.h> +# include <proj.h> +# include "projects.h" PROJ_HEAD(healpix, "HEALPix") "\n\tSph., Ellps."; PROJ_HEAD(rhealpix, "rHEALPix") "\n\tSph., Ellps.\n\tnorth_square= south_square="; @@ -649,10 +650,12 @@ PJ *PROJECTION(rhealpix) { /* Check for valid north_square and south_square inputs. */ if (Q->north_square < 0 || Q->north_square > 3) { - E_ERROR(-47); + proj_errno_set(P, PJD_ERR_AXIS); + return freeup_new(P); } if (Q->south_square < 0 || Q->south_square > 3) { - E_ERROR(-47); + proj_errno_set(P, PJD_ERR_AXIS); + return freeup_new(P); } if (P->es != 0.0) { Q->apa = pj_authset(P->es); /* For auth_lat(). */ diff --git a/src/PJ_igh.c b/src/PJ_igh.c index 5a019a9a..9b5f7075 100644 --- a/src/PJ_igh.c +++ b/src/PJ_igh.c @@ -1,5 +1,5 @@ #define PJ_LIB__ -#include <projects.h> +#include "projects.h" PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph."; @@ -175,8 +175,8 @@ static void freeup (PJ *P) { */ #define SETUP(n, proj, x_0, y_0, lon_0) \ - if (!(Q->pj[n-1] = pj_##proj(0))) E_ERROR_0; \ - if (!(Q->pj[n-1] = pj_##proj(Q->pj[n-1]))) E_ERROR_0; \ + if (!(Q->pj[n-1] = pj_##proj(0))) return freeup_new(P); \ + if (!(Q->pj[n-1] = pj_##proj(Q->pj[n-1]))) return freeup_new(P); \ Q->pj[n-1]->ctx = P->ctx; \ Q->pj[n-1]->x0 = x_0; \ Q->pj[n-1]->y0 = y_0; \ diff --git a/src/PJ_imw_p.c b/src/PJ_imw_p.c index 0820dd49..0b2a6602 100644 --- a/src/PJ_imw_p.c +++ b/src/PJ_imw_p.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(imw_p, "International Map of the World Polyconic") "\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]"; @@ -27,7 +28,7 @@ static int phi12(PJ *P, double *del, double *sig) { Q->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f; *del = 0.5 * (Q->phi_2 - Q->phi_1); *sig = 0.5 * (Q->phi_2 + Q->phi_1); - err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0; + err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? PJD_ERR_ABS_LAT1_EQ_ABS_LAT2 : 0; } return err; } @@ -155,9 +156,11 @@ PJ *PROJECTION(imw_p) { return freeup_new (P); P->opaque = Q; - if (!(Q->en = pj_enfn(P->es))) E_ERROR_0; - if( (i = phi12(P, &del, &sig)) != 0) - E_ERROR(i); + if (!(Q->en = pj_enfn(P->es))) return freeup_new(P); + if( (i = phi12(P, &del, &sig)) != 0) { + proj_errno_set(P, i); + return freeup_new(P); + } if (Q->phi_2 < Q->phi_1) { /* make sure P->phi_1 most southerly */ del = Q->phi_1; Q->phi_1 = Q->phi_2; diff --git a/src/PJ_isea.c b/src/PJ_isea.c index 6bc6a66c..deaecb0f 100644 --- a/src/PJ_isea.c +++ b/src/PJ_isea.c @@ -1034,7 +1034,8 @@ isea_forward(struct isea_dgg *g, struct isea_geo *in) */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph"; @@ -1099,7 +1100,8 @@ PJ *PROJECTION(isea) { } else if (!strcmp(opt, "pole")) { isea_orient_pole(&Q->dgg); } else { - E_ERROR(-34); + proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + return freeup_new(P); } } @@ -1138,7 +1140,8 @@ PJ *PROJECTION(isea) { } else { /* TODO verify error code. Possibly eliminate magic */ - E_ERROR(-34); + proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + return freeup_new(P); } } diff --git a/src/PJ_laea.c b/src/PJ_laea.c index f4976cab..0a638fb0 100644 --- a/src/PJ_laea.c +++ b/src/PJ_laea.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell"; @@ -55,7 +56,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ q = Q->qp + q; break; } - if (fabs(b) < EPS10) F_ERROR; + if (fabs(b) < EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } switch (Q->mode) { case OBLIQ: @@ -98,7 +102,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ case OBLIQ: xy.y = 1. + Q->sinb1 * sinphi + Q->cosb1 * cosphi * coslam; oblcon: - if (xy.y <= EPS10) F_ERROR; + if (xy.y <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.y = sqrt(2. / xy.y); xy.x = xy.y * cosphi * sin(lp.lam); xy.y *= Q->mode == EQUIT ? sinphi : @@ -108,7 +115,10 @@ oblcon: coslam = -coslam; /*-fallthrough*/ case S_POLE: - if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR; + if (fabs(lp.phi + P->phi0) < EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.y = M_FORTPI - lp.phi * .5; xy.y = 2. * (Q->mode == S_POLE ? cos(xy.y) : sin(xy.y)); xy.x = xy.y * sin(lp.lam); @@ -174,7 +184,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ double cosz=0.0, rh, sinz=0.0; rh = hypot(xy.x, xy.y); - if ((lp.phi = rh * .5 ) > 1.) I_ERROR; + if ((lp.phi = rh * .5 ) > 1.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } lp.phi = 2. * asin(lp.phi); if (Q->mode == OBLIQ || Q->mode == EQUIT) { sinz = sin(lp.phi); diff --git a/src/PJ_lagrng.c b/src/PJ_lagrng.c index 142c28ab..98500900 100644 --- a/src/PJ_lagrng.c +++ b/src/PJ_lagrng.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(lagrng, "Lagrange") "\n\tMisc Sph, no inv.\n\tW="; @@ -23,8 +24,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ } else { lp.phi = sin(lp.phi); v = Q->a1 * pow((1. + lp.phi)/(1. - lp.phi), Q->hrw); - if ((c = 0.5 * (v + 1./v) + cos(lp.lam *= Q->rw)) < TOL) - F_ERROR; + if ((c = 0.5 * (v + 1./v) + cos(lp.lam *= Q->rw)) < TOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = 2. * sin(lp.lam) / c; xy.y = (v - 1./v) / c; } @@ -56,12 +59,18 @@ PJ *PROJECTION(lagrng) { P->opaque = Q; Q->rw = pj_param(P->ctx, P->params, "dW").f; - if (Q->rw <= 0) E_ERROR(-27); + if (Q->rw <= 0) { + proj_errno_set(P, PJD_ERR_W_OR_M_ZERO_OR_LESS); + return freeup_new(P); + } Q->rw = 1. / Q->rw; Q->hrw = 0.5 * Q->rw; phi1 = sin(pj_param(P->ctx, P->params, "rlat_1").f); - if (fabs(fabs(phi1) - 1.) < TOL) E_ERROR(-22); + if (fabs(fabs(phi1) - 1.) < TOL) { + proj_errno_set(P, PJD_ERR_LAT_LARGER_THAN_90); + return freeup_new(P); + } Q->a1 = pow((1. - phi1)/(1. + phi1), Q->hrw); diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c index f1f0b324..905f7086 100644 --- a/src/PJ_lcc.c +++ b/src/PJ_lcc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(lcc, "Lambert Conformal Conic") "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0"; @@ -22,7 +23,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ double rho; if (fabs(fabs(lp.phi) - M_HALFPI) < EPS10) { - if ((lp.phi * Q->n) <= 0.) F_ERROR; + if ((lp.phi * Q->n) <= 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } rho = 0.; } else { rho = Q->c * (Q->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi), @@ -53,8 +57,11 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ } if (Q->ellips) { lp.phi = pj_phi2(P->ctx, pow(rho / Q->c, 1./Q->n), P->e); - if (lp.phi == HUGE_VAL) - I_ERROR; + if (lp.phi == HUGE_VAL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + } else lp.phi = 2. * atan(pow(Q->c / rho, 1./Q->n)) - M_HALFPI; lp.lam = atan2(xy.x, xy.y) / Q->n; @@ -115,7 +122,10 @@ PJ *PROJECTION(lcc) { if (!pj_param(P->ctx, P->params, "tlat_0").i) P->phi0 = Q->phi1; } - if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21); + if (fabs(Q->phi1 + Q->phi2) < EPS10) { + proj_errno_set(P, PJD_ERR_CONIC_LAT_EQUAL); + return freeup_new(P); + } Q->n = sinphi = sin(Q->phi1); cosphi = cos(Q->phi1); secant = fabs(Q->phi1 - Q->phi2) >= EPS10; diff --git a/src/PJ_lcca.c b/src/PJ_lcca.c index 4a586308..77704375 100644 --- a/src/PJ_lcca.c +++ b/src/PJ_lcca.c @@ -1,7 +1,8 @@ /* PROJ.4 Cartographic Projection System */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative") "\n\tConic, Sph&Ell\n\tlat_0="; @@ -58,7 +59,10 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ S -= (dif = (fS(S, Q->C) - dr) / fSp(S, Q->C)); if (fabs(dif) < DEL_TOL) break; } - if (!i) I_ERROR + if (!i) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } lp.phi = pj_inv_mlfn(P->ctx, S + Q->M0, P->es, Q->en); return lp; @@ -90,9 +94,13 @@ PJ *PROJECTION(lcca) { P->opaque = Q; (Q->en = pj_enfn(P->es)); - if (!Q->en) E_ERROR_0; - if (!pj_param(P->ctx, P->params, "tlat_0").i) E_ERROR(50); - if (P->phi0 == 0.) E_ERROR(51); + if (!Q->en) + return freeup_new(P); + + if (P->phi0 == 0.) { + proj_errno_set(P, PJD_ERR_LAT_0_IS_ZERO); + return freeup_new(P); + } Q->l = sin(P->phi0); Q->M0 = pj_mlfn(P->phi0, Q->l, cos(P->phi0), Q->en); s2p0 = Q->l * Q->l; diff --git a/src/PJ_loxim.c b/src/PJ_loxim.c index 57f07048..6cd56eef 100644 --- a/src/PJ_loxim.c +++ b/src/PJ_loxim.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(loxim, "Loximuthal") "\n\tPCyl Sph"; @@ -72,8 +73,10 @@ PJ *PROJECTION(loxim) { Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f; Q->cosphi1 = cos(Q->phi1); - if (Q->cosphi1 < EPS) - E_ERROR(-22); + if (Q->cosphi1 < EPS) { + proj_errno_set(P, PJD_ERR_LAT_LARGER_THAN_90); + return freeup_new(P); + } Q->tanphi1 = tan(M_FORTPI + 0.5 * Q->phi1); diff --git a/src/PJ_lsat.c b/src/PJ_lsat.c index a141d568..15009d61 100644 --- a/src/PJ_lsat.c +++ b/src/PJ_lsat.c @@ -1,6 +1,7 @@ /* based upon Snyder and Linck, USGS-NMD */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(lsat, "Space oblique for LANDSAT") "\n\tCyl, Sph&Ell\n\tlsat= path="; @@ -172,9 +173,15 @@ PJ *PROJECTION(lsat) { P->opaque = Q; land = pj_param(P->ctx, P->params, "ilsat").i; - if (land <= 0 || land > 5) E_ERROR(-28); + if (land <= 0 || land > 5) { + proj_errno_set(P, PJD_ERR_LSAT_NOT_IN_RANGE); + return freeup_new(P); + } path = pj_param(P->ctx, P->params, "ipath").i; - if (path <= 0 || path > (land <= 3 ? 251 : 233)) E_ERROR(-29); + if (path <= 0 || path > (land <= 3 ? 251 : 233)) { + proj_errno_set(P, PJD_ERR_PATH_NOT_IN_RANGE); + return freeup_new(P); + } if (land <= 3) { P->lam0 = DEG_TO_RAD * 128.87 - M_TWOPI / 251. * path; Q->p22 = 103.2669323; diff --git a/src/PJ_mbtfpp.c b/src/PJ_mbtfpp.c index 172be236..f8a5e807 100644 --- a/src/PJ_mbtfpp.c +++ b/src/PJ_mbtfpp.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(mbtfpp, "McBride-Thomas Flat-Polar Parabolic") "\n\tCyl., Sph."; @@ -27,19 +28,23 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ lp.phi = xy.y / FYC; if (fabs(lp.phi) >= 1.) { - if (fabs(lp.phi) > ONEEPS) - I_ERROR - else + if (fabs(lp.phi) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else { lp.phi = (lp.phi < 0.) ? -M_HALFPI : M_HALFPI; + } } else lp.phi = asin(lp.phi); lp.lam = xy.x / ( FXC * (2. * cos(C23 * (lp.phi *= 3.)) - 1.) ); if (fabs(lp.phi = sin(lp.phi) / CS) >= 1.) { - if (fabs(lp.phi) > ONEEPS) - I_ERROR - else + if (fabs(lp.phi) > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else { lp.phi = (lp.phi < 0.) ? -M_HALFPI : M_HALFPI; + } } else lp.phi = asin(lp.phi); diff --git a/src/PJ_mbtfpq.c b/src/PJ_mbtfpq.c index 7a436721..b6910a45 100644 --- a/src/PJ_mbtfpq.c +++ b/src/PJ_mbtfpq.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic") "\n\tCyl., Sph."; @@ -38,7 +39,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ lp.phi = RYC * xy.y; if (fabs(lp.phi) > 1.) { - if (fabs(lp.phi) > ONETOL) I_ERROR + if (fabs(lp.phi) > ONETOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else if (lp.phi < 0.) { t = -1.; lp.phi = -M_PI; } else { t = 1.; lp.phi = M_PI; } } else @@ -46,7 +50,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ lp.lam = RXC * xy.x / (1. + 2. * cos(lp.phi)/cos(0.5 * lp.phi)); lp.phi = RC * (t + sin(lp.phi)); if (fabs(lp.phi) > 1.) - if (fabs(lp.phi) > ONETOL) I_ERROR + if (fabs(lp.phi) > ONETOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI; else lp.phi = asin(lp.phi); diff --git a/src/PJ_merc.c b/src/PJ_merc.c index aafe7b85..04e5dc21 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; @@ -7,8 +8,10 @@ PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ XY xy = {0.0,0.0}; - if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) - F_ERROR; + if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = P->k0 * lp.lam; xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e)); return xy; @@ -17,8 +20,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; - if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) - F_ERROR; + if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; +} xy.x = P->k0 * lp.lam; xy.y = P->k0 * log(tan(M_FORTPI + .5 * lp.phi)); return xy; @@ -27,8 +32,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ LP lp = {0.0,0.0}; - if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) - I_ERROR; + if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; +} lp.lam = xy.x / P->k0; return lp; } @@ -53,7 +60,11 @@ PJ *PROJECTION(merc) { if( (is_phits = pj_param(P->ctx, P->params, "tlat_ts").i) ) { phits = fabs(pj_param(P->ctx, P->params, "rlat_ts").f); - if (phits >= M_HALFPI) E_ERROR(-24); + if (phits >= M_HALFPI) { + proj_errno_set(P, PJD_ERR_LAT_TS_LARGER_THAN_90); + freeup(P); + return 0; + } } if (P->es != 0.0) { /* ellipsoid */ diff --git a/src/PJ_misrsom.c b/src/PJ_misrsom.c index 6f1c9b29..49be1669 100644 --- a/src/PJ_misrsom.c +++ b/src/PJ_misrsom.c @@ -21,7 +21,8 @@ *****************************************************************************/ /* based upon Snyder and Linck, USGS-NMD */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(misrsom, "Space oblique for MISR") "\n\tCyl, Sph&Ell\n\tpath="; @@ -189,7 +190,10 @@ PJ *PROJECTION(misrsom) { P->opaque = Q; path = pj_param(P->ctx, P->params, "ipath").i; - if (path <= 0 || path > 233) E_ERROR(-29); + if (path <= 0 || path > 233) { + proj_errno_set(P, PJD_ERR_PATH_NOT_IN_RANGE); + return freeup_new(P); + } P->lam0 = DEG_TO_RAD * 129.3056 - M_TWOPI / 233. * path; alf = 98.30382 * DEG_TO_RAD; Q->p22 = 98.88 / 1440.0; diff --git a/src/PJ_nsper.c b/src/PJ_nsper.c index 32059bcf..589b6203 100644 --- a/src/PJ_nsper.c +++ b/src/PJ_nsper.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double height; @@ -50,7 +51,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ xy.y = sinphi; break; } - if (xy.y < Q->rp) F_ERROR; + if (xy.y < Q->rp) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.y = Q->pn1 / (Q->p - xy.y); xy.x = xy.y * cosphi * sin(lp.lam); switch (Q->mode) { @@ -95,7 +99,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ xy.y = bq * Q->cg - bm * Q->sg; } rh = hypot(xy.x, xy.y); - if ((sinz = 1. - rh * rh * Q->pfact) < 0.) I_ERROR; + if ((sinz = 1. - rh * rh * Q->pfact) < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } sinz = (Q->p - sqrt(sinz)) / (Q->pn1 / rh + rh / Q->pn1); cosz = sqrt(1. - sinz * sinz); if (fabs(rh) <= EPS10) { @@ -146,7 +153,10 @@ static void freeup (PJ *P) { static PJ *setup(PJ *P) { struct pj_opaque *Q = P->opaque; - if ((Q->height = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30); + if ((Q->height = pj_param(P->ctx, P->params, "dh").f) <= 0.) { + proj_errno_set(P, PJD_ERR_H_LESS_THAN_ZERO); + return freeup_new(P); + } if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) Q->mode = P->phi0 < 0. ? S_POLE : N_POLE; else if (fabs(P->phi0) < EPS10) diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c index d789b3f5..238dba67 100644 --- a/src/PJ_ob_tran.c +++ b/src/PJ_ob_tran.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" #include <string.h> struct pj_opaque { @@ -110,11 +111,20 @@ PJ *PROJECTION(ob_tran) { P->opaque = Q; /* get name of projection to be translated */ - if (!(name = pj_param(P->ctx, P->params, "so_proj").s)) E_ERROR(-26); + if (!(name = pj_param(P->ctx, P->params, "so_proj").s)) { + proj_errno_set(P, PJD_ERR_NO_ROTATION_PROJ); + return freeup_new(P); + } /* avoid endless recursion */ - if( strcmp(name, "ob_tran") == 0 ) E_ERROR(-37); + if( strcmp(name, "ob_tran") == 0 ) { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ; - if (!s || !(Q->link = (*pj_list[i].proj)(0))) E_ERROR(-37); + if (!s || !(Q->link = (*pj_list[i].proj)(0))) { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } /* copy existing header into new */ P->es = 0.; /* force to spherical */ Q->link->params = P->params; @@ -149,8 +159,10 @@ PJ *PROJECTION(ob_tran) { fabs(fabs(phic) - HALFPI) <= TOL || fabs(fabs(alpha) - HALFPI) <= TOL) */ - if (fabs(fabs(phic) - M_HALFPI) <= TOL) - E_ERROR(-32); + if (fabs(fabs(phic) - M_HALFPI) <= TOL) { + proj_errno_set(P, PJD_ERR_LAT_0_OR_ALPHA_EQ_90); + return freeup_new(P); + } Q->lamp = lamc + aatan2(-cos(alpha), -sin(alpha) * sin(phic)); phip = aasin(P->ctx,cos(phic) * sin(alpha)); } else if (pj_param(P->ctx, P->params, "to_lat_p").i) { /* specified new pole */ @@ -163,10 +175,11 @@ PJ *PROJECTION(ob_tran) { phi1 = pj_param(P->ctx, P->params, "ro_lat_1").f; lam2 = pj_param(P->ctx, P->params, "ro_lon_2").f; phi2 = pj_param(P->ctx, P->params, "ro_lat_2").f; - if (fabs(phi1 - phi2) <= TOL || - (con = fabs(phi1)) <= TOL || - fabs(con - M_HALFPI) <= TOL || - fabs(fabs(phi2) - M_HALFPI) <= TOL) E_ERROR(-33); + if (fabs(phi1 - phi2) <= TOL || (con = fabs(phi1)) <= TOL || + fabs(con - M_HALFPI) <= TOL || fabs(fabs(phi2) - M_HALFPI) <= TOL) { + proj_errno_set(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90); + return freeup_new(P); + } Q->lamp = atan2(cos(phi1) * sin(phi2) * cos(lam1) - sin(phi1) * cos(phi2) * cos(lam2), sin(phi1) * cos(phi2) * sin(lam2) - diff --git a/src/PJ_oea.c b/src/PJ_oea.c index 11efb194..2ba7917f 100644 --- a/src/PJ_oea.c +++ b/src/PJ_oea.c @@ -1,53 +1,54 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(oea, "Oblated Equal Area") "\n\tMisc Sph\n\tn= m= theta="; struct pj_opaque { - double theta; - double m, n; - double two_r_m, two_r_n, rm, rn, hm, hn; - double cp0, sp0; + double theta; + double m, n; + double two_r_m, two_r_n, rm, rn, hm, hn; + double cp0, sp0; }; static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double Az, M, N, cp, sp, cl, shz; - - cp = cos(lp.phi); - sp = sin(lp.phi); - cl = cos(lp.lam); - Az = aatan2(cp * sin(lp.lam), Q->cp0 * sp - Q->sp0 * cp * cl) + Q->theta; - shz = sin(0.5 * aacos(P->ctx, Q->sp0 * sp + Q->cp0 * cp * cl)); - M = aasin(P->ctx, shz * sin(Az)); - N = aasin(P->ctx, shz * cos(Az) * cos(M) / cos(M * Q->two_r_m)); - xy.y = Q->n * sin(N * Q->two_r_n); - xy.x = Q->m * sin(M * Q->two_r_m) * cos(N) / cos(N * Q->two_r_n); - - return xy; + double Az, M, N, cp, sp, cl, shz; + + cp = cos(lp.phi); + sp = sin(lp.phi); + cl = cos(lp.lam); + Az = aatan2(cp * sin(lp.lam), Q->cp0 * sp - Q->sp0 * cp * cl) + Q->theta; + shz = sin(0.5 * aacos(P->ctx, Q->sp0 * sp + Q->cp0 * cp * cl)); + M = aasin(P->ctx, shz * sin(Az)); + N = aasin(P->ctx, shz * cos(Az) * cos(M) / cos(M * Q->two_r_m)); + xy.y = Q->n * sin(N * Q->two_r_n); + xy.x = Q->m * sin(M * Q->two_r_m) * cos(N) / cos(N * Q->two_r_n); + + return xy; } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double N, M, xp, yp, z, Az, cz, sz, cAz; - - N = Q->hn * aasin(P->ctx,xy.y * Q->rn); - M = Q->hm * aasin(P->ctx,xy.x * Q->rm * cos(N * Q->two_r_n) / cos(N)); - xp = 2. * sin(M); - yp = 2. * sin(N) * cos(M * Q->two_r_m) / cos(M); - cAz = cos(Az = aatan2(xp, yp) - Q->theta); - z = 2. * aasin(P->ctx, 0.5 * hypot(xp, yp)); - sz = sin(z); - cz = cos(z); - lp.phi = aasin(P->ctx, Q->sp0 * cz + Q->cp0 * sz * cAz); - lp.lam = aatan2(sz * sin(Az), - Q->cp0 * cz - Q->sp0 * sz * cAz); - - return lp; + double N, M, xp, yp, z, Az, cz, sz, cAz; + + N = Q->hn * aasin(P->ctx,xy.y * Q->rn); + M = Q->hm * aasin(P->ctx,xy.x * Q->rm * cos(N * Q->two_r_n) / cos(N)); + xp = 2. * sin(M); + yp = 2. * sin(N) * cos(M * Q->two_r_m) / cos(M); + cAz = cos(Az = aatan2(xp, yp) - Q->theta); + z = 2. * aasin(P->ctx, 0.5 * hypot(xp, yp)); + sz = sin(z); + cz = cos(z); + lp.phi = aasin(P->ctx, Q->sp0 * cz + Q->cp0 * sz * cAz); + lp.lam = aatan2(sz * sin(Az), + Q->cp0 * cz - Q->sp0 * sz * cAz); + + return lp; } @@ -73,25 +74,24 @@ PJ *PROJECTION(oea) { return freeup_new (P); P->opaque = Q; - if (((Q->n = pj_param(P->ctx, P->params, "dn").f) <= 0.) || - ((Q->m = pj_param(P->ctx, P->params, "dm").f) <= 0.)) - { - E_ERROR(-39) - } - else { - Q->theta = pj_param(P->ctx, P->params, "rtheta").f; - Q->sp0 = sin(P->phi0); - Q->cp0 = cos(P->phi0); - Q->rn = 1./ Q->n; - Q->rm = 1./ Q->m; - Q->two_r_n = 2. * Q->rn; - Q->two_r_m = 2. * Q->rm; - Q->hm = 0.5 * Q->m; - Q->hn = 0.5 * Q->n; - P->fwd = s_forward; - P->inv = s_inverse; - P->es = 0.; - } + if (((Q->n = pj_param(P->ctx, P->params, "dn").f) <= 0.) || + ((Q->m = pj_param(P->ctx, P->params, "dm").f) <= 0.)) { + proj_errno_set(P, PJD_ERR_INVALID_M_OR_N); + return freeup_new(P); + } else { + Q->theta = pj_param(P->ctx, P->params, "rtheta").f; + Q->sp0 = sin(P->phi0); + Q->cp0 = cos(P->phi0); + Q->rn = 1./ Q->n; + Q->rm = 1./ Q->m; + Q->two_r_n = 2. * Q->rn; + Q->two_r_m = 2. * Q->rm; + Q->hm = 0.5 * Q->m; + Q->hn = 0.5 * Q->n; + P->fwd = s_forward; + P->inv = s_inverse; + P->es = 0.; + } return P; } diff --git a/src/PJ_omerc.c b/src/PJ_omerc.c index 8874784d..0d86e460 100644 --- a/src/PJ_omerc.c +++ b/src/PJ_omerc.c @@ -22,7 +22,8 @@ ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(omerc, "Oblique Mercator") "\n\tCyl, Sph&Ell no_rot\n\t" @@ -50,8 +51,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ T = .5 * (W + temp); V = sin(Q->B * lp.lam); U = (S * Q->singam - V * Q->cosgam) / T; - if (fabs(fabs(U) - 1.0) < EPS) - F_ERROR; + if (fabs(fabs(U) - 1.0) < EPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } v = 0.5 * Q->ArB * log((1. - U)/(1. + U)); temp = cos(Q->B * lp.lam); if(fabs(temp) < TOL) { @@ -97,8 +100,10 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ lp.phi = Up < 0. ? -M_HALFPI : M_HALFPI; } else { lp.phi = Q->E / sqrt((1. + Up) / (1. - Up)); - if ((lp.phi = pj_phi2(P->ctx, pow(lp.phi, 1. / Q->B), P->e)) == HUGE_VAL) - I_ERROR; + if ((lp.phi = pj_phi2(P->ctx, pow(lp.phi, 1. / Q->B), P->e)) == HUGE_VAL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } lp.lam = - Q->rB * atan2((Sp * Q->cosgam - Vp * Q->singam), cos(Q->BrA * u)); } @@ -159,7 +164,10 @@ PJ *PROJECTION(omerc) { (con = fabs(phi1)) <= TOL || fabs(con - M_HALFPI) <= TOL || fabs(fabs(P->phi0) - M_HALFPI) <= TOL || - fabs(fabs(phi2) - M_HALFPI) <= TOL) E_ERROR(-33); + fabs(fabs(phi2) - M_HALFPI) <= TOL) { + proj_errno_set(P, PJD_ERR_LAT_0_OR_ALPHA_EQ_90); + return freeup_new(P); + } } com = sqrt(P->one_es); if (fabs(P->phi0) > EPS) { diff --git a/src/PJ_ortho.c b/src/PJ_ortho.c index d0897a72..bb150f05 100644 --- a/src/PJ_ortho.c +++ b/src/PJ_ortho.c @@ -1,48 +1,57 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph."; struct pj_opaque { - double sinph0; - double cosph0; - int mode; + double sinph0; + double cosph0; + int mode; }; #define EPS10 1.e-10 -#define N_POLE 0 +#define N_POLE 0 #define S_POLE 1 -#define EQUIT 2 -#define OBLIQ 3 +#define EQUIT 2 +#define OBLIQ 3 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double coslam, cosphi, sinphi; - - cosphi = cos(lp.phi); - coslam = cos(lp.lam); - switch (Q->mode) { - case EQUIT: - if (cosphi * coslam < - EPS10) F_ERROR; - xy.y = sin(lp.phi); - break; - case OBLIQ: - if (Q->sinph0 * (sinphi = sin(lp.phi)) + - Q->cosph0 * cosphi * coslam < - EPS10) F_ERROR; - xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam; - break; - case N_POLE: - coslam = - coslam; + double coslam, cosphi, sinphi; + + cosphi = cos(lp.phi); + coslam = cos(lp.lam); + switch (Q->mode) { + case EQUIT: + if (cosphi * coslam < - EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + xy.y = sin(lp.phi); + break; + case OBLIQ: + if (Q->sinph0 * (sinphi = sin(lp.phi)) + Q->cosph0 * cosphi * coslam < - EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam; + break; + case N_POLE: + coslam = - coslam; /*-fallthrough*/ - case S_POLE: - if (fabs(lp.phi - P->phi0) - EPS10 > M_HALFPI) F_ERROR; - xy.y = cosphi * coslam; - break; - } - xy.x = cosphi * sin(lp.lam); - return xy; + case S_POLE: + if (fabs(lp.phi - P->phi0) - EPS10 > M_HALFPI) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + xy.y = cosphi * coslam; + break; + } + xy.x = cosphi * sin(lp.lam); + return xy; } @@ -52,7 +61,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ double rh, cosc, sinc; if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) { - if ((sinc - 1.) > EPS10) I_ERROR; + if ((sinc - 1.) > EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } sinc = 1.; } cosc = sqrt(1. - sinc * sinc); /* in this range OK */ @@ -114,17 +126,17 @@ PJ *PROJECTION(ortho) { return freeup_new (P); P->opaque = Q; - if (fabs(fabs(P->phi0) - M_HALFPI) <= EPS10) - Q->mode = P->phi0 < 0. ? S_POLE : N_POLE; - else if (fabs(P->phi0) > EPS10) { - Q->mode = OBLIQ; - Q->sinph0 = sin(P->phi0); - Q->cosph0 = cos(P->phi0); - } else - Q->mode = EQUIT; - P->inv = s_inverse; - P->fwd = s_forward; - P->es = 0.; + if (fabs(fabs(P->phi0) - M_HALFPI) <= EPS10) + Q->mode = P->phi0 < 0. ? S_POLE : N_POLE; + else if (fabs(P->phi0) > EPS10) { + Q->mode = OBLIQ; + Q->sinph0 = sin(P->phi0); + Q->cosph0 = cos(P->phi0); + } else + Q->mode = EQUIT; + P->inv = s_inverse; + P->fwd = s_forward; + P->es = 0.; return P; } diff --git a/src/PJ_poly.c b/src/PJ_poly.c index 8fb671af..f5af36ac 100644 --- a/src/PJ_poly.c +++ b/src/PJ_poly.c @@ -1,17 +1,18 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(poly, "Polyconic (American)") - "\n\tConic, Sph&Ell"; + "\n\tConic, Sph&Ell"; struct pj_opaque { - double ml0; \ - double *en; + double ml0; \ + double *en; }; -#define TOL 1e-10 -#define CONV 1e-10 -#define N_ITER 10 +#define TOL 1e-10 +#define CONV 1e-10 +#define N_ITER 10 #define I_ITER 20 #define ITOL 1.e-12 @@ -19,37 +20,37 @@ struct pj_opaque { static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double ms, sp, cp; + double ms, sp, cp; - if (fabs(lp.phi) <= TOL) { + if (fabs(lp.phi) <= TOL) { xy.x = lp.lam; xy.y = -Q->ml0; } else { - sp = sin(lp.phi); - ms = fabs(cp = cos(lp.phi)) > TOL ? pj_msfn(sp, cp, P->es) / sp : 0.; - xy.x = ms * sin(lp.lam *= sp); - xy.y = (pj_mlfn(lp.phi, sp, cp, Q->en) - Q->ml0) + ms * (1. - cos(lp.lam)); - } + sp = sin(lp.phi); + ms = fabs(cp = cos(lp.phi)) > TOL ? pj_msfn(sp, cp, P->es) / sp : 0.; + xy.x = ms * sin(lp.lam *= sp); + xy.y = (pj_mlfn(lp.phi, sp, cp, Q->en) - Q->ml0) + ms * (1. - cos(lp.lam)); + } - return xy; + 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 cot, E; + double cot, E; - if (fabs(lp.phi) <= TOL) { + if (fabs(lp.phi) <= TOL) { xy.x = lp.lam; xy.y = Q->ml0; } else { - cot = 1. / tan(lp.phi); - xy.x = sin(E = lp.lam * sin(lp.phi)) * cot; - xy.y = lp.phi - P->phi0 + cot * (1. - cos(E)); - } + cot = 1. / tan(lp.phi); + xy.x = sin(E = lp.lam * sin(lp.phi)) * cot; + xy.y = lp.phi - P->phi0 + cot * (1. - cos(E)); + } - return xy; + return xy; } @@ -57,64 +58,71 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - xy.y += Q->ml0; - if (fabs(xy.y) <= TOL) { + xy.y += Q->ml0; + if (fabs(xy.y) <= TOL) { lp.lam = xy.x; lp.phi = 0.; } else { - double r, c, sp, cp, s2ph, ml, mlb, mlp, dPhi; - int i; - - r = xy.y * xy.y + xy.x * xy.x; - for (lp.phi = xy.y, i = I_ITER; i ; --i) { - sp = sin(lp.phi); - s2ph = sp * ( cp = cos(lp.phi)); - if (fabs(cp) < ITOL) - I_ERROR; - c = sp * (mlp = sqrt(1. - P->es * sp * sp)) / cp; - ml = pj_mlfn(lp.phi, sp, cp, Q->en); - mlb = ml * ml + r; - mlp = P->one_es / (mlp * mlp * mlp); - lp.phi += ( dPhi = - ( ml + ml + c * mlb - 2. * xy.y * (c * ml + 1.) ) / ( - P->es * s2ph * (mlb - 2. * xy.y * ml) / c + - 2.* (xy.y - ml) * (c * mlp - 1. / s2ph) - mlp - mlp )); - if (fabs(dPhi) <= ITOL) - break; - } - if (!i) - I_ERROR; - c = sin(lp.phi); - lp.lam = asin(xy.x * tan(lp.phi) * sqrt(1. - P->es * c * c)) / sin(lp.phi); - } - - return lp; + double r, c, sp, cp, s2ph, ml, mlb, mlp, dPhi; + int i; + + r = xy.y * xy.y + xy.x * xy.x; + for (lp.phi = xy.y, i = I_ITER; i ; --i) { + sp = sin(lp.phi); + s2ph = sp * ( cp = cos(lp.phi)); + if (fabs(cp) < ITOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + c = sp * (mlp = sqrt(1. - P->es * sp * sp)) / cp; + ml = pj_mlfn(lp.phi, sp, cp, Q->en); + mlb = ml * ml + r; + mlp = P->one_es / (mlp * mlp * mlp); + lp.phi += ( dPhi = + ( ml + ml + c * mlb - 2. * xy.y * (c * ml + 1.) ) / ( + P->es * s2ph * (mlb - 2. * xy.y * ml) / c + + 2.* (xy.y - ml) * (c * mlp - 1. / s2ph) - mlp - mlp )); + if (fabs(dPhi) <= ITOL) + break; + } + if (!i) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + c = sin(lp.phi); + lp.lam = asin(xy.x * tan(lp.phi) * sqrt(1. - P->es * c * c)) / sin(lp.phi); + } + + return lp; } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; - double B, dphi, tp; - int i; + double B, dphi, tp; + int i; - if (fabs(xy.y = P->phi0 + xy.y) <= TOL) { + if (fabs(xy.y = P->phi0 + xy.y) <= TOL) { lp.lam = xy.x; lp.phi = 0.; } else { - lp.phi = xy.y; - B = xy.x * xy.x + xy.y * xy.y; - i = N_ITER; - do { - tp = tan(lp.phi); - lp.phi -= (dphi = (xy.y * (lp.phi * tp + 1.) - lp.phi - - .5 * ( lp.phi * lp.phi + B) * tp) / - ((lp.phi - xy.y) / tp - 1.)); - } while (fabs(dphi) > CONV && --i); - if (! i) I_ERROR; - lp.lam = asin(xy.x * tan(lp.phi)) / sin(lp.phi); - } - - return lp; + lp.phi = xy.y; + B = xy.x * xy.x + xy.y * xy.y; + i = N_ITER; + do { + tp = tan(lp.phi); + lp.phi -= (dphi = (xy.y * (lp.phi * tp + 1.) - lp.phi - + .5 * ( lp.phi * lp.phi + B) * tp) / + ((lp.phi - xy.y) / tp - 1.)); + } while (fabs(dphi) > CONV && --i); + if (! i) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + lp.lam = asin(xy.x * tan(lp.phi)) / sin(lp.phi); + } + + return lp; } @@ -142,16 +150,16 @@ PJ *PROJECTION(poly) { return freeup_new (P); P->opaque = Q; - if (P->es != 0.0) { - if (!(Q->en = pj_enfn(P->es))) E_ERROR_0; - Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en); - P->inv = e_inverse; - P->fwd = e_forward; - } else { - Q->ml0 = -P->phi0; - P->inv = s_inverse; - P->fwd = s_forward; - } + if (P->es != 0.0) { + if (!(Q->en = pj_enfn(P->es))) return freeup_new(P); + Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en); + P->inv = e_inverse; + P->fwd = e_forward; + } else { + Q->ml0 = -P->phi0; + P->inv = s_inverse; + P->fwd = s_forward; + } return P; } diff --git a/src/PJ_robin.c b/src/PJ_robin.c index 4eab552a..60a52324 100644 --- a/src/PJ_robin.c +++ b/src/PJ_robin.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(robin, "Robinson") "\n\tPCyl., Sph."; @@ -80,8 +81,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ (void) P; i = (int)floor((dphi = fabs(lp.phi)) * C1); - if( i < 0 ) - F_ERROR; + if( i < 0 ){ + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } if (i >= NODES) i = NODES - 1; dphi = RAD_TO_DEG * (dphi - RC1 * i); xy.x = V(X[i], dphi) * FXC * lp.lam; @@ -102,7 +105,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ lp.lam = xy.x / FXC; lp.phi = fabs(xy.y / FYC); if (lp.phi >= 1.) { /* simple pathologic cases */ - if (lp.phi > ONEEPS) I_ERROR + if (lp.phi > ONEEPS) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } else { lp.phi = xy.y < 0. ? -M_HALFPI : M_HALFPI; lp.lam /= X[NODES].c0; @@ -110,8 +116,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ } else { /* general problem */ /* in Y space, reduce to table interval */ i = (int)floor(lp.phi * NODES); - if( i < 0 || i >= NODES ) - I_ERROR; + if( i < 0 || i >= NODES ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } for (;;) { if (Y[i].c0 > lp.phi) --i; else if (Y[i+1].c0 <= lp.phi) ++i; diff --git a/src/PJ_sch.c b/src/PJ_sch.c index 98c0d2a1..6c97a4f3 100644 --- a/src/PJ_sch.c +++ b/src/PJ_sch.c @@ -33,7 +33,8 @@ ****************************************************************************/ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" #include "geocent.h" struct pj_opaque { @@ -61,9 +62,10 @@ static LPZ inverse3d(XYZ xyz, PJ *P) { pxyz[1] = xyz.x * P->a / Q->rcurv; pxyz[2] = xyz.z; - if( pj_Convert_Geodetic_To_Geocentric( &(Q->sph), pxyz[0], pxyz[1], pxyz[2], - temp, temp+1, temp+2) != 0) - I3_ERROR; + if( pj_Convert_Geodetic_To_Geocentric( &(Q->sph), pxyz[0], pxyz[1], pxyz[2], temp, temp+1, temp+2) != 0) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lpz; + } /* Apply rotation */ pxyz[0] = Q->transMat[0] * temp[0] + Q->transMat[1] * temp[1] + Q->transMat[2] * temp[2]; @@ -100,9 +102,10 @@ static XYZ forward3d(LPZ lpz, PJ *P) { /* Convert lat lon to geocentric coordinates */ - if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), lpz.phi, lpz.lam, lpz.z, - temp, temp+1, temp+2 ) != 0 ) - F3_ERROR; + if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), lpz.phi, lpz.lam, lpz.z, temp, temp+1, temp+2 ) != 0 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xyz; + } /* Adjust for offset */ @@ -162,8 +165,10 @@ static PJ *setup(PJ *P) { /* general initialization */ temp = P->a * sqrt(1.0 - P->es); /* Setup original geocentric system */ - if ( pj_Set_Geocentric_Parameters(&(Q->elp_0), P->a, temp) != 0) - E_ERROR(-37); + if ( pj_Set_Geocentric_Parameters(&(Q->elp_0), P->a, temp) != 0) { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } clt = cos(Q->plat); slt = sin(Q->plat); @@ -187,9 +192,10 @@ static PJ *setup(PJ *P) { /* general initialization */ #endif /* Set up local sphere at the given peg point */ - if ( pj_Set_Geocentric_Parameters(&(Q->sph), Q->rcurv, Q->rcurv) != 0) - E_ERROR(-37); - + if ( pj_Set_Geocentric_Parameters(&(Q->sph), Q->rcurv, Q->rcurv) != 0) { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } /* Set up the transformation matrices */ Q->transMat[0] = clt * clo; Q->transMat[1] = -shdg*slo - slt*clo * chdg; @@ -205,7 +211,8 @@ static PJ *setup(PJ *P) { /* general initialization */ if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), Q->plat, Q->plon, Q->h0, pxyz, pxyz+1, pxyz+2 ) != 0 ) { - E_ERROR(-14) + proj_errno_set(P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT); + return freeup_new(P); } @@ -234,20 +241,26 @@ PJ *PROJECTION(sch) { /* Check if peg latitude was defined */ if (pj_param(P->ctx, P->params, "tplat_0").i) Q->plat = pj_param(P->ctx, P->params, "rplat_0").f; - else - E_ERROR(-37); + else { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } /* Check if peg longitude was defined */ if (pj_param(P->ctx, P->params, "tplon_0").i) Q->plon = pj_param(P->ctx, P->params, "rplon_0").f; - else - E_ERROR(-37); + else { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } /* Check if peg latitude is defined */ if (pj_param(P->ctx, P->params, "tphdg_0").i) Q->phdg = pj_param(P->ctx, P->params, "rphdg_0").f; - else - E_ERROR(-37); + else { + proj_errno_set(P, PJD_ERR_FAILED_TO_FIND_PROJ); + return freeup_new(P); + } /* Check if average height was defined - If so read it in */ diff --git a/src/PJ_sconics.c b/src/PJ_sconics.c index df7fb661..8759c66e 100644 --- a/src/PJ_sconics.c +++ b/src/PJ_sconics.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { @@ -46,7 +47,7 @@ static int phi12(PJ *P, double *del) { p2 = pj_param(P->ctx, P->params, "rlat_2").f; *del = 0.5 * (p2 - p1); P->opaque->sig = 0.5 * (p2 + p1); - err = (fabs(*del) < EPS || fabs(P->opaque->sig) < EPS) ? -42 : 0; + err = (fabs(*del) < EPS || fabs(P->opaque->sig) < EPS) ? PJD_ERR_ABS_LAT1_EQ_ABS_LAT2 : 0; *del = *del; } return err; @@ -129,8 +130,10 @@ static PJ *setup(PJ *P, int type) { Q->type = type; i = phi12 (P, &del); - if(i) - E_ERROR(i); + if(i) { + proj_errno_set(P, i); + return freeup_new(P); + } switch (Q->type) { case TISSOT: @@ -169,8 +172,10 @@ static PJ *setup(PJ *P, int type) { Q->n = sin (Q->sig); Q->c2 = cos (del); Q->c1 = 1./tan (Q->sig); - if (fabs (del = P->phi0 - Q->sig) - EPS10 >= M_HALFPI) - E_ERROR(-43); + if (fabs (del = P->phi0 - Q->sig) - EPS10 >= M_HALFPI) { + proj_errno_set(P, PJD_ERR_LAT_0_HALF_PI_FROM_MEAN); + return freeup_new(P); + } Q->rho_0 = Q->c2 * (Q->c1 - tan (del)); break; diff --git a/src/PJ_somerc.c b/src/PJ_somerc.c index f7438a92..751d0c77 100644 --- a/src/PJ_somerc.c +++ b/src/PJ_somerc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(somerc, "Swiss. Obl. Mercator") "\n\tCyl, Ell\n\tFor CH1903"; @@ -54,8 +55,10 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ if (i) { lp.phi = phip; lp.lam = lamp / Q->c; - } else - I_ERROR + } else { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } return (lp); } diff --git a/src/PJ_stere.c b/src/PJ_stere.c index f338d16d..befca308 100644 --- a/src/PJ_stere.c +++ b/src/PJ_stere.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts="; PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth"; @@ -96,7 +97,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ case OBLIQ: xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam; oblcon: - if (xy.y <= EPS10) F_ERROR; + if (xy.y <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = (xy.y = Q->akm1 / xy.y) * cosphi * sinlam; xy.y *= (Q->mode == EQUIT) ? sinphi : cosph0 * sinphi - sinph0 * cosphi * coslam; @@ -106,7 +110,10 @@ oblcon: lp.phi = - lp.phi; /*-fallthrough*/ case S_POLE: - if (fabs (lp.phi - M_HALFPI) < TOL) F_ERROR; + if (fabs (lp.phi - M_HALFPI) < TOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = sinlam * ( xy.y = Q->akm1 * tan (M_FORTPI + .5 * lp.phi) ); xy.y *= coslam; break; @@ -159,7 +166,9 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ return lp; } } - I_ERROR; + + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; } @@ -298,14 +307,17 @@ PJ *PROJECTION(ups) { return freeup_new (P); P->opaque = Q; - /* International Ellipsoid */ - P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI; - if (P->es == 0.0) E_ERROR(-34); - P->k0 = .994; - P->x0 = 2000000.; - P->y0 = 2000000.; - Q->phits = M_HALFPI; - P->lam0 = 0.; + /* International Ellipsoid */ + P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI; + if (P->es == 0.0) { + proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + return freeup_new(P); + } + P->k0 = .994; + P->x0 = 2000000.; + P->y0 = 2000000.; + Q->phits = M_HALFPI; + P->lam0 = 0.; return setup(P); } diff --git a/src/PJ_tcc.c b/src/PJ_tcc.c index 316d5e50..80bb6bc8 100644 --- a/src/PJ_tcc.c +++ b/src/PJ_tcc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(tcc, "Transverse Central Cylindrical") "\n\tCyl, Sph, no inv."; @@ -11,7 +12,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ double b, bt; b = cos (lp.phi) * sin (lp.lam); - if ((bt = 1. - b * b) < EPS10) F_ERROR; + if ((bt = 1. - b * b) < EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = b / sqrt(bt); xy.y = atan2 (tan (lp.phi) , cos (lp.lam)); return xy; diff --git a/src/PJ_tmerc.c b/src/PJ_tmerc.c index 5cbb2b5d..04afc64c 100644 --- a/src/PJ_tmerc.c +++ b/src/PJ_tmerc.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell"; @@ -83,16 +84,20 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ cosphi = cos(lp.phi); b = cosphi * sin (lp.lam); - if (fabs (fabs (b) - 1.) <= EPS10) - F_ERROR; + if (fabs (fabs (b) - 1.) <= EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } xy.x = P->opaque->ml0 * log ((1. + b) / (1. - b)); xy.y = cosphi * cos (lp.lam) / sqrt (1. - b * b); b = fabs ( xy.y ); if (b >= 1.) { - if ((b - 1.) > EPS10) - F_ERROR + if ((b - 1.) > EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } else xy.y = 0.; } else xy.y = acos (xy.y); @@ -174,7 +179,7 @@ static PJ *setup(PJ *P) { /* general initialization */ struct pj_opaque *Q = P->opaque; if (P->es != 0.0) { if (!(Q->en = pj_enfn(P->es))) - E_ERROR_0; + return freeup_new(P); Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en); Q->esp = P->es / (1. - P->es); P->inv = e_inverse; diff --git a/src/PJ_tpeqd.c b/src/PJ_tpeqd.c index 1b289a0b..90dd7568 100644 --- a/src/PJ_tpeqd.c +++ b/src/PJ_tpeqd.c @@ -1,56 +1,57 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(tpeqd, "Two Point Equidistant") - "\n\tMisc Sph\n\tlat_1= lon_1= lat_2= lon_2="; + "\n\tMisc Sph\n\tlat_1= lon_1= lat_2= lon_2="; struct pj_opaque { - double cp1, sp1, cp2, sp2, ccs, cs, sc, r2z0, z02, dlam2; \ - double hz0, thz0, rhshz0, ca, sa, lp, lamc; + double cp1, sp1, cp2, sp2, ccs, cs, sc, r2z0, z02, dlam2; \ + double hz0, thz0, rhshz0, ca, sa, lp, lamc; }; static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0, 0.0}; struct pj_opaque *Q = P->opaque; - double t, z1, z2, dl1, dl2, sp, cp; - - sp = sin(lp.phi); - cp = cos(lp.phi); - z1 = aacos(P->ctx, Q->sp1 * sp + Q->cp1 * cp * cos (dl1 = lp.lam + Q->dlam2)); - z2 = aacos(P->ctx, Q->sp2 * sp + Q->cp2 * cp * cos (dl2 = lp.lam - Q->dlam2)); - z1 *= z1; - z2 *= z2; - - xy.x = Q->r2z0 * (t = z1 - z2); - t = Q->z02 - t; - xy.y = Q->r2z0 * asqrt (4. * Q->z02 * z2 - t * t); - if ((Q->ccs * sp - cp * (Q->cs * sin(dl1) - Q->sc * sin(dl2))) < 0.) - xy.y = -xy.y; - return xy; + double t, z1, z2, dl1, dl2, sp, cp; + + sp = sin(lp.phi); + cp = cos(lp.phi); + z1 = aacos(P->ctx, Q->sp1 * sp + Q->cp1 * cp * cos (dl1 = lp.lam + Q->dlam2)); + z2 = aacos(P->ctx, Q->sp2 * sp + Q->cp2 * cp * cos (dl2 = lp.lam - Q->dlam2)); + z1 *= z1; + z2 *= z2; + + xy.x = Q->r2z0 * (t = z1 - z2); + t = Q->z02 - t; + xy.y = Q->r2z0 * asqrt (4. * Q->z02 * z2 - t * t); + if ((Q->ccs * sp - cp * (Q->cs * sin(dl1) - Q->sc * sin(dl2))) < 0.) + xy.y = -xy.y; + return xy; } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double cz1, cz2, s, d, cp, sp; - - cz1 = cos (hypot(xy.y, xy.x + Q->hz0)); - cz2 = cos (hypot(xy.y, xy.x - Q->hz0)); - s = cz1 + cz2; - d = cz1 - cz2; - lp.lam = - atan2(d, (s * Q->thz0)); - lp.phi = aacos(P->ctx, hypot (Q->thz0 * s, d) * Q->rhshz0); - if ( xy.y < 0. ) - lp.phi = - lp.phi; - /* lam--phi now in system relative to P1--P2 base equator */ - sp = sin (lp.phi); - cp = cos (lp.phi); - lp.phi = aasin (P->ctx, Q->sa * sp + Q->ca * cp * (s = cos(lp.lam -= Q->lp))); - lp.lam = atan2 (cp * sin(lp.lam), Q->sa * cp * s - Q->ca * sp) + Q->lamc; - return lp; + double cz1, cz2, s, d, cp, sp; + + cz1 = cos (hypot(xy.y, xy.x + Q->hz0)); + cz2 = cos (hypot(xy.y, xy.x - Q->hz0)); + s = cz1 + cz2; + d = cz1 - cz2; + lp.lam = - atan2(d, (s * Q->thz0)); + lp.phi = aacos(P->ctx, hypot (Q->thz0 * s, d) * Q->rhshz0); + if ( xy.y < 0. ) + lp.phi = - lp.phi; + /* lam--phi now in system relative to P1--P2 base equator */ + sp = sin (lp.phi); + cp = cos (lp.phi); + lp.phi = aasin (P->ctx, Q->sa * sp + Q->ca * cp * (s = cos(lp.lam -= Q->lp))); + lp.lam = atan2 (cp * sin(lp.lam), Q->sa * cp * s - Q->ca * sp) + Q->lamc; + return lp; } @@ -70,48 +71,50 @@ static void freeup (PJ *P) { PJ *PROJECTION(tpeqd) { - double lam_1, lam_2, phi_1, phi_2, A12, pp; + double lam_1, lam_2, phi_1, phi_2, A12, pp; struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque)); if (0==Q) return freeup_new (P); P->opaque = Q; - /* get control point locations */ - phi_1 = pj_param(P->ctx, P->params, "rlat_1").f; - lam_1 = pj_param(P->ctx, P->params, "rlon_1").f; - phi_2 = pj_param(P->ctx, P->params, "rlat_2").f; - lam_2 = pj_param(P->ctx, P->params, "rlon_2").f; - - if (phi_1 == phi_2 && lam_1 == lam_2) - E_ERROR(-25); - P->lam0 = adjlon (0.5 * (lam_1 + lam_2)); - Q->dlam2 = adjlon (lam_2 - lam_1); - - Q->cp1 = cos (phi_1); - Q->cp2 = cos (phi_2); - Q->sp1 = sin (phi_1); - Q->sp2 = sin (phi_2); - Q->cs = Q->cp1 * Q->sp2; - Q->sc = Q->sp1 * Q->cp2; - Q->ccs = Q->cp1 * Q->cp2 * sin(Q->dlam2); - Q->z02 = aacos(P->ctx, Q->sp1 * Q->sp2 + Q->cp1 * Q->cp2 * cos (Q->dlam2)); - Q->hz0 = .5 * Q->z02; - A12 = atan2(Q->cp2 * sin (Q->dlam2), - Q->cp1 * Q->sp2 - Q->sp1 * Q->cp2 * cos (Q->dlam2)); - Q->ca = cos(pp = aasin(P->ctx, Q->cp1 * sin(A12))); - Q->sa = sin(pp); - Q->lp = adjlon ( atan2 (Q->cp1 * cos(A12), Q->sp1) - Q->hz0); - Q->dlam2 *= .5; - Q->lamc = M_HALFPI - atan2(sin(A12) * Q->sp1, cos(A12)) - Q->dlam2; - Q->thz0 = tan (Q->hz0); - Q->rhshz0 = .5 / sin (Q->hz0); - Q->r2z0 = 0.5 / Q->z02; - Q->z02 *= Q->z02; + /* get control point locations */ + phi_1 = pj_param(P->ctx, P->params, "rlat_1").f; + lam_1 = pj_param(P->ctx, P->params, "rlon_1").f; + phi_2 = pj_param(P->ctx, P->params, "rlat_2").f; + lam_2 = pj_param(P->ctx, P->params, "rlon_2").f; + + if (phi_1 == phi_2 && lam_1 == lam_2) { + proj_errno_set(P, PJD_ERR_CONTROL_POINT_NO_DIST); + return freeup_new(P); + } + P->lam0 = adjlon (0.5 * (lam_1 + lam_2)); + Q->dlam2 = adjlon (lam_2 - lam_1); + + Q->cp1 = cos (phi_1); + Q->cp2 = cos (phi_2); + Q->sp1 = sin (phi_1); + Q->sp2 = sin (phi_2); + Q->cs = Q->cp1 * Q->sp2; + Q->sc = Q->sp1 * Q->cp2; + Q->ccs = Q->cp1 * Q->cp2 * sin(Q->dlam2); + Q->z02 = aacos(P->ctx, Q->sp1 * Q->sp2 + Q->cp1 * Q->cp2 * cos (Q->dlam2)); + Q->hz0 = .5 * Q->z02; + A12 = atan2(Q->cp2 * sin (Q->dlam2), + Q->cp1 * Q->sp2 - Q->sp1 * Q->cp2 * cos (Q->dlam2)); + Q->ca = cos(pp = aasin(P->ctx, Q->cp1 * sin(A12))); + Q->sa = sin(pp); + Q->lp = adjlon ( atan2 (Q->cp1 * cos(A12), Q->sp1) - Q->hz0); + Q->dlam2 *= .5; + Q->lamc = M_HALFPI - atan2(sin(A12) * Q->sp1, cos(A12)) - Q->dlam2; + Q->thz0 = tan (Q->hz0); + Q->rhshz0 = .5 / sin (Q->hz0); + Q->r2z0 = 0.5 / Q->z02; + Q->z02 *= Q->z02; P->inv = s_inverse; P->fwd = s_forward; - P->es = 0.; + P->es = 0.; return P; } diff --git a/src/PJ_urm5.c b/src/PJ_urm5.c index ab853e0e..9d904fc7 100644 --- a/src/PJ_urm5.c +++ b/src/PJ_urm5.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(urm5, "Urmaev V") "\n\tPCyl., Sph., no inv.\n\tn= q= alpha="; @@ -44,10 +45,14 @@ PJ *PROJECTION(urm5) { if (pj_param(P->ctx, P->params, "tn").i) { Q->n = pj_param(P->ctx, P->params, "dn").f; - if (Q->n <= 0. || Q->n > 1.) - E_ERROR(-40) - } else - E_ERROR(-40) + if (Q->n <= 0. || Q->n > 1.) { + proj_errno_set(P, PJD_ERR_N_OUT_OF_RANGE); + return freeup_new(0); + } + } else { + proj_errno_set(P, PJD_ERR_N_OUT_OF_RANGE); + return freeup_new(0); + } Q->q3 = pj_param(P->ctx, P->params, "dq").f / 3.; alpha = pj_param(P->ctx, P->params, "ralpha").f; t = Q->n * sin (alpha); diff --git a/src/PJ_urmfps.c b/src/PJ_urmfps.c index ff775c25..fcc7b853 100644 --- a/src/PJ_urmfps.c +++ b/src/PJ_urmfps.c @@ -1,11 +1,12 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(urmfps, "Urmaev Flat-Polar Sinusoidal") "\n\tPCyl, Sph.\n\tn="; PROJ_HEAD(wag1, "Wagner I (Kavraisky VI)") "\n\tPCyl, Sph."; struct pj_opaque { - double n, C_y; + double n, C_y; }; #define C_x 0.8773826753 @@ -14,19 +15,19 @@ struct pj_opaque { static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ XY xy = {0.0, 0.0}; - lp.phi = aasin (P->ctx,P->opaque->n * sin (lp.phi)); - xy.x = C_x * lp.lam * cos (lp.phi); - xy.y = P->opaque->C_y * lp.phi; - return xy; + lp.phi = aasin (P->ctx,P->opaque->n * sin (lp.phi)); + xy.x = C_x * lp.lam * cos (lp.phi); + xy.y = P->opaque->C_y * lp.phi; + return xy; } static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0, 0.0}; - xy.y /= P->opaque->C_y; - lp.phi = aasin(P->ctx, sin (xy.y) / P->opaque->n); - lp.lam = xy.x / (C_x * cos (xy.y)); - return lp; + xy.y /= P->opaque->C_y; + lp.phi = aasin(P->ctx, sin (xy.y) / P->opaque->n); + lp.lam = xy.x / (C_x * cos (xy.y)); + return lp; } @@ -45,11 +46,11 @@ static void freeup (PJ *P) { } static PJ *setup(PJ *P) { - P->opaque->C_y = Cy / P->opaque->n; - P->es = 0.; - P->inv = s_inverse; - P->fwd = s_forward; - return P; + P->opaque->C_y = Cy / P->opaque->n; + P->es = 0.; + P->inv = s_inverse; + P->fwd = s_forward; + return P; } @@ -59,12 +60,16 @@ PJ *PROJECTION(urmfps) { return freeup_new (P); P->opaque = Q; - if (pj_param(P->ctx, P->params, "tn").i) { - P->opaque->n = pj_param(P->ctx, P->params, "dn").f; - if (P->opaque->n <= 0. || P->opaque->n > 1.) - E_ERROR(-40) - } else - E_ERROR(-40) + if (pj_param(P->ctx, P->params, "tn").i) { + P->opaque->n = pj_param(P->ctx, P->params, "dn").f; + if (P->opaque->n <= 0. || P->opaque->n > 1.) { + proj_errno_set(P, PJD_ERR_N_OUT_OF_RANGE); + return freeup_new(P); + } + } else { + proj_errno_set(P, PJD_ERR_N_OUT_OF_RANGE); + return freeup_new(P); + } return setup(P); } @@ -76,7 +81,7 @@ PJ *PROJECTION(wag1) { return freeup_new (P); P->opaque = Q; - P->opaque->n = 0.8660254037844386467637231707; + P->opaque->n = 0.8660254037844386467637231707; return setup(P); } diff --git a/src/PJ_vandg.c b/src/PJ_vandg.c index 18d75f12..282eb661 100644 --- a/src/PJ_vandg.c +++ b/src/PJ_vandg.c @@ -1,5 +1,6 @@ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" PROJ_HEAD(vandg, "van der Grinten (I)") "\n\tMisc Sph"; @@ -18,7 +19,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ double al, al2, g, g2, p2; p2 = fabs(lp.phi / M_HALFPI); - if ((p2 - TOL) > 1.) F_ERROR; + if ((p2 - TOL) > 1.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } if (p2 > 1.) p2 = 1.; if (fabs(lp.phi) <= TOL) { @@ -41,7 +45,10 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ if (lp.lam < 0.) xy.x = -xy.x; xy.y = fabs(xy.x / M_PI); xy.y = 1. - xy.y * (xy.y + 2. * al); - if (xy.y < -TOL) F_ERROR; + if (xy.y < -TOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } if (xy.y < 0.) xy.y = 0.; else @@ -81,8 +88,10 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ t = r2 + TPISQ * (x2 - y2 + HPISQ); lp.lam = fabs(xy.x) <= TOL ? 0. : .5 * (r - PISQ + (t <= 0. ? 0. : sqrt(t))) / xy.x; - } else - I_ERROR; + } else { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } return lp; } diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c index bd94551d..74a7eea0 100644 --- a/src/pj_strerrno.c +++ b/src/pj_strerrno.c @@ -60,6 +60,7 @@ pj_err_list[] = { "invalid scale", /* -52 */ "non-convergent computation", /* -53 */ "missing required arguments", /* -54 */ + "lat_0 = 0", /* -55 */ }; char *pj_strerrno(int err) { diff --git a/src/proj_etmerc.c b/src/proj_etmerc.c index cbca9bbb..ff466ea8 100644 --- a/src/proj_etmerc.c +++ b/src/proj_etmerc.c @@ -42,7 +42,8 @@ #define PROJ_LIB__ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { @@ -57,7 +58,7 @@ struct pj_opaque { PROJ_HEAD(etmerc, "Extended Transverse Mercator") "\n\tCyl, Sph\n\tlat_ts=(0)\nlat_0=(0)"; PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") - "\n\tCyl, Sph\n\tzone= south"; + "\n\tCyl, Sph\n\tzone= south"; #define PROJ_ETMERC_ORDER 6 @@ -250,8 +251,10 @@ static PJ *setup(PJ *P) { /* general initialization */ double f, n, np, Z; struct pj_opaque *Q = P->opaque; - if (P->es <= 0) - E_ERROR(-34); + if (P->es <= 0) { + proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + return freeup_new(P); + } /* flattening */ f = P->es / (1 + sqrt (1 - P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */ @@ -329,7 +332,7 @@ static PJ *setup(PJ *P) { /* general initialization */ Q->Zb = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z)); P->inv = e_inverse; P->fwd = e_forward; - return P; + return P; } @@ -405,35 +408,39 @@ int pj_etmerc_selftest (void) { PJ *PROJECTION(utm) { - int zone; + int zone; struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque)); if (0==Q) return freeup_new (P); P->opaque = Q; - if (P->es == 0.0) - E_ERROR(-34); - P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.; - P->x0 = 500000.; - if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */ - { - zone = pj_param(P->ctx, P->params, "izone").i; - if (zone > 0 && zone <= 60) - --zone; - else - E_ERROR(-35) - } - else /* nearest central meridian input */ - { - zone = (int)(floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)); - if (zone < 0) - zone = 0; - else if (zone >= 60) - zone = 59; - } - P->lam0 = (zone + .5) * M_PI / 30. - M_PI; - P->k0 = 0.9996; - P->phi0 = 0.; + if (P->es == 0.0) { + proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + return freeup_new(P); + } + P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.; + P->x0 = 500000.; + if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */ + { + zone = pj_param(P->ctx, P->params, "izone").i; + if (zone > 0 && zone <= 60) + --zone; + else { + proj_errno_set(P, PJD_ERR_INVALID_UTM_ZONE); + return freeup_new(P); + } + } + else /* nearest central meridian input */ + { + zone = (int)(floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)); + if (zone < 0) + zone = 0; + else if (zone >= 60) + zone = 59; + } + P->lam0 = (zone + .5) * M_PI / 30. - M_PI; + P->k0 = 0.9996; + P->phi0 = 0.; return setup (P); } diff --git a/src/proj_rouss.c b/src/proj_rouss.c index 264b1b77..b33e7926 100644 --- a/src/proj_rouss.c +++ b/src/proj_rouss.c @@ -24,7 +24,8 @@ ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #define PJ_LIB__ -#include <projects.h> +#include <proj.h> +#include "projects.h" struct pj_opaque { double s0; @@ -105,7 +106,7 @@ PJ *PROJECTION(rouss) { P->opaque = Q; if (!((Q->en = proj_mdist_ini(P->es)))) - E_ERROR_0; + return freeup_new(P); es2 = sin(P->phi0); Q->s0 = proj_mdist(P->phi0, es2, cos(P->phi0), Q->en); t = 1. - (es2 = P->es * es2 * es2); diff --git a/src/projects.h b/src/projects.h index a2dc10b9..4be77d9b 100644 --- a/src/projects.h +++ b/src/projects.h @@ -523,15 +523,37 @@ struct FACTORS { #define PJD_WGS84 4 /* WGS84 (or anything considered equivalent) */ /* library errors */ -#define PJD_ERR_NO_ARGS -1 -#define PJD_ERR_INVALID_M_OR_N -39 -#define PJD_ERR_GEOCENTRIC -45 -#define PJD_ERR_AXIS -47 -#define PJD_ERR_GRID_AREA -48 -#define PJD_ERR_CATALOG -49 -#define PJD_ERR_INVALID_SCALE -52 -#define PJD_ERR_NON_CONVERGENT -53 -#define PJD_ERR_MISSING_ARGS -54 +#define PJD_ERR_NO_ARGS -1 +#define PJD_ERR_LAT_OR_LON_EXCEED_LIMIT -14 +#define PJD_ERR_TOLERANCE_CONDITION -20 +#define PJD_ERR_CONIC_LAT_EQUAL -21 +#define PJD_ERR_LAT_LARGER_THAN_90 -22 +#define PJD_ERR_LAT1_IS_ZERO -23 +#define PJD_ERR_LAT_TS_LARGER_THAN_90 -24 +#define PJD_ERR_CONTROL_POINT_NO_DIST -25 +#define PJD_ERR_NO_ROTATION_PROJ -26 +#define PJD_ERR_W_OR_M_ZERO_OR_LESS -27 +#define PJD_ERR_LSAT_NOT_IN_RANGE -28 +#define PJD_ERR_PATH_NOT_IN_RANGE -29 +#define PJD_ERR_H_LESS_THAN_ZERO -30 +#define PJD_ERR_LAT_1_OR_2_ZERO_OR_90 -32 +#define PJD_ERR_LAT_0_OR_ALPHA_EQ_90 -33 +#define PJD_ERR_ELLIPSOID_USE_REQUIRED -34 +#define PJD_ERR_INVALID_UTM_ZONE -35 +#define PJD_ERR_FAILED_TO_FIND_PROJ -37 +#define PJD_ERR_INVALID_M_OR_N -39 +#define PJD_ERR_N_OUT_OF_RANGE -40 +#define PJD_ERR_ABS_LAT1_EQ_ABS_LAT2 -42 +#define PJD_ERR_LAT_0_HALF_PI_FROM_MEAN -43 +#define PJD_ERR_GEOCENTRIC -45 +#define PJD_ERR_UNKNOWN_PRIME_MERIDIAN -46 +#define PJD_ERR_AXIS -47 +#define PJD_ERR_GRID_AREA -48 +#define PJD_ERR_INVALID_SWEEP_AXIS -49 +#define PJD_ERR_INVALID_SCALE -52 +#define PJD_ERR_NON_CONVERGENT -53 +#define PJD_ERR_MISSING_ARGS -54 +#define PJD_ERR_LAT_0_IS_ZERO -55 struct projFileAPI_t; @@ -582,12 +604,6 @@ extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[]; #ifdef PJ_LIB__ #define PROJ_HEAD(id, name) static const char des_##id [] = name -#define E_ERROR(err) { pj_ctx_set_errno( P->ctx, err); freeup(P); return(0); } -#define E_ERROR_0 { freeup(P); return(0); } -#define F_ERROR { pj_ctx_set_errno( P->ctx, -20); return(xy); } -#define F3_ERROR { pj_ctx_set_errno( P->ctx, -20); return(xyz); } -#define I_ERROR { pj_ctx_set_errno( P->ctx, -20); return(lp); } -#define I3_ERROR { pj_ctx_set_errno( P->ctx, -20); return(lpz); } #define PROJECTION(name) \ pj_projection_specific_setup_##name (PJ *P); \ |