/Proj4/PJ_cass.c
C | 82 lines | 81 code | 1 blank | 0 comment | 5 complexity | ae3dae96ab70eb78789fe671ff78d8dd MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_cass.c 4.1 94/02/15 GIE REL"; 3#endif 4#define PROJ_PARMS__ \ 5 double m0; \ 6 double n; \ 7 double t; \ 8 double a1; \ 9 double c; \ 10 double r; \ 11 double dd; \ 12 double d2; \ 13 double a2; \ 14 double tn; \ 15 double *en; 16#define PJ_LIB__ 17# include "projects.h" 18PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell"; 19# define EPS10 1e-10 20# define C1 .16666666666666666666 21# define C2 .00833333333333333333 22# define C3 .04166666666666666666 23# define C4 .33333333333333333333 24# define C5 .06666666666666666666 25FORWARD(e_forward); /* ellipsoid */ 26 xy.y = pj_mlfn(lp.phi, P->n = sin(lp.phi), P->c = cos(lp.phi), P->en); 27 P->n = 1./sqrt(1. - P->es * P->n * P->n); 28 P->tn = tan(lp.phi); P->t = P->tn * P->tn; 29 P->a1 = lp.lam * P->c; 30 P->c *= P->es * P->c / (1 - P->es); 31 P->a2 = P->a1 * P->a1; 32 xy.x = P->n * P->a1 * (1. - P->a2 * P->t * 33 (C1 - (8. - P->t + 8. * P->c) * P->a2 * C2)); 34 xy.y -= P->m0 - P->n * P->tn * P->a2 * 35 (.5 + (5. - P->t + 6. * P->c) * P->a2 * C3); 36 return (xy); 37} 38FORWARD(s_forward); /* spheroid */ 39 xy.x = asin(cos(lp.phi) * sin(lp.lam)); 40 xy.y = atan2(tan(lp.phi) , cos(lp.lam)) - P->phi0; 41 return (xy); 42} 43INVERSE(e_inverse); /* ellipsoid */ 44 double ph1; 45 46 ph1 = pj_inv_mlfn(P->m0 + xy.y, P->es, P->en); 47 P->tn = tan(ph1); P->t = P->tn * P->tn; 48 P->n = sin(ph1); 49 P->r = 1. / (1. - P->es * P->n * P->n); 50 P->n = sqrt(P->r); 51 P->r *= (1. - P->es) * P->n; 52 P->dd = xy.x / P->n; 53 P->d2 = P->dd * P->dd; 54 lp.phi = ph1 - (P->n * P->tn / P->r) * P->d2 * 55 (.5 - (1. + 3. * P->t) * P->d2 * C3); 56 lp.lam = P->dd * (1. + P->t * P->d2 * 57 (-C4 + (1. + 3. * P->t) * P->d2 * C5)) / cos(ph1); 58 return (lp); 59} 60INVERSE(s_inverse); /* spheroid */ 61 lp.phi = asin(sin(P->dd = xy.y + P->phi0) * cos(xy.x)); 62 lp.lam = atan2(tan(xy.x), cos(P->dd)); 63 return (lp); 64} 65FREEUP; 66 if (P) { 67 if (P->en) 68 pj_dalloc(P->en); 69 pj_dalloc(P); 70 } 71} 72ENTRY1(cass, en) 73 if (P->es) { 74 if (!(P->en = pj_enfn(P->es))) E_ERROR_0; 75 P->m0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en); 76 P->inv = e_inverse; 77 P->fwd = e_forward; 78 } else { 79 P->inv = s_inverse; 80 P->fwd = s_forward; 81 } 82ENDENTRY(P)