rnmo.c

seismic software,very useful

    float *temp1, *temp2 ;
    float **trnmo_buf ;
    float *atnz ;
    
    int   *ind_ex ;
    int   ip ;
    float *extrema ;
    float clip_semb ;

    xmax = MAX(ABS(xx[0]), ABS(xx[nx-1])) ;
    xmax = xmax*xmax ;
    amin1 = amin/xmax ;
    amax1 = amax/xmax ;

    df = 1./dt ;
    da = (amax1-amin1)/(float)na ;
    a0 = amin1 ;

    /* allocate temporary arrays */
    xx2 = ealloc1float(nx) ;
    tau = ealloc1float(nt) ;
    temp1 = ealloc1float(nt) ;
    temp2 = ealloc1float(nt) ;
    trnmo_buf = ealloc2float(nt, nx) ;
    atnz = ealloc1float(nt) ;

    ind_ex = ealloc1int(nt) ;
    extrema = ealloc1float(nt) ;
    
    for( ix=0; ix<nx; ix++ ) {
       xx[ix] = ABS(xx[ix]) ;
       xx2[ix] = xx[ix]*xx[ix] ;
    }

    for( it=0; it<nt; it++ ) {
       tau[it] = t0 + it*dt ;
    }
    
    /* initialize arrays for sorting semblance */
    for( it=0; it<nt; it++ ) {
       ind_ex[it] = it ;
       extrema[it] = -1.0 + EPS*franuni() ; 
    }

    /* 
     * Scan and search for optimum parameter of residual moveout equation
     * a near offset approximation
     *    Trnmo = Tnmo + aa*X^2
     */

    /* Loop over NMO time */
    for( it=0; it<nt; it++ ) {
       tnmo = tau[it] ;
       sum1 = 0 ;  
       sum2 = 0. ;
       for( ix=0; ix<nx; ix++ ) {
          a1 = cmp[ix][it] ;
          sum1 += a1 ;
          sum2 += a1*a1 ;
       }
       sembmax[it] = (sum1*sum1)/(sum2+EPS) ;
       aaopt[it] = 0. ;

       /* Loop over residual moveouts */
       for( ia=0; ia<=na; ia++ ) {

          aa = ia*da + a0 ;

          sum1 = 0. ;
          sum2 = 0. ;

          /* Loop over offset */
          for( ix=0; ix<nx; ix++ ) {
             trnmo = tnmo + aa*xx2[ix] ;
             it2 = (trnmo-t0)/dt + 0.5 ;
             if( it2 >= 1 && it2 < nt-1 ) {
               it1 = it2-1 ;
               it3 = it2+1 ;
               a1 = cmp[ix][it1] ; 
               a2 = 2.*cmp[ix][it2] ; 
               a3 = cmp[ix][it3] ;

               amp1 = 0.25*(a1 + a2 + a3) ;
               amp2 = amp1*amp1 ;
               sum1 += amp1 ;
               sum2 += amp2 ;
             }
          } /* end of ix loop */
          
          semb = (sum1*sum1)/(sum2+EPS) ;

          if( semb > sembmax[it] && ABS(sembmax[it]) > 1000.0*EPS ) {
            sembmax[it] = semb ;
            aaopt[it] = aa ;
          }
       } /* end of ia loop */
    } /* end of it loop */


    /* 
     * Smoothing by triangle integration with semblance weight
     */ 

    for( it=0; it<nt; it++ ) {
       temp1[it] = sembmax[it]*aaopt[it] ;
    }

    for( it=0; it<nt; it++ ) {
       temp2[it] = sembmax[it] ;
    }

    /* Smooth aaopt by triangle smoother */ 

    if( ntri > 1 ) {
      triangle( ntri, nt, temp1, temp1 ) ;
      triangle( ntri, nt, temp2, temp2 ) ;
    }

    for( it=0; it<nt; it++ ) {
       aaopt[it] = temp1[it]/(temp2[it]+2.0*EPS) ;
    }

    /* 
     * compute the semblance using the optimal moveout 
     */

    /* Loop over NMO time */
    for( it=0; it<nt; it++ ) {

       tnmo = tau[it] ;
       aa = aaopt[it] ;
       sembmax[it] = 0. ;

       sum1 = 0. ;
       sum2 = 0. ;

       /* Loop over offset */
       for( ix=0; ix<nx; ix++ ) {
          trnmo = tnmo + aa*xx2[ix] ;
          it2 = (trnmo-t0)/dt + 0.5 ;
          if( it2 >= 1 && it2 < nt-1 ) {
            it1 = it2-1 ;
            it3 = it2+1 ;
            a1 = cmp[ix][it1] ; 
            a2 = 2.*cmp[ix][it2] ;
            a3 = cmp[ix][it3] ;

            amp1 = 0.25*(a1 + a2 + a3) ;
            amp2 = amp1*amp1 ;
            sum1 += amp1 ;
            sum2 += amp2 ;
          }
       } /* end of ix loop */
       sembmax[it] = (sum1*sum1)/(sum2+EPS) ;

    } /* end of it loop */

    /* smooth sembmax using 9-point triangle smoother */
    triangle( 9, nt, sembmax, sembmax ) ;
    triangle( 3, nt, sembmax, sembmax ) ;

    /* 
     * sort out the first nrmo largest semblance values 
     */

    /* search for extrema along the semblance curve */
    for( it=3; it<nt-3; it++ ) {
       if( sembmax[it] > sembmax[it-1] &&
           sembmax[it] > sembmax[it-2] &&
           sembmax[it] > sembmax[it-3] &&
           sembmax[it] > sembmax[it+1] &&
           sembmax[it] > sembmax[it+2] &&
           sembmax[it] > sembmax[it+3] ) {
         extrema[it] = sembmax[it] ;
       }
    }
    
    /* sort an index array based on semblance */
    qkisort( nt, extrema, ind_ex ) ;
    

    /* save the nrmo largest extrema, set the rest to -1.0 */
    for( it=0; it<nt-nrmo; it++ ) {
       extrema[ind_ex[it]] = -1.0 ;
    }

    /* find a clip value for sembalance */
    clip_semb = extrema[ind_ex[NINT(0.98*nt-1.0)]] ;

    /* save the selected optimal residual movemouts */ 
    ip = 0 ;
    for( it=0; it<nt-1; it++ ) {
       if( extrema[it] > 0.10*clip_semb ) {
/*
         sembmax[it] = -1.0 ;
*/
         tnout[ip] = tau[it] ;
         rnout[ip] = aaopt[it] ;
         ip++ ;
       }
    }
    *nrmo1 = ip ;

    /* 
     * Apply optimal estimates of RNMO
     */

    if( rmo_option == 2 ) {

      vresamp(ip, tnout, rnout, nt, t0, dt, aaopt) ;
/*
    for( it=0; it<nt; it++ ) {
       sembmax[it] = aaopt[it] ;
    }
*/
      /* Zero output storage */
      for( ix=0; ix<nx; ix++ )
         for( it=0; it<nt; it++)
            flat[ix][it] = 0. ;

      /* Compute the residual moveouts */
      for( it=0; it<nt; it++ ) {
     
         tnmo = tau[it] ;
         aa = aaopt[it] ;

         for( ix=0; ix<nx; ix++ ) {
            if( ABS(cmp[ix][it]) > EPS ) {
              trnmo_buf[ix][it] = (tnmo + aa*xx2[ix])*df ;
            } else {
              trnmo_buf[ix][it] = tnmo*df ;
            }
         }
      }
   
      /* Perform residual moveout correction */
      for( ix=0; ix<nx; ix++ ) {

         /* Do rnmo via 8-point sinc interploation */
         ints8r(nt, 1.0, t0*df, cmp[ix], 0.0, 0.0, nt, trnmo_buf[ix], flat[ix]);
         if( sscale ) {
           /* Compute inverse of stretch facter */
           atnz[0] = trnmo_buf[ix][1] - trnmo_buf[ix][0] ;
           for( it=1; it<nt; it++ ) {
              atnz[it] = trnmo_buf[ix][it] - trnmo_buf[ix][it-1] ;
           }

           /* Scale by the NMO stretch factor */
           for( it=0; it<nt; it++ ) {
              flat[ix][it] *= atnz[it] ;
           }
         }
      }
    }

    /* free temp arrays */
    free1float(xx2) ;
    free1float(tau) ;
    free1float(temp1) ;
    free1float(temp2) ;
    free2float(trnmo_buf) ;
    free1float(atnz) ;
    free1int(ind_ex) ;
    free1float(extrema) ;
}

