"Fossies" - the Fresh Open Source Software Archive

Member "proj-6.2.1/src/projections/nsper.cpp" (28 Oct 2019, 5058 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 "nsper.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 namespace { // anonymous namespace
    8 enum Mode {
    9     N_POLE = 0,
   10     S_POLE = 1,
   11     EQUIT  = 2,
   12     OBLIQ  = 3
   13 };
   14 } // anonymous namespace
   15 
   16 namespace { // anonymous namespace
   17 struct pj_opaque {
   18     double  height;
   19     double  sinph0;
   20     double  cosph0;
   21     double  p;
   22     double  rp;
   23     double  pn1;
   24     double  pfact;
   25     double  h;
   26     double  cg;
   27     double  sg;
   28     double  sw;
   29     double  cw;
   30     enum Mode mode;
   31     int     tilt;
   32 };
   33 } // anonymous namespace
   34 
   35 PROJ_HEAD(nsper, "Near-sided perspective") "\n\tAzi, Sph\n\th=";
   36 PROJ_HEAD(tpers, "Tilted perspective") "\n\tAzi, Sph\n\ttilt= azi= h=";
   37 
   38 # define EPS10 1.e-10
   39 
   40 
   41 static PJ_XY nsper_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
   42     PJ_XY xy = {0.0,0.0};
   43     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
   44     double  coslam, cosphi, sinphi;
   45 
   46     sinphi = sin(lp.phi);
   47     cosphi = cos(lp.phi);
   48     coslam = cos(lp.lam);
   49     switch (Q->mode) {
   50     case OBLIQ:
   51         xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam;
   52         break;
   53     case EQUIT:
   54         xy.y = cosphi * coslam;
   55         break;
   56     case S_POLE:
   57         xy.y = - sinphi;
   58         break;
   59     case N_POLE:
   60         xy.y = sinphi;
   61         break;
   62     }
   63     if (xy.y < Q->rp) {
   64         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
   65         return xy;
   66     }
   67     xy.y = Q->pn1 / (Q->p - xy.y);
   68     xy.x = xy.y * cosphi * sin(lp.lam);
   69     switch (Q->mode) {
   70     case OBLIQ:
   71         xy.y *= (Q->cosph0 * sinphi -
   72            Q->sinph0 * cosphi * coslam);
   73         break;
   74     case EQUIT:
   75         xy.y *= sinphi;
   76         break;
   77     case N_POLE:
   78         coslam = - coslam;
   79         /*-fallthrough*/
   80     case S_POLE:
   81         xy.y *= cosphi * coslam;
   82         break;
   83     }
   84     if (Q->tilt) {
   85         double yt, ba;
   86 
   87         yt = xy.y * Q->cg + xy.x * Q->sg;
   88         ba = 1. / (yt * Q->sw * Q->h + Q->cw);
   89         xy.x = (xy.x * Q->cg - xy.y * Q->sg) * Q->cw * ba;
   90         xy.y = yt * ba;
   91     }
   92     return xy;
   93 }
   94 
   95 
   96 static PJ_LP nsper_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
   97     PJ_LP lp = {0.0,0.0};
   98     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
   99     double  rh;
  100 
  101     if (Q->tilt) {
  102         double bm, bq, yt;
  103 
  104         yt = 1./(Q->pn1 - xy.y * Q->sw);
  105         bm = Q->pn1 * xy.x * yt;
  106         bq = Q->pn1 * xy.y * Q->cw * yt;
  107         xy.x = bm * Q->cg + bq * Q->sg;
  108         xy.y = bq * Q->cg - bm * Q->sg;
  109     }
  110     rh = hypot(xy.x, xy.y);
  111     if (fabs(rh) <= EPS10) {
  112         lp.lam = 0.;
  113         lp.phi = P->phi0;
  114     } else {
  115         double cosz, sinz;
  116         sinz = 1. - rh * rh * Q->pfact;
  117         if (sinz < 0.) {
  118             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
  119             return lp;
  120         }
  121         sinz = (Q->p - sqrt(sinz)) / (Q->pn1 / rh + rh / Q->pn1);
  122         cosz = sqrt(1. - sinz * sinz);
  123         switch (Q->mode) {
  124         case OBLIQ:
  125             lp.phi = asin(cosz * Q->sinph0 + xy.y * sinz * Q->cosph0 / rh);
  126             xy.y = (cosz - Q->sinph0 * sin(lp.phi)) * rh;
  127             xy.x *= sinz * Q->cosph0;
  128             break;
  129         case EQUIT:
  130             lp.phi = asin(xy.y * sinz / rh);
  131             xy.y = cosz * rh;
  132             xy.x *= sinz;
  133             break;
  134         case N_POLE:
  135             lp.phi = asin(cosz);
  136             xy.y = -xy.y;
  137             break;
  138         case S_POLE:
  139             lp.phi = - asin(cosz);
  140             break;
  141         }
  142         lp.lam = atan2(xy.x, xy.y);
  143     }
  144     return lp;
  145 }
  146 
  147 
  148 static PJ *setup(PJ *P) {
  149     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
  150 
  151     Q->height = pj_param(P->ctx, P->params, "dh").f;
  152 
  153     if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10)
  154         Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
  155     else if (fabs(P->phi0) < EPS10)
  156         Q->mode = EQUIT;
  157     else {
  158         Q->mode = OBLIQ;
  159         Q->sinph0 = sin(P->phi0);
  160         Q->cosph0 = cos(P->phi0);
  161     }
  162     Q->pn1 = Q->height / P->a; /* normalize by radius */
  163     if ( Q->pn1 <= 0 || Q->pn1 > 1e10 )
  164         return pj_default_destructor (P, PJD_ERR_INVALID_H);
  165     Q->p = 1. + Q->pn1;
  166     Q->rp = 1. / Q->p;
  167     Q->h = 1. / Q->pn1;
  168     Q->pfact = (Q->p + 1.) * Q->h;
  169     P->inv = nsper_s_inverse;
  170     P->fwd = nsper_s_forward;
  171     P->es = 0.;
  172 
  173     return P;
  174 }
  175 
  176 
  177 PJ *PROJECTION(nsper) {
  178     struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
  179     if (nullptr==Q)
  180         return pj_default_destructor (P, ENOMEM);
  181     P->opaque = Q;
  182 
  183     Q->tilt = 0;
  184 
  185    return setup(P);
  186 }
  187 
  188 
  189 PJ *PROJECTION(tpers) {
  190     double omega, gamma;
  191 
  192     struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
  193     if (nullptr==Q)
  194         return pj_default_destructor (P, ENOMEM);
  195     P->opaque = Q;
  196 
  197     omega = pj_param(P->ctx, P->params, "rtilt").f;
  198     gamma = pj_param(P->ctx, P->params, "razi").f;
  199     Q->tilt = 1;
  200     Q->cg = cos(gamma); Q->sg = sin(gamma);
  201     Q->cw = cos(omega); Q->sw = sin(omega);
  202 
  203     return setup(P);
  204 }
  205