-
Notifications
You must be signed in to change notification settings - Fork 749
/
factors.cpp
114 lines (90 loc) · 3.06 KB
/
factors.cpp
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
/* projection scale factors */
#include "proj.h"
#include "proj_internal.h"
#include <math.h>
#include <errno.h>
#ifndef DEFAULT_H
#define DEFAULT_H 1e-5 /* radian default for numeric h */
#endif
#define EPS 1.0e-12
int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) {
double cosphi, t, n, r;
int err;
PJ_COORD coo = {{0, 0, 0, 0}};
coo.lp = lp;
/* Failing the 3 initial checks will most likely be due to */
/* earlier errors, so we leave errno alone */
if (nullptr == fac)
return 1;
if (nullptr == P)
return 1;
if (HUGE_VAL == lp.lam)
return 1;
/* But from here, we're ready to make our own mistakes */
err = proj_errno_reset(P);
/* Indicate that all factors are numerical approximations */
fac->code = 0;
/* Check for latitude or longitude overange */
if ((fabs(lp.phi) - M_HALFPI) > EPS) {
proj_log_error(P, _("Invalid latitude"));
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
return 1;
}
if (fabs(lp.lam) > 10.) {
proj_log_error(P, _("Invalid longitude"));
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
return 1;
}
/* Set a reasonable step size for the numerical derivatives */
h = fabs(h);
if (h < EPS)
h = DEFAULT_H;
/* If input latitudes are geocentric, convert to geographic */
if (P->geoc)
lp = pj_geocentric_latitude(P, PJ_INV, coo).lp;
/* If latitude + one step overshoots the pole, move it slightly inside, */
/* so the numerical derivative still exists */
if (fabs(lp.phi) > (M_HALFPI - h))
lp.phi = lp.phi < 0. ? -(M_HALFPI - h) : (M_HALFPI - h);
/* Longitudinal distance from central meridian */
lp.lam -= P->lam0;
if (!P->over)
lp.lam = adjlon(lp.lam);
/* Derivatives */
if (pj_deriv(lp, h, P, &(fac->der))) {
proj_log_error(P, _("Invalid latitude or longitude"));
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
return 1;
}
/* Scale factors */
cosphi = cos(lp.phi);
fac->h = hypot(fac->der.x_p, fac->der.y_p);
fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
if (P->es != 0.0) {
t = sin(lp.phi);
t = 1. - P->es * t * t;
n = sqrt(t);
fac->h *= t * n / P->one_es;
fac->k *= n;
r = t * t / P->one_es;
} else
r = 1.;
/* Convergence */
fac->conv = -atan2(fac->der.x_p, fac->der.y_p);
/* Areal scale factor */
fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) * r /
cosphi;
/* Meridian-parallel angle (theta prime) */
fac->thetap = aasin(P->ctx, fac->s / (fac->h * fac->k));
/* Tissot ellipse axis */
t = fac->k * fac->k + fac->h * fac->h;
fac->a = sqrt(t + 2. * fac->s);
t = t - 2. * fac->s;
t = t > 0 ? sqrt(t) : 0;
fac->b = 0.5 * (fac->a - t);
fac->a = 0.5 * (fac->a + t);
/* Angular distortion */
fac->omega = 2. * aasin(P->ctx, (fac->a - fac->b) / (fac->a + fac->b));
proj_errno_restore(P, err);
return 0;
}