/Proj4/PJ_ocea.c
C | 71 lines | 65 code | 3 blank | 3 comment | 4 complexity | 9845d1090fd8992c0d2fa94d056c08dc MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_ocea.c 4.1 94/02/15 GIE REL"; 3#endif 4#define PROJ_PARMS__ \ 5 double rok; \ 6 double rtk; \ 7 double sinphi; \ 8 double cosphi; \ 9 double singam; \ 10 double cosgam; 11#define PJ_LIB__ 12#include "projects.h" 13PROJ_HEAD(ocea, "Oblique Cylindrical Equal Area") "\n\tCyl, Sph" 14 "lonc= alpha= or\n\tlat_1= lat_2= lon_1= lon_2="; 15FORWARD(s_forward); /* spheroid */ 16 double t; 17 18 xy.y = sin(lp.lam); 19/* 20 xy.x = atan2((tan(lp.phi) * P->cosphi + P->sinphi * xy.y) , cos(lp.lam)); 21*/ 22 t = cos(lp.lam); 23 xy.x = atan((tan(lp.phi) * P->cosphi + P->sinphi * xy.y) / t); 24 if (t < 0.) 25 xy.x += PI; 26 xy.x *= P->rtk; 27 xy.y = P->rok * (P->sinphi * sin(lp.phi) - P->cosphi * cos(lp.phi) * xy.y); 28 return (xy); 29} 30INVERSE(s_inverse); /* spheroid */ 31 double t, s; 32 33 xy.y /= P->rok; 34 xy.x /= P->rtk; 35 t = sqrt(1. - xy.y * xy.y); 36 lp.phi = asin(xy.y * P->sinphi + t * P->cosphi * (s = sin(xy.x))); 37 lp.lam = atan2(t * P->sinphi * s - xy.y * P->cosphi, 38 t * cos(xy.x)); 39 return (lp); 40} 41FREEUP; if (P) pj_dalloc(P); } 42ENTRY0(ocea) 43 double phi_0=0.0, phi_1, phi_2, lam_1, lam_2, lonz, alpha; 44 45 P->rok = P->a / P->k0; 46 P->rtk = P->a * P->k0; 47 if ( pj_param(P->params, "talpha").i) { 48 alpha = pj_param(P->params, "ralpha").f; 49 lonz = pj_param(P->params, "rlonc").f; 50 P->singam = atan(-cos(alpha)/(-sin(phi_0) * sin(alpha))) + lonz; 51 P->sinphi = asin(cos(phi_0) * sin(alpha)); 52 } else { 53 phi_1 = pj_param(P->params, "rlat_1").f; 54 phi_2 = pj_param(P->params, "rlat_2").f; 55 lam_1 = pj_param(P->params, "rlon_1").f; 56 lam_2 = pj_param(P->params, "rlon_2").f; 57 P->singam = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) - 58 sin(phi_1) * cos(phi_2) * cos(lam_2), 59 sin(phi_1) * cos(phi_2) * sin(lam_2) - 60 cos(phi_1) * sin(phi_2) * sin(lam_1) ); 61 P->sinphi = atan(-cos(P->singam - lam_1) / tan(phi_1)); 62 } 63 P->lam0 = P->singam + HALFPI; 64 P->cosphi = cos(P->sinphi); 65 P->sinphi = sin(P->sinphi); 66 P->cosgam = cos(P->singam); 67 P->singam = sin(P->singam); 68 P->inv = s_inverse; 69 P->fwd = s_forward; 70 P->es = 0.; 71ENDENTRY(P)