/Proj4/PJ_lsat.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_lsat.c	4.1	94/02/15	GIE	REL";
  3#endif
  4/* based upon Snyder and Linck, USGS-NMD */
  5#define PROJ_PARMS__ \
  6    double a2, a4, b, c1, c3; \
  7    double q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
  8#define PJ_LIB__
  9#include	"projects.h"
 10PROJ_HEAD(lsat, "Space oblique for LANDSAT")
 11	"\n\tCyl, Sph&Ell\n\tlsat= path=";
 12#define TOL 1e-7
 13#define PI_HALFPI 4.71238898038468985766
 14#define TWOPI_HALFPI 7.85398163397448309610
 15	static void
 16seraz0(double lam, double mult, PJ *P) {
 17    double sdsq, h, s, fc, sd, sq, d__1;
 18
 19    lam *= DEG_TO_RAD;
 20    sd = sin(lam);
 21    sdsq = sd * sd;
 22    s = P->p22 * P->sa * cos(lam) * sqrt((1. + P->t * sdsq) / ((
 23	    1. + P->w * sdsq) * (1. + P->q * sdsq)));
 24    d__1 = 1. + P->q * sdsq;
 25    h = sqrt((1. + P->q * sdsq) / (1. + P->w * sdsq)) * ((1. + 
 26	    P->w * sdsq) / (d__1 * d__1) - P->p22 * P->ca);
 27    sq = sqrt(P->xj * P->xj + s * s);
 28    P->b += fc = mult * (h * P->xj - s * s) / sq;
 29    P->a2 += fc * cos(lam + lam);
 30    P->a4 += fc * cos(lam * 4.);
 31    fc = mult * s * (h + P->xj) / sq;
 32    P->c1 += fc * cos(lam);
 33    P->c3 += fc * cos(lam * 3.);
 34}
 35FORWARD(e_forward); /* ellipsoid */
 36    int l, nn;
 37    double lamt, xlam, sdsq, c, d, s, lamdp, phidp, lampp, tanph,
 38		lamtp, cl, sd, sp, fac, sav, tanphi;
 39
 40	if (lp.phi > HALFPI)
 41	    lp.phi = HALFPI;
 42	else if (lp.phi < -HALFPI)
 43	    lp.phi = -HALFPI;
 44	lampp = lp.phi >= 0. ? HALFPI : PI_HALFPI;
 45	tanphi = tan(lp.phi);
 46	for (nn = 0;;) {
 47		sav = lampp;
 48		lamtp = lp.lam + P->p22 * lampp;
 49		cl = cos(lamtp);
 50		if (fabs(cl) < TOL)
 51		    lamtp -= TOL;
 52		fac = lampp - sin(lampp) * (cl < 0. ? -HALFPI : HALFPI);
 53		for (l = 50; l; --l) {
 54			lamt = lp.lam + P->p22 * sav;
 55			if (fabs(c = cos(lamt)) < TOL)
 56			    lamt -= TOL;
 57			xlam = (P->one_es * tanphi * P->sa + sin(lamt) * P->ca) / c;
 58			lamdp = atan(xlam) + fac;
 59			if (fabs(fabs(sav) - fabs(lamdp)) < TOL)
 60			    break;
 61			sav = lamdp;
 62		}
 63		if (!l || ++nn >= 3 || (lamdp > P->rlm && lamdp < P->rlm2))
 64			break;
 65		if (lamdp <= P->rlm)
 66		    lampp = TWOPI_HALFPI;
 67		else if (lamdp >= P->rlm2)
 68		    lampp = HALFPI;
 69	}
 70	if (l) {
 71		sp = sin(lp.phi);
 72		phidp = aasin((P->one_es * P->ca * sp - P->sa * cos(lp.phi) * 
 73			sin(lamt)) / sqrt(1. - P->es * sp * sp));
 74		tanph = log(tan(FORTPI + .5 * phidp));
 75		sd = sin(lamdp);
 76		sdsq = sd * sd;
 77		s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)
 78			 / ((1. + P->w * sdsq) * (1. + P->q * sdsq)));
 79		d = sqrt(P->xj * P->xj + s * s);
 80		xy.x = P->b * lamdp + P->a2 * sin(2. * lamdp) + P->a4 *
 81			sin(lamdp * 4.) - tanph * s / d;
 82		xy.y = P->c1 * sd + P->c3 * sin(lamdp * 3.) + tanph * P->xj / d;
 83	} else
 84		xy.x = xy.y = HUGE_VAL;
 85	return xy;
 86}
 87INVERSE(e_inverse); /* ellipsoid */
 88    int nn;
 89    double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
 90
 91	lamdp = xy.x / P->b;
 92	nn = 50;
 93	do {
 94		sav = lamdp;
 95		sd = sin(lamdp);
 96		sdsq = sd * sd;
 97		s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)
 98			 / ((1. + P->w * sdsq) * (1. + P->q * sdsq)));
 99		lamdp = xy.x + xy.y * s / P->xj - P->a2 * sin(
100			2. * lamdp) - P->a4 * sin(lamdp * 4.) - s / P->xj * (
101			P->c1 * sin(lamdp) + P->c3 * sin(lamdp * 3.));
102		lamdp /= P->b;
103	} while (fabs(lamdp - sav) >= TOL && --nn);
104	sl = sin(lamdp);
105	fac = exp(sqrt(1. + s * s / P->xj / P->xj) * (xy.y - 
106		P->c1 * sl - P->c3 * sin(lamdp * 3.)));
107	phidp = 2. * (atan(fac) - FORTPI);
108	dd = sl * sl;
109	if (fabs(cos(lamdp)) < TOL)
110	    lamdp -= TOL;
111	spp = sin(phidp);
112	sppsq = spp * spp;
113	lamt = atan(((1. - sppsq * P->rone_es) * tan(lamdp) * 
114		P->ca - spp * P->sa * sqrt((1. + P->q * dd) * (
115		1. - sppsq) - sppsq * P->u) / cos(lamdp)) / (1. - sppsq 
116		* (1. + P->u)));
117	sl = lamt >= 0. ? 1. : -1.;
118	scl = cos(lamdp) >= 0. ? 1. : -1;
119	lamt -= HALFPI * (1. - scl) * sl;
120	lp.lam = lamt - P->p22 * lamdp;
121	if (fabs(P->sa) < TOL)
122	    lp.phi = aasin(spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
123	else
124		lp.phi = atan((tan(lamdp) * cos(lamt) - P->ca * sin(lamt)) /
125			(P->one_es * P->sa));
126	return lp;
127}
128FREEUP; if (P) pj_dalloc(P); }
129ENTRY0(lsat)
130    int land, path;
131    double lam, alf, esc, ess;
132
133	land = pj_param(P->params, "ilsat").i;
134	if (land <= 0 || land > 5) E_ERROR(-28);
135	path = pj_param(P->params, "ipath").i;
136	if (path <= 0 || path > (land <= 3 ? 251 : 233)) E_ERROR(-29);
137	if (land <= 3) {
138		P->lam0 = DEG_TO_RAD * 128.87 - TWOPI / 251. * path;
139	    P->p22 = 103.2669323;
140	    alf = DEG_TO_RAD * 99.092;
141	} else {
142		P->lam0 = DEG_TO_RAD * 129.3 - TWOPI / 233. * path;
143	    P->p22 = 98.8841202;
144	    alf = DEG_TO_RAD * 98.2;
145	}
146	P->p22 /= 1440.;
147	P->sa = sin(alf);
148	P->ca = cos(alf);
149	if (fabs(P->ca) < 1e-9)
150	    P->ca = 1e-9;
151	esc = P->es * P->ca * P->ca;
152	ess = P->es * P->sa * P->sa;
153	P->w = (1. - esc) * P->rone_es;
154	P->w = P->w * P->w - 1.;
155	P->q = ess * P->rone_es;
156	P->t = ess * (2. - P->es) * P->rone_es * P->rone_es;
157	P->u = esc * P->rone_es;
158	P->xj = P->one_es * P->one_es * P->one_es;
159	P->rlm = PI * (1. / 248. + .5161290322580645);
160	P->rlm2 = P->rlm + TWOPI;
161    P->a2 = P->a4 = P->b = P->c1 = P->c3 = 0.;
162	seraz0(0., 1., P);
163	for (lam = 9.; lam <= 81.0001; lam += 18.)
164	    seraz0(lam, 4., P);
165	for (lam = 18; lam <= 72.0001; lam += 18.)
166	    seraz0(lam, 2., P);
167	seraz0(90., 1., P);
168	P->a2 /= 30.;
169	P->a4 /= 60.;
170	P->b /= 30.;
171	P->c1 /= 15.;
172	P->c3 /= 45.;
173	P->inv = e_inverse; P->fwd = e_forward;
174ENDENTRY(P)