PJ_cass.c
2.07 KB
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
#ifndef lint
static const char SCCSID[]="@(#)PJ_cass.c 4.1 94/02/15 GIE REL";
#endif
#define PROJ_PARMS__ \
double m0; \
double n; \
double t; \
double a1; \
double c; \
double r; \
double dd; \
double d2; \
double a2; \
double tn; \
double *en;
#define PJ_LIB__
# include "projects.h"
PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";
# define EPS10 1e-10
# define C1 .16666666666666666666
# define C2 .00833333333333333333
# define C3 .04166666666666666666
# define C4 .33333333333333333333
# define C5 .06666666666666666666
FORWARD(e_forward); /* ellipsoid */
xy.y = pj_mlfn(lp.phi, P->n = sin(lp.phi), P->c = cos(lp.phi), P->en);
P->n = 1./sqrt(1. - P->es * P->n * P->n);
P->tn = tan(lp.phi); P->t = P->tn * P->tn;
P->a1 = lp.lam * P->c;
P->c *= P->es * P->c / (1 - P->es);
P->a2 = P->a1 * P->a1;
xy.x = P->n * P->a1 * (1. - P->a2 * P->t *
(C1 - (8. - P->t + 8. * P->c) * P->a2 * C2));
xy.y -= P->m0 - P->n * P->tn * P->a2 *
(.5 + (5. - P->t + 6. * P->c) * P->a2 * C3);
return (xy);
}
FORWARD(s_forward); /* spheroid */
xy.x = asin(cos(lp.phi) * sin(lp.lam));
xy.y = atan2(tan(lp.phi) , cos(lp.lam)) - P->phi0;
return (xy);
}
INVERSE(e_inverse); /* ellipsoid */
double ph1;
ph1 = pj_inv_mlfn(P->m0 + xy.y, P->es, P->en);
P->tn = tan(ph1); P->t = P->tn * P->tn;
P->n = sin(ph1);
P->r = 1. / (1. - P->es * P->n * P->n);
P->n = sqrt(P->r);
P->r *= (1. - P->es) * P->n;
P->dd = xy.x / P->n;
P->d2 = P->dd * P->dd;
lp.phi = ph1 - (P->n * P->tn / P->r) * P->d2 *
(.5 - (1. + 3. * P->t) * P->d2 * C3);
lp.lam = P->dd * (1. + P->t * P->d2 *
(-C4 + (1. + 3. * P->t) * P->d2 * C5)) / cos(ph1);
return (lp);
}
INVERSE(s_inverse); /* spheroid */
lp.phi = asin(sin(P->dd = xy.y + P->phi0) * cos(xy.x));
lp.lam = atan2(tan(xy.x), cos(P->dd));
return (lp);
}
FREEUP;
if (P) {
if (P->en)
pj_dalloc(P->en);
pj_dalloc(P);
}
}
ENTRY1(cass, en)
if (P->es) {
if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
P->m0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en);
P->inv = e_inverse;
P->fwd = e_forward;
} else {
P->inv = s_inverse;
P->fwd = s_forward;
}
ENDENTRY(P)