get_isohyet.c

Go to the documentation of this file.

#include "prototypes_new.h"
 
float get_isohyet(int newx,int newy,int isom)
 
{
 
extern struct isoh *isoh;
extern struct dval dval;
int kk,jj,ix,iy;
float distance,value,dist1,dist2,dist,isohyet;
double conv=.0174,lat,lon,c,lonf;
double t1,t2,t3,t4,t5,lat1,lat2;
 
t1=.0001;
 
lonf= atan(((double)newy-dval.yo)/(t1))/conv-90.;
 
t1=(double)newx-dval.xo;
 
lon= atan(((double)newy-dval.yo)/(t1))/conv;
 
if(lon > 0)
          lon=-180+lon-2*lonf;
 
lon=lon+90+dval.lo;
 
c=dval.a*cos((lon-dval.lo-90)*conv)/((double)newx-dval.xo);
c=pow(c,2);
 
t1=1/(1+c);
t2=1-c;
t3=pow(2*t1,2);
t4=4*t2*t1;
t5=pow(t3-t4,.5)/2;
 
if(fabs(-t1+t5) > 1.0)
          lat1=-9999;
else
          lat1=asin(-t1+t5)/conv;
 
if(fabs(-t1-t5) > 1.0)
          lat2=-9999;
else
          lat2=asin(-t1-t5)/conv;
 
if(lat1 > 0 && lat1 < 90)
                 lat=lat1;
 
else
                 lat=lat2;
 
iy=(isoh->max_lat-lat)/isoh->delta_lat;
ix=(isoh->max_lon-lon)/isoh->delta_lon;
 
if(ix < 0 || iy < 0 || ix > isoh->maxi-1 || iy> isoh->maxj-1) 
                            return(-1.0);
 
distance=0.0;
value=0.0;
 
for(kk=ix;kk<=ix+1;kk++) {
      
for(jj=iy;jj<=iy+1;jj++) {
 
           dist1=(float)newx-(float)isoh->coord[kk][jj].x;
           dist2=(float)newy-(float)isoh->coord[kk][jj].y;
 
           dist=pow(dist1,2)+pow(dist2,2);
 
           if(dist==0.0)
                      dist=1.0;
 
           dist=1/dist;
 
           value=value+dist*(float)isoh->value[isom][kk][jj];
           distance=distance+dist;
 
         }
             }
           
isohyet=value/distance;
 
return(isohyet);
 
 
}