/*
** 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.
*/
#define PROJ_PARMS__ \
  double  Az, kRg, p0s, A, C, Ca, Cb, Cc, Cd; \
  int    rot;
#define PROJ_LIB__
#include  <lib_proj.h>
PROJ_HEAD(labrd, "Laborde") "\n\tCyl, Sph\n\tSpecial for Madagascar";
#define EPS  1.e-10
FORWARD(e_forward);
  double V1, V2, ps, sinps, cosps, sinps2, cosps2, I1, I2, I3, I4, I5, I6,
    x2, y2, t;

  V1 = P->A * log( tan(FORTPI + .5 * lp.phi) );
  t = P->e * sin(lp.phi);
  V2 = .5 * P->e * P->A * log ((1. + t)/(1. - t));
  ps = 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);
  I1 = ps - P->p0s;
  cosps = cos(ps);  cosps2 = cosps * cosps;
  sinps = sin(ps);  sinps2 = sinps * sinps;
  I4 = P->A * cosps;
  I2 = .5 * P->A * I4 * sinps;
  I3 = I2 * P->A * P->A * (5. * cosps2 - sinps2) / 12.;
  I6 = I4 * P->A * P->A;
  I5 = I6 * (cosps2 - sinps2) / 6.;
  I6 *= P->A * P->A *
    (5. * cosps2 * cosps2 + sinps2 * (sinps2 - 18. * cosps2)) / 120.;
  t = lp.lam * lp.lam;
  xy.x = P->kRg * lp.lam * (I4 + t * (I5 + t * I6));
  xy.y = P->kRg * (I1 + t * (I2 + t * I3));
  x2 = xy.x * xy.x;
  y2 = xy.y * xy.y;
  V1 = 3. * xy.x * y2 - xy.x * x2;
  V2 = xy.y * y2 - 3. * x2 * xy.y;
  xy.x += P->Ca * V1 + P->Cb * V2;
  xy.y += P->Ca * V2 - P->Cb * V1;
  return (xy);
}
INVERSE(e_inverse); /* ellipsoid & spheroid */
  double x2, y2, V1, V2, V3, V4, t, t2, ps, pe, tpe, s,
    I7, I8, I9, I10, I11, d, Re;
  int i;

  x2 = xy.x * xy.x;
  y2 = xy.y * xy.y;
  V1 = 3. * xy.x * y2 - xy.x * x2;
  V2 = xy.y * y2 - 3. * x2 * xy.y;
  V3 = xy.x * (5. * y2 * y2 + x2 * (-10. * y2 + x2 ));
  V4 = xy.y * (5. * x2 * x2 + y2 * (-10. * x2 + y2 ));
  xy.x += - P->Ca * V1 - P->Cb * V2 + P->Cc * V3 + P->Cd * V4;
  xy.y +=   P->Cb * V1 - P->Ca * V2 - P->Cd * V3 + P->Cc * V4;
  ps = P->p0s + xy.y / P->kRg;
  pe = ps + P->phi0 - P->p0s;
  for ( i = 20; i; --i) {
    V1 = P->A * log(tan(FORTPI + .5 * pe));
    tpe = P->e * sin(pe);
    V2 = .5 * P->e * P->A * log((1. + tpe)/(1. - tpe));
    t = ps - 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);
    pe += t;
    if (fabs(t) < EPS)
      break;
  }
/*
  if (!i) {
  } else {
  }
*/
  t = P->e * sin(pe);
  t = 1. - t * t;
  Re = P->one_es / ( t * sqrt(t) );
  t = tan(ps);
  t2 = t * t;
  s = P->kRg * P->kRg;
  d = Re * P->k0 * P->kRg;
  I7 = t / (2. * d);
  I8 = t * (5. + 3. * t2) / (24. * d * s);
  d = cos(ps) * P->kRg * P->A;
  I9 = 1. / d;
  d *= s;
  I10 = (1. + 2. * t2) / (6. * d);
  I11 = (5. + t2 * (28. + 24. * t2)) / (120. * d * s);
  x2 = xy.x * xy.x;
  lp.phi = pe + x2 * (-I7 + I8 * x2);
  lp.lam = xy.x * (I9 + x2 * (-I10 + x2 * I11));
  return (lp);
}
FREEUP; if (P) free(P); }
ENTRY0(labrd)
  double Az, sinp, R, N, t;

  P->rot  = proj_param(P->params, "bno_rot").i == 0;
  Az = proj_param(P->params, "razi").f;
  sinp = sin(P->phi0);
  t = 1. - P->es * sinp * sinp;
  N = 1. / sqrt(t);
  R = P->one_es * N / t;
  P->kRg = P->k0 * sqrt( N * R );
  P->p0s = atan( sqrt(R / N) * tan(P->phi0) );
  P->A = sinp / sin(P->p0s);
  t = P->e * sinp;
  P->C = .5 * P->e * P->A * log((1. + t)/(1. - t)) +
    - P->A * log( tan(FORTPI + .5 * P->phi0))
    + log( tan(FORTPI + .5 * P->p0s));
  t = Az + Az;
  P->Ca = (1. - cos(t)) * ( P->Cb = 1. / (12. * P->kRg * P->kRg) );
  P->Cb *= sin(t);
  P->Cc = 3. * (P->Ca * P->Ca - P->Cb * P->Cb);
  P->Cd = 6. * P->Ca * P->Cb;
  P->inv = e_inverse;
  P->fwd = e_forward;
ENDENTRY(P)
/*
** Log: proj_labrd.c
** Revision 3.1  2006/01/11 01:38:18  gie
** Initial
**
*/