susedrate.c

seismic software,very useful

static const char ID[] = "Source file: \
   $RCSfile: susedrate.c,v $ \
   $Revision: 1.12 $ \
   $Date: 2005/04/27 12:43:42 $ ";


char *sdoc =
"SUSEDRATE - 2D arbline sedimentation rate computation \n"
"\n"
"susedrate	 <in.su hznfile= horizon_name=age  >out.su \n"
"\n"
"Required Parameters:\n"
"\n"
"\nhznfile=            Emerald City format horizon file"
"\nhzns=               number of horizons in file"
"\nhorizon_name=age    names and ages of horizons"
"\n"
"Optional parameters:\n"
"\n"
"\nsalt=        name of file containing salt polygons in SeisWorks"
"\n             fault polygon format"
"\n"
"\nzscale=1e-3  scale factor to apply to trace header dt value"
"\n"
"Decompaction to depositional volume using an exponential compaction\n"
"trend of the form:\n"
"    rho = (rho_min*rho_max)/((exp(-k_rho*(z-wd))+rho_min)\n"
"\n"
"All the following parameters must be supplied for this option\n"
"\n"
"\nrho_min=       minimum value for density"
"\nrho_max=       maximum value for density"
"\nk_rho=         density trend exponent"
"\n"
"\n The water bottom horizon must be supplied with an age of 0"
"\n"
"\nNotes: "
"\n"
"\n The horizons must extend the entire length of the input arbline "
"\n which cannot cross itself."
"\n"
"\nExample:"
"\n"
"\n susedrate <input.su >output.su \\"
"\nrho_min=1.8                     \\"
"\nrho_max=2.65                    \\"
"\nk_rho=1.2e-4                    \\"
"\nsalt=salt.ply                   \\"
"\nhzns=5                          \\"
"\nhzn_wb=0                        \\"
"\nhzn_a1=5.5                      \\"
"\nhzn_a2=7.5                      \\"
"\nhzn_a3=9.5                      \\"
"\nhzn_a4=15.5                     \\"
"\n"
"\n"
"\n";


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "su.h"
#include "par.h"

int InPolygon( float X0 ,float Y0 ,float* X ,float* Y ,int n );
int agecmp( const void* e1 ,const void* e2 );

segy tr;
segy out;

typedef struct {
   float* depth;
   float  age;
}HRZN;

#define MAX_HRZ 100

