1 /*
2 ** libproj -- library of cartographic projections
3 **
4 ** Copyright (c) 2003, 2006 Gerald I. Evenden
5 */
6 /*
7 ** Permission is hereby granted, free of charge, to any person obtaining
8 ** a copy of this software and associated documentation files (the
9 ** "Software"), to deal in the Software without restriction, including
10 ** without limitation the rights to use, copy, modify, merge, publish,
11 ** distribute, sublicense, and/or sell copies of the Software, and to
12 ** permit persons to whom the Software is furnished to do so, subject to
13 ** the following conditions:
14 **
15 ** The above copyright notice and this permission notice shall be
16 ** included in all copies or substantial portions of the Software.
17 **
18 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
19 ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20 ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
21 ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
22 ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25 */
26 /* Computes distance from equator along the meridian to latitude phi
27 ** and inverse on unit ellipsoid.
28 ** Precision commensurate with double precision.
29 */
30 #define PROJ_LIB__
31 #include <projects.h>
32 #define MAX_ITER 20
33 #define TOL 1e-14
34
35 struct MDIST {
36 int nb;
37 double es;
38 double E;
39 double b[1];
40 };
41 #define B ((struct MDIST *)b)
42 void *
proj_mdist_ini(double es)43 proj_mdist_ini(double es) {
44 double numf, numfi, twon1, denf, denfi, ens, T, twon;
45 double den, El, Es;
46 double E[MAX_ITER];
47 struct MDIST *b;
48 int i, j;
49
50 /* generate E(e^2) and its terms E[] */
51 ens = es;
52 numf = twon1 = denfi = 1.;
53 denf = 1.;
54 twon = 4.;
55 Es = El = E[0] = 1.;
56 for (i = 1; i < MAX_ITER ; ++i) {
57 numf *= (twon1 * twon1);
58 den = twon * denf * denf * twon1;
59 T = numf/den;
60 Es -= (E[i] = T * ens);
61 ens *= es;
62 twon *= 4.;
63 denf *= ++denfi;
64 twon1 += 2.;
65 if (Es == El) /* jump out if no change */
66 break;
67 El = Es;
68 }
69 if ((b = (struct MDIST *)malloc(sizeof(struct MDIST)+
70 (i*sizeof(double)))) == NULL)
71 return(NULL);
72 b->nb = i - 1;
73 b->es = es;
74 b->E = Es;
75 /* generate b_n coefficients--note: collapse with prefix ratios */
76 b->b[0] = Es = 1. - Es;
77 numf = denf = 1.;
78 numfi = 2.;
79 denfi = 3.;
80 for (j = 1; j < i; ++j) {
81 Es -= E[j];
82 numf *= numfi;
83 denf *= denfi;
84 b->b[j] = Es * numf / denf;
85 numfi += 2.;
86 denfi += 2.;
87 }
88 return (b);
89 }
90 double
proj_mdist(double phi,double sphi,double cphi,const void * b)91 proj_mdist(double phi, double sphi, double cphi, const void *b) {
92 double sc, sum, sphi2, D;
93 int i;
94
95 sc = sphi * cphi;
96 sphi2 = sphi * sphi;
97 D = phi * B->E - B->es * sc / sqrt(1. - B->es * sphi2);
98 sum = B->b[i = B->nb];
99 while (i) sum = B->b[--i] + sphi2 * sum;
100 return(D + sc * sum);
101 }
102 double
proj_inv_mdist(projCtx ctx,double dist,const void * b)103 proj_inv_mdist(projCtx ctx, double dist, const void *b) {
104 double s, t, phi, k;
105 int i;
106
107 k = 1./(1.- B->es);
108 i = MAX_ITER;
109 phi = dist;
110 while ( i-- ) {
111 s = sin(phi);
112 t = 1. - B->es * s * s;
113 phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
114 (t * sqrt(t)) * k;
115 if (fabs(t) < TOL) /* that is no change */
116 return phi;
117 }
118 /* convergence failed */
119 pj_ctx_set_errno(ctx, -17);
120 return phi;
121 }
122