/Proj4/PJ_lsat.c
C | 174 lines | 169 code | 4 blank | 1 comment | 30 complexity | 97e211823d8d045e8a79bbf58880a50c MD5 | raw file
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)