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