/*
 * Apply residual NMO correction 
 */

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

    int   it, ix ;
    float df, trnmo, tnmo ;
    float aa ;
    
    float *xx2, *tau ;
    float **trnmo_buf ;
    float *atnz ;

    df = 1./dt ;

    /* allocate temporary arrays */
    xx2 = ealloc1float(nx) ;
    tau = ealloc1float(nt) ;
    trnmo_buf = ealloc2float(nt, nx) ;
    atnz = ealloc1float(nt) ;
    
    for( ix=0; ix<nx; ix++ ) {
       xx[ix] = ABS(xx[ix]) ;
       xx2[ix] = xx[ix]*xx[ix] ;
    }

    for( it=0; it<nt; it++ ) {
       tau[it] = t0 + it*dt ;
    }

    /* 
     * Apply optimal estimates of RNMO
     */

    /* Zero output storage */
    for( ix=0; ix<nx; ix++ )
       for( it=0; it<nt; it++)
          flat[ix][it] = 0. ;

    /* Compute the residual moveouts */
    for( it=0; it<nt; it++ ) {
     
       tnmo = tau[it] ;
       aa = aaopt[it] ;

       for( ix=0; ix<nx; ix++ ) {
          if( ABS(cmp[ix][it]) > EPS ) { 
            trnmo_buf[ix][it] = (tnmo + aa*xx2[ix])*df ;
          } else {
            trnmo_buf[ix][it] = tnmo*df ;
          }
       }
    }
   
    /* Perform residual moveout correction */
    for( ix=0; ix<nx; ix++ ) {

       /* Do rnmo via 8-point sinc interploation */
       ints8r(nt, 1.0, t0*df, cmp[ix], 0.0, 0.0, nt, trnmo_buf[ix], flat[ix]) ;

       if( sscale ) {
         /* Compute inverse of stretch facter */
         atnz[0] = trnmo_buf[ix][1] - trnmo_buf[ix][0] ;
         for( it=1; it<nt; it++ ) {
            atnz[it] = trnmo_buf[ix][it] - trnmo_buf[ix][it-1] ;
         }

         /* Scale by the NMO stretch factor */
         for( it=0; it<nt; it++ ) {
            flat[ix][it] *= atnz[it] ;
         }
       }
    }

    /* free temp arrays */
    free1float(xx2) ;
    free1float(tau) ;
    free2float(trnmo_buf) ;
    free1float(atnz) ;
} 




