/Proj4/PJ_lcca.c
C | 76 lines | 64 code | 6 blank | 6 comment | 9 complexity | 343d16310027e7fdf17c1c91c36adeb3 MD5 | raw file
1static const char RCS_ID[] = "$Id: PJ_lcca.c,v 1.1 2003/03/04 02:59:41 warmerda Exp $"; 2/* PROJ.4 Cartographic Projection System -- Revision Log: 3**$Log: PJ_lcca.c,v $ 4**Revision 1.1 2003/03/04 02:59:41 warmerda 5**New 6** 7*/ 8#define MAX_ITER 10 9#define DEL_TOL 1e-12 10#define PROJ_PARMS__ \ 11 double *en; \ 12 double r0, l, M0; \ 13 double C; 14#define PJ_LIB__ 15#include "projects.h" 16 17PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative") 18 "\n\tConic, Sph&Ell\n\tlat_0="; 19 20 static double /* func to compute dr */ 21fS(double S, double C) { 22 return(S * ( 1. + S * S * C)); 23} 24 static double /* deriv of fs */ 25fSp(double S, double C) { 26 return(1. + 3.* S * S * C); 27} 28FORWARD(e_forward); /* ellipsoid */ 29 double S, S3, r, dr; 30 31 S = pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), P->en) - P->M0; 32 dr = fS(S, P->C); 33 r = P->r0 - dr; 34 xy.x = P->k0 * (r * sin( lp.lam *= P->l ) ); 35 xy.y = P->k0 * (P->r0 - r * cos(lp.lam) ); 36 return (xy); 37} 38INVERSE(e_inverse); /* ellipsoid & spheroid */ 39 double theta, dr, S, dif; 40 int i; 41 42 xy.x /= P->k0; 43 xy.y /= P->k0; 44 theta = atan2(xy.x , P->r0 - xy.y); 45 dr = xy.y - xy.x * tan(0.5 * theta); 46 lp.lam = theta / P->l; 47 S = dr; 48 for (i = MAX_ITER; i ; --i) { 49 S -= (dif = (fS(S, P->C) - dr) / fSp(S, P->C)); 50 if (fabs(dif) < DEL_TOL) break; 51 } 52 if (!i) I_ERROR 53 lp.phi = pj_inv_mlfn(S + P->M0, P->es, P->en); 54 return (lp); 55} 56FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } } 57ENTRY0(lcca) 58 double s2p0, N0, R0, tan0, tan20; 59 60 if (!(P->en = pj_enfn(P->es))) E_ERROR_0; 61 if (!pj_param(P->params, "tlat_0").i) E_ERROR(50); 62 if (P->phi0 == 0.) E_ERROR(51); 63 P->l = sin(P->phi0); 64 P->M0 = pj_mlfn(P->phi0, P->l, cos(P->phi0), P->en); 65 s2p0 = P->l * P->l; 66 R0 = 1. / (1. - P->es * s2p0); 67 N0 = sqrt(R0); 68 R0 *= P->one_es * N0; 69 tan0 = tan(P->phi0); 70 tan20 = tan0 * tan0; 71 P->r0 = N0 / tan0; 72 P->C = 1. / (6. * R0 * N0); 73 P->inv = e_inverse; 74 P->fwd = e_forward; 75ENDENTRY(P) 76