view release on metacpan or  search on metacpan

gctpc/robinv.c  view on Meta::CPAN

         pr[13]=0.6769; 
         xlr[13]=0.835;  
         pr[14]=0.7346;
         xlr[14]=0.7986;
         pr[15]=0.7903;  
         xlr[15]=0.7597;  
         pr[16]=0.8435;
         xlr[16]=0.7186;
         pr[17]=0.8936; 
         xlr[17]=0.6732; 
         pr[18]=0.9394; 
         xlr[18]=0.6213;
         pr[19]=0.9761;  
         xlr[19]=0.5722;  
         pr[20]=1.0;  
         xlr[20]=0.5322;

         for (i = 0; i < 21; i++)
            xlr[i] *= 0.9858;   

/* Report parameters to the user
  -----------------------------*/
ptitle("ROBINSON"); 
radius(r);
cenlon(center_long);
offsetp(false_easting,false_northing);
return(OK);
}

/* Robinson inverse equations--mapping x,y to lat/long
  ------------------------------------------------------------*/
long robinv(x, y, lon, lat)
double x;			/* (O) X projection coordinate */
double y;			/* (O) Y projection coordinate */
double *lon;			/* (I) Longitude */
double *lat;			/* (I) Latitude */

{
double yy;
double p2;
double u,v,t,c;
double phid;
double temp;
double y1;
long ip1;
long i;


/* Inverse equations
  -----------------*/
x -= false_easting;
y -= false_northing;

yy = 2.0 * y / PI / R;
phid = yy * 90.0;
p2 = fabs(phid / 5.0);
ip1 = (long) (p2 - EPSLN);
if (ip1 == 0)
   ip1 = 1;

/* Stirling's interpolation formula as used in forward transformation is 
   reversed for first estimation of LAT. from rectangular coordinates. LAT.
   is then adjusted by iteration until use of forward series provides correct 
   value of Y within tolerance.
---------------------------------------------------------------------------*/
for (i = 0;;)
   {
   u = pr[ip1 + 3] - pr[ip1 + 1];
   v = pr[ip1 + 3] - 2.0 * pr[ip1 + 2] + pr[ip1 + 1];
   t = 2.0 * (fabs(yy) - pr[ip1 + 2]) / u;
   c = v / u;
   p2 = t * (1.0 - c * t * (1.0 - 2.0 * c * t));

   if ((p2 >= 0.0) || (ip1 == 1))
      {
      if (y >= 0)
         phid = (p2 + (double) ip1 ) * 5.0;
      else
         phid = -(p2 + (double) ip1 ) * 5.0;
      
      do
        {
        p2 = fabs(phid / 5.0);
        ip1 = (long) (p2 - EPSLN);
        p2 -= (double) ip1;
 
        if (y >= 0)
           y1 = R * (pr[ip1 +2] + p2 *(pr[ip1 + 3] - pr[ip1 +1]) / 2.0 + p2 
                  * p2 * (pr[ip1 + 3] - 2.0 * pr[ip1 + 2] + pr[ip1 + 1])/2.0) 
                  * PI / 2.0; 
        else
           y1 = -R * (pr[ip1 +2] + p2 *(pr[ip1 + 3] - pr[ip1 +1]) / 2.0 + p2 
                   * p2 * (pr[ip1 + 3] - 2.0 * pr[ip1 + 2] + pr[ip1 + 1])/2.0) 
                   * PI / 2.0; 
        phid += -180.0 * (y1 - y) / PI / R;
        i++;
        if (i > 75)
           {
           p_error("Too many iterations in inverse","robinv-conv");
           return(234);
           }
        }
      while (fabs(y1 - y) > .00001);
      break;
      }
   else
      {
      ip1 -= 1;
      if (ip1 < 0)
           {
           p_error("Too many iterations in inverse","robinv-conv");
           return(234);
           }
      }
   }
*lat  = phid * .01745329252;

/* calculate  LONG. using final LAT. with transposed forward Stirling's 
   interpolation formula.
---------------------------------------------------------------------*/
*lon = lon_center + x / R / (xlr[ip1 + 2] + p2 * (xlr[ip1 + 3] - xlr[ip1 + 1])
                      / 2.0 + p2 * p2 * (xlr[ip1 + 3] - 2.0 * xlr[ip1 + 2] + 
                      xlr[ip1 + 1]) / 2.0);
*lon = adjust_lon(*lon);

return(OK);
}