gmo.c

seismic software,very useful

/* GMO apply gamma moveout to traces*/

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


/*********************** self documentation **********************/
string sdoc = 
"GMO - apply gamma moveout to traces 		 	\n"
"\n"
"gmo [parameters] < input.data >output.data 				\n" 
"\n"
"Required parameters:							\n"
"ggfile=	         file of gamma grid file (generated by ggrid)	\n"
"Optional parameters:							\n"
"fzg=from-ggfile         minimum depth of input gamma grid 		\n"
"dzg=from-ggfile         depth interval of input gamma grid 		\n"
"nzg=from-ggfile         number of depth intervals of input gamma grid 	\n"
"cdpming=from-ggfile     first cdp number of input gamma grid		\n" 
"dcdpg=from-ggfile       cdp number increment of input gamma grid	\n" 
"ncdpg=from-ggfile       number of cdp's of input gamma grid		\n" 
"smute=1.5	  samples with GMO stretch exceeding smute are zeroed	\n"  
"\n"
"NOTE: \n"
"  1.   If the ggfile is a standard grid file, (i.e., with the header),	\n"
"       fzg, dzg, nzg, cdpming, dcdpg, and ncdpg are default from the  	\n"
"       header of the ggfile when they are not specified.		\n"
"       If the ggfile is not standard, (i.e., without grid header), 	\n"
"       these six parameters  will be default to: 			\n"
"           fzg = minimum depth of migrated trace			\n"
"           dzg = depth interval of migrated trace			\n"
"           nzg = number of depth samples of migrated trace		\n"
"           cdpming = cdp of first input migrated trace			\n"
"           dcdpg = 1							\n"
"           ncdpg = computed from the size of the ggfile 		\n"
"       If the ggfile is standard, and user also supplies these six	\n"
"       parameters, a check will be done to make sure they are 		\n"
"       consistant.							\n"  
"  2.  gamma is defined as Vnew/Vmig					\n"
"\n"
"AUTHOR:		Zhiming Li,       ,	9/12/91   \n"		    
" \n"
;
/**************** end self doc ***********************************/

segy tr;
ghed gh;