/* 
 * Convole with triangle
 */

void triangle(int nr, int n12, float *uu, float *vv)  
{
/* 
 * input:  nr -- rectangle width (points) (triangle base twice as wide)
 *         uu[i2], i2=0, n12-1 -- a vector of data
 * output: vv[i2], i2=0, n12-1 -- may be on top of uu 
 */
    int   i, np, nq ;
    float *pp, *qq, *tt ;
    
    /* allocate memory */
    pp = ealloc1float(n12+nr-1) ; 
    qq = ealloc1float(n12+nr+nr-2) ; 
    tt = ealloc1float(n12) ; 

    for( i=0; i<n12; i++ ) qq[i] = uu[i] ;
    if( n12 == 1 ) {
      for( i=0; i<n12; i++ ) tt[i] = qq[i] ;
    } else {
      boxconv( nr, n12, qq, pp) ;
      np = nr + n12 - 1 ;
      boxconv( nr, np, pp, qq ) ;
      nq = nr + np - 1 ;
      for( i=0; i<n12; i++ ) tt[i] = qq[i+nr-1] ;

      /* fold back near end */
      for( i=0; i<nr-1; i++ ) tt[i] = tt[i] + qq[nr-i-2] ;

      /* fold back far end */
      for( i=0; i<nr-1; i++ ) tt[n12-i-1] = tt[n12-i-1] + qq[n12+(nr-1)+i] ;
    }

    for( i=0; i<n12; i++ ) vv[i] = tt[i] ;

    /* free memory */
    free1float(pp) ;
    free1float(qq) ;
    free1float(tt) ;
}

