"Fossies" - the Fresh Open Source Software Archive

Member "proj-6.2.1/src/projections/ortho.cpp" (28 Oct 2019, 3902 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 "ortho.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(ortho, "Orthographic") "\n\tAzi, Sph";
    8 
    9 namespace { // anonymous namespace
   10 enum Mode {
   11     N_POLE = 0,
   12     S_POLE = 1,
   13     EQUIT  = 2,
   14     OBLIQ  = 3
   15 };
   16 } // anonymous namespace
   17 
   18 namespace { // anonymous namespace
   19 struct pj_opaque {
   20     double  sinph0;
   21     double  cosph0;
   22     enum Mode mode;
   23 };
   24 } // anonymous namespace
   25 
   26 #define EPS10 1.e-10
   27 
   28 static PJ_XY forward_error(PJ *P, PJ_LP lp, PJ_XY xy) {
   29     proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
   30     proj_log_trace(P, "Coordinate (%.3f, %.3f) is on the unprojected hemisphere",
   31                    proj_todeg(lp.lam), proj_todeg(lp.phi));
   32     return xy;
   33 }
   34 
   35 static PJ_XY ortho_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
   36     PJ_XY xy;
   37     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
   38     double  coslam, cosphi, sinphi;
   39 
   40     xy.x = HUGE_VAL; xy.y = HUGE_VAL;
   41 
   42     cosphi = cos(lp.phi);
   43     coslam = cos(lp.lam);
   44     switch (Q->mode) {
   45     case EQUIT:
   46         if (cosphi * coslam < - EPS10)
   47             return forward_error(P, lp, xy);
   48         xy.y = sin(lp.phi);
   49         break;
   50     case OBLIQ:
   51         if (Q->sinph0 * (sinphi = sin(lp.phi)) + Q->cosph0 * cosphi * coslam < - EPS10)
   52             return forward_error(P, lp, xy);
   53         xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam;
   54         break;
   55     case N_POLE:
   56         coslam = - coslam;
   57                 /*-fallthrough*/
   58     case S_POLE:
   59         if (fabs(lp.phi - P->phi0) - EPS10 > M_HALFPI)
   60             return forward_error(P, lp, xy);
   61         xy.y = cosphi * coslam;
   62         break;
   63     }
   64     xy.x = cosphi * sin(lp.lam);
   65     return xy;
   66 }
   67 
   68 
   69 static PJ_LP ortho_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
   70     PJ_LP lp;
   71     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
   72     double  rh, cosc, sinc;
   73 
   74     lp.lam = HUGE_VAL; lp.phi = HUGE_VAL;
   75 
   76     if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) {
   77         if ((sinc - 1.) > EPS10) {
   78             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
   79             proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary");
   80             return lp;
   81         }
   82         sinc = 1.;
   83     }
   84     cosc = sqrt(1. - sinc * sinc); /* in this range OK */
   85     if (fabs(rh) <= EPS10) {
   86         lp.phi = P->phi0;
   87         lp.lam = 0.0;
   88     } else {
   89         switch (Q->mode) {
   90         case N_POLE:
   91             xy.y = -xy.y;
   92             lp.phi = acos(sinc);
   93             break;
   94         case S_POLE:
   95             lp.phi = - acos(sinc);
   96             break;
   97         case EQUIT:
   98             lp.phi = xy.y * sinc / rh;
   99             xy.x *= sinc;
  100             xy.y = cosc * rh;
  101             goto sinchk;
  102         case OBLIQ:
  103             lp.phi = cosc * Q->sinph0 + xy.y * sinc * Q->cosph0 /rh;
  104             xy.y = (cosc - Q->sinph0 * lp.phi) * rh;
  105             xy.x *= sinc * Q->cosph0;
  106         sinchk:
  107             if (fabs(lp.phi) >= 1.)
  108                 lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
  109             else
  110                 lp.phi = asin(lp.phi);
  111             break;
  112         }
  113         lp.lam = (xy.y == 0. && (Q->mode == OBLIQ || Q->mode == EQUIT))
  114              ? (xy.x == 0. ? 0. : xy.x < 0. ? -M_HALFPI : M_HALFPI)
  115                            : atan2(xy.x, xy.y);
  116     }
  117     return lp;
  118 }
  119 
  120 
  121 
  122 PJ *PROJECTION(ortho) {
  123     struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
  124     if (nullptr==Q)
  125         return pj_default_destructor(P, ENOMEM);
  126     P->opaque = Q;
  127 
  128     if (fabs(fabs(P->phi0) - M_HALFPI) <= EPS10)
  129         Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
  130     else if (fabs(P->phi0) > EPS10) {
  131         Q->mode = OBLIQ;
  132         Q->sinph0 = sin(P->phi0);
  133         Q->cosph0 = cos(P->phi0);
  134     } else
  135         Q->mode = EQUIT;
  136     P->inv = ortho_s_inverse;
  137     P->fwd = ortho_s_forward;
  138     P->es = 0.;
  139 
  140     return P;
  141 }
  142