/Proj4/PJ_somerc.c
C | 69 lines | 66 code | 3 blank | 0 comment | 4 complexity | 68f924fecb7cc74952401f35c15b0ca8 MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_somerc.c 4.1 95/08/09 GIE REL"; 3#endif 4#define PROJ_PARMS__ \ 5 double K, c, hlf_e, kR, cosp0, sinp0; 6#define PJ_LIB__ 7#include "projects.h" 8PROJ_HEAD(somerc, "Swiss. Obl. Mercator") "\n\tCyl, Ell\n\tFor CH1903"; 9#define EPS 1.e-10 10#define NITER 6 11FORWARD(e_forward); 12 double phip, lamp, phipp, lampp, sp, cp; 13 14 sp = P->e * sin(lp.phi); 15 phip = 2.* atan( exp( P->c * ( 16 log(tan(FORTPI + 0.5 * lp.phi)) - P->hlf_e * log((1. + sp)/(1. - sp))) 17 + P->K)) - HALFPI; 18 lamp = P->c * lp.lam; 19 cp = cos(phip); 20 phipp = aasin(P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp)); 21 lampp = aasin(cp * sin(lamp) / cos(phipp)); 22 xy.x = P->kR * lampp; 23 xy.y = P->kR * log(tan(FORTPI + 0.5 * phipp)); 24 return (xy); 25} 26INVERSE(e_inverse); /* ellipsoid & spheroid */ 27 double phip, lamp, phipp, lampp, cp, esp, con, delp; 28 int i; 29 30 phipp = 2. * (atan(exp(xy.y / P->kR)) - FORTPI); 31 lampp = xy.x / P->kR; 32 cp = cos(phipp); 33 phip = aasin(P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp)); 34 lamp = aasin(cp * sin(lampp) / cos(phip)); 35 con = (P->K - log(tan(FORTPI + 0.5 * phip)))/P->c; 36 for (i = NITER; i ; --i) { 37 esp = P->e * sin(phip); 38 delp = (con + log(tan(FORTPI + 0.5 * phip)) - P->hlf_e * 39 log((1. + esp)/(1. - esp)) ) * 40 (1. - esp * esp) * cos(phip) * P->rone_es; 41 phip -= delp; 42 if (fabs(delp) < EPS) 43 break; 44 } 45 if (i) { 46 lp.phi = phip; 47 lp.lam = lamp / P->c; 48 } else 49 I_ERROR 50 return (lp); 51} 52FREEUP; if (P) pj_dalloc(P); } 53ENTRY0(somerc) 54 double cp, phip0, sp; 55 56 P->hlf_e = 0.5 * P->e; 57 cp = cos(P->phi0); 58 cp *= cp; 59 P->c = sqrt(1 + P->es * cp * cp * P->rone_es); 60 sp = sin(P->phi0); 61 P->cosp0 = cos( phip0 = aasin(P->sinp0 = sp / P->c) ); 62 sp *= P->e; 63 P->K = log(tan(FORTPI + 0.5 * phip0)) - P->c * ( 64 log(tan(FORTPI + 0.5 * P->phi0)) - P->hlf_e * 65 log((1. + sp) / (1. - sp))); 66 P->kR = P->k0 * sqrt(P->one_es) / (1. - sp * sp); 67 P->inv = e_inverse; 68 P->fwd = e_forward; 69ENDENTRY(P)