PJ_somerc.c
1.81 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
#define PROJ_PARMS__ \
double K, c, hlf_e, kR, cosp0, sinp0;
#define PJ_LIB__
#include <projects.h>
PROJ_HEAD(somerc, "Swiss. Obl. Mercator") "\n\tCyl, Ell\n\tFor CH1903";
#define EPS 1.e-10
#define NITER 6
FORWARD(e_forward);
double phip, lamp, phipp, lampp, sp, cp;
sp = P->e * sin(lp.phi);
phip = 2.* atan( exp( P->c * (
log(tan(FORTPI + 0.5 * lp.phi)) - P->hlf_e * log((1. + sp)/(1. - sp)))
+ P->K)) - HALFPI;
lamp = P->c * lp.lam;
cp = cos(phip);
phipp = aasin(P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp));
lampp = aasin(cp * sin(lamp) / cos(phipp));
xy.x = P->kR * lampp;
xy.y = P->kR * log(tan(FORTPI + 0.5 * phipp));
return (xy);
}
INVERSE(e_inverse); /* ellipsoid & spheroid */
double phip, lamp, phipp, lampp, cp, esp, con, delp;
int i;
phipp = 2. * (atan(exp(xy.y / P->kR)) - FORTPI);
lampp = xy.x / P->kR;
cp = cos(phipp);
phip = aasin(P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp));
lamp = aasin(cp * sin(lampp) / cos(phip));
con = (P->K - log(tan(FORTPI + 0.5 * phip)))/P->c;
for (i = NITER; i ; --i) {
esp = P->e * sin(phip);
delp = (con + log(tan(FORTPI + 0.5 * phip)) - P->hlf_e *
log((1. + esp)/(1. - esp)) ) *
(1. - esp * esp) * cos(phip) * P->rone_es;
phip -= delp;
if (fabs(delp) < EPS)
break;
}
if (i) {
lp.phi = phip;
lp.lam = lamp / P->c;
} else
I_ERROR
return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
ENTRY0(somerc)
double cp, phip0, sp;
P->hlf_e = 0.5 * P->e;
cp = cos(P->phi0);
cp *= cp;
P->c = sqrt(1 + P->es * cp * cp * P->rone_es);
sp = sin(P->phi0);
P->cosp0 = cos( phip0 = aasin(P->sinp0 = sp / P->c) );
sp *= P->e;
P->K = log(tan(FORTPI + 0.5 * phip0)) - P->c * (
log(tan(FORTPI + 0.5 * P->phi0)) - P->hlf_e *
log((1. + sp) / (1. - sp)));
P->kR = P->k0 * sqrt(P->one_es) / (1. - sp * sp);
P->inv = e_inverse;
P->fwd = e_forward;
ENDENTRY(P)