| 50 | |
| 51 | |
| 52 | PJ *PROJECTION(ocea) { |
| 53 | double phi_1, phi_2, lam_1, lam_2, lonz, alpha; |
| 54 | |
| 55 | struct pj_opaque *Q = static_cast<struct pj_opaque*>(calloc (1, sizeof (struct pj_opaque))); |
| 56 | if (nullptr==Q) |
| 57 | return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/); |
| 58 | P->opaque = Q; |
| 59 | |
| 60 | Q->rok = 1. / P->k0; |
| 61 | Q->rtk = P->k0; |
| 62 | double lam_p, phi_p; |
| 63 | /*If the keyword "alpha" is found in the sentence then use 1point+1azimuth*/ |
| 64 | if ( pj_param(P->ctx, P->params, "talpha").i) { |
| 65 | /*Define Pole of oblique transformation from 1 point & 1 azimuth*/ |
| 66 | // ERO: I've added M_PI so that the alpha is the angle from point 1 to point 2 |
| 67 | // from the North in a clockwise direction |
| 68 | // (to be consistent with omerc behavior) |
| 69 | alpha = M_PI + pj_param(P->ctx, P->params, "ralpha").f; |
| 70 | lonz = pj_param(P->ctx, P->params, "rlonc").f; |
| 71 | /*Equation 9-8 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ |
| 72 | // Actually slightliy modified to use atan2(), as it is suggested by |
| 73 | // Snyder for equation 9-1, but this is not mentioned here |
| 74 | lam_p = atan2(-cos(alpha) , -sin(P->phi0) * sin(alpha)) + lonz; |
| 75 | /*Equation 9-7 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ |
| 76 | phi_p = asin(cos(P->phi0) * sin(alpha)); |
| 77 | /*If the keyword "alpha" is NOT found in the sentence then use 2points*/ |
| 78 | } else { |
| 79 | /*Define Pole of oblique transformation from 2 points*/ |
| 80 | phi_1 = pj_param(P->ctx, P->params, "rlat_1").f; |
| 81 | phi_2 = pj_param(P->ctx, P->params, "rlat_2").f; |
| 82 | lam_1 = pj_param(P->ctx, P->params, "rlon_1").f; |
| 83 | lam_2 = pj_param(P->ctx, P->params, "rlon_2").f; |
| 84 | /*Equation 9-1 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ |
| 85 | lam_p = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) - |
| 86 | sin(phi_1) * cos(phi_2) * cos(lam_2), |
| 87 | sin(phi_1) * cos(phi_2) * sin(lam_2) - |
| 88 | cos(phi_1) * sin(phi_2) * sin(lam_1) ); |
| 89 | |
| 90 | /* take care of P->lam0 wrap-around when +lam_1=-90*/ |
| 91 | if (lam_1 == -M_HALFPI) |
| 92 | lam_p = -lam_p; |
| 93 | |
| 94 | /*Equation 9-2 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ |
| 95 | double cos_lamp_m_minus_lam_1 = cos(lam_p - lam_1); |
| 96 | double tan_phi_1 = tan(phi_1); |
| 97 | if( tan_phi_1 == 0.0 ) { |
| 98 | // Not sure if we want to support this case, but at least this avoids |
| 99 | // a division by zero, and gives the same result as the below atan() |
| 100 | phi_p = (cos_lamp_m_minus_lam_1 >= 0.0 ) ? -M_HALFPI : M_HALFPI; |
| 101 | } |
| 102 | else { |
| 103 | phi_p = atan(- cos_lamp_m_minus_lam_1 / tan_phi_1); |
| 104 | } |
| 105 | } |
| 106 | P->lam0 = lam_p + M_HALFPI; |
| 107 | Q->cosphi = cos(phi_p); |
| 108 | Q->sinphi = sin(phi_p); |
| 109 | P->inv = ocea_s_inverse; |