/Proj4/PJ_cass.c

 1#ifndef lint
 2static 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"
18PROJ_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
25FORWARD(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}
38FORWARD(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}
43INVERSE(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}
60INVERSE(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}
65FREEUP;
66	if (P) {
67		if (P->en)
68			pj_dalloc(P->en);
69		pj_dalloc(P);
70	}
71}
72ENTRY1(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	}
82ENDENTRY(P)