svn-gvsig-desktop / tags / v1_10_0_Build_1257 / libraries / libjni-proj4 / src / PJ_tmerc.c @ 42039
History | View | Annotate | Download (3.87 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_tmerc.c 4.2 94/06/02 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double esp; \
|
6 |
double ml0; \
|
7 |
double *en;
|
8 |
#define PJ_LIB__
|
9 |
#include <projects.h> |
10 |
PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell"; |
11 |
PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)")
|
12 |
"\n\tCyl, Sph\n\tzone= south";
|
13 |
#define EPS10 1.e-10 |
14 |
#define aks0 P->esp
|
15 |
#define aks5 P->ml0
|
16 |
#define FC1 1. |
17 |
#define FC2 .5 |
18 |
#define FC3 .16666666666666666666 |
19 |
#define FC4 .08333333333333333333 |
20 |
#define FC5 .05 |
21 |
#define FC6 .03333333333333333333 |
22 |
#define FC7 .02380952380952380952 |
23 |
#define FC8 .01785714285714285714 |
24 |
FORWARD(e_forward); /* ellipse */
|
25 |
double al, als, n, cosphi, sinphi, t;
|
26 |
|
27 |
sinphi = sin(lp.phi); cosphi = cos(lp.phi); |
28 |
t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; |
29 |
t *= t; |
30 |
al = cosphi * lp.lam; |
31 |
als = al * al; |
32 |
al /= sqrt(1. - P->es * sinphi * sinphi);
|
33 |
n = P->esp * cosphi * cosphi; |
34 |
xy.x = P->k0 * al * (FC1 + |
35 |
FC3 * als * (1. - t + n +
|
36 |
FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t) |
37 |
+ FC7 * als * (61. + t * ( t * (179. - t) - 479. ) ) |
38 |
))); |
39 |
xy.y = P->k0 * (pj_mlfn(lp.phi, sinphi, cosphi, P->en) - P->ml0 + |
40 |
sinphi * al * lp.lam * FC2 * ( 1. +
|
41 |
FC4 * als * (5. - t + n * (9. + 4. * n) + |
42 |
FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t) |
43 |
+ FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) ) |
44 |
)))); |
45 |
return (xy);
|
46 |
} |
47 |
FORWARD(s_forward); /* sphere */
|
48 |
double b, cosphi;
|
49 |
|
50 |
b = (cosphi = cos(lp.phi)) * sin(lp.lam); |
51 |
if (fabs(fabs(b) - 1.) <= EPS10) F_ERROR; |
52 |
xy.x = aks5 * log((1. + b) / (1. - b)); |
53 |
if ((b = fabs( xy.y = cosphi * cos(lp.lam) / sqrt(1. - b * b) )) >= 1.) { |
54 |
if ((b - 1.) > EPS10) F_ERROR |
55 |
else xy.y = 0.; |
56 |
} else
|
57 |
xy.y = acos(xy.y); |
58 |
if (lp.phi < 0.) xy.y = -xy.y; |
59 |
xy.y = aks0 * (xy.y - P->phi0); |
60 |
return (xy);
|
61 |
} |
62 |
INVERSE(e_inverse); /* ellipsoid */
|
63 |
double n, con, cosphi, d, ds, sinphi, t;
|
64 |
|
65 |
lp.phi = pj_inv_mlfn(P->ml0 + xy.y / P->k0, P->es, P->en); |
66 |
if (fabs(lp.phi) >= HALFPI) {
|
67 |
lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
|
68 |
lp.lam = 0.;
|
69 |
} else {
|
70 |
sinphi = sin(lp.phi); |
71 |
cosphi = cos(lp.phi); |
72 |
t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; |
73 |
n = P->esp * cosphi * cosphi; |
74 |
d = xy.x * sqrt(con = 1. - P->es * sinphi * sinphi) / P->k0;
|
75 |
con *= t; |
76 |
t *= t; |
77 |
ds = d * d; |
78 |
lp.phi -= (con * ds / (1.-P->es)) * FC2 * (1. - |
79 |
ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) - |
80 |
ds * FC6 * (61. + t * (90. - 252. * n + |
81 |
45. * t) + 46. * n |
82 |
- ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) |
83 |
))); |
84 |
lp.lam = d*(FC1 - |
85 |
ds*FC3*( 1. + 2.*t + n - |
86 |
ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n |
87 |
- ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) ) |
88 |
))) / cosphi; |
89 |
} |
90 |
return (lp);
|
91 |
} |
92 |
INVERSE(s_inverse); /* sphere */
|
93 |
double h, g;
|
94 |
|
95 |
h = exp(xy.x / aks0); |
96 |
g = .5 * (h - 1. / h); |
97 |
h = cos(P->phi0 + xy.y / aks0); |
98 |
lp.phi = asin(sqrt((1. - h * h) / (1. + g * g))); |
99 |
if (xy.y < 0.) lp.phi = -lp.phi; |
100 |
lp.lam = (g || h) ? atan2(g, h) : 0.;
|
101 |
return (lp);
|
102 |
} |
103 |
FREEUP; |
104 |
if (P) {
|
105 |
if (P->en)
|
106 |
pj_dalloc(P->en); |
107 |
pj_dalloc(P); |
108 |
} |
109 |
} |
110 |
static PJ *
|
111 |
setup(PJ *P) { /* general initialization */
|
112 |
if (P->es) {
|
113 |
if (!(P->en = pj_enfn(P->es)))
|
114 |
E_ERROR_0; |
115 |
P->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en); |
116 |
P->esp = P->es / (1. - P->es);
|
117 |
P->inv = e_inverse; |
118 |
P->fwd = e_forward; |
119 |
} else {
|
120 |
aks0 = P->k0; |
121 |
aks5 = .5 * aks0;
|
122 |
P->inv = s_inverse; |
123 |
P->fwd = s_forward; |
124 |
} |
125 |
return P;
|
126 |
} |
127 |
ENTRY1(tmerc, en) |
128 |
ENDENTRY(setup(P)) |
129 |
ENTRY1(utm, en) |
130 |
int zone;
|
131 |
|
132 |
if (!P->es) E_ERROR(-34); |
133 |
P->y0 = pj_param(P->params, "bsouth").i ? 10000000. : 0.; |
134 |
P->x0 = 500000.; |
135 |
if (pj_param(P->params, "tzone").i) /* zone input ? */ |
136 |
if ((zone = pj_param(P->params, "izone").i) > 0 && zone <= 60) |
137 |
--zone; |
138 |
else
|
139 |
E_ERROR(-35)
|
140 |
else /* nearest central meridian input */ |
141 |
if ((zone = floor((adjlon(P->lam0) + PI) * 30. / PI)) < 0) |
142 |
zone = 0;
|
143 |
else if (zone >= 60) |
144 |
zone = 59;
|
145 |
P->lam0 = (zone + .5) * PI / 30. - PI; |
146 |
P->k0 = 0.9996; |
147 |
P->phi0 = 0.;
|
148 |
ENDENTRY(setup(P)) |