/Proj4/PJ_ortho.c
C | 98 lines | 94 code | 4 blank | 0 comment | 20 complexity | ed1e498babb0b308f366671db80941a9 MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_ortho.c 4.1 94/02/15 GIE REL"; 3#endif 4#define PROJ_PARMS__ \ 5 double sinph0; \ 6 double cosph0; \ 7 int mode; 8#define PJ_LIB__ 9#include "projects.h" 10PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph."; 11#define EPS10 1.e-10 12#define N_POLE 0 13#define S_POLE 1 14#define EQUIT 2 15#define OBLIQ 3 16FORWARD(s_forward); /* spheroid */ 17 double coslam, cosphi, sinphi; 18 19 cosphi = cos(lp.phi); 20 coslam = cos(lp.lam); 21 switch (P->mode) { 22 case EQUIT: 23 if (cosphi * coslam < - EPS10) F_ERROR; 24 xy.y = sin(lp.phi); 25 break; 26 case OBLIQ: 27 if (P->sinph0 * (sinphi = sin(lp.phi)) + 28 P->cosph0 * cosphi * coslam < - EPS10) F_ERROR; 29 xy.y = P->cosph0 * sinphi - P->sinph0 * cosphi * coslam; 30 break; 31 case N_POLE: 32 coslam = - coslam; 33 case S_POLE: 34 if (fabs(lp.phi - P->phi0) - EPS10 > HALFPI) F_ERROR; 35 xy.y = cosphi * coslam; 36 break; 37 } 38 xy.x = cosphi * sin(lp.lam); 39 return (xy); 40} 41 42INVERSE(s_inverse); /* spheroid */ 43 double rh, cosc, sinc; 44 45 if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) { 46 if ((sinc - 1.) > EPS10) I_ERROR; 47 sinc = 1.; 48 } 49 cosc = sqrt(1. - sinc * sinc); /* in this range OK */ 50 if (fabs(rh) <= EPS10) { 51 lp.phi = P->phi0; 52 lp.lam = 0.0; 53 } else { 54 switch (P->mode) { 55 case N_POLE: 56 xy.y = -xy.y; 57 lp.phi = acos(sinc); 58 break; 59 case S_POLE: 60 lp.phi = - acos(sinc); 61 break; 62 case EQUIT: 63 lp.phi = xy.y * sinc / rh; 64 xy.x *= sinc; 65 xy.y = cosc * rh; 66 goto sinchk; 67 case OBLIQ: 68 lp.phi = cosc * P->sinph0 + xy.y * sinc * P->cosph0 /rh; 69 xy.y = (cosc - P->sinph0 * lp.phi) * rh; 70 xy.x *= sinc * P->cosph0; 71 sinchk: 72 if (fabs(lp.phi) >= 1.) 73 lp.phi = lp.phi < 0. ? -HALFPI : HALFPI; 74 else 75 lp.phi = asin(lp.phi); 76 break; 77 } 78 lp.lam = (xy.y == 0. && (P->mode == OBLIQ || P->mode == EQUIT)) 79 ? (xy.x == 0. ? 0. : xy.x < 0. ? -HALFPI : HALFPI) 80 : atan2(xy.x, xy.y); 81 } 82 return (lp); 83} 84 85FREEUP; if (P) pj_dalloc(P); } 86ENTRY0(ortho) 87 if (fabs(fabs(P->phi0) - HALFPI) <= EPS10) 88 P->mode = P->phi0 < 0. ? S_POLE : N_POLE; 89 else if (fabs(P->phi0) > EPS10) { 90 P->mode = OBLIQ; 91 P->sinph0 = sin(P->phi0); 92 P->cosph0 = cos(P->phi0); 93 } else 94 P->mode = EQUIT; 95 P->inv = s_inverse; 96 P->fwd = s_forward; 97 P->es = 0.; 98ENDENTRY(P)