surnmo.c

seismic software,very useful

#include "su.h"
#include "cwp.h"
#include "segy.h"
#include "ghdr.h"
#include "gridhd.h"
#include "grid.h"

/*********************** self documentation ******************************/
string sdoc =
"                                                                       \n"
" NAME                                                                  \n"
"       RNmo -- Residual NMO transform                                  \n"
"               Input is an NMO corrected gather,                       \n"
"               Output is a Residual NMO corrected gather.              \n"
"                                                                       \n"
" SYNOPSIS                                                              \n"
"       RNmo < stdin > stdout [input parameters]                        \n"
"                                                                       \n"
"                                                                       \n"
" PARAMETERS:                                                           \n"
"                                                                       \n"
"  rmo_option=1    options for running RNmo                             \n"
"                  1 - scan and output rmo corrections in HANDVEL format \n"
"                  2 - scan and output rmo corrections in HANDVEL format \n"
"                       and rmo corrected gathers (for testing purpose)  \n" 
"                  3 - apply rmo corrections (input grid in vgrid format)\n"
"                       and output rmo corrected gathers                 \n"
"                                                                       \n"
"  amin=-150       negative max moveout time/depth (ms or meters)       \n"
"  amax=150        positive max moveout time/depth (ms or meters)       \n"
"  na=30           number of moveouts to try from a=amin to a=amax      \n"
"  ntri=25         halfwidth of triangle smoothing filter (in points)   \n"
"  sscale=1        if sscale=0, don't apply NMO stretch scaling         \n"
"                                                                       \n"
"  tracekey=tracr  segy key word defining trace number inline (vgrid)   \n"
"  linekey=fldr    segy key word defining line number (vgrid)           \n"
"  gatherkey=cdp   segy key word defining gather assemble               \n"
"                                                                       \n"
"  f_rmo=rmo.hvel  file for residual moveout corrections at maximum     \n"
"                   offset in HANDVEL format                            \n" 
"  f_semb=semb.su  file for optimal semblance traces                    \n"
"  maxoffset=      max. offset (in meters or feet) for computing        \n"
"                   residual moveout correction                         \n" 
"  nrmo=40         max. number of residual moveout corrections of each  \n"
"                   gather for output                                   \n" 
"                                                                       \n"
"                                                                       \n"
"  rmogrid=      file name of residual moveout correction grid          \n"
 "                (stored as (t,x,y) order with grid header             \n"
"                                                                       \n"
"  ocdp2=from_rmogrid   inline trace number of first rmo grid trace     \n"
"  dcdp2=from_rmogrid   inline trace number increment of rmo grid       \n"
"  oline3=from_rmogrid  line number of first rmo grid trace             \n"
"  dline3=from_rmogrid  line number increment of rmo grid               \n"
"                                                                       \n"
"  ntrgrid=from_rmogrid no. of time or depth samples per trace in rmogrid\n"
"  dtrgrid=from_rmogrid time or depth interval (in ms or m) of rmogrid  \n"
"  ftrgrid=from_rmogrid minimum time or depth (in ms or m) of rmogrid   \n"
"                                                                       \n"
"                                                                       \n"
" Keyword Residual NMO                                                  \n"
"                                                                       \n"
"                                                                       \n"
" AUTHOR:                                                               \n"
"       Wenlong Xu                                                      \n"
"                                                                       \n"
" DATE: Oct.9, 2000                                                     \n"
"                                                                       \n"
" REFERENCE:                                                            \n"
"                                                                       \n"
"\n";
/**************** end self doc *******************************************/


segychdr ch ;
segybhdr bh ;
segytrace tr, tr_out ;


void rnmo(float **cmp, int nt, int nx, float t0, float dt, float amin,
          float amax, int na, int ntri, int sscale, float *xx, float **flat,
          float *aaopt, float *sembmax, int rmo_option, int nrmo, int *nrmo1, 
          float *tnout, float *rnout) ;

void rnmo_apply(float **cmp, int nt, int nx, float t0, float dt, int sscale, 
          float *xx, float **flat, float *aaopt ) ;

/* Note: copied from supstack */
void readv(FILE *vgfp, int ntvgrid, int nx, int ny, int ix, int iy,
    float dtvgrid, float ftvgrid, int nt, float dt, float ft, float *vel);

