/Proj4/PJ_mbt_fps.c
C | 41 lines | 39 code | 2 blank | 0 comment | 3 complexity | c85925fedd16e409b2cb925dbee37078 MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_mbt_fps.c 4.1 94/02/15 GIE REL"; 3#endif 4#define PJ_LIB__ 5#include "projects.h" 6PROJ_HEAD(mbt_fps, "McBryde-Thomas Flat-Pole Sine (No. 2)") "\n\tCyl., Sph."; 7#define MAX_ITER 10 8#define LOOP_TOL 1e-7 9#define C1 0.45503 10#define C2 1.36509 11#define C3 1.41546 12#define C_x 0.22248 13#define C_y 1.44492 14#define C1_2 0.33333333333333333333333333 15FORWARD(s_forward); /* spheroid */ 16 double k, V, t; 17 int i; 18 19 k = C3 * sin(lp.phi); 20 for (i = MAX_ITER; i ; --i) { 21 t = lp.phi / C2; 22 lp.phi -= V = (C1 * sin(t) + sin(lp.phi) - k) / 23 (C1_2 * cos(t) + cos(lp.phi)); 24 if (fabs(V) < LOOP_TOL) 25 break; 26 } 27 t = lp.phi / C2; 28 xy.x = C_x * lp.lam * (1. + 3. * cos(lp.phi)/cos(t) ); 29 xy.y = C_y * sin(t); 30 return (xy); 31} 32INVERSE(s_inverse); /* spheroid */ 33 double t, s; 34 35 lp.phi = C2 * (t = aasin(xy.y / C_y)); 36 lp.lam = xy.x / (C_x * (1. + 3. * cos(lp.phi)/cos(t))); 37 lp.phi = aasin((C1 * sin(t) + sin(lp.phi)) / C3); 38 return (lp); 39} 40FREEUP; if (P) pj_dalloc(P); } 41ENTRY0(mbt_fps) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)