view release on metacpan or  search on metacpan

gctpc/tminv.c  view on Meta::CPAN

double r_maj;			/* major axis			*/
double r_min;			/* minor axis			*/
double scale_fact;		/* scale factor			*/
double center_lon;		/* center longitude		*/
double center_lat;		/* center latitude		*/
double false_east;		/* x offset in meters		*/
double false_north;		/* y offset in meters		*/
{
double temp;			/* temporary variable		*/

/* Place parameters in static storage for common use
  -------------------------------------------------*/
r_major = r_maj;
r_minor = r_min;
scale_factor = scale_fact;
lon_center = center_lon;
lat_origin = center_lat;
false_northing = false_north;
false_easting = false_east;

temp = r_minor / r_major;
es   = 1.0 - SQUARE(temp);
e    = sqrt(es);
e0   = e0fn(es);
e1   = e1fn(es);
e2   = e2fn(es);
e3   = e3fn(es);
ml0  = r_major * mlfn(e0, e1, e2, e3, lat_origin);
esp  = es / (1.0 - es);

if (es < .00001)
   ind = 1;

/* Report parameters to the user
  -----------------------------*/
ptitle("TRANSVERSE MERCATOR (TM)"); 
radius2(r_major, r_minor);
genrpt(scale_factor,"Scale Factor at C. Meridian:    ");
cenlonmer(lon_center);
origin(lat_origin);
offsetp(false_easting,false_northing);
return(OK);
}

/* Transverse Mercator inverse equations--mapping x,y to lat,long 
   Note:  The algorithm for UTM is exactly the same as TM and therefore
          if a change is implemented, also make the change to UTMINV.c
  --------------------------------------------------------------*/
long tminv(x, y, lon, lat)
double x;		/* (I) X projection coordinate 			*/
double y;		/* (I) Y projection coordinate 			*/
double *lon;		/* (O) Longitude 				*/
double *lat;		/* (O) Latitude 				*/
{
double con,phi;		/* temporary angles				*/
double delta_phi;	/* difference between longitudes		*/
long i;			/* counter variable				*/
double sin_phi, cos_phi, tan_phi;	/* sin cos and tangent values	*/
double c, cs, t, ts, n, r, d, ds;	/* temporary variables		*/
double f, h, g, temp;			/* temporary variables		*/
long max_iter = 6;			/* maximun number of iterations	*/

/* fortran code for spherical form 
--------------------------------*/
if (ind != 0)
   {
   f = exp(x/(r_major * scale_factor));
   g = .5 * (f - 1/f);
   temp = lat_origin + y/(r_major * scale_factor);
   h = cos(temp);
   con = sqrt((1.0 - h * h)/(1.0 + g * g));
   *lat = asinz(con);
   if (temp < 0)
     *lat = -*lat;
   if ((g == 0) && (h == 0))
     {
     *lon = lon_center;
     return(OK);
     }
   else
     {
     *lon = adjust_lon(atan2(g,h) + lon_center);
     return(OK);
     }
   }

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

con = (ml0 + y / scale_factor) / r_major;
phi = con;
for (i=0;;i++)
   {
   delta_phi = ((con + e1 * sin(2.0*phi) - e2 * sin(4.0*phi) + e3 * sin(6.0*phi))
               / e0) - phi;
/*
   delta_phi = ((con + e1 * sin(2.0*phi) - e2 * sin(4.0*phi)) / e0) - phi;
*/
   phi += delta_phi;
   if (fabs(delta_phi) <= EPSLN) break;
   if (i >= max_iter) 
      { 
      p_error("Latitude failed to converge","TM-INVERSE"); 
      return(95);
      }
   }
if (fabs(phi) < HALF_PI)
   {
   sincos(phi, &sin_phi, &cos_phi);
   tan_phi = tan(phi);
   c    = esp * SQUARE(cos_phi);
   cs   = SQUARE(c);
   t    = SQUARE(tan_phi);
   ts   = SQUARE(t);
   con  = 1.0 - es * SQUARE(sin_phi); 
   n    = r_major / sqrt(con);
   r    = n * (1.0 - es) / con;
   d    = x / (n * scale_factor);
   ds   = SQUARE(d);