zstretch.c

seismic software,very useful

/* ZSTRETCH apply depth stretch to input traces via gamma-z relation*/

#include "su.h"
#include "segy.h"
#include "header.h"


/*********************** self documentation **********************/
string sdoc = 
"ZSTRETCH - apply depth stretch to traces 		 	\n"
"\n"
"zstretch ggfile= [parameters] < input.data >output.data 		\n" 
"\n"
"Required parameters:							\n"
"ggfile=	file of gamma (Vnew/Vold) grid file (generated by ggrid)\n"
"Optional parameters:							\n"
"fzg=fz         minimum depth of input gamma grid (fz is the minimum    \n"
"               depth from trace header of input depth migrated trace)  \n"
"dzg=dz         depth interval of input gamma grid (dz is the depth     \n"
"               interval from trace header of input depth migrated trace)\n"
"nzg=nt         number of depth intervals of input gamma grid (nz is the \n"
"               number of samples per trace from trace header of input  \n"
"               depth migrated trace)                                   \n"
"cdpming=       first cdp number of input gamma grid                    \n"
"               (default=the same as the cdp of first input trace       \n"
"dcdpg=1        cdp number increment of input gamma grid                \n"
"ncdpg=         number of cdp's of input gamma grid                     \n"
"               (default=computed internally from ggfile size)          \n"
"\n"
"NOTE: \n"
"       By default, fzg, dzg and nzg in input gamma grid file must be   \n"
"       the same as those in input depth migrated trace. The first cdp  \n"
"       number of the input gamma grid file is default to the cdp number\n"
"       of the first input trace. The cdp number increment in the gamma \n"
"       grid file defaults to 1. When input cdp is outside the range    \n"
"       [cdpming,cdpming+(ncdpg-1)*dcdpg], gamma function at edge(s) of \n"
"       the gamma grid will be used.                                    \n"
"AUTHOR:		Zhiming Li,       ,	9/12/91   \n"		    
" \n"
;
/**************** end self doc ***********************************/

segytrace tr;

main(int argc, char **argv)
{
    string ggfile;
    FILE *infp=stdin,*outfp=stdout,*gfp;
    int iz, nz;
    float dz, fz, odz; 	
    float *gamma, *znew, *zold, *datanew;
    int indx;
    int cdpold, idz;
    float res, zz;
    float fzg, dzg, tmp;
    int nzg, cdpming, dcdpg, ncdpg, jz;


    /* get parameters */
    initargs(argc,argv);
    askdoc(1);
    if (!getparstring("ggfile",&ggfile)) err("can't get ggfile name !\n");


    /* read in first trace */
    if (!fgettr(infp,&tr)) err("can't get first trace");
    idz = tr.dt;
    if(idz>1000) idz=idz/1000;
    dz = idz;
    nz = tr.ns;
    fz = tr.delrt;
    if(fz>1000.) fz=fz/1000.;
    if(tr.dz!=0.) {
	fz=tr.fz;
	dz=tr.dz;
    }
    odz = 1./dz;

    /* read in optional parameters */
    if ( !getparint("cdpming",&cdpming) ) cdpming = tr.cdp;
    if ( !getparint("dcdpg",&dcdpg) ) dcdpg = 1;
    if ( !getparfloat("fzg",&fzg) ) fzg = fz;
    if ( !getparfloat("dzg",&dzg) ) dzg = dz;
    if ( !getparint("nzg",&nzg) ) nzg = nz;

    /* open input gamma grid file */
    gfp = fopen(ggfile,"r");
    if ( !getparint("ncdpg",&ncdpg) ) {
        fseek(gfp,0L,2);
        ncdpg = ftell(gfp)/sizeof(float)/nzg;
        fseek(gfp,0L,0);
    }

    /* memory allocation */
    znew = (float*)malloc(nz*sizeof(float));
    zold = (float*)malloc(nz*sizeof(float));
    gamma = (float*)malloc(nzg*sizeof(float));
    datanew = (float*)malloc(nz*sizeof(float));
    for(iz=0;iz<nz;iz++) znew[iz] = fz + iz*dz;

    cdpold = tr.cdp - 1;
    /* main loop over input traces */
    do {
       if ( cdpold != tr.cdp ) {
          cdpold = tr.cdp;
	  indx = (tr.cdp-cdpming)/dcdpg;
          if(indx <0) {
             indx = 0;
          }
          else if (indx >ncdpg-1) {
             indx = ncdpg - 1;
          }
          fseek(gfp,indx*nzg*sizeof(float),0);
          fread(gamma,sizeof(float),nzg,gfp);
          for(iz=0;iz<nz;iz++) {
                /* interpolate input ggrid */
                tmp = (znew[iz]-fzg)/dzg;
                jz = tmp;
                if(jz<1) {
                        jz = 1;
                } else if (jz>nzg-1) {
                        jz = nzg - 1;
                }

		zold[iz] = ( znew[iz] / gamma[jz] - fz ) * odz;
	  }
       }
       for(iz=0;iz<nz;iz++) {
	  zz = zold[iz]; 
	  indx = zz;
	  res = zz - indx;
	  if ( indx>=0 && indx<nz-1 ){
	     datanew[iz]=tr.data[indx]+res*(tr.data[indx+1]-tr.data[indx]); 
	  }
	  else { 
	     datanew[iz] = 0.;
	  }
       }
       for(iz=0;iz<nz;iz++) tr.data[iz]=datanew[iz];		

       fputtr(outfp,&tr);

    }while(fgettr(infp,&tr));

    return EXIT_SUCCESS;
}