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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
#ifndef lint
static const char SCCSID[]="@(#)PJ_imw_p.c 4.1 94/05/22 GIE REL";
#endif
#define PROJ_PARMS__ \
double P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
double phi_1, phi_2, lam_1; \
double *en; \
int mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
#define PJ_LIB__
#include <projects.h>
PROJ_HEAD(imw_p, "International Map of the World Polyconic")
"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
#define TOL 1e-10
#define EPS 1e-10
static int
phi12(PJ *P, double *del, double *sig) {
int err = 0;
if (!pj_param(P->params, "tlat_1").i ||
!pj_param(P->params, "tlat_2").i) {
err = -41;
} else {
P->phi_1 = pj_param(P->params, "rlat_1").f;
P->phi_2 = pj_param(P->params, "rlat_2").f;
*del = 0.5 * (P->phi_2 - P->phi_1);
*sig = 0.5 * (P->phi_2 + P->phi_1);
err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
}
return err;
}
static XY
loc_for(LP lp, PJ *P, double *yc) {
XY xy;
if (! lp.phi) {
xy.x = lp.lam;
xy.y = 0.;
} else {
double xa, ya, xb, yb, xc, yc, D, B, m, sp, t, R, C;
sp = sin(lp.phi);
m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en);
xa = P->Pp + P->Qp * m;
ya = P->P + P->Q * m;
R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
C = sqrt(R * R - xa * xa);
if (lp.phi < 0.) C = - C;
C += ya - R;
if (P->mode < 0) {
xb = lp.lam;
yb = P->C2;
} else {
t = lp.lam * P->sphi_2;
xb = P->R_2 * sin(t);
yb = P->C2 + P->R_2 * (1. - cos(t));
}
if (P->mode > 0) {
xc = lp.lam;
yc = 0.;
} else {
t = lp.lam * P->sphi_1;
xc = P->R_1 * sin(t);
yc = P->R_1 * (1. - cos(t));
}
D = (xb - xc)/(yb - yc);
B = xc + D * (C + R - yc);
xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
if (lp.phi > 0)
xy.x = - xy.x;
xy.x = (B + xy.x) / (1. + D * D);
xy.y = sqrt(R * R - xy.x * xy.x);
if (lp.phi > 0)
xy.y = - xy.y;
xy.y += C + R;
}
return (xy);
}
FORWARD(e_forward); /* ellipsoid */
double yc;
xy = loc_for(lp, P, &yc);
return (xy);
}
INVERSE(e_inverse); /* ellipsoid */
XY t;
double yc;
lp.phi = P->phi_2;
lp.lam = xy.x / cos(lp.phi);
do {
t = loc_for(lp, P, &yc);
lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1;
lp.lam = lp.lam * xy.x / t.x;
} while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
return (lp);
}
static void
xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
double t, F;
*sp = sin(phi);
*R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
F = P->lam_1 * *sp;
*y = *R * (1 - cos(F));
*x = *R * sin(F);
}
FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
ENTRY1(imw_p, en)
double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
int i;
if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
if (i = phi12(P, &del, &sig))
E_ERROR(i);
if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */
del = P->phi_1;
P->phi_1 = P->phi_2;
P->phi_2 = del;
}
if (pj_param(P->params, "tlon_1").i)
P->lam_1 = pj_param(P->params, "rlon_1").f;
else { /* use predefined based upon latitude */
sig = fabs(sig * RAD_TO_DEG);
if (sig <= 60) sig = 2.;
else if (sig <= 76) sig = 4.;
else sig = 8.;
P->lam_1 = sig * DEG_TO_RAD;
}
P->mode = 0;
if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
else {
P->mode = 1;
y1 = 0.;
x1 = P->lam_1;
}
if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
else {
P->mode = -1;
T2 = 0.;
x2 = P->lam_1;
}
m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en);
m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en);
t = m2 - m1;
s = x2 - x1;
y2 = sqrt(t * t - s * s) + y1;
P->C2 = y2 - T2;
t = 1. / t;
P->P = (m2 * y1 - m1 * y2) * t;
P->Q = (y2 - y1) * t;
P->Pp = (m2 * x1 - m1 * x2) * t;
P->Qp = (x2 - x1) * t;
P->fwd = e_forward;
P->inv = e_inverse;
ENDENTRY(P)
|
