svn-gvsig-desktop / tags / v1_9_Build_1241 / libraries / libjni-proj4 / src / PJ_lcc.c @ 34066
History | View | Annotate | Download (2.93 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_lcc.c 4.2 94/03/18 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double phi1; \
|
6 |
double phi2; \
|
7 |
double n; \
|
8 |
double rho; \
|
9 |
double rho0; \
|
10 |
double c; \
|
11 |
int ellips;
|
12 |
#define PJ_LIB__
|
13 |
#include <projects.h> |
14 |
PROJ_HEAD(lcc, "Lambert Conformal Conic")
|
15 |
"\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0";
|
16 |
# define EPS10 1.e-10 |
17 |
FORWARD(e_forward); /* ellipsoid & spheroid */
|
18 |
if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
|
19 |
if ((lp.phi * P->n) <= 0.) F_ERROR; |
20 |
P->rho = 0.;
|
21 |
} |
22 |
else
|
23 |
P->rho = P->c * (P->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi), |
24 |
P->e), P->n) : pow(tan(FORTPI + .5 * lp.phi), -P->n));
|
25 |
xy.x = P->k0 * (P->rho * sin( lp.lam *= P->n ) ); |
26 |
xy.y = P->k0 * (P->rho0 - P->rho * cos(lp.lam) ); |
27 |
return (xy);
|
28 |
} |
29 |
INVERSE(e_inverse); /* ellipsoid & spheroid */
|
30 |
xy.x /= P->k0; |
31 |
xy.y /= P->k0; |
32 |
if( (P->rho = hypot(xy.x, xy.y = P->rho0 - xy.y)) != 0.0) { |
33 |
if (P->n < 0.) { |
34 |
P->rho = -P->rho; |
35 |
xy.x = -xy.x; |
36 |
xy.y = -xy.y; |
37 |
} |
38 |
if (P->ellips) {
|
39 |
if ((lp.phi = pj_phi2(pow(P->rho / P->c, 1./P->n), P->e)) |
40 |
== HUGE_VAL) |
41 |
I_ERROR; |
42 |
} else
|
43 |
lp.phi = 2. * atan(pow(P->c / P->rho, 1./P->n)) - HALFPI; |
44 |
lp.lam = atan2(xy.x, xy.y) / P->n; |
45 |
} else {
|
46 |
lp.lam = 0.;
|
47 |
lp.phi = P->n > 0. ? HALFPI : - HALFPI;
|
48 |
} |
49 |
return (lp);
|
50 |
} |
51 |
SPECIAL(fac) { |
52 |
if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
|
53 |
if ((lp.phi * P->n) <= 0.) return; |
54 |
P->rho = 0.;
|
55 |
} else
|
56 |
P->rho = P->c * (P->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi), |
57 |
P->e), P->n) : pow(tan(FORTPI + .5 * lp.phi), -P->n));
|
58 |
fac->code |= IS_ANAL_HK + IS_ANAL_CONV; |
59 |
fac->k = fac->h = P->k0 * P->n * P->rho / |
60 |
pj_msfn(sin(lp.phi), cos(lp.phi), P->es); |
61 |
fac->conv = - P->n * lp.lam; |
62 |
} |
63 |
FREEUP; if (P) pj_dalloc(P); }
|
64 |
ENTRY0(lcc) |
65 |
double cosphi, sinphi;
|
66 |
int secant;
|
67 |
|
68 |
P->phi1 = pj_param(P->params, "rlat_1").f;
|
69 |
if (pj_param(P->params, "tlat_2").i) |
70 |
P->phi2 = pj_param(P->params, "rlat_2").f;
|
71 |
else {
|
72 |
P->phi2 = P->phi1; |
73 |
if (!pj_param(P->params, "tlat_0").i) |
74 |
P->phi0 = P->phi1; |
75 |
} |
76 |
if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21); |
77 |
P->n = sinphi = sin(P->phi1); |
78 |
cosphi = cos(P->phi1); |
79 |
secant = fabs(P->phi1 - P->phi2) >= EPS10; |
80 |
if( (P->ellips = (P->es != 0.)) ) { |
81 |
double ml1, m1;
|
82 |
|
83 |
P->e = sqrt(P->es); |
84 |
m1 = pj_msfn(sinphi, cosphi, P->es); |
85 |
ml1 = pj_tsfn(P->phi1, sinphi, P->e); |
86 |
if (secant) { /* secant cone */ |
87 |
P->n = log(m1 / |
88 |
pj_msfn(sinphi = sin(P->phi2), cos(P->phi2), P->es)); |
89 |
P->n /= log(ml1 / pj_tsfn(P->phi2, sinphi, P->e)); |
90 |
} |
91 |
P->c = (P->rho0 = m1 * pow(ml1, -P->n) / P->n); |
92 |
P->rho0 *= (fabs(fabs(P->phi0) - HALFPI) < EPS10) ? 0. :
|
93 |
pow(pj_tsfn(P->phi0, sin(P->phi0), P->e), P->n); |
94 |
} else {
|
95 |
if (secant)
|
96 |
P->n = log(cosphi / cos(P->phi2)) / |
97 |
log(tan(FORTPI + .5 * P->phi2) /
|
98 |
tan(FORTPI + .5 * P->phi1));
|
99 |
P->c = cosphi * pow(tan(FORTPI + .5 * P->phi1), P->n) / P->n;
|
100 |
P->rho0 = (fabs(fabs(P->phi0) - HALFPI) < EPS10) ? 0. :
|
101 |
P->c * pow(tan(FORTPI + .5 * P->phi0), -P->n);
|
102 |
} |
103 |
P->inv = e_inverse; |
104 |
P->fwd = e_forward; |
105 |
P->spc = fac; |
106 |
ENDENTRY(P) |