susedrate.c

seismic software,very useful

      saltfp=efopen(salt,"r");
      fprintf( stderr ,"Processing %s\n" ,salt );

      fgets( buf ,sizeof(buf) ,saltfp );
      fseek( saltfp ,0 ,SEEK_END );
      nlines = ftell(saltfp) / strlen(buf);

      saltpoly_x  = (float*)ealloc1( 2*nlines ,sizeof(float) );
      saltpoly_y  = (float*)ealloc1( 2*nlines ,sizeof(float) );
      saltpoly_z  = (float*)ealloc1( 2*nlines ,sizeof(float) );
      saltpoly_i  = (float*)ealloc1( 2*nlines ,sizeof(float) );
      saltpoly_m  = (float*)ealloc1( 2*nlines ,sizeof(float) );

      salt_n  = (int*)  ealloc1(   nlines ,sizeof(int)   );

      salt_x   = (float**)ealloc1( nlines ,sizeof(float*) );
      salt_y   = (float**)ealloc1( nlines ,sizeof(float*) );
      salt_z   = (float**)ealloc1( nlines ,sizeof(float*) );
      salt_idx = (float**)ealloc1( nlines ,sizeof(float*) );
      salt_min = (float**)ealloc1( nlines ,sizeof(float*) );

      salt_x[0]   = saltpoly_x;
      salt_y[0]   = saltpoly_y;
      salt_z[0]   = saltpoly_z;
      salt_idx[0] = saltpoly_i;
      salt_min[0] = saltpoly_m;

      j = 0;
      k = 0;

      xptr = salt_x[0];
      yptr = salt_y[0];
      zptr = salt_z[0];

      iptr = salt_idx[0];
      mptr = salt_min[0];

      fseek( saltfp ,0 ,SEEK_SET );

      while( !feof( saltfp ) && fgets( buf ,sizeof(buf) ,saltfp ) ){

         sscanf( &(buf[0])  ,"%12lf" ,&x );
         sscanf( &(buf[12]) ,"%12lf" ,&y );
         sscanf( &(buf[24]) ,"%12lf" ,&z );

/*--------------------------------------------------------------------*\
         Check the value in card column 43 to see if we've reached the end
         of a polygon.  If so close the polygon by repeating the first 
         point.  
\*--------------------------------------------------------------------*/

         if( !(strncmp( &(buf[42]) ,"1" ,1 )) ){

            /*------------------*/
            /* start of polygon */
            /*------------------*/
          
            *xptr     = x;
            *yptr     = y;
            *zptr     = z;
            *mptr     = 1.0e30;

            first_x   = x;
            first_y   = y;
            first_z   = z;

            salt_x[j]   = xptr++;
            salt_y[j]   = yptr++;
            salt_z[j]   = zptr++;
            salt_idx[j] = iptr++;
            salt_min[j] = mptr++;

            k=1;

         }else if( !(strncmp( &(buf[42]) ,"3" ,1 )) ){

            /*----------------*/
            /* end of polygon */
            /*----------------*/

            *xptr++   = x;
            *yptr++   = y;
            *zptr++   = z;
            *mptr++   = 1.0e30;
             iptr++;

            *xptr++   = first_x;
            *yptr++   = first_y;
            *zptr++   = first_z;
            *mptr++   = 1.0e30;
             iptr++;

            salt_n[j] = k+2;
            j++;

         }else{

            /*------------------------------*/
            /* interior of polygon polyline */
            /*------------------------------*/

            *xptr++   = x;
            *yptr++   = y;
            *zptr++   = z;
            *mptr++   = 1.0e30;

            iptr++;
            k++;

         }

      }

      npolygons = j;

/*--------------------------------------------------------------------*\
      Determine the trace index to use for the X coordinate of the salt
      polygon vertices. The strategy is to find the minimum distance
      from the point to a trace.
\*--------------------------------------------------------------------*/

      ptr = (char*)trace_buf;

      for( i=0; i<n_traces; i++ ){

         sx = ((segy*)ptr)->sx*scaleco;      
         sy = ((segy*)ptr)->sy*scaleco;      

         for( j=0; j<npolygons; j++ ){

            for( k=0; k<salt_n[j]; k++ ){
           
               dx = sqrt( pow( (sx - salt_x[j][k]) ,2.0 )
                        + pow( (sy - salt_y[j][k]) ,2.0 ) );
               if( dx < salt_min[j][k] ){
                   salt_min[j][k] = dx;
                   salt_idx[j][k] = i;
               }

            }
         }

         ptr += 240+tr.ns*sizeof(float);

      }

/*--------------------------------------------------------------------*\
      dump out the salt polygons in X-Y and trace index coordinates
      so the coordinate transform can be checked.
\*--------------------------------------------------------------------*/

      if( (ptr=getenv( "POLY")) && (dbgfp=fopen( ptr ,"w")) ){

         for( i=0; i<npolygons; i++ ){

            for( j=0; j<salt_n[i]; j++ ){

               fprintf( dbgfp ,"%10.1f " ,salt_x[i][j] );
               fprintf( dbgfp ,"%10.1f " ,salt_y[i][j] );
               fprintf( dbgfp ,"%10.1f " ,salt_z[i][j] );
               fprintf( dbgfp ,"%10.1f " ,salt_min[i][j] );
               fprintf( dbgfp ,"%10.1f " ,salt_idx[i][j] );

               fprintf( dbgfp ,"%10.1f " 
                              ,trace[(int)salt_idx[i][j]]->sx*scaleco );
               fprintf( dbgfp ,"%10.1f " 
                              ,trace[(int)salt_idx[i][j]]->sy*scaleco );

               fprintf( dbgfp ,"\n" );

            }
            fprintf( dbgfp ,"\n" );
         }

         fclose( dbgfp );
      }

   }

/*--------------------------------------------------------------------*\
   Create an output attribute trace for each input trace header.  The
   output must be sampled exactly as the input was sampled so that
   SeisWorks can properly make the overlays.
\*--------------------------------------------------------------------*/

   ptr = (char*)trace_buf;

   for( j=0; j<n_traces; j++ ){

      idx = j;

      memcpy( (char*)&out ,ptr ,240 );
      memset( (char*)out.data ,0 ,out.ns*sizeof(float) );

      i = 1;

      start = (hrz[0].depth[j]-out.delrt)/(out.dt*zscale);

      if( start < 0 ){
         start = 0;
      }

      for( k=start; k<out.ns; k++ ){

         z = out.delrt+k*out.dt*zscale;

         while( i < nhrz-1 && z > hrz[i].depth[j] ){
            i++;
         }

         /*-------------------------------*/
         /* check the horizon depth order */
         /*-------------------------------*/

         if( hrz[i].depth[j] < hrz[i-1].depth[j] ){

            continue;
         }

         /*----------------------------------*/
         /* calculate the sedimentation rate */
         /*----------------------------------*/

         if( rho_min > 0.0 ){

            rho = (rho_max*rho_min) / ( (rho_max - rho_min)
                 * exp( -k_rho*(z-hrz[0].depth[j]) ) + rho_min );

            factor = rho /rho_min;

         }

         out.data[k] = factor * (hrz[i].depth[j]-hrz[i-1].depth[j])
                      / (hrz[i].age-hrz[i-1].age);

         /*--------------------------------*/
         /* apply the salt polygon masking */
         /*--------------------------------*/

         for( p=0; p<npolygons; p++ ){

            if( InPolygon( idx ,z ,salt_idx[p] ,salt_z[p] ,salt_n[p]) ){
               out.data[k] = 0.0;
            }
         }
      }

      /*------------------------*/
      /* write the output trace */
      /*------------------------*/

      fputtr( stdout ,&out );
      ptr += 240+tr.ns*sizeof(float);

   }

   return(0);

}

/*--------------------------------------------------------------------*\
   InPolygon() determines if a point is within a polygon denoted by
   a polyline with identical first and last points.

   The basic algorithm is from "Computational Geometry in C", by
   Joseph O'Rourke but has been completely rewritten to make it 
   suitable for industrial use.

   The function returns integer values indicating the location of
   the input point relative to the bounding polygon:

   0 - point is exterior to the polygon
   1 - point is a vertex of the polygon
   2 - point lies on a line segment
   3 - point is strictly interior to the polygon

   Reginald H. Beardsley                            [e-mail]
\*--------------------------------------------------------------------*/

#include <stdio.h>

int InPolygon( float X0 ,float Y0 ,float* X ,float* Y ,int n ){

   int i;            /* vertex          */
   int j;            /* adjacent vertex */

   int R = 0;        /* number of right crossings */
   int L = 0;        /* number of left crossings  */

   float x;

   for( i=0; i<n; i++ ){

      if( X[i] == X0 && Y[i] == Y0 ){
         return( 1 );

      }

      j = (i + n - 1) % n;

      if(   ((Y[i] > Y0) && (Y[j] <= Y0)) 
          ||((Y[j] > Y0) && (Y[i] <= Y0)) ){
         
         x = ((X[i]-X0) * (Y[j]-Y0) - (X[j]-X0) * (Y[i]-Y0))
            / ((Y[j]-Y0) - (Y[i]-Y0));

         if( x > 0 ){
            R++;
         }
      }

      if(   ((Y[i] > Y0) && (Y[j] <= Y0)) 
          ||((Y[j] > Y0) && (Y[i] <= Y0)) ){
         
         x = ((X[i]-X0) * (Y[j]-Y0) - (X[j]-X0) * (Y[i]-Y0))
            / ((Y[j]-Y0) - (Y[i]-Y0));

         if( x < 0 ){
            L++;
         }
      }

   }

   if( (R % 2) != (L % 2) ){
      return( 2 ); 
   }

   if( (R % 2) == 1 ){
      return( 3 );

   }else{
      return( 0 );

   }
}

int agecmp( const void* e1 ,const void* e2 ){

   if( ((HRZN*)e1)->age < ((HRZN*)e2)->age ){
      return( -1 );

   }else if( ((HRZN*)e1)->age > ((HRZN*)e2)->age ){
      return( 1 );

   }else{
      return( 0 );

   }
}