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
128
129
130
131
|
#ifndef lint
static const char SCCSID[]="@(#)PJ_aea.c 4.2 94/03/18 GIE REL";
#endif
# define EPS10 1.e-10
# define TOL7 1.e-7
#define PROJ_PARMS__ \
double ec; \
double n; \
double c; \
double dd; \
double n2; \
double rho0; \
double rho; \
double phi1; \
double phi2; \
double *en; \
int ellips;
#define PJ_LIB__
# include <projects.h>
PROJ_HEAD(aea, "Albers Equal Area")
"\n\tConic Sph&Ell\n\tlat_1= lat_2=";
PROJ_HEAD(leac, "Lambert Equal Area Conic")
"\n\tConic, Sph&Ell\n\tlat_1= south";
/* determine latitude angle phi-1 */
# define N_ITER 15
# define EPSILON 1.0e-7
# define TOL 1.0e-10
static double
phi1_(double qs, double Te, double Tone_es) {
int i;
double Phi, sinpi, cospi, con, com, dphi;
Phi = asin (.5 * qs);
if (Te < EPSILON)
return( Phi );
i = N_ITER;
do {
sinpi = sin (Phi);
cospi = cos (Phi);
con = Te * sinpi;
com = 1. - con * con;
dphi = .5 * com * com / cospi * (qs / Tone_es -
sinpi / com + .5 / Te * log ((1. - con) /
(1. + con)));
Phi += dphi;
} while (fabs(dphi) > TOL && --i);
return( i ? Phi : HUGE_VAL );
}
FORWARD(e_forward); /* ellipsoid & spheroid */
if ((P->rho = P->c - (P->ellips ? P->n * pj_qsfn(sin(lp.phi),
P->e, P->one_es) : P->n2 * sin(lp.phi))) < 0.) F_ERROR
P->rho = P->dd * sqrt(P->rho);
xy.x = P->rho * sin( lp.lam *= P->n );
xy.y = P->rho0 - P->rho * cos(lp.lam);
return (xy);
}
INVERSE(e_inverse) /* ellipsoid & spheroid */;
if( (P->rho = hypot(xy.x, xy.y = P->rho0 - xy.y)) != 0.0 ) {
if (P->n < 0.) {
P->rho = -P->rho;
xy.x = -xy.x;
xy.y = -xy.y;
}
lp.phi = P->rho / P->dd;
if (P->ellips) {
lp.phi = (P->c - lp.phi * lp.phi) / P->n;
if (fabs(P->ec - fabs(lp.phi)) > TOL7) {
if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
I_ERROR
} else
lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
} else if (fabs(lp.phi = (P->c - lp.phi * lp.phi) / P->n2) <= 1.)
lp.phi = asin(lp.phi);
else
lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
lp.lam = atan2(xy.x, xy.y) / P->n;
} else {
lp.lam = 0.;
lp.phi = P->n > 0. ? HALFPI : - HALFPI;
}
return (lp);
}
FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
static PJ *
setup(PJ *P) {
double cosphi, sinphi;
int secant;
if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21);
P->n = sinphi = sin(P->phi1);
cosphi = cos(P->phi1);
secant = fabs(P->phi1 - P->phi2) >= EPS10;
if( (P->ellips = (P->es > 0.))) {
double ml1, m1;
if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
m1 = pj_msfn(sinphi, cosphi, P->es);
ml1 = pj_qsfn(sinphi, P->e, P->one_es);
if (secant) { /* secant cone */
double ml2, m2;
sinphi = sin(P->phi2);
cosphi = cos(P->phi2);
m2 = pj_msfn(sinphi, cosphi, P->es);
ml2 = pj_qsfn(sinphi, P->e, P->one_es);
P->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
}
P->ec = 1. - .5 * P->one_es * log((1. - P->e) /
(1. + P->e)) / P->e;
P->c = m1 * m1 + P->n * ml1;
P->dd = 1. / P->n;
P->rho0 = P->dd * sqrt(P->c - P->n * pj_qsfn(sin(P->phi0),
P->e, P->one_es));
} else {
if (secant) P->n = .5 * (P->n + sin(P->phi2));
P->n2 = P->n + P->n;
P->c = cosphi * cosphi + P->n2 * sinphi;
P->dd = 1. / P->n;
P->rho0 = P->dd * sqrt(P->c - P->n2 * sin(P->phi0));
}
P->inv = e_inverse; P->fwd = e_forward;
return P;
}
ENTRY0(aea)
P->phi1 = pj_param(P->params, "rlat_1").f;
P->phi2 = pj_param(P->params, "rlat_2").f;
ENDENTRY(setup(P))
ENTRY0(leac)
P->phi2 = pj_param(P->params, "rlat_1").f;
P->phi1 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI;
ENDENTRY(setup(P))
|
