svn-gvsig-desktop / trunk / libraries / libjni-proj4 / src / PJ_cass.c @ 7098
History | View | Annotate | Download (2.07 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_cass.c 4.1 94/02/15 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double m0; \
|
6 |
double n; \
|
7 |
double t; \
|
8 |
double a1; \
|
9 |
double c; \
|
10 |
double r; \
|
11 |
double dd; \
|
12 |
double d2; \
|
13 |
double a2; \
|
14 |
double tn; \
|
15 |
double *en;
|
16 |
#define PJ_LIB__
|
17 |
# include <projects.h> |
18 |
PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell"; |
19 |
# define EPS10 1e-10 |
20 |
# define C1 .16666666666666666666 |
21 |
# define C2 .00833333333333333333 |
22 |
# define C3 .04166666666666666666 |
23 |
# define C4 .33333333333333333333 |
24 |
# define C5 .06666666666666666666 |
25 |
FORWARD(e_forward); /* ellipsoid */
|
26 |
xy.y = pj_mlfn(lp.phi, P->n = sin(lp.phi), P->c = cos(lp.phi), P->en); |
27 |
P->n = 1./sqrt(1. - P->es * P->n * P->n); |
28 |
P->tn = tan(lp.phi); P->t = P->tn * P->tn; |
29 |
P->a1 = lp.lam * P->c; |
30 |
P->c *= P->es * P->c / (1 - P->es);
|
31 |
P->a2 = P->a1 * P->a1; |
32 |
xy.x = P->n * P->a1 * (1. - P->a2 * P->t *
|
33 |
(C1 - (8. - P->t + 8. * P->c) * P->a2 * C2)); |
34 |
xy.y -= P->m0 - P->n * P->tn * P->a2 * |
35 |
(.5 + (5. - P->t + 6. * P->c) * P->a2 * C3); |
36 |
return (xy);
|
37 |
} |
38 |
FORWARD(s_forward); /* spheroid */
|
39 |
xy.x = asin(cos(lp.phi) * sin(lp.lam)); |
40 |
xy.y = atan2(tan(lp.phi) , cos(lp.lam)) - P->phi0; |
41 |
return (xy);
|
42 |
} |
43 |
INVERSE(e_inverse); /* ellipsoid */
|
44 |
double ph1;
|
45 |
|
46 |
ph1 = pj_inv_mlfn(P->m0 + xy.y, P->es, P->en); |
47 |
P->tn = tan(ph1); P->t = P->tn * P->tn; |
48 |
P->n = sin(ph1); |
49 |
P->r = 1. / (1. - P->es * P->n * P->n); |
50 |
P->n = sqrt(P->r); |
51 |
P->r *= (1. - P->es) * P->n;
|
52 |
P->dd = xy.x / P->n; |
53 |
P->d2 = P->dd * P->dd; |
54 |
lp.phi = ph1 - (P->n * P->tn / P->r) * P->d2 * |
55 |
(.5 - (1. + 3. * P->t) * P->d2 * C3); |
56 |
lp.lam = P->dd * (1. + P->t * P->d2 *
|
57 |
(-C4 + (1. + 3. * P->t) * P->d2 * C5)) / cos(ph1); |
58 |
return (lp);
|
59 |
} |
60 |
INVERSE(s_inverse); /* spheroid */
|
61 |
lp.phi = asin(sin(P->dd = xy.y + P->phi0) * cos(xy.x)); |
62 |
lp.lam = atan2(tan(xy.x), cos(P->dd)); |
63 |
return (lp);
|
64 |
} |
65 |
FREEUP; |
66 |
if (P) {
|
67 |
if (P->en)
|
68 |
pj_dalloc(P->en); |
69 |
pj_dalloc(P); |
70 |
} |
71 |
} |
72 |
ENTRY1(cass, en) |
73 |
if (P->es) {
|
74 |
if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
|
75 |
P->m0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en); |
76 |
P->inv = e_inverse; |
77 |
P->fwd = e_forward; |
78 |
} else {
|
79 |
P->inv = s_inverse; |
80 |
P->fwd = s_forward; |
81 |
} |
82 |
ENDENTRY(P) |