"Fossies" - the Fresh Open Source Software Archive

Member "proj-6.2.1/src/projections/stere.cpp" (28 Oct 2019, 9130 Bytes) of package /linux/privat/proj-6.2.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "stere.cpp" see the Fossies "Dox" file reference documentation.

    1 #define PJ_LIB__
    2 #include <errno.h>
    3 #include "proj.h"
    4 #include "proj_internal.h"
    5 #include "proj_math.h"
    6 
    7 PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
    8 PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";
    9 
   10 
   11 namespace { // anonymous namespace
   12 enum Mode {
   13     S_POLE = 0,
   14     N_POLE = 1,
   15     OBLIQ  = 2,
   16     EQUIT  = 3
   17 };
   18 } // anonymous namespace
   19 
   20 namespace { // anonymous namespace
   21 struct pj_opaque {
   22     double phits;
   23     double sinX1;
   24     double cosX1;
   25     double akm1;
   26     enum Mode mode;
   27 };
   28 } // anonymous namespace
   29 
   30 #define sinph0  static_cast<struct pj_opaque*>(P->opaque)->sinX1
   31 #define cosph0  static_cast<struct pj_opaque*>(P->opaque)->cosX1
   32 #define EPS10   1.e-10
   33 #define TOL 1.e-8
   34 #define NITER   8
   35 #define CONV    1.e-10
   36 
   37 static double ssfn_ (double phit, double sinphi, double eccen) {
   38     sinphi *= eccen;
   39     return (tan (.5 * (M_HALFPI + phit)) *
   40        pow ((1. - sinphi) / (1. + sinphi), .5 * eccen));
   41 }
   42 
   43 
   44 static PJ_XY stere_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
   45     PJ_XY xy = {0.0,0.0};
   46     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
   47     double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A = 0.0, sinphi;
   48 
   49     coslam = cos (lp.lam);
   50     sinlam = sin (lp.lam);
   51     sinphi = sin (lp.phi);
   52     if (Q->mode == OBLIQ || Q->mode == EQUIT) {
   53         sinX = sin (X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - M_HALFPI);
   54         cosX = cos (X);
   55     }
   56 
   57     switch (Q->mode) {
   58     case OBLIQ: {
   59         const double denom = Q->cosX1 * (1. + Q->sinX1 * sinX +
   60            Q->cosX1 * cosX * coslam);
   61         if( denom == 0 ) {
   62             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
   63             return proj_coord_error().xy;
   64         }
   65         A = Q->akm1 / denom;
   66         xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam);
   67         xy.x = A * cosX;
   68         break;
   69     }
   70 
   71     case EQUIT:
   72         /* avoid zero division */
   73         if (1. + cosX * coslam == 0.0) {
   74             xy.y = HUGE_VAL;
   75         } else {
   76             A = Q->akm1 / (1. + cosX * coslam);
   77             xy.y = A * sinX;
   78         }
   79         xy.x = A * cosX;
   80         break;
   81 
   82     case S_POLE:
   83         lp.phi = -lp.phi;
   84         coslam = - coslam;
   85         sinphi = -sinphi;
   86         /*-fallthrough*/
   87     case N_POLE:
   88         xy.x = Q->akm1 * pj_tsfn (lp.phi, sinphi, P->e);
   89         xy.y = - xy.x * coslam;
   90         break;
   91     }
   92 
   93     xy.x = xy.x * sinlam;
   94     return xy;
   95 }
   96 
   97 
   98 static PJ_XY stere_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
   99     PJ_XY xy = {0.0,0.0};
  100     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
  101     double  sinphi, cosphi, coslam, sinlam;
  102 
  103     sinphi = sin(lp.phi);
  104     cosphi = cos(lp.phi);
  105     coslam = cos(lp.lam);
  106     sinlam = sin(lp.lam);
  107 
  108     switch (Q->mode) {
  109     case EQUIT:
  110         xy.y = 1. + cosphi * coslam;
  111         goto oblcon;
  112     case OBLIQ:
  113         xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
  114 oblcon:
  115         if (xy.y <= EPS10) {
  116             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
  117             return xy;
  118         }
  119         xy.x = (xy.y = Q->akm1 / xy.y) * cosphi * sinlam;
  120         xy.y *= (Q->mode == EQUIT) ? sinphi :
  121            cosph0 * sinphi - sinph0 * cosphi * coslam;
  122         break;
  123     case N_POLE:
  124         coslam = - coslam;
  125         lp.phi = - lp.phi;
  126         /*-fallthrough*/
  127     case S_POLE:
  128         if (fabs (lp.phi - M_HALFPI) < TOL) {
  129             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
  130             return xy;
  131         }
  132         xy.x = sinlam * ( xy.y = Q->akm1 * tan (M_FORTPI + .5 * lp.phi) );
  133         xy.y *= coslam;
  134         break;
  135     }
  136     return xy;
  137 }
  138 
  139 
  140 static PJ_LP stere_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
  141     PJ_LP lp = {0.0,0.0};
  142     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
  143     double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
  144     int i;
  145 
  146     rho = hypot (xy.x, xy.y);
  147 
  148     switch (Q->mode) {
  149     case OBLIQ:
  150     case EQUIT:
  151         cosphi = cos ( tp = 2. * atan2 (rho * Q->cosX1 , Q->akm1) );
  152         sinphi = sin (tp);
  153                 if ( rho == 0.0 )
  154             phi_l = asin (cosphi * Q->sinX1);
  155                 else
  156             phi_l = asin (cosphi * Q->sinX1 + (xy.y * sinphi * Q->cosX1 / rho));
  157 
  158         tp = tan (.5 * (M_HALFPI + phi_l));
  159         xy.x *= sinphi;
  160         xy.y = rho * Q->cosX1 * cosphi - xy.y * Q->sinX1* sinphi;
  161         halfpi = M_HALFPI;
  162         halfe = .5 * P->e;
  163         break;
  164     case N_POLE:
  165         xy.y = -xy.y;
  166         /*-fallthrough*/
  167     case S_POLE:
  168         phi_l = M_HALFPI - 2. * atan (tp = - rho / Q->akm1);
  169         halfpi = -M_HALFPI;
  170         halfe = -.5 * P->e;
  171         break;
  172     }
  173 
  174     for (i = NITER;  i--; ) {
  175         sinphi = P->e * sin(phi_l);
  176         lp.phi = 2. * atan (tp * pow ((1.+sinphi)/(1.-sinphi), halfe)) - halfpi;
  177         if (fabs (phi_l - lp.phi) < CONV) {
  178             if (Q->mode == S_POLE)
  179                 lp.phi = -lp.phi;
  180             lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y);
  181             return lp;
  182         }
  183         phi_l = lp.phi;
  184     }
  185 
  186     proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
  187     return lp;
  188 }
  189 
  190 
  191 static PJ_LP stere_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
  192     PJ_LP lp = {0.0,0.0};
  193     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
  194     double  c, rh, sinc, cosc;
  195 
  196     sinc = sin (c = 2. * atan ((rh = hypot (xy.x, xy.y)) / Q->akm1));
  197     cosc = cos (c);
  198     lp.lam = 0.;
  199 
  200     switch (Q->mode) {
  201     case EQUIT:
  202         if (fabs (rh) <= EPS10)
  203             lp.phi = 0.;
  204         else
  205             lp.phi = asin (xy.y * sinc / rh);
  206         if (cosc != 0. || xy.x != 0.)
  207             lp.lam = atan2 (xy.x * sinc, cosc * rh);
  208         break;
  209     case OBLIQ:
  210         if (fabs (rh) <= EPS10)
  211             lp.phi = P->phi0;
  212         else
  213             lp.phi = asin (cosc * sinph0 + xy.y * sinc * cosph0 / rh);
  214         if ((c = cosc - sinph0 * sin (lp.phi)) != 0. || xy.x != 0.)
  215             lp.lam = atan2 (xy.x * sinc * cosph0, c * rh);
  216         break;
  217     case N_POLE:
  218         xy.y = -xy.y;
  219         /*-fallthrough*/
  220     case S_POLE:
  221         if (fabs (rh) <= EPS10)
  222             lp.phi = P->phi0;
  223         else
  224             lp.phi = asin (Q->mode == S_POLE ? - cosc : cosc);
  225         lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y);
  226         break;
  227     }
  228     return lp;
  229 }
  230 
  231 
  232 static PJ *setup(PJ *P) {                   /* general initialization */
  233     double t;
  234     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
  235 
  236     if (fabs ((t = fabs (P->phi0)) - M_HALFPI) < EPS10)
  237         Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
  238     else
  239         Q->mode = t > EPS10 ? OBLIQ : EQUIT;
  240     Q->phits = fabs (Q->phits);
  241 
  242     if (P->es != 0.0) {
  243         double X;
  244 
  245         switch (Q->mode) {
  246         case N_POLE:
  247         case S_POLE:
  248             if (fabs (Q->phits - M_HALFPI) < EPS10)
  249                 Q->akm1 = 2. * P->k0 /
  250                    sqrt (pow (1+P->e,1+P->e) * pow (1-P->e,1-P->e));
  251             else {
  252                 Q->akm1 = cos (Q->phits) /
  253                    pj_tsfn (Q->phits, t = sin (Q->phits), P->e);
  254                 t *= P->e;
  255                 Q->akm1 /= sqrt(1. - t * t);
  256             }
  257             break;
  258         case EQUIT:
  259         case OBLIQ:
  260             t = sin (P->phi0);
  261             X = 2. * atan (ssfn_(P->phi0, t, P->e)) - M_HALFPI;
  262             t *= P->e;
  263             Q->akm1 = 2. * P->k0 * cos (P->phi0) / sqrt(1. - t * t);
  264             Q->sinX1 = sin (X);
  265             Q->cosX1 = cos (X);
  266             break;
  267         }
  268         P->inv = stere_e_inverse;
  269         P->fwd = stere_e_forward;
  270     } else {
  271         switch (Q->mode) {
  272         case OBLIQ:
  273             sinph0 = sin (P->phi0);
  274             cosph0 = cos (P->phi0);
  275             /*-fallthrough*/
  276         case EQUIT:
  277             Q->akm1 = 2. * P->k0;
  278             break;
  279         case S_POLE:
  280         case N_POLE:
  281             Q->akm1 = fabs (Q->phits - M_HALFPI) >= EPS10 ?
  282                cos (Q->phits) / tan (M_FORTPI - .5 * Q->phits) :
  283                2. * P->k0 ;
  284             break;
  285         }
  286 
  287         P->inv = stere_s_inverse;
  288         P->fwd = stere_s_forward;
  289     }
  290     return P;
  291 }
  292 
  293 
  294 PJ *PROJECTION(stere) {
  295     struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
  296     if (nullptr==Q)
  297         return pj_default_destructor (P, ENOMEM);
  298     P->opaque = Q;
  299 
  300     Q->phits = pj_param (P->ctx, P->params, "tlat_ts").i ?
  301                pj_param (P->ctx, P->params, "rlat_ts").f : M_HALFPI;
  302 
  303     return setup(P);
  304 }
  305 
  306 
  307 PJ *PROJECTION(ups) {
  308     struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
  309     if (nullptr==Q)
  310         return pj_default_destructor (P, ENOMEM);
  311     P->opaque = Q;
  312 
  313     /* International Ellipsoid */
  314     P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI;
  315     if (P->es == 0.0) {
  316         proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
  317         return pj_default_destructor (P, ENOMEM);
  318     }
  319     P->k0 = .994;
  320     P->x0 = 2000000.;
  321     P->y0 = 2000000.;
  322     Q->phits = M_HALFPI;
  323     P->lam0 = 0.;
  324 
  325     return setup(P);
  326 }
  327