.nfsee02

seismic software,very useful

static const char ID[] = "Source file: \
   $RCSfile: supressure.c,v $ \
   $Revision: 1.29 $ \
   $Date: 2005/02/15 16:58:03 $ ";


char *sdoc =
"SUSEDRATE - 2D arbline sedimentation rate computation \n"
"\n"
"susedrate	 <in.su hrz= [hrz=] age= [age=] >out.su \n"
"\n"
"Required Parameters:\n"
"hrz=           2D Landmark depth horizon filenames  \n"
"               (ordered from shallowest to deepest) \n"
"age=           ages of horizon (MA)\n"
"\n"
"\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"
"rho_min=       minimum value for density\n"
"rho_max=       maximum value for density\n"
"k_rho=         density trend exponent\n"
"wd=            name of water depth horizon grid\n"
"\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                    \\"
"\nwd=water_depth.hrz              \\"
"\nsalt=salt.ply                   \\"
"\nhrz=hz1.hrz   age=5.6           \\"
"\nhrz=hz2.hrz   age=7.5           \\"
"\nhrz=hz3.hrz   age=8.7           \\"
"\nhrz=hz4.hrz   age=11.7          \\"
"\nhrz=hz5.hrz   age=12.6          \\"
"\nhrz=hz6.hrz   age=12.9          \\"
"\nhrz=hz7.hrz   age=13.6          \\"
"\nhrz=hz8.hrz   age=14.6          \\"
"\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 );

segy tr;
segy out;

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

   /*-------------*/
   /* water depth */
   /*-------------*/

   char*   wdname;
   FILE*   wdfp;
   float*  wd_z  = 0;

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

   char*   hrzname;
   FILE*   hrzfp;
   int     nhrz = 0;
   float** hrz_z = 0;
   float   line;
   float   xline;

   /*------*/
   /* 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;
   int     nlines;

   /*------*/
   /* ages */
   /*------*/

   float* age   = 0;
   int    nage = 0;

   /* decompaction */

   float  rho;
   float  rho_min;
   float  rho_max;
   float  k_rho;

   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;

   char* ptr;

   int debug = 0;

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

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

   /*-----------------------------*/
   /* 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( getparstring( "wd" ,&wdname )){
      i++;
   }

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

   /*----------------------*/
   /* check parameter list */
   /*----------------------*/

   if( (nhrz = countparname("hrz")) == 1 ){
      err(" at least 2 horizons are needed \n");
   }

   if( (nage = countparname("age")) != nhrz ){
      err(" %d ages do not match %d horizons \n",nage ,nhrz );
   }

   /*-------------------------*/
   /* read and check the ages */
   /*-------------------------*/

   age = (float*) emalloc( nage*sizeof(float) );

   for( i=0; i<nage; i++ ){
      getnparfloat( i+1 ,"age" ,&(age[i]) );
   }

   for( i=1; i<nage; i++ ){

      if( age[i] <= age[i-1] ){
         fprintf( stderr ,"Age[%d] < age[%d]\n" ,i ,i-1 );
         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);
   }

/*--------------------------------------------------------------------*\
   Post the horizon depths to the list for each input trace
\*--------------------------------------------------------------------*/

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

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

      getnparstring(i+1,"hrz",&hrzname);
      hrzfp=efopen(hrzname,"r");
      fprintf( stderr ,"Processing %s\n" ,hrzname );

      while( !feof( hrzfp ) ){
         fgets( buf ,sizeof(buf) ,hrzfp );
         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[j][i] = depth;
            }
            ptr += 240+tr.ns*sizeof(float);
         }
      }
   }

/*--------------------------------------------------------------------*\
   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.
\*--------------------------------------------------------------------*/

   getparstring( "salt" ,&salt );
   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) );