Skip to content
Newer
Older
100644 108 lines (100 sloc) 2.93 KB
1a37e1c @josephg Readded source
josephg authored
1 /*
2 ** libproj -- library of cartographic projections
3 **
4 ** Copyright (c) 2003 Gerald I. Evenden
5 */
6 static const char
7 LIBPROJ_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
33 struct 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
42 srat(double esinp, double exp) {
43 return(pow((1.-esinp)/(1.+esinp), exp));
44 }
45
46 void *
47 pj_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
67 pj_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
77 pj_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
Something went wrong with that request. Please try again.