diff options
Diffstat (limited to 'src/projections/stere.cpp')
| -rw-r--r-- | src/projections/stere.cpp | 35 |
1 files changed, 21 insertions, 14 deletions
diff --git a/src/projections/stere.cpp b/src/projections/stere.cpp index 9b24a596..683d484c 100644 --- a/src/projections/stere.cpp +++ b/src/projections/stere.cpp @@ -41,7 +41,7 @@ static double ssfn_ (double phit, double sinphi, double eccen) { } -static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ +static PJ_XY stere_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ PJ_XY xy = {0.0,0.0}; struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A = 0.0, sinphi; @@ -55,11 +55,18 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ } switch (Q->mode) { - case OBLIQ: - A = Q->akm1 / (Q->cosX1 * (1. + Q->sinX1 * sinX + - Q->cosX1 * cosX * coslam)); + case OBLIQ: { + const double denom = Q->cosX1 * (1. + Q->sinX1 * sinX + + Q->cosX1 * cosX * coslam); + if( denom == 0 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return proj_coord_error().xy; + } + A = Q->akm1 / denom; xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam); - goto xmul; /* but why not just xy.x = A * cosX; break; ? */ + xy.x = A * cosX; + break; + } case EQUIT: /* avoid zero division */ @@ -69,7 +76,6 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ A = Q->akm1 / (1. + cosX * coslam); xy.y = A * sinX; } -xmul: xy.x = A * cosX; break; @@ -89,7 +95,7 @@ xmul: } -static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ +static PJ_XY stere_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ PJ_XY xy = {0.0,0.0}; struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); double sinphi, cosphi, coslam, sinlam; @@ -131,7 +137,7 @@ oblcon: } -static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ +static PJ_LP stere_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ PJ_LP lp = {0.0,0.0}; struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0; @@ -165,7 +171,7 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ break; } - for (i = NITER; i--; phi_l = lp.phi) { + for (i = NITER; i--; ) { sinphi = P->e * sin(phi_l); lp.phi = 2. * atan (tp * pow ((1.+sinphi)/(1.-sinphi), halfe)) - halfpi; if (fabs (phi_l - lp.phi) < CONV) { @@ -174,6 +180,7 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y); return lp; } + phi_l = lp.phi; } proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); @@ -181,7 +188,7 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ } -static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ +static PJ_LP stere_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ PJ_LP lp = {0.0,0.0}; struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); double c, rh, sinc, cosc; @@ -258,8 +265,8 @@ static PJ *setup(PJ *P) { /* general initialization */ Q->cosX1 = cos (X); break; } - P->inv = e_inverse; - P->fwd = e_forward; + P->inv = stere_e_inverse; + P->fwd = stere_e_forward; } else { switch (Q->mode) { case OBLIQ: @@ -277,8 +284,8 @@ static PJ *setup(PJ *P) { /* general initialization */ break; } - P->inv = s_inverse; - P->fwd = s_forward; + P->inv = stere_s_inverse; + P->fwd = stere_s_forward; } return P; } |