svn-gvsig-desktop / tags / v1_12_0_Build_1417 / libraries / libjni-proj4 / src / PJ_imw_p.c @ 40245
History | View | Annotate | Download (3.72 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_imw_p.c 4.1 94/05/22 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
|
6 |
double phi_1, phi_2, lam_1; \
|
7 |
double *en; \
|
8 |
int mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */ |
9 |
#define PJ_LIB__
|
10 |
#include <projects.h> |
11 |
PROJ_HEAD(imw_p, "International Map of the World Polyconic")
|
12 |
"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
|
13 |
#define TOL 1e-10 |
14 |
#define EPS 1e-10 |
15 |
static int |
16 |
phi12(PJ *P, double *del, double *sig) { |
17 |
int err = 0; |
18 |
|
19 |
if (!pj_param(P->params, "tlat_1").i || |
20 |
!pj_param(P->params, "tlat_2").i) {
|
21 |
err = -41;
|
22 |
} else {
|
23 |
P->phi_1 = pj_param(P->params, "rlat_1").f;
|
24 |
P->phi_2 = pj_param(P->params, "rlat_2").f;
|
25 |
*del = 0.5 * (P->phi_2 - P->phi_1); |
26 |
*sig = 0.5 * (P->phi_2 + P->phi_1); |
27 |
err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0; |
28 |
} |
29 |
return err;
|
30 |
} |
31 |
static XY
|
32 |
loc_for(LP lp, PJ *P, double *yc) {
|
33 |
XY xy; |
34 |
|
35 |
if (! lp.phi) {
|
36 |
xy.x = lp.lam; |
37 |
xy.y = 0.;
|
38 |
} else {
|
39 |
double xa, ya, xb, yb, xc, yc, D, B, m, sp, t, R, C;
|
40 |
|
41 |
sp = sin(lp.phi); |
42 |
m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en); |
43 |
xa = P->Pp + P->Qp * m; |
44 |
ya = P->P + P->Q * m; |
45 |
R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp)); |
46 |
C = sqrt(R * R - xa * xa); |
47 |
if (lp.phi < 0.) C = - C; |
48 |
C += ya - R; |
49 |
if (P->mode < 0) { |
50 |
xb = lp.lam; |
51 |
yb = P->C2; |
52 |
} else {
|
53 |
t = lp.lam * P->sphi_2; |
54 |
xb = P->R_2 * sin(t); |
55 |
yb = P->C2 + P->R_2 * (1. - cos(t));
|
56 |
} |
57 |
if (P->mode > 0) { |
58 |
xc = lp.lam; |
59 |
yc = 0.;
|
60 |
} else {
|
61 |
t = lp.lam * P->sphi_1; |
62 |
xc = P->R_1 * sin(t); |
63 |
yc = P->R_1 * (1. - cos(t));
|
64 |
} |
65 |
D = (xb - xc)/(yb - yc); |
66 |
B = xc + D * (C + R - yc); |
67 |
xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
|
68 |
if (lp.phi > 0) |
69 |
xy.x = - xy.x; |
70 |
xy.x = (B + xy.x) / (1. + D * D);
|
71 |
xy.y = sqrt(R * R - xy.x * xy.x); |
72 |
if (lp.phi > 0) |
73 |
xy.y = - xy.y; |
74 |
xy.y += C + R; |
75 |
} |
76 |
return (xy);
|
77 |
} |
78 |
FORWARD(e_forward); /* ellipsoid */
|
79 |
double yc;
|
80 |
xy = loc_for(lp, P, &yc); |
81 |
return (xy);
|
82 |
} |
83 |
INVERSE(e_inverse); /* ellipsoid */
|
84 |
XY t; |
85 |
double yc;
|
86 |
|
87 |
lp.phi = P->phi_2; |
88 |
lp.lam = xy.x / cos(lp.phi); |
89 |
do {
|
90 |
t = loc_for(lp, P, &yc); |
91 |
lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1; |
92 |
lp.lam = lp.lam * xy.x / t.x; |
93 |
} while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
|
94 |
return (lp);
|
95 |
} |
96 |
static void |
97 |
xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) { |
98 |
double F;
|
99 |
|
100 |
*sp = sin(phi); |
101 |
*R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp )); |
102 |
F = P->lam_1 * *sp; |
103 |
*y = *R * (1 - cos(F));
|
104 |
*x = *R * sin(F); |
105 |
} |
106 |
FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } } |
107 |
ENTRY1(imw_p, en) |
108 |
double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
|
109 |
int i;
|
110 |
|
111 |
if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
|
112 |
if( (i = phi12(P, &del, &sig)) != 0) |
113 |
E_ERROR(i); |
114 |
if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */ |
115 |
del = P->phi_1; |
116 |
P->phi_1 = P->phi_2; |
117 |
P->phi_2 = del; |
118 |
} |
119 |
if (pj_param(P->params, "tlon_1").i) |
120 |
P->lam_1 = pj_param(P->params, "rlon_1").f;
|
121 |
else { /* use predefined based upon latitude */ |
122 |
sig = fabs(sig * RAD_TO_DEG); |
123 |
if (sig <= 60) sig = 2.; |
124 |
else if (sig <= 76) sig = 4.; |
125 |
else sig = 8.; |
126 |
P->lam_1 = sig * DEG_TO_RAD; |
127 |
} |
128 |
P->mode = 0;
|
129 |
if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
|
130 |
else {
|
131 |
P->mode = 1;
|
132 |
y1 = 0.;
|
133 |
x1 = P->lam_1; |
134 |
} |
135 |
if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
|
136 |
else {
|
137 |
P->mode = -1;
|
138 |
T2 = 0.;
|
139 |
x2 = P->lam_1; |
140 |
} |
141 |
m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en); |
142 |
m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en); |
143 |
t = m2 - m1; |
144 |
s = x2 - x1; |
145 |
y2 = sqrt(t * t - s * s) + y1; |
146 |
P->C2 = y2 - T2; |
147 |
t = 1. / t;
|
148 |
P->P = (m2 * y1 - m1 * y2) * t; |
149 |
P->Q = (y2 - y1) * t; |
150 |
P->Pp = (m2 * x1 - m1 * x2) * t; |
151 |
P->Qp = (x2 - x1) * t; |
152 |
P->fwd = e_forward; |
153 |
P->inv = e_inverse; |
154 |
ENDENTRY(P) |