svn-gvsig-desktop / trunk / libraries / libjni-proj4 / src / pj_factors.c @ 13395
History | View | Annotate | Download (2.4 KB)
1 |
/* projection scale factors */
|
---|---|
2 |
#ifndef lint
|
3 |
static const char SCCSID[]="@(#)pj_factors.c 4.9 94/03/17 GIE REL"; |
4 |
#endif
|
5 |
#define PJ_LIB__
|
6 |
#include <projects.h> |
7 |
#include <errno.h> |
8 |
#ifndef DEFAULT_H
|
9 |
#define DEFAULT_H 1e-5 /* radian default for numeric h */ |
10 |
#endif
|
11 |
#define EPS 1.0e-12 |
12 |
int
|
13 |
pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) { |
14 |
struct DERIVS der;
|
15 |
double cosphi, t, n, r;
|
16 |
|
17 |
/* check for forward and latitude or longitude overange */
|
18 |
if ((t = fabs(lp.phi)-HALFPI) > EPS || fabs(lp.lam) > 10.) { |
19 |
pj_errno = -14;
|
20 |
return 1; |
21 |
} else { /* proceed */ |
22 |
errno = pj_errno = 0;
|
23 |
if (fabs(t) <= EPS) /* adjust to pi/2 */ |
24 |
lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
|
25 |
else if (P->geoc) |
26 |
lp.phi = atan(P->rone_es * tan(lp.phi)); |
27 |
lp.lam -= P->lam0; /* compute del lp.lam */
|
28 |
if (!P->over)
|
29 |
lp.lam = adjlon(lp.lam); /* adjust del longitude */
|
30 |
if (h <= 0.) |
31 |
h = DEFAULT_H; |
32 |
if (P->spc) /* get what projection analytic values */ |
33 |
P->spc(lp, P, fac); |
34 |
if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
|
35 |
(IS_ANAL_XL_YL+IS_ANAL_XP_YP)) && |
36 |
pj_deriv(lp, h, P, &der)) |
37 |
return 1; |
38 |
if (!(fac->code & IS_ANAL_XL_YL)) {
|
39 |
fac->der.x_l = der.x_l; |
40 |
fac->der.y_l = der.y_l; |
41 |
} |
42 |
if (!(fac->code & IS_ANAL_XP_YP)) {
|
43 |
fac->der.x_p = der.x_p; |
44 |
fac->der.y_p = der.y_p; |
45 |
} |
46 |
cosphi = cos(lp.phi); |
47 |
if (!(fac->code & IS_ANAL_HK)) {
|
48 |
fac->h = hypot(fac->der.x_p, fac->der.y_p); |
49 |
fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi; |
50 |
if (P->es) {
|
51 |
t = sin(lp.phi); |
52 |
t = 1. - P->es * t * t;
|
53 |
n = sqrt(t); |
54 |
fac->h *= t * n / P->one_es; |
55 |
fac->k *= n; |
56 |
r = t * t / P->one_es; |
57 |
} else
|
58 |
r = 1.;
|
59 |
} else if (P->es) { |
60 |
r = sin(lp.phi); |
61 |
r = 1. - P->es * r * r;
|
62 |
r = r * r / P->one_es; |
63 |
} else
|
64 |
r = 1.;
|
65 |
/* convergence */
|
66 |
if (!(fac->code & IS_ANAL_CONV)) {
|
67 |
fac->conv = - atan2(fac->der.y_l, fac->der.x_l); |
68 |
if (fac->code & IS_ANAL_XL_YL)
|
69 |
fac->code |= IS_ANAL_CONV; |
70 |
} |
71 |
/* areal scale factor */
|
72 |
fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) * |
73 |
r / cosphi; |
74 |
/* meridian-parallel angle theta prime */
|
75 |
fac->thetap = aasin(fac->s / (fac->h * fac->k)); |
76 |
/* Tissot ellips axis */
|
77 |
t = fac->k * fac->k + fac->h * fac->h; |
78 |
fac->a = sqrt(t + 2. * fac->s);
|
79 |
t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t); |
80 |
fac->b = 0.5 * (fac->a - t); |
81 |
fac->a = 0.5 * (fac->a + t); |
82 |
/* omega */
|
83 |
fac->omega = 2. * aasin((fac->a - fac->b)/(fac->a + fac->b));
|
84 |
} |
85 |
return 0; |
86 |
} |