pj_ell_set.c
3.14 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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
/* set ellipsoid parameters a and es */
#include <projects.h>
#include <string.h>
#define SIXTH .1666666666666666667 /* 1/6 */
#define RA4 .04722222222222222222 /* 17/360 */
#define RA6 .02215608465608465608 /* 67/3024 */
#define RV4 .06944444444444444444 /* 5/72 */
#define RV6 .04243827160493827160 /* 55/1296 */
int /* initialize geographic shape parameters */
pj_ell_set(paralist *pl, double *a, double *es) {
int i;
double b=0.0, e;
char *name;
paralist *start = 0, *curr;
/* check for varying forms of ellipsoid input */
*a = *es = 0.;
/* R takes precedence */
if (pj_param(pl, "tR").i)
*a = pj_param(pl, "dR").f;
else { /* probable elliptical figure */
/* check if ellps present and temporarily append its values to pl */
if ((name = pj_param(pl, "sellps").s)) {
char *s;
for (start = pl; start && start->next ; start = start->next) ;
curr = start;
for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i) ;
if (!s) { pj_errno = -9; return 1; }
curr = curr->next = pj_mkparam(pj_ellps[i].major);
curr = curr->next = pj_mkparam(pj_ellps[i].ell);
}
*a = pj_param(pl, "da").f;
if (pj_param(pl, "tes").i) /* eccentricity squared */
*es = pj_param(pl, "des").f;
else if (pj_param(pl, "te").i) { /* eccentricity */
e = pj_param(pl, "de").f;
*es = e * e;
} else if (pj_param(pl, "trf").i) { /* recip flattening */
*es = pj_param(pl, "drf").f;
if (!*es) {
pj_errno = -10;
goto bomb;
}
*es = 1./ *es;
*es = *es * (2. - *es);
} else if (pj_param(pl, "tf").i) { /* flattening */
*es = pj_param(pl, "df").f;
*es = *es * (2. - *es);
} else if (pj_param(pl, "tb").i) { /* minor axis */
b = pj_param(pl, "db").f;
*es = 1. - (b * b) / (*a * *a);
} /* else *es == 0. and sphere of radius *a */
if (!b)
b = *a * sqrt(1. - *es);
/* following options turn ellipsoid into equivalent sphere */
if (pj_param(pl, "bR_A").i) { /* sphere--area of ellipsoid */
*a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
*es = 0.;
} else if (pj_param(pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
*a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));
*es = 0.;
} else if (pj_param(pl, "bR_a").i) { /* sphere--arithmetic mean */
*a = .5 * (*a + b);
*es = 0.;
} else if (pj_param(pl, "bR_g").i) { /* sphere--geometric mean */
*a = sqrt(*a * b);
*es = 0.;
} else if (pj_param(pl, "bR_h").i) { /* sphere--harmonic mean */
*a = 2. * *a * b / (*a + b);
*es = 0.;
} else if ((i = pj_param(pl, "tR_lat_a").i) || /* sphere--arith. */
pj_param(pl, "tR_lat_g").i) { /* or geom. mean at latitude */
double tmp;
tmp = sin(pj_param(pl, i ? "rR_lat_a" : "rR_lat_g").f);
if (fabs(tmp) > HALFPI) {
pj_errno = -11;
goto bomb;
}
tmp = 1. - *es * tmp * tmp;
*a *= i ? .5 * (1. - *es + tmp) / ( tmp * sqrt(tmp)) :
sqrt(1. - *es) / tmp;
*es = 0.;
}
bomb:
if (start) { /* clean up temporary extension of list */
pj_dalloc(start->next->next);
pj_dalloc(start->next);
start->next = 0;
}
if (pj_errno)
return 1;
}
/* some remaining checks */
if (*es < 0.)
{ pj_errno = -12; return 1; }
if (*a <= 0.)
{ pj_errno = -13; return 1; }
return 0;
}