diff options
Diffstat (limited to 'src/projections/ocea.cpp')
| -rw-r--r-- | src/projections/ocea.cpp | 47 |
1 files changed, 29 insertions, 18 deletions
diff --git a/src/projections/ocea.cpp b/src/projections/ocea.cpp index 75aa6666..646b8638 100644 --- a/src/projections/ocea.cpp +++ b/src/projections/ocea.cpp @@ -15,13 +15,11 @@ struct pj_opaque { double rtk; double sinphi; double cosphi; - double singam; - double cosgam; }; } // anonymous namespace -static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ +static PJ_XY ocea_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 t; @@ -36,7 +34,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ } -static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ +static PJ_LP ocea_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 t, s; @@ -52,7 +50,7 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ PJ *PROJECTION(ocea) { - double phi_0=0.0, phi_1, phi_2, lam_1, lam_2, lonz, alpha; + double phi_1, phi_2, lam_1, lam_2, lonz, alpha; struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); if (nullptr==Q) @@ -61,15 +59,21 @@ PJ *PROJECTION(ocea) { Q->rok = 1. / P->k0; Q->rtk = P->k0; + double lam_p, phi_p; /*If the keyword "alpha" is found in the sentence then use 1point+1azimuth*/ if ( pj_param(P->ctx, P->params, "talpha").i) { /*Define Pole of oblique transformation from 1 point & 1 azimuth*/ - alpha = pj_param(P->ctx, P->params, "ralpha").f; + // ERO: I've added M_PI so that the alpha is the angle from point 1 to point 2 + // from the North in a clockwise direction + // (to be consistent with omerc behaviour) + alpha = M_PI + pj_param(P->ctx, P->params, "ralpha").f; lonz = pj_param(P->ctx, P->params, "rlonc").f; /*Equation 9-8 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ - Q->singam = atan(-cos(alpha)/(-sin(phi_0) * sin(alpha))) + lonz; + // Actually slightliy modified to use atan2(), as it is suggested by + // Snyder for equation 9-1, but this is not mentionned here + lam_p = atan2(-cos(alpha) , -sin(P->phi0) * sin(alpha)) + lonz; /*Equation 9-7 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ - Q->sinphi = asin(cos(phi_0) * sin(alpha)); + phi_p = asin(cos(P->phi0) * sin(alpha)); /*If the keyword "alpha" is NOT found in the sentence then use 2points*/ } else { /*Define Pole of oblique transformation from 2 points*/ @@ -78,25 +82,32 @@ PJ *PROJECTION(ocea) { lam_1 = pj_param(P->ctx, P->params, "rlon_1").f; lam_2 = pj_param(P->ctx, P->params, "rlon_2").f; /*Equation 9-1 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ - Q->singam = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) - + lam_p = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) - sin(phi_1) * cos(phi_2) * cos(lam_2), sin(phi_1) * cos(phi_2) * sin(lam_2) - cos(phi_1) * sin(phi_2) * sin(lam_1) ); /* take care of P->lam0 wrap-around when +lam_1=-90*/ if (lam_1 == -M_HALFPI) - Q->singam = -Q->singam; + lam_p = -lam_p; /*Equation 9-2 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ - Q->sinphi = atan(-cos(Q->singam - lam_1) / tan(phi_1)); + double cos_lamp_m_minus_lam_1 = cos(lam_p - lam_1); + double tan_phi_1 = tan(phi_1); + if( tan_phi_1 == 0.0 ) { + // Not sure if we want to support this case, but at least this avoids + // a division by zero, and gives the same result as the below atan() + phi_p = (cos_lamp_m_minus_lam_1 >= 0.0 ) ? -M_HALFPI : M_HALFPI; + } + else { + phi_p = atan(- cos_lamp_m_minus_lam_1 / tan_phi_1); + } } - P->lam0 = Q->singam + M_HALFPI; - Q->cosphi = cos(Q->sinphi); - Q->sinphi = sin(Q->sinphi); - Q->cosgam = cos(Q->singam); - Q->singam = sin(Q->singam); - P->inv = s_inverse; - P->fwd = s_forward; + P->lam0 = lam_p + M_HALFPI; + Q->cosphi = cos(phi_p); + Q->sinphi = sin(phi_p); + P->inv = ocea_s_inverse; + P->fwd = ocea_s_forward; P->es = 0.; return P; |