main(int argc, char **argv)
{
    string ggfile;
    FILE *infp=stdin,*outfp=stdout,*gfp;
    int idz, nz, rsampz;
    int nzg, jz;
    float dz, dzint,z; 	
    float fzg, dzg; 
    float *gamgrid, *g2m1, *mout2, *dataint, *zint, *zx , *mute;
    float *g2;
    int oldcdp, oldoffset,ncdpg;
    float cdpming, dcdpg;
    float smute, osmute, odz, odz2, z0, fz;
    int mutes, newgamma ;
    int indz, indx, nzint, iz;
    float temp, res, zzint, odzint, str;
    int ierr;


    /* get parameters */
    initargs(argc,argv);
    askdoc(1);
/* 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;
    fz = tr.delrt/dz;
    if ( fz >1000. ) fz=fz/1000.;

    /* if dz defined in segy trace header, get it */
    if ( tr.dz != 0. ) {
	dz = tr.dz;
	fz = tr.fz;
    }

    nz = tr.ns;
    if ( !getparstring("ggfile",&ggfile)) err("ggfile must be given ! \n");
    if ( !getparfloat("smute",&smute)) smute=1.5;
    /* open ggfile and read header */
    gfp = efopen(ggfile,"r");
    ierr = fgetghdr(gfp,&gh);

	if(ierr!=0) warn(" grid header error \n");

    if ( !getparfloat("cdpming",&cdpming) ) {
	if(ierr==0) {
		getgval(&gh,"ocdp2",&cdpming);
		if(cdpming==0.) 
	      err(" grid file gh.ocdp2=0.; update gh.ocdp2 or spicify cdpming");
	} else {
		cdpming = tr.cdp;
	}
    }
    if ( !getparfloat("dcdpg",&dcdpg) ) {
	if(ierr==0) {
		getgval(&gh,"dcdp2",&dcdpg);
		if(dcdpg==0.) 
		err(" grid file gh.dcdp2=0.; update gh.dcdp2 or spicify dcdpg");
	} else {
		dcdpg = 1;
	}
    }
    if ( !getparfloat("fzg",&fzg) ) {
	if(ierr==0) {
		getgval(&gh,"o1",&fzg);
	} else {
		fzg = fz;
	}
    }
    if ( !getparfloat("dzg",&dzg) ) {
	if(ierr==0) {
		getgval(&gh,"d1",&dzg);
	} else {
		dzg = dz;
	} 
    }	

    if ( !getparint("nzg",&nzg) ) {
	if(ierr==0) {
		nzg = (int) (float)( gh.n1/gh.scale + 0.00001);
	} else {
		nzg = nz;
	}
    }
    if ( !getparint("ncdpg",&ncdpg) ) {
	if(ierr==0) {
		ncdpg = (int) (float)( gh.n2/gh.scale + 0.00001);
	} else {
    		fseek(gfp,0L,2);
    		ncdpg = ftell(gfp)/sizeof(float)/nzg;
        	fseek(gfp,0L,0);
	}
    }
    /* error checking */
    if (ierr==0) {
		if(ncdpg!=(int) (float)( gh.n2/gh.scale + 0.00001))
			err(" check gh.n2 ");
		if(nzg!=(int) (float)( gh.n1/gh.scale + 0.00001))
			err(" check gh.n1 ");
		if(dzg!=gh.d1/gh.scale)
			warn("dzg=%g not the same as gh.d1 = %g ", dzg, gh.d1/gh.scale);
		if(fzg!=gh.o1/gh.scale)
			warn("fzg not the same as gh.o1");
		if(dcdpg!=gh.dcdp2/gh.scale)
			warn("dcdpg not the same as gh.dcdp2");
		if(cdpming!=gh.ocdp2/gh.scale)
			warn("cdpming not the same as gh.ocdp2");
    } else {
                warn(" input ggfile non-standard, defaults used ");
    }

    if(dcdpg==0.) err(" dcdpg must not be 0"); 

    gamgrid = (float*)malloc(nzg*sizeof(float));

    /* sinc interpolate input traces before gmo for better accuracy */
    rsampz = 4.;
    nzint = nz * rsampz;
    dzint = 1./rsampz;
    odzint = rsampz;

    /* memory allocation */
    g2m1 = (float*)malloc(nz*sizeof(float));
    g2 = (float*)malloc(nz*sizeof(float));
    dataint = (float*)malloc(nzint*sizeof(float));
    zint = (float*)malloc(nzint*sizeof(float));
    mout2 = (float*)malloc(nz*sizeof(float));
    zx = (float*)malloc(nz*sizeof(float));
    mute = (float*)malloc(nz*sizeof(float));

    /* main loop over input traces */
    oldcdp = tr.cdp - 1;
    oldoffset=tr.offset - 1;
    odz = 1./dz;
    odz2 = odz*odz;

    for(iz=0;iz<nzint;iz++) zint[iz] = fz + iz*dzint;
    osmute = 1./smute;


    do { 

	/* find cdp index in gamma grid */
	if(tr.cdp!=oldcdp) {
	   newgamma = 1;
	   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(gamgrid,sizeof(float),nzg,gfp);
	   /* compute (gamma^2-1) */
	   for(iz=0;iz<nz;iz++) { 
		z = fz + iz * dz;
		/* interpolate input ggrid */
		z = (z-fzg)/dzg;
		jz = z;
		if(jz<1) {
			jz = 1;
		} else if (jz>nzg-1) {
			jz = nzg - 1; 
		}
		/* store as 1./g2 */
		g2[iz] = 1./(gamgrid[jz]*gamgrid[jz]);
		g2m1[iz] = g2[iz]-1.;
	   }		
	} else {
	   newgamma = 0;
	}
	if (newgamma || tr.offset!=oldoffset) {
	   /* compute 0.25*offset^2/dz^2*(gamma^2-1)  */
	   temp = (tr.offset*tr.offset)*odz2*0.25;
	   for(iz=0;iz<nz;iz++) mout2[iz] = temp*g2m1[iz];
	}

	/* sinc interpolate the input trace */
	if(nzint>nz) {
	   ints8r(nz,1.0,fz,tr.data,0.0,0.0,nzint,zint,dataint);
	} else {
	   for(iz=0;iz<nz;iz++) dataint[iz]=tr.data[iz];
	}

	/* compute depth z(z0) (in sample) */
	for(iz=0,z0=fz; iz<nz; ++iz, z0+=1.0) { 
	   temp = z0*z0 + mout2[iz];
	   if (temp >= 0.) {
	      zx[iz] = sqrt(temp);
	      mute[iz] = 1.;
	   }
	   else {
	      zx[iz] = fz - 1.;
	      mute[iz] = 0.;
	   }
	}
	/* gmo */
	/* determine if stretch mute is needed */
	str = zx[1]-zx[0];
	if (g2m1[0]>=0.&& str < smute) {
	   mutes = 0;
	} else if(g2m1[0]<0.&& str > osmute) {
	   mutes = 0;
	} else {
	   mutes = 1;
	}
	if(mutes!=1){
	   zzint = (zx[0] - fz)*odzint;
	   indz = zzint;
	   res = zzint - indz;
	   if(indz>0 && indz<nzint-1) {
	      tr.data[0] = dataint[indz]+res*(dataint[indz+1]-dataint[indz]);	
	   } else {
	      tr.data[0] = 0.;	
	   }		
	} else {
	      tr.data[0] = 0.;
	}
	for(iz=1;iz<nz;iz++) {
	   str = zx[iz]-zx[iz-1];
	   if (g2m1[iz]>=0.&& str < smute) {
	      mutes = 0;
	   } else if(g2m1[iz]<0.&& str > osmute) {
	      mutes = 0;
	   } else {
	      mutes = 1;
	   }
	   if(mutes!=1) {
	      zzint = (zx[iz] - fz)*odzint;
	      indz = zzint;
	      res = zzint - indz;
	      if(indz>0 && indz<nzint-1) {
	         tr.data[iz]=dataint[indz]+res*(dataint[indz+1]-dataint[indz]);	
	      } else {
	         tr.data[iz]=0.;	
	      }
	   } else {
	         tr.data[iz]=0.;	
	   }
	}
	/* zero out those with negative inside sqrt */
	for(iz=0;iz<nz;iz++)tr.data[iz]=tr.data[iz]*mute[iz];
	   
	/* output */
	fputtr(outfp,&tr);
	/* update offset and cdp */
	oldcdp = tr.cdp;
	oldoffset = tr.offset;

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

    free(gamgrid);
    free(g2m1);
    free(mout2);
    free(dataint);
    free(mute);
    free(zx);
    free(zint);

    fclose(outfp);
    

    return EXIT_SUCCESS;
}