main(int argc, char **argv)
{
    FILE   *rmofp ;
    FILE   *rmogrdfp ;
    FILE   *sembfp ;
    string f_rmo ;
    string rmogrid ;
    string f_semb ;

    int    rmo_option ;

    /* variable for handling the trace header parameters */
    String tracekey, linekey, gatherkey ;
    String trktype,  lnktype, gatherktype ;
    Value  trkval,   lnkval,  gatherkval ;
    int    indxtrk,  indxlnk, indxgatherk ;
    float  fx, dx, fy, dy, x, y ;
    int    i2, i3 ;

    int    traceno, lineno, cdplbl ;
    float  *tnout, *rnout ;
    int    ip, nrmo, nrmo1 ;

    /* variables for rnmo grid files */
    int    ntrgrid ;
    float  dtrgrid ;
    float  ftrgrid ;
    int    n1,n2,n3,n4,n5 ;
    float  o1,o2,o3,o4,o5,d1,d2,d3,d4,d5 ;
    float  scale, ocdp2, oline3, dcdp2, dline3 ;
    float  rmin, rmax ;
    int    dtype,ierr,orient,gtype ;
    ghed   gh ;

    int   nt, nx, maxnx ;
    int   it, ix ;
    float t0, dt ;

    float **cmp ;
    float **flat ;
    float *xx ;
    float *aaopt ;
    float *sembmax ;

    float maxoffset ;
    int   maxoffset_status ;  /* 0 = read from header of 1st gather, */ 
                              /* 1 = provided by user */
    int   icdp ;

    float amin, amax ;
    int   na ;
    int   ntri ;
    int   sscale;   /* if non-zero, apply NMO stretch scaling */

    char  *headers ;  /* trace header info. */
    int   cdp;        /* cdp from current input trace header */
    int   cdpprev;    /* cdp from previous input trace header */
    int   gottrace;   /* =1 if an input trace was read */
    int   done;       /* =1 if done with everything */
    
    
    /* hook up getpar */
    initargs(argc, argv) ;
    askdoc(0) ;
   
    /* read headers */
    gethdr(&ch, &bh) ;

    /* optional parameters */

    if( !getparint  ("rmo_option", &rmo_option) ) rmo_option = 1 ;

    if( !getparstring("gatherkey",&gatherkey) ) gatherkey = "cdp" ;
    if( !getparstring("tracekey",&tracekey) ) tracekey ="tracf" ;
    if( !getparstring("linekey",&linekey) ) linekey = "fldr" ;

    if( !getparint  ("maxnx", &maxnx) ) maxnx = 120 ;

    if( !getparfloat("maxoffset", &maxoffset) ) {
      maxoffset_status = 0 ;
    } else {
      maxoffset_status = 1 ;
      maxoffset *= 0.001 ;
    }
   
    if( !getparfloat("amin", &amin) ) amin =-150. ;
    if( !getparfloat("amax", &amax) ) amax = 150. ;

    amin = amin*0.001 ;
    amax = amax*0.001 ;

    if( !getparint  ("na", &na) ) na = 30 ;
    if( !getparint  ("ntri", &ntri) ) ntri = 25 ;
    if( !getparint("sscale", &sscale) ) sscale = 1 ;

    if( !getparint("nrmo", &nrmo) ) nrmo = 40 ;
    if( !getparstring("f_rmo",&f_rmo) ) f_rmo = "rmo.hvel" ;
    if( !getparstring("f_semb",&f_semb) ) f_semb = "semb.su" ;

    if( rmo_option == 1 || rmo_option == 2 ) {
      if( (rmofp=fopen(f_rmo, "w")) == NULL ) {
        fprintf (stderr, "%s: can not open %s \n", argv[0], f_rmo);
        exit(-1) ;
      }
/* uncomment for handle file size larger than 2GB */
      file2g(rmofp) ;

      if( (sembfp=fopen(f_semb, "w")) == NULL ) {
        fprintf (stderr, "%s: can not open %s \n", argv[0], f_semb);
        exit(-1) ;
      }
/* uncomment for handle file size larger than 2GB */
      file2g(sembfp) ;
    }

    if( rmo_option == 3 ) {
      if (!getparstring("rmogrid",&rmogrid) )
        err(" File rmogrid must be specified ") ;
      if( (rmogrdfp = fopen(rmogrid,"r"))==NULL )
        err("Input rmogrid file %s not found \n",rmogrid) ;
 
      /* obtain grid header info */
/* uncomment for handle file size larger than 2GB */
      file2g(rmogrdfp);
      ierr = fgetghdr(rmogrdfp, &gh);
      if(ierr==0) fromghdr(&gh,&scale,&dtype,&n1,&n2,&n3,&n4,&n5,
                          &d1,&d2,&d3,&d4,&d5,&o1,&o2,&o3,&o4,&o5,
                          &dcdp2,&dline3,&ocdp2,&oline3,
                          &rmin,&rmax,&orient,&gtype) ;
    }
    
    /* get parameters from the first trace */
    if( !gettr(&tr) ) err( "can't get first trace" );
    nt = tr.ns ;
    dt = (float)tr.dt / 1000000.0 ;
    t0 = tr.delrt / 1000.0 ;


    /* handle the trace keys */
    gatherktype = hdtype(gatherkey);
    trktype = hdtype(tracekey);
    lnktype = hdtype(linekey);

    indxgatherk = getindex(gatherkey);
    indxtrk = getindex(tracekey);
    indxlnk = getindex(linekey);

    if( rmo_option == 3 ) {
      if(!getparfloat("ocdp2",&fx) && ierr==0) fx = ocdp2;
      if(!getparfloat("dcdp2",&dx) && ierr==0) dx = dcdp2;
      if(!getparfloat("oline3",&fy) && ierr==0) fy = oline3;
      if(!getparfloat("dline3",&dy) && ierr==0) dy = dline3;

      if( !getparint("ntrgrid",&ntrgrid) ) {
        if(ierr==0) {
          ntrgrid = n1;
        } else {
          ntrgrid = nt;
        }
      }
      if( !getparfloat("dtrgrid",&dtrgrid) ) {
        if(ierr==0) {
          dtrgrid = d1 * 0.001;
        } else {
          dtrgrid = dt;
        }
      } else {
        dtrgrid = dtrgrid*0.001;
      }
      if( !getparfloat("ftrgrid",&ftrgrid) ) {
        if(ierr==0) {
          ftrgrid = o1 * 0.001;
        } else {
          ftrgrid = t0;
        }
      } else {
        ftrgrid = ftrgrid*0.001;
      }
    }


    /* trace number and line number */ 
    gethval(&tr,indxtrk,&trkval);
    traceno = vtoi(trktype,trkval);
    gethval(&tr,indxlnk,&lnkval);
    lineno = vtoi(lnktype,lnkval);

    if( rmo_option == 3 ) {
      x = (traceno - fx)/dx + 0.5;
      i2 = x;
      y = (lineno - fy)/dy + 0.5;
      i3 = y;
  
      if(i2<0) i2=0;
      if(i2>n2-1) i2=n2-1;
      if(i3<0) i3=0;
      if(i3>n3-1) i3=n3-1;
    }


    /* gather assemble identifier */
    gethval(&tr,indxgatherk,&gatherkval);
    cdp = vtoi(gatherktype,gatherkval);

    /* allocate memory */
    headers = (char*) emalloc(maxnx*HDRBYTES);
    cmp  = ealloc2float(nt, maxnx) ;
    flat = ealloc2float(nt, maxnx) ;
    aaopt = ealloc1float(nt) ;
    sembmax = ealloc1float(nt) ;
    xx   = ealloc1float(maxnx) ;

    tnout = ealloc1float(nrmo) ;
    rnout = ealloc1float(nrmo) ;


    /* output binary header */
    if( rmo_option == 1 || rmo_option == 2 ) {
      // fputhdr(sembfp, &ch, &bh) ;
    }
    if( rmo_option == 2 || rmo_option == 3 ) {
      // puthdr(&ch, &bh) ;
    }
  

    /* initialize flags */
    gottrace = 1 ;
    done = 0 ;

    /* remember previous cdp */
    cdpprev = cdp ;
    nx = 0 ;

    icdp = 0 ;
    /* loop all cdps */
    do {
       /* if got a trace */
       if( gottrace ) {

         /* determine offset and cdp */
         gethval(&tr,indxgatherk,&gatherkval);
         cdp = vtoi(gatherktype,gatherkval);

         if( cdp == cdpprev ) {

           /* get trace samples */

           /* save the trace header info. */
           bcopy( (char*)&tr, headers+nx*HDRBYTES, HDRBYTES ) ;
           /* for( it=0; it<nt; it++ ) cmp[nx][it] = tr.data[it] ; */
           bcopy( tr.data, cmp[nx], nt*sizeof(float) ) ;
           xx[nx] = tr.offset/1000.0 ;
           nx++ ;
         }
       }

       /* if cdp has changed or no more input traces */
       if( cdp != cdpprev || !gottrace ) {

         if( maxoffset_status == 0 ) {
           maxoffset = xx[nx-1] ;
           maxoffset_status = 1 ;  /* reset the status */
           fprintf(stderr,"Warning:\n") ;
           fprintf(stderr,"Max.offset of the 1st gather is used.  %.1f\n", 
                   maxoffset*1000.0) ;
         }

         /* perform Residual NMO correction */

         if( rmo_option == 1 || rmo_option == 2 ) {

           nrmo1 = nrmo ;
           for( ip=0; ip<nrmo; ip++ ) {
              tnout[ip] = -9.999 ; 
              rnout[ip] = 0.0 ;
           }

           rnmo(cmp, nt, nx, t0, dt, amin, amax, na, ntri, sscale, xx, flat,
                aaopt, sembmax, rmo_option, nrmo, &nrmo1, tnout, rnout) ;

           /* output the optimal semblance trace */
           bcopy( headers, (char*)&tr_out, HDRBYTES ) ;
           bcopy( sembmax, tr_out.data, nt*sizeof(float) ) ;
           tr_out.tracl = icdp + 1 ;
           fputtr(sembfp, &tr_out) ;


         } else if( rmo_option == 3 ) {

           /* prepare rmo correction function at the current cdp location */
           readv(rmogrdfp,ntrgrid,n2,n3,i2,i3,dtrgrid,ftrgrid,nt,dt,t0,aaopt);


           for( it=0; it<nt; it++ ) {
              aaopt[it] /= (maxoffset*maxoffset*1000.0) ; 
           }

           rnmo_apply(cmp, nt, nx, t0, dt, sscale, xx, flat, aaopt) ;
         }

         icdp++ ;

         /* output the NMO corrected gather */
         if( rmo_option == 2 || rmo_option == 3 ) {
           for( ix=0; ix<nx; ix++ ) {
              /* put trace header back */
              bcopy( headers+ix*HDRBYTES, (char*)&tr_out, HDRBYTES ) ;
              /* for( it=0; it<nt; it++ ) tr_out.data[it] = flat[ix][it] ; */
              bcopy( flat[ix], tr_out.data, nt*sizeof(float) ) ;
              puttr(&tr_out) ;
           }
         }
         
         /* output the residual correction at maximum offset in HANDVEL */

         if( rmo_option == 1 || rmo_option == 2 ) {
           cdplbl = lineno*10000 + traceno ;

           for( ip=0; ip<nrmo1; ip++ ) {
              tnout[ip] *= 1000.0 ;
              rnout[ip] *= maxoffset*maxoffset*1000.0 ;
           }

           printhvel(cdplbl, nrmo1, tnout, rnout, rmofp) ;
         }

         /* save the first trace of the next cdp */
         if( gottrace ) {
           /* start a new cdp gather, reset the offset counter */
           nx = 0 ;
           /* save the trace header info. */
           bcopy( (char*)&tr, headers+nx*HDRBYTES, HDRBYTES ) ;
           /* for( it=0; it<nt; it++ ) cmp[nx][it] = tr.data[it] ; */
           bcopy( tr.data, cmp[nx], nt*sizeof(float) ) ;
           xx[nx] = tr.offset/1000.0 ;
           nx++ ;

           /* update trace number and line number */
           gethval(&tr,indxtrk,&trkval);
           traceno = vtoi(trktype,trkval);
           gethval(&tr,indxlnk,&lnkval);
           lineno = vtoi(lnktype,lnkval);

           if( rmo_option == 3 ) {
             x = (traceno - fx)/dx + 0.5;
             i2 = x;
             y = (lineno - fy)/dy + 0.5;
             i3 = y;
  
             if(i2<0) i2=0;
             if(i2>n2-1) i2=n2-1;
             if(i3<0) i3=0;
             if(i3>n3-1) i3=n3-1;
           }
 
         }

         /* if no more input traces, break input trace loop */
         if( !gottrace )  break ;

         /* remember previous cdp and line number */
         cdpprev = cdp ;
       }

       /* get next trace (if there is one) */
       if( !gettr(&tr) ) gottrace = 0 ;

    } while( !done ) ;
   
    if( rmo_option == 1 || rmo_option == 2 ) {
      fclose(rmofp) ;
      fclose(sembfp) ;
    }

    /* free memory */
    free2float(cmp) ;
    free2float(flat) ;
    free1float(aaopt) ;
    free1float(sembmax) ;
    free1float(xx) ;
    free1float(tnout) ;
    free1float(rnout) ;
    free(headers) ;


    return EXIT_SUCCESS ;
}

 

/* 
 * Residual NMO scan and correction 
 */

#define  EPS 1.0e-12

void boxconv(int nb, int nx, float *xx, float *yy) ;
void triangle(int nr, int n12, float *uu, float *vv) ;

void vresamp(int np, float *ts, float *vs, int nt, float t0, float dt,
             float *v ) ;
 
void rnmo(float **cmp, int nt, int nx, float t0, float dt, float amin, 
          float amax, int na, int ntri, int sscale, float *xx, float **flat,
          float *aaopt, float *sembmax, int rmo_option, int nrmo, int *nrmo1, 
          float *tnout, float *rnout)
{

    int   it, ix, ia, it1, it2, it3 ;
    float df, trnmo, tnmo ;
    float sum1, sum2, semb ;
    float amp1, amp2 ;
    float aa, a0, da ;
    float a1, a2, a3 ; 
    float t1, slope ;

    float xmax, amin1, amax1 ;
    
    float *xx2, *tau ;