view release on metacpan or  search on metacpan

gctpc/alconinv.c  view on Meta::CPAN

/* Inverse equations
  -----------------*/
x = (x - false_easting) / r_major;
y = (y - false_northing) / r_major;
xp = x;
yp = y;
nn = 0;

/* Use Knuth algorithm for summing complex terms, to convert Modified-
   Stereographic conformal to Oblique Stereographic coordinates.
--------------------------------------------------------------------*/
do
  {
  r = xp + xp;
  s = xp * xp + yp * yp;
  ar = acoef[n];
  ai = bcoef[n];
  br = acoef[n -1];
  bi = bcoef[n - 1];
  cr = (double) (n) * ar;
  ci = (double) (n) * ai;
  dr = (double) (n -1) * br;
  di = (double) (n -1) * bi;

  for (j = 2; j <= n; j++)
      {
      arn = br + r * ar;
      ain = bi + r * ai;
      if (j < n)
        {
        br = acoef[n -j] - s * ar;
        bi = bcoef[n - j] - s * ai;
        ar = arn;
        ai = ain;
        crn = dr  + r * cr;
        cin = di  + r * ci;
        dr = (double) (n - j) * acoef[n -j] - s * cr;
        di = (double) (n - j) * bcoef[n -j] - s * ci;
        cr = crn;
        ci = cin;
        }
      }
  br = -s * ar;
  bi = -s * ai;
  ar = arn;
  ai = ain;
  fxyr = xp * ar - yp * ai + br - x;
  fxyi = yp * ar + xp * ai + bi - y;
  fpxyr = xp * cr - yp * ci + dr;
  fpxyi = yp * cr + xp * ci + ci;
  den = fpxyr * fpxyr + fpxyi * fpxyi;
  dxp = -(fxyr * fpxyr + fxyi * fpxyi) / den;
  dyp = -(fxyi * fpxyr - fxyr * fpxyi) / den;
  xp = xp + dxp;
  yp = yp + dyp;
  ds = fabs(dxp) + fabs(dyp);
  nn++;
  if (nn > 20)
     {
     p_error("Too many iterations in inverse","alcon-inv");
     return(235);
     }
  }
while (ds > EPSLN);

/* convert Oblique Stereographic coordinates to LAT/LONG
------------------------------------------------------*/
rh = sqrt(xp * xp + yp * yp);
z = 2.0 * atan(rh / 2.0);
sincos(z,&sinz,&cosz);
*lon = lon_center;
if (fabs(rh) <= EPSLN)
   {
   *lat = lat_center;
   return(OK);
   }
chi = asinz(cosz * sin_p26 + (yp * sinz * cos_p26) / rh);
nn = 0;
phi = chi;
do
  {
  esphi = e * sin(phi);
  dphi = 2.0 * atan(tan((HALF_PI + chi) / 2.0) * 
         pow(((1.0 + esphi) / (1.0 - esphi)),(e / 2.0))) - HALF_PI - phi;
  phi += dphi;
  nn++;
  if (nn > 20)
     {
     p_error("Too many iterations in inverse","alcon-inv");
     return(236);
     }
  }
while(fabs(dphi) > EPSLN);

*lat = phi;
*lon = adjust_lon (lon_center + atan2((xp * sinz), (rh * cos_p26 * cosz - yp *
                   sin_p26 * sinz)));
     

return(OK);
}