void boxconv(int nb, int nx, float *xx, float *yy) 
{
/* 
 * input:  nx, xx[i], i=0, nx-1    the data
 *         nb -- the box length
 * output: yy[i], i=0, nx+nb-2     smoothed data
 */
    int   ny, i ;
    float *bb ;
    
    /* allocate memory */
    bb = ealloc1float(nx+nb) ; 

    if( nb < 1 || nb > nx ) err("error in boxconv's input") ;
    ny = nx + nb - 1 ;
    for( i=0; i<ny; i++ ) bb[i] = 0. ;

    /* make B(Z) = X(Z)/(1-Z) */ 
    bb[0] = xx[0] ;
    for( i=1; i<nx; i++ ) bb[i] = bb[i-1] + xx[i] ; 
    for( i=nx; i<ny; i++ ) bb[i] = bb[i-1] ;
    for( i=0; i<nb; i++ ) yy[i] = bb[i] ;  
    /* make Y(Z) = B(Z)*(1-Z**nb) */ 
    for( i=nb; i<ny; i++ ) yy[i] = bb[i] - bb[i-nb] ;
    for( i=0; i<ny; i++ ) yy[i] = yy[i] / nb ;
    
    /* free memory */
    free1float(bb) ;
}

/* 
 * read and interpolate velocity function from an Vnmo grid file 
 */

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) {

    float *vread, ratio, t, ftr;
    int icdp, it, itread;
    float tmp;
    long lpos;

    vread = (float *) malloc(ntvgrid*sizeof(float));

    icdp = ix + iy*nx;

    lpos = icdp;
    lpos = lpos*ntvgrid*sizeof(float);

    fseek(vgfp,lpos,0);
    fread(vread,sizeof(float),ntvgrid,vgfp);

    ratio = dt/dtvgrid;
    ftr = (ft - ftvgrid)/dtvgrid;

    for(it=0;it<nt;it++) {
        t = it*ratio + ftr;
        itread = t;
        if(itread < 0) {
            tmp = vread[0];
        } else if(itread >= ntvgrid-1) {
            tmp = vread[ntvgrid-1];
        } else {
            tmp = vread[itread] + (t-itread)*
                (vread[itread+1]-vread[itread]);
        }
        vel[it] = tmp;
    }

    free(vread);
}

/* 
 * Re-interpolate a (t,v) function to every sample point
 */ 

void vresamp(int np, float *ts, float *vs, int nt, float t0, float dt,
             float *v )
{
    int it, j, j1 ;
    float tmax, tout, fraction ;
    
    tmax = (nt-1)*dt + t0 ;

    j1 = 0 ;
    it = 0 ;
    tout = t0 ; 

    do {
       for( j=j1; j<np; j++ ) {
          if( ts[j] >= tout) {
            if( ts[j] == tout) {
              v[it] = vs[j] ;
            } else {
              if( tout < ts[0] ) {
                v[it] = vs[0] ;
              } else {
                fraction = (ts[j] - tout) / (ts[j] - ts[j-1]) ;
                v[it] = fraction*vs[j-1] + (1.-fraction)*vs[j] ;
              }
            }
            j1 = j ;
            break ;
          } else {
            fraction = (tmax - tout) / (tmax - ts[np-1]) ;
            v[it] = fraction*vs[np-1] ;
          }
       }
       it++ ;
       tout = tout + dt ;
    } while( tout <= tmax ) ; 
}