int main( int argc ,char* argv[] ){

   /*---------*/
   /* horizon */
   /*---------*/

   char    hrzname[1024];
   char*   hznfile;
   FILE*   hrzfp;
   int     nhrz = 0;
   float** hrz_z = 0;
   float   line;
   float   xline;

   HRZN    hrz[MAX_HRZ];

   /*------*/
   /* salt */
   /*------*/

   char*   salt;
   FILE*   saltfp;

   float   idx;
   float   first_x;
   float   first_y;
   float   first_z;

   /* polygon buffers */

   float*  saltpoly_i;
   float*  saltpoly_m;
   float*  saltpoly_x;
   float*  saltpoly_y;
   float*  saltpoly_z;

   /* input pointers */

   float*  iptr;
   float*  mptr;
   float*  xptr;
   float*  yptr;
   float*  zptr;

   /* polygon arrays */

   float** salt_x;
   float** salt_y;
   float** salt_z;
   float** salt_idx;
   float** salt_min;
   int*    salt_n;

   /* polygon counts */

   int     npolygons = 0;
   int     nlines;

   /* decompaction */

   float  rho;
   float  rho_min=0;
   float  rho_max=0;
   float  k_rho=0;

   float factor = 1.0;

   /*---------------------*/
   /* seismic trace input */
   /*---------------------*/

   segy** trace;
   segy*  trace_buf;
   int    n_traces;

   /*---------------*/
   /* miscellaneous */
   /*---------------*/

   char   buf[1024];

   float  scaleco;
   float  zscale = 1.0e-3;
   float  depth;

   double sx;
   double sy;
   double dx;

   double x;
   double y;
   double z;

   int i;
   int j;
   int k;
   int p;
   int start;

   char* ptr;

   FILE* dbgfp=0;

   int debug = 0;

   /*----------------*/
   /* hook up getpar */
   /*----------------*/

   initargs(argc,argv);
   askdoc(0);

   if( !getparint( "hrzns" ,&nhrz )  ){
      err( "number of horizons must be specified" );

   }else if( nhrz > MAX_HRZ ){
      fprintf( stderr ,"maximum horizons is %d\n" ,MAX_HRZ );
      exit(-1);

   }

   /*-----------------------------*/
   /* check compaction parameters */
   /*-----------------------------*/

   i = 0;

   if( getparfloat( "rho_min" ,&rho_min )){
      i++;
   }

   if( getparfloat( "rho_max" ,&rho_max )){
      i++;
   }

   if( getparfloat( "k_rho" ,&k_rho )){
      i++;
   }

   if( i > 0 && i < 3 ){
      fprintf( stderr ,"All compaction arguments must be specified\n" );
      exit(-1);
   }

/*--------------------------------------------------------------------*\
   Read trace headers. This is complicated by the use of a staticly
   defined data member in the SU trace structure . First find out
   how much space we need. Then allocate the space and read all the
   traces using a pointer to the appropriate location.
\*--------------------------------------------------------------------*/

   if( !(n_traces=fgettra( stdin ,&tr ,0 )) ){
      err( "Unable to get first trace!" );
   }

   if( tr.scalco < 0 ){
      scaleco = -1.0 / tr.scalco;

   }else if( tr.scalco > 0 ){
      scaleco = tr.scalco;

   }else{
      scaleco = 1.0;
   }

   if( !(trace=calloc( n_traces ,sizeof(segy*) )) ){
      err( "Unable to allocate trace table" );
   }
   if( !(trace_buf=calloc( n_traces ,240+tr.ns*sizeof(float) )) ){
      err( "Unable to allocate trace buffer" );
   }

   ptr = (char*)trace_buf;
   for( i=0; i<n_traces; i++ ){
      fgettra( stdin ,(segy*)ptr ,i );
      trace[i] = (segy*)ptr;
      ptr += 240+tr.ns*sizeof(float);
   }

/*--------------------------------------------------------------------*\
   Allocate the horizon buffer space
\*--------------------------------------------------------------------*/

   hrz_z=(float**)ealloc2( n_traces ,nhrz+1 ,sizeof(float));

   memset( (char*)hrz ,0 ,sizeof(hrz) );

   if( !getparstring( "hznfile" ,&hznfile ) ){
      err( "missing hznfile=" );
   }

   hrzfp=efopen( hznfile ,"r" );
   fprintf( stderr ,"Processing %s\n" ,hznfile );

   i=0;

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

      if( buf[0] == '#' || buf[0] == '@' ){
         continue;

      }else if( !strcmp( buf ,"EOD\n" ) ){
         i++;

      }else if( strcspn( buf ," \t" ) == strlen( buf ) ){
         sscanf( buf ,"%s" ,hrzname );

         if( !getparfloat( hrzname ,&hrz[i].age ) ){
            fprintf( stderr ,"missing age for %s" ,hrzname );
            exit(-1);
         }
 
         hrz[i].depth=hrz_z[i];

      }else{
         sscanf( buf ,"%f%f%f%f%f" ,&line ,&xline ,&x ,&y ,&depth );

         ptr = (char*)trace_buf;

         for( j=0; j<n_traces; j++ ){
            if(   ((segy*)ptr)->fldr == (int)rint(line)
               && ((segy*)ptr)->cdp  == (int)rint(xline) ){
               hrz_z[i][j] = depth;
            }
            ptr += 240+tr.ns*sizeof(float);
         }
      }

   }

   nhrz=i;

   qsort( hrz ,nhrz ,sizeof(HRZN) ,agecmp );

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

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

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

            fprintf( dbgfp ,"%8d "    ,j               );
            fprintf( dbgfp ,"%8d "    ,i               );
            fprintf( dbgfp ,"%10.2f " ,hrz[i].age      );
            fprintf( dbgfp ,"%6.2f "  ,hrz[i].depth[j] );
            fprintf( dbgfp ,"%6.2f "  ,hrz_z[i][j]     );

            fprintf( dbgfp ,"\n" );
         }

         fprintf( dbgfp ,"\n" );

      }

      fclose( dbgfp );

   }

/*--------------------------------------------------------------------*\
   Read the salt polygons. This is complicated by the use of X-Y-Z
   coordinates for the polygons rather than the LINE-XLINE-Z used
   for the seismic trace loading.

   Determine approximate total number of lines in polygon file and
   allocate space for twice as many points. This wastes a bit of
   space, but not enough to worry about.  Note that we don't
   know how many points are in a polygon until we read the last point
   in the polygon.  The method of determining the file size is dumb,
   but so widely used in SU there's little point in doing better.
\*--------------------------------------------------------------------*/

   if( getparstring( "salt" ,&salt ) ){