/Proj4/PJ_nsper.c
C | 152 lines | 147 code | 5 blank | 0 comment | 15 complexity | 48a66bc194b7d17201ee4e63be916f63 MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_nsper.c 4.1 94/02/15 GIE REL"; 3#endif 4#define PROJ_PARMS__ \ 5 double height; \ 6 double sinph0; \ 7 double cosph0; \ 8 double p; \ 9 double rp; \ 10 double pn1; \ 11 double pfact; \ 12 double h; \ 13 double cg; \ 14 double sg; \ 15 double sw; \ 16 double cw; \ 17 int mode; \ 18 int tilt; 19#define PJ_LIB__ 20#include "projects.h" 21PROJ_HEAD(nsper, "Near-sided perspective") "\n\tAzi, Sph\n\th="; 22PROJ_HEAD(tpers, "Tilted perspective") "\n\tAzi, Sph\n\ttilt= azi= h="; 23# define EPS10 1.e-10 24# define N_POLE 0 25# define S_POLE 1 26# define EQUIT 2 27# define OBLIQ 3 28FORWARD(s_forward); /* spheroid */ 29 double coslam, cosphi, sinphi; 30 31 sinphi = sin(lp.phi); 32 cosphi = cos(lp.phi); 33 coslam = cos(lp.lam); 34 switch (P->mode) { 35 case OBLIQ: 36 xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam; 37 break; 38 case EQUIT: 39 xy.y = cosphi * coslam; 40 break; 41 case S_POLE: 42 xy.y = - sinphi; 43 break; 44 case N_POLE: 45 xy.y = sinphi; 46 break; 47 } 48 if (xy.y < P->rp) F_ERROR; 49 xy.y = P->pn1 / (P->p - xy.y); 50 xy.x = xy.y * cosphi * sin(lp.lam); 51 switch (P->mode) { 52 case OBLIQ: 53 xy.y *= (P->cosph0 * sinphi - 54 P->sinph0 * cosphi * coslam); 55 break; 56 case EQUIT: 57 xy.y *= sinphi; 58 break; 59 case N_POLE: 60 coslam = - coslam; 61 case S_POLE: 62 xy.y *= cosphi * coslam; 63 break; 64 } 65 if (P->tilt) { 66 double yt, ba; 67 68 yt = xy.y * P->cg + xy.x * P->sg; 69 ba = 1. / (yt * P->sw * P->h + P->cw); 70 xy.x = (xy.x * P->cg - xy.y * P->sg) * P->cw * ba; 71 xy.y = yt * ba; 72 } 73 return (xy); 74} 75INVERSE(s_inverse); /* spheroid */ 76 double rh, cosz, sinz; 77 78 if (P->tilt) { 79 double bm, bq, yt; 80 81 yt = 1./(P->pn1 - xy.y * P->sw); 82 bm = P->pn1 * xy.x * yt; 83 bq = P->pn1 * xy.y * P->cw * yt; 84 xy.x = bm * P->cg + bq * P->sg; 85 xy.y = bq * P->cg - bm * P->sg; 86 } 87 rh = hypot(xy.x, xy.y); 88 if ((sinz = 1. - rh * rh * P->pfact) < 0.) I_ERROR; 89 sinz = (P->p - sqrt(sinz)) / (P->pn1 / rh + rh / P->pn1); 90 cosz = sqrt(1. - sinz * sinz); 91 if (fabs(rh) <= EPS10) { 92 lp.lam = 0.; 93 lp.phi = P->phi0; 94 } else { 95 switch (P->mode) { 96 case OBLIQ: 97 lp.phi = asin(cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh); 98 xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh; 99 xy.x *= sinz * P->cosph0; 100 break; 101 case EQUIT: 102 lp.phi = asin(xy.y * sinz / rh); 103 xy.y = cosz * rh; 104 xy.x *= sinz; 105 break; 106 case N_POLE: 107 lp.phi = asin(cosz); 108 xy.y = -xy.y; 109 break; 110 case S_POLE: 111 lp.phi = - asin(cosz); 112 break; 113 } 114 lp.lam = atan2(xy.x, xy.y); 115 } 116 return (lp); 117} 118FREEUP; if (P) pj_dalloc(P); } 119 static PJ * 120setup(PJ *P) { 121 if ((P->height = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30); 122 if (fabs(fabs(P->phi0) - HALFPI) < EPS10) 123 P->mode = P->phi0 < 0. ? S_POLE : N_POLE; 124 else if (fabs(P->phi0) < EPS10) 125 P->mode = EQUIT; 126 else { 127 P->mode = OBLIQ; 128 P->sinph0 = sin(P->phi0); 129 P->cosph0 = cos(P->phi0); 130 } 131 P->pn1 = P->height / P->a; /* normalize by radius */ 132 P->p = 1. + P->pn1; 133 P->rp = 1. / P->p; 134 P->h = 1. / P->pn1; 135 P->pfact = (P->p + 1.) * P->h; 136 P->inv = s_inverse; 137 P->fwd = s_forward; 138 P->es = 0.; 139 return P; 140} 141ENTRY0(nsper) 142 P->tilt = 0; 143ENDENTRY(setup(P)) 144ENTRY0(tpers) 145 double omega, gamma; 146 147 omega = pj_param(P->params, "dtilt").f * DEG_TO_RAD; 148 gamma = pj_param(P->params, "dazi").f * DEG_TO_RAD; 149 P->tilt = 1; 150 P->cg = cos(gamma); P->sg = sin(gamma); 151 P->cw = cos(omega); P->sw = sin(omega); 152ENDENTRY(setup(P))