/Proj4/pj_gauss.c
C | 107 lines | 62 code | 7 blank | 38 comment | 5 complexity | 6664bdff4c492a935f3b4b40da7ef1b1 MD5 | raw file
1/* 2** libproj -- library of cartographic projections 3** 4** Copyright (c) 2003 Gerald I. Evenden 5*/ 6static const char 7LIBPROJ_ID[] = "$Id: pj_gauss.c,v 1.1 2004/10/20 17:04:00 fwarmerdam Exp $"; 8/* 9** Permission is hereby granted, free of charge, to any person obtaining 10** a copy of this software and associated documentation files (the 11** "Software"), to deal in the Software without restriction, including 12** without limitation the rights to use, copy, modify, merge, publish, 13** distribute, sublicense, and/or sell copies of the Software, and to 14** permit persons to whom the Software is furnished to do so, subject to 15** the following conditions: 16** 17** The above copyright notice and this permission notice shall be 18** included in all copies or substantial portions of the Software. 19** 20** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 21** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 22** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 23** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY 24** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, 25** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 26** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 27*/ 28#define PJ_LIB__ 29#include "projects.h" 30 31#define MAX_ITER 20 32 33struct GAUSS { 34 double C; 35 double K; 36 double e; 37 double ratexp; 38}; 39#define EN ((struct GAUSS *)en) 40#define DEL_TOL 1e-14 41 static double 42srat(double esinp, double exp) { 43 return(pow((1.-esinp)/(1.+esinp), exp)); 44} 45 46 void * 47pj_gauss_ini(double e, double phi0, double *chi, double *rc) { 48 double sphi, cphi, es; 49 struct GAUSS *en; 50 51 if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL) 52 return (NULL); 53 es = e * e; 54 EN->e = e; 55 sphi = sin(phi0); 56 cphi = cos(phi0); cphi *= cphi; 57 *rc = sqrt(1. - es) / (1. - es * sphi * sphi); 58 EN->C = sqrt(1. + es * cphi * cphi / (1. - es)); 59 *chi = asin(sphi / EN->C); 60 EN->ratexp = 0.5 * EN->C * e; 61 EN->K = tan(.5 * *chi + FORTPI) / ( 62 pow(tan(.5 * phi0 + FORTPI), EN->C) * 63 srat(EN->e * sphi, EN->ratexp) ); 64 return ((void *)en); 65} 66 LP 67pj_gauss(LP elp, const void *en) { 68 LP slp; 69 70 slp.phi = 2. * atan( EN->K * 71 pow(tan(.5 * elp.phi + FORTPI), EN->C) * 72 srat(EN->e * sin(elp.phi), EN->ratexp) ) - HALFPI; 73 slp.lam = EN->C * (elp.lam); 74 return(slp); 75} 76 LP 77pj_inv_gauss(LP slp, const void *en) { 78 LP elp; 79 double num; 80 int i; 81 82 elp.lam = slp.lam / EN->C; 83 num = pow(tan(.5 * slp.phi + FORTPI)/EN->K, 1./EN->C); 84 for (i = MAX_ITER; i; --i) { 85 elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e)) 86 - HALFPI; 87 if (fabs(elp.phi - slp.phi) < DEL_TOL) break; 88 slp.phi = elp.phi; 89 } 90 /* convergence failed */ 91 if (!i) 92 pj_errno = -17; 93 return (elp); 94} 95/* Revision Log: 96** $Log: pj_gauss.c,v $ 97** Revision 1.1 2004/10/20 17:04:00 fwarmerdam 98** New 99** 100** Revision 2.2 2004/03/15 16:07:42 gie 101** removed es from init structure 102** 103** Revision 2.1 2003/03/28 01:44:30 gie 104** Initial 105** 106*/ 107