PageRenderTime 34ms CodeModel.GetById 21ms RepoModel.GetById 0ms app.codeStats 1ms

/Proj4/PJ_lsat.c

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