/*
** libproj -- library of cartographic projections
**
** Copyright (c) 2003, 2006   Gerald I. Evenden
*/
static const char
LIBPROJ_ID[] = "Id";
/*
** Permission is hereby granted, free of charge, to any person obtaining
** a copy of this software and associated documentation files (the
** "Software"), to deal in the Software without restriction, including
** without limitation the rights to use, copy, modify, merge, publish,
** distribute, sublicense, and/or sell copies of the Software, and to
** permit persons to whom the Software is furnished to do so, subject to
** the following conditions:
**
** The above copyright notice and this permission notice shall be
** included in all copies or substantial portions of the Software.
**
** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
/* based upon Snyder and Linck, USGS-NMD */
#define PROJ_PARMS__ \
    double a2, a4, b, c1, c3; \
    double q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
#define PROJ_LIB__
#include  <lib_proj.h>
PROJ_HEAD(lsat, "Space oblique for LANDSAT")
  "\n\tCyl, Sph&Ell\n\tlsat= path=";
#define TOL 1e-7
#define PI_HALFPI 4.71238898038468985766
#define TWOPI_HALFPI 7.85398163397448309610
  static void
seraz0(double lam, double mult, PROJ *P) {
    double sdsq, h, s, fc, sd, sq, d__1;

    lam *= DEG_TO_RAD;
    sd = sin(lam);
    sdsq = sd * sd;
    s = P->p22 * P->sa * cos(lam) * sqrt((1. + P->t * sdsq) / ((
      1. + P->w * sdsq) * (1. + P->q * sdsq)));
    d__1 = 1. + P->q * sdsq;
    h = sqrt((1. + P->q * sdsq) / (1. + P->w * sdsq)) * ((1. + 
      P->w * sdsq) / (d__1 * d__1) - P->p22 * P->ca);
    sq = sqrt(P->xj * P->xj + s * s);
    P->b += fc = mult * (h * P->xj - s * s) / sq;
    P->a2 += fc * cos(lam + lam);
    P->a4 += fc * cos(lam * 4.);
    fc = mult * s * (h + P->xj) / sq;
    P->c1 += fc * cos(lam);
    P->c3 += fc * cos(lam * 3.);
}
FORWARD(e_forward); /* ellipsoid */
    int l, nn;
    double lamt = 0., xlam, sdsq, c, d, s, lamdp = 0., phidp,
    lampp, tanph, lamtp, cl, sd, sp, fac, sav, tanphi;

  if (lp.phi > HALFPI)
      lp.phi = HALFPI;
  else if (lp.phi < -HALFPI)
      lp.phi = -HALFPI;
  lampp = lp.phi >= 0. ? HALFPI : PI_HALFPI;
  tanphi = tan(lp.phi);
  for (nn = 0;;) {
    sav = lampp;
    lamtp = lp.lam + P->p22 * lampp;
    cl = cos(lamtp);
    if (fabs(cl) < TOL)
        lamtp -= TOL;
    fac = lampp - sin(lampp) * (cl < 0. ? -HALFPI : HALFPI);
    for (l = 50; l; --l) {
      lamt = lp.lam + P->p22 * sav;
      if (fabs(c = cos(lamt)) < TOL)
          lamt -= TOL;
      xlam = (P->one_es * tanphi * P->sa + sin(lamt) * P->ca) / c;
      lamdp = atan(xlam) + fac;
      if (fabs(fabs(sav) - fabs(lamdp)) < TOL)
          break;
      sav = lamdp;
    }
    if (!l || ++nn >= 3 || (lamdp > P->rlm && lamdp < P->rlm2))
      break;
    if (lamdp <= P->rlm)
        lampp = TWOPI_HALFPI;
    else if (lamdp >= P->rlm2)
        lampp = HALFPI;
  }
  if (l) {
    sp = sin(lp.phi);
    phidp = proj_asin((P->one_es * P->ca * sp - P->sa * cos(lp.phi) * 
      sin(lamt)) / sqrt(1. - P->es * sp * sp));
    tanph = log(tan(FORTPI + .5 * phidp));
    sd = sin(lamdp);
    sdsq = sd * sd;
    s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)
       / ((1. + P->w * sdsq) * (1. + P->q * sdsq)));
    d = sqrt(P->xj * P->xj + s * s);
    xy.x = P->b * lamdp + P->a2 * sin(2. * lamdp) + P->a4 *
      sin(lamdp * 4.) - tanph * s / d;
    xy.y = P->c1 * sd + P->c3 * sin(lamdp * 3.) + tanph * P->xj / d;
  } else
    xy.x = xy.y = HUGE_VAL;
  return xy;
}
INVERSE(e_inverse); /* ellipsoid */
    int nn;
    double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;

  lamdp = xy.x / P->b;
  nn = 50;
  do {
    sav = lamdp;
    sd = sin(lamdp);
    sdsq = sd * sd;
    s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)
       / ((1. + P->w * sdsq) * (1. + P->q * sdsq)));
    lamdp = xy.x + xy.y * s / P->xj - P->a2 * sin(
      2. * lamdp) - P->a4 * sin(lamdp * 4.) - s / P->xj * (
      P->c1 * sin(lamdp) + P->c3 * sin(lamdp * 3.));
    lamdp /= P->b;
  } while (fabs(lamdp - sav) >= TOL && --nn);
  sl = sin(lamdp);
  fac = exp(sqrt(1. + s * s / P->xj / P->xj) * (xy.y - 
    P->c1 * sl - P->c3 * sin(lamdp * 3.)));
  phidp = 2. * (atan(fac) - FORTPI);
  dd = sl * sl;
  if (fabs(cos(lamdp)) < TOL)
      lamdp -= TOL;
  spp = sin(phidp);
  sppsq = spp * spp;
  lamt = atan(((1. - sppsq * P->rone_es) * tan(lamdp) * 
    P->ca - spp * P->sa * sqrt((1. + P->q * dd) * (
    1. - sppsq) - sppsq * P->u) / cos(lamdp)) / (1. - sppsq 
    * (1. + P->u)));
  sl = lamt >= 0. ? 1. : -1.;
  scl = cos(lamdp) >= 0. ? 1. : -1;
  lamt -= HALFPI * (1. - scl) * sl;
  lp.lam = lamt - P->p22 * lamdp;
  if (fabs(P->sa) < TOL)
      lp.phi = proj_asin(spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
  else
    lp.phi = atan((tan(lamdp) * cos(lamt) - P->ca * sin(lamt)) /
      (P->one_es * P->sa));
  return lp;
}
FREEUP; if (P) free(P); }
ENTRY0(lsat)
    int land, path;
    double lam, alf, esc, ess;

  land = proj_param(P->params, "ilsat").i;
  if (land <= 0 || land > 5) E_ERROR(-28);
  path = proj_param(P->params, "ipath").i;
  if (path <= 0 || path > (land <= 3 ? 251 : 233)) E_ERROR(-29);
  if (land <= 3) {
    P->lam0 = DEG_TO_RAD * 128.87 - TWOPI / 251. * path;
      P->p22 = 103.2669323;
      alf = DEG_TO_RAD * 99.092;
  } else {
    P->lam0 = DEG_TO_RAD * 129.3 - TWOPI / 233. * path;
      P->p22 = 98.8841202;
      alf = DEG_TO_RAD * 98.2;
  }
  P->p22 /= 1440.;
  P->sa = sin(alf);
  P->ca = cos(alf);
  if (fabs(P->ca) < 1e-9)
      P->ca = 1e-9;
  esc = P->es * P->ca * P->ca;
  ess = P->es * P->sa * P->sa;
  P->w = (1. - esc) * P->rone_es;
  P->w = P->w * P->w - 1.;
  P->q = ess * P->rone_es;
  P->t = ess * (2. - P->es) * P->rone_es * P->rone_es;
  P->u = esc * P->rone_es;
  P->xj = P->one_es * P->one_es * P->one_es;
  P->rlm = PI * (1. / 248. + .5161290322580645);
  P->rlm2 = P->rlm + TWOPI;
    P->a2 = P->a4 = P->b = P->c1 = P->c3 = 0.;
  seraz0(0., 1., P);
  for (lam = 9.; lam <= 81.0001; lam += 18.)
      seraz0(lam, 4., P);
  for (lam = 18; lam <= 72.0001; lam += 18.)
      seraz0(lam, 2., P);
  seraz0(90., 1., P);
  P->a2 /= 30.;
  P->a4 /= 60.;
  P->b /= 30.;
  P->c1 /= 15.;
  P->c3 /= 45.;
  P->inv = e_inverse; P->fwd = e_forward;
ENDENTRY(P)
/*
** Log: proj_lsat.c
** Revision 1.1  2008-11-07 16:41:14  jeff
** ENH: Adding a 2D geoview. Adding the geographic projection library libproj4
** to Utilities. Updating the architecture of the geospatial views. All
** multi-resolution sources are now subclasses of vtkGeoSource. Each source
** has its own worker thread for fetching refined images or geometry.
** On the 3D side, vtkGeoGlobeSource is an appropriate source for vtkGeoTerrain,
** and vtkGeoAlignedImageSource is an appropriate source for
** vtkGeoAlignedImageRepresentation. On the 2D side, vtkGeoProjectionSource is an
** appropriate source for vtkGeoTerrain2D, and the image source is the same.
**
** Revision 3.1  2006/01/11 01:38:18  gie
** Initial
**
*/