/Proj4/PJ_somerc.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)PJ_somerc.c	4.1	95/08/09	GIE	REL";
 3#endif
 4#define PROJ_PARMS__ \
 5	double	K, c, hlf_e, kR, cosp0, sinp0;
 6#define PJ_LIB__
 7#include	"projects.h"
 8PROJ_HEAD(somerc, "Swiss. Obl. Mercator") "\n\tCyl, Ell\n\tFor CH1903";
 9#define EPS	1.e-10
10#define NITER 6
11FORWARD(e_forward);
12	double phip, lamp, phipp, lampp, sp, cp;
13
14	sp = P->e * sin(lp.phi);
15	phip = 2.* atan( exp( P->c * (
16		log(tan(FORTPI + 0.5 * lp.phi)) - P->hlf_e * log((1. + sp)/(1. - sp)))
17		+ P->K)) - HALFPI;
18	lamp = P->c * lp.lam;
19	cp = cos(phip);
20	phipp = aasin(P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp));
21	lampp = aasin(cp * sin(lamp) / cos(phipp));
22	xy.x = P->kR * lampp;
23	xy.y = P->kR * log(tan(FORTPI + 0.5 * phipp));
24	return (xy);
25}
26INVERSE(e_inverse); /* ellipsoid & spheroid */
27	double phip, lamp, phipp, lampp, cp, esp, con, delp;
28	int i;
29
30	phipp = 2. * (atan(exp(xy.y / P->kR)) - FORTPI);
31	lampp = xy.x / P->kR;
32	cp = cos(phipp);
33	phip = aasin(P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp));
34	lamp = aasin(cp * sin(lampp) / cos(phip));
35	con = (P->K - log(tan(FORTPI + 0.5 * phip)))/P->c;
36	for (i = NITER; i ; --i) {
37		esp = P->e * sin(phip);
38		delp = (con + log(tan(FORTPI + 0.5 * phip)) - P->hlf_e *
39			log((1. + esp)/(1. - esp)) ) *
40			(1. - esp * esp) * cos(phip) * P->rone_es;
41		phip -= delp;
42		if (fabs(delp) < EPS)
43			break;
44	}
45	if (i) {
46		lp.phi = phip;
47		lp.lam = lamp / P->c;
48	} else
49		I_ERROR
50	return (lp);
51}
52FREEUP; if (P) pj_dalloc(P); }
53ENTRY0(somerc)
54	double cp, phip0, sp;
55
56	P->hlf_e = 0.5 * P->e;
57	cp = cos(P->phi0);
58	cp *= cp;
59	P->c = sqrt(1 + P->es * cp * cp * P->rone_es);
60	sp = sin(P->phi0);
61	P->cosp0 = cos( phip0 = aasin(P->sinp0 = sp / P->c) );
62	sp *= P->e;
63	P->K = log(tan(FORTPI + 0.5 * phip0)) - P->c * (
64		log(tan(FORTPI + 0.5 * P->phi0)) - P->hlf_e *
65		log((1. + sp) / (1. - sp)));
66	P->kR = P->k0 * sqrt(P->one_es) / (1. - sp * sp);
67	P->inv = e_inverse;
68	P->fwd = e_forward;
69ENDENTRY(P)