view release on metacpan or  search on metacpan

gctpc/somfor.c  view on Meta::CPAN

double lat;		/* (I) Latitude 		*/
double *x;		/* (O) X projection coordinate 	*/
double *y;		/* (O) Y projection coordinate 	*/
{
long n,i,l;
double delta_lon;
double rlm,tabs,tlam,xlam,c,xlamt,ab2,sc1,ab1,xlamp,sav;
double d,sdsq,sd,tanlg,xtan,tphi,dp,dd,ds,rlm2;
double scl,tlamp,conv,delta_lat,radlt,radln;
double temp;
char   errorbuf[80];

/* Forward equations
  -----------------*/
conv=1.e-7;
delta_lat=lat;
delta_lon= lon-lon_center;

/* Test for latitude and longitude approaching 90 degrees
   ----------------------------------------------------*/
if (delta_lat>1.570796) delta_lat=1.570796;
if (delta_lat<-1.570796) delta_lat= -1.570796;
radlt=delta_lat;
radln=delta_lon;
if(delta_lat>=0.0)tlamp=PI/2.0; 
if(start!= 0.0)tlamp=2.5*PI;
if(delta_lat<0.0) tlamp=1.5*PI;
n=0;

L230:  sav=tlamp;
       l=0;
       xlamp=radln+p21*tlamp;
       ab1=cos(xlamp);
       if(fabs(ab1)<conv) xlamp=xlamp-1.e-7;
       if(ab1>=0.0) scl=1.0;
       if(ab1<0.0) scl= -1.0;
       ab2=tlamp-(scl)*sin(tlamp)*HALF_PI;
L240:  xlamt=radln+p21*sav;
       c=cos(xlamt);
       if (fabs(c)<1.e-7) xlamt=xlamt-1.e-7;
       xlam=(((1.0-es)*tan(radlt)*sa)+sin(xlamt)*ca)/c;
       tlam=atan(xlam);
       tlam=tlam+ab2;
       tabs=fabs(sav)-fabs(tlam);
       if(fabs(tabs)<conv) goto L250;
       l=l+1;
       if (l > 50) goto L260;
       sav=tlam;
       goto L240;

/* Adjust for confusion at beginning and end of landsat orbits
  -----------------------------------------------------------*/
L250:  rlm=PI*LANDSAT_RATIO;
       rlm2=rlm+2.0*PI;
       n++;
       if(n>=3) goto L300;
       if(tlam>rlm&&tlam<rlm2) goto L300;
       if(tlam<rlm)tlamp=2.50*PI;
       if(tlam>=rlm2) tlamp=HALF_PI;
       goto L230;
L260:  sprintf(errorbuf,"50 iterations without conv\n");
       p_error(errorbuf,"som-forward");
       return(214);

/* tlam computed - now compute tphi
  --------------------------------*/
L300: ds=sin(tlam);
      dd=ds*ds;
      dp=sin(radlt);
      tphi=asin(((1.0-es)*ca*dp-sa*cos(radlt)*sin(xlamt))/sqrt(1.0-es*dp*dp));

/* compute x and y
  ---------------*/
xtan = (PI/4.0) + (tphi/2.0);
tanlg = log(tan(xtan));
sd=sin(tlam);
sdsq=sd*sd;
s=p21*sa*cos(tlam)*sqrt((1.0+t*sdsq)/((1.0+w*sdsq)*(1.0+q*sdsq)));
d=sqrt(xj*xj+s*s);
*x=b*tlam+a2*sin(2.0*tlam)+a4*sin(4.0*tlam)-tanlg*s/d;
*x = a* *x;
*y=c1*sd+c3*sin(3.0*tlam)+tanlg*xj/d;
*y = a* *y;

/* Negate x & swap x,y
  -------------------*/
temp=  *x;
*x= *y + false_easting;
*y=temp + false_northing;;
return(OK);
}
/* Series to calculate a,b,c coefficients to convert from transform
   latitude,longitude to Space Oblique Mercator (SOM) rectangular coordinates

   Mathematical analysis by John Snyder 6/82
  --------------------------------------------------------------------------*/
static double som_series(fb,fa2,fa4,fc1,fc3,dlam)
double *fb,*fa2,*fa4,*fc1,*fc3,*dlam;
{
double sd,sdsq,h,sq,fc;

*dlam= *dlam*0.0174532925;               /* Convert dlam to radians */
sd=sin(*dlam);
sdsq=sd*sd;
s=p21*sa*cos(*dlam)*sqrt((1.0+t*sdsq)/((1.0+w*sdsq)*(1.0+q*sdsq)));
h=sqrt((1.0+q*sdsq)/(1.0+w*sdsq))*(((1.0+w*sdsq)/((1.0+q*sdsq)*(1.0+
     q*sdsq)))-p21*ca);
sq=sqrt(xj*xj+s*s);
*fb=(h*xj-s*s)/sq;
*fa2= *fb*cos(2.0* *dlam);
*fa4= *fb*cos(4.0* *dlam);
fc=s*(h+xj)/sq;
*fc1=fc*cos(*dlam);
*fc3=fc*cos(3.0* *dlam);
}