view release on metacpan or  search on metacpan

gctpc/robfor.c  view on Meta::CPAN

         pr[7]=0.31;     
         xlr[7]=0.973;    
         pr[8]=0.372;      
         xlr[8]=0.96;
         pr[9]=0.434; 
         xlr[9]=0.9427;
         pr[10]=0.4958;
         xlr[10]=0.9216;
         pr[11]=0.5571;  
         xlr[11]=0.8962;  
         pr[12]=0.6176;
         xlr[12]=0.8679;
         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 forward equations--mapping lat,long to x,y
  ------------------------------------------------------------*/
long robfor(lon, lat, x, y)
double lon;			/* (I) Longitude */
double lat;			/* (I) Latitude */
double *x;			/* (O) X projection coordinate */
double *y;			/* (O) Y projection coordinate */

{
double dlon;
double p2;
long ip1;

/* Forward equations
  -----------------*/
dlon = adjust_lon(lon - lon_center);
p2 = fabs(lat / 5.0 / .01745329252);
ip1 = (long) (p2 - EPSLN);

/* Stirling's interpolation formula (using 2nd Diff.)
---------------------------------------------------*/
p2 -= (double) ip1;
*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) * 
          dlon + false_easting;

if (lat >= 0)
   *y = 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 +
	    false_northing;
else
   *y = -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 +
	     false_northing;

return(OK);
}