PageRenderTime 34ms CodeModel.GetById 12ms app.highlight 18ms RepoModel.GetById 1ms app.codeStats 0ms

/Proj4/PJ_mod_ster.c

http://github.com/route-me/route-me
C | 214 lines | 204 code | 9 blank | 1 comment | 13 complexity | 0691829146fffd6d10b33c1bc645a0d9 MD5 | raw file
  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_mod_ster.c	4.1	94/02/15	GIE	REL";
  3#endif
  4/* based upon Snyder and Linck, USGS-NMD */
  5#define PROJ_PARMS__ \
  6    COMPLEX	*zcoeff; \
  7	double	cchio, schio; \
  8	int		n;
  9#define PJ_LIB__
 10#include	"projects.h"
 11PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)";
 12PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)";
 13PROJ_HEAD(gs48, "Mod. Stererographics of 48 U.S.") "\n\tAzi(mod)";
 14PROJ_HEAD(alsk, "Mod. Stererographics of Alaska") "\n\tAzi(mod)";
 15PROJ_HEAD(gs50, "Mod. Stererographics of 50 U.S.") "\n\tAzi(mod)";
 16#define EPSLN 1e-10
 17
 18FORWARD(e_forward); /* ellipsoid */
 19	double sinlon, coslon, esphi, chi, schi, cchi, s;
 20	COMPLEX p;
 21
 22	sinlon = sin(lp.lam);
 23	coslon = cos(lp.lam);
 24	esphi = P->e * sin(lp.phi);
 25	chi = 2. * atan(tan((HALFPI + lp.phi) * .5) *
 26		pow((1. - esphi) / (1. + esphi), P->e * .5)) - HALFPI;
 27	schi = sin(chi);
 28	cchi = cos(chi);
 29	s = 2. / (1. + P->schio * schi + P->cchio * cchi * coslon);
 30	p.r = s * cchi * sinlon;
 31	p.i = s * (P->cchio * schi - P->schio * cchi * coslon);
 32	p = pj_zpoly1(p, P->zcoeff, P->n);
 33	xy.x = p.r;
 34	xy.y = p.i;
 35	return xy;
 36}
 37INVERSE(e_inverse); /* ellipsoid */
 38	int nn;
 39	COMPLEX p, fxy, fpxy, dp;
 40	double den, rh, z, sinz, cosz, chi, phi, dphi, esphi;
 41
 42	p.r = xy.x;
 43	p.i = xy.y;
 44	for (nn = 20; nn ;--nn) {
 45		fxy = pj_zpolyd1(p, P->zcoeff, P->n, &fpxy);
 46		fxy.r -= xy.x;
 47		fxy.i -= xy.y;
 48		den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
 49		dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
 50		dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
 51		p.r += dp.r;
 52		p.i += dp.i;
 53		if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
 54			break;
 55	}
 56	if (nn) {
 57		rh = hypot(p.r, p.i);
 58		z = 2. * atan(.5 * rh);
 59		sinz = sin(z);
 60		cosz = cos(z);
 61		lp.lam = P->lam0;
 62		if (fabs(rh) <= EPSLN) {
 63			lp.phi = P->phi0;
 64			return lp;
 65		}
 66		chi = aasin(cosz * P->schio + p.i * sinz * P->cchio / rh);
 67		phi = chi;
 68		for (nn = 20; nn ;--nn) {
 69			esphi = P->e * sin(phi);
 70			dphi = 2. * atan(tan((HALFPI + chi) * .5) *
 71				pow((1. + esphi) / (1. - esphi), P->e * .5)) - HALFPI - phi;
 72			phi += dphi;
 73			if (fabs(dphi) <= EPSLN)
 74				break;
 75		}
 76	}
 77	if (nn) {
 78		lp.phi = phi;
 79		lp.lam = atan2(p.r * sinz, rh * P->cchio * cosz - p.i * 
 80			P->schio * sinz);
 81    } else
 82		lp.lam = lp.phi = HUGE_VAL;
 83	return lp;
 84}
 85FREEUP; if (P) pj_dalloc(P); }
 86	static PJ *
 87setup(PJ *P) { /* general initialization */
 88	double esphi, chio;
 89
 90	if (P->es) {
 91		esphi = P->e * sin(P->phi0);
 92		chio = 2. * atan(tan((HALFPI + P->phi0) * .5) *
 93			pow((1. - esphi) / (1. + esphi), P->e * .5)) - HALFPI;
 94	} else
 95		chio = P->phi0;
 96	P->schio = sin(chio);
 97	P->cchio = cos(chio);
 98	P->inv = e_inverse; P->fwd = e_forward;
 99	return P;
100}
101ENTRY0(mil_os)
102	static COMPLEX /* Miller Oblated Stereographic */
103AB[] = {
104	{0.924500,	0.},
105	{0.,			0.},
106	{0.019430,	0.}
107};
108
109	P->n = 2;
110	P->lam0 = DEG_TO_RAD * 20.;
111	P->phi0 = DEG_TO_RAD * 18.;
112	P->zcoeff = AB;
113	P->es = 0.;
114ENDENTRY(setup(P))
115ENTRY0(lee_os)
116	static COMPLEX /* Lee Oblated Stereographic */
117AB[] = {
118	{0.721316,	0.},
119	{0.,			0.},
120        {-0.0088162,	 -0.00617325}
121};
122
123	P->n = 2;
124	P->lam0 = DEG_TO_RAD * -165.;
125	P->phi0 = DEG_TO_RAD * -10.;
126	P->zcoeff = AB;
127	P->es = 0.;
128ENDENTRY(setup(P))
129ENTRY0(gs48)
130	static COMPLEX /* 48 United States */
131AB[] = {
132	{0.98879,	0.},
133	{0.,		0.},
134	{-0.050909,	0.},
135	{0.,		0.},
136        {0.075528,	0.}
137};
138
139	P->n = 4;
140	P->lam0 = DEG_TO_RAD * -96.;
141	P->phi0 = DEG_TO_RAD * -39.;
142	P->zcoeff = AB;
143	P->es = 0.;
144	P->a = 6370997.;
145ENDENTRY(setup(P))
146ENTRY0(alsk)
147	static COMPLEX
148ABe[] = { /* Alaska ellipsoid */
149	{.9945303,	0.},
150	{.0052083,	-.0027404},
151	{.0072721,	.0048181},
152	{-.0151089,	-.1932526},
153	{.0642675,	-.1381226},
154	{.3582802,	-.2884586}},
155ABs[] = { /* Alaska sphere */
156	{.9972523,	0.},
157	{.0052513,	-.0041175},
158	{.0074606,	.0048125},
159	{-.0153783,	-.1968253},
160	{.0636871,	-.1408027},
161        {.3660976,	-.2937382}
162};
163
164	P->n = 5;
165	P->lam0 = DEG_TO_RAD * -152.;
166	P->phi0 = DEG_TO_RAD * 64.;
167	if (P->es) { /* fixed ellipsoid/sphere */
168		P->zcoeff = ABe;
169		P->a = 6378206.4;
170		P->e = sqrt(P->es = 0.00676866);
171	} else {
172		P->zcoeff = ABs;
173		P->a = 6370997.;
174	}
175ENDENTRY(setup(P))
176ENTRY0(gs50)
177	static COMPLEX
178ABe[] = { /* GS50 ellipsoid */
179	{.9827497,	0.},
180	{.0210669,	.0053804},
181	{-.1031415,	-.0571664},
182	{-.0323337,	-.0322847},
183	{.0502303,	.1211983},
184	{.0251805,	.0895678},
185	{-.0012315,	-.1416121},
186	{.0072202,	-.1317091},
187	{-.0194029,	.0759677},
188        {-.0210072,	.0834037}
189},
190ABs[] = { /* GS50 sphere */
191	{.9842990,	0.},
192	{.0211642,	.0037608},
193	{-.1036018,	-.0575102},
194	{-.0329095,	-.0320119},
195	{.0499471,	.1223335},
196	{.0260460,	.0899805},
197	{.0007388,	-.1435792},
198	{.0075848,	-.1334108},
199	{-.0216473,	.0776645},
200        {-.0225161,	.0853673}
201};
202
203	P->n = 9;
204	P->lam0 = DEG_TO_RAD * -120.;
205	P->phi0 = DEG_TO_RAD * 45.;
206	if (P->es) { /* fixed ellipsoid/sphere */
207		P->zcoeff = ABe;
208		P->a = 6378206.4;
209		P->e = sqrt(P->es = 0.00676866);
210	} else {
211		P->zcoeff = ABs;
212		P->a = 6370997.;
213	}
214ENDENTRY(setup(P))