sudmofk.c

seismic software,very useful

/* SUDMOFK: $Revision: 1.6 $ ; $Date: 90/06/14 14:38:34 $		*/

/*----------------------------------------------------------------------
 * Copyright (c) Colorado School of Mines, 1990.
 * All rights reserved.
 *
 * This code is part of SU.  SU stands for Seismic Unix, a processing line
 * developed at the Colorado School of Mines, partially based on Stanford
 * Exploration Project (SEP) software.  Inquiries should be addressed to:
 *
 *  Jack K. Cohen, Center for Wave Phenomena, Colorado School of Mines,
 *  Golden, CO 80401  ([e-mail])
 *----------------------------------------------------------------------
 */

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

/*********************** self documentation ******************************/
string sdoc =
"\n"
"SUDMOFK - DMO via F-K domain (log-stretch) method for common-offset gathers\n"
"\n"
"sudmofk <stdin >stdout cdpmin= cdpmax= dxcdp= noffmix= [optional parms]\n"
"\n"
"Required Parameters:\n"
"cdpmin                  minimum cdp (integer number) for which to apply DMO\n"
"cdpmax                  maximum cdp (integer number) for which to apply DMO\n"
"dxcdp                   distance between successive cdps\n"
"noffmix                 number of offsets to mix (see notes)\n"
"\n"
"Optional Parameters:\n"
"vmin=1500.0             minimum velocity used to determine maximum slope\n"
"tdmo=0.0                times corresponding to rms velocities in vdmo\n"
"vdmo=vmin               rms velocities corresponding to times in tdmo\n"
"verbose=0               =1 for diagnostic print\n"
"\n"
"Notes:\n"
"Input traces should be sorted into common-offset gathers.  One common-\n"
"offset gather ends and another begins when the offset field of the trace\n"
"headers changes.\n"
"\n"
"The cdp field of the input trace headers must be the cdp bin NUMBER, NOT\n"
"the cdp location expressed in units of meters or feet.\n"
"\n"
"The number of offsets to mix (noffmix) should typically equal the ratio of\n"
"the shotpoint spacing to the cdp spacing.  This choice ensures that every\n"
"cdp will be represented in each offset mix.  Traces in each mix will\n"
"contribute through DMO to other traces in adjacent cdps within that mix.\n"
"\n"
"The tdmo and vdmo arrays specify a velocity function of time that is\n"
"used to implement a first-order correction for depth-variable velocity.\n"
"The times in tdmo must be monotonically increasing.\n"
"\n"
"For each offset, the minimum time at which a non-zero sample exists is\n"
"used to determine a mute time.  Output samples for times earlier than this\n" 
"mute time will be zeroed.  Computation time may be significantly reduced\n"
"if the input traces are zeroed (muted) for early times at large offsets.\n"
"\n"
"Trace header fields accessed:  nt, dt, delrt, offset, cdp.\n"
"\n";
/**************** end self doc *******************************************/

/* Credits:
 *	CWP: Dave
 *
 * Technical Reference:
 *      Dip-Moveout Processing - SEG Course Notes
 *      Dave Hale, 1988
 */

void dmooff (float offset, int nxfft, float dx, 
	int nt, float dt, float ft,
	float vmin, int ntvdmo, float *tdmo, float *vdmo,
	float **ptx);
void maketu (float offset, int itmin, float vmin, 
	int ntvdmo, float *tdmo, float *vdmo,
	int nt, float dt, float ft, float **uoftp,
	int *nup, float *dup, float *fup, float **tofup, float *tconp);

segy tr,tro;

main(int argc, char **argv)
{
	int nt;		/* number of time samples per trace */
	float dt;	/* time sampling interval */
	float ft;	/* time of first sample */
	int it;		/* time sample index */
	int cdpmin;	/* minimum cdp to process */
	int cdpmax;	/* maximum cdp to process */
	float dx;	/* cdp sampling interval */
	int nx;	        /* number of cdps to process */
	int nxfft;	/* number of cdps after zero padding for fft */
	int nxpad;	/* minimum number of cdps for zero padding */
	int ix;		/* cdp index, starting with ix=0 */
	int noffmix;	/* number of offsets to mix */
	float vmin;	/* minimum velocity */
	float *tdmo;	/* times at which rms velocities are specified */
	float *vdmo;	/* rms velocities at times specified in tdmo */
	int ntdmo;	/* number tnmo values specified */
	int itdmo;	/* index into tnmo array */
	int nvdmo;	/* number vnmo values specified */
	float **p;	/* traces for one offset - common-offset gather */
	float **q;	/* DMO-corrected and mixed traces to be output */
	float offset;	/* source-receiver offset of current trace */
	float oldoffset;/* offset of previous trace */
	int noff;	/* number of offsets processed in current mix */
	int ntrace;	/* number of traces processed in current mix */
	int itrace;	/* trace index */
	int gottrace;	/* non-zero if an input trace was read */
	int done;	/* non-zero if done */
	int verbose;	/* =1 for diagnostic print */
	FILE *hfp;	/* file pointer for temporary header file */

	/* hook up getpar */
	initargs(argc, argv);
	askdoc(1);

	/* get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = (float)tr.dt/1000000.0;
	ft = (float)tr.delrt/1000.0;
	offset = tr.offset;

	/* get parameters */
	if (!getparint("cdpmin",&cdpmin)) err("must specify cdpmin");
	if (!getparint("cdpmax",&cdpmax)) err("must specify cdpmax");
	if (cdpmin>cdpmax) err("cdpmin must not be greater than cdpmax");
	if (!getparfloat("dxcdp",&dx)) err("must specify dxcdp");
	if (!getparint("noffmix",&noffmix)) err("must specify noffmix");
	if (!getparfloat("vmin",&vmin)) vmin = 1500.0;
	ntdmo = countparval("tdmo");
	if (ntdmo==0) ntdmo = 1;
	tdmo = ealloc1float(ntdmo);
	if (!getparfloat("tdmo",tdmo)) tdmo[0] = 0.0;
	nvdmo = countparval("vdmo");
	if (nvdmo==0) nvdmo = 1;
	if (nvdmo!=ntdmo) err("number of tdmo and vdmo must be equal");
	vdmo = ealloc1float(nvdmo);
	if (!getparfloat("vdmo",vdmo)) vdmo[0] = vmin;
	for (itdmo=1; itdmo<ntdmo; ++itdmo)
		if (tdmo[itdmo]<=tdmo[itdmo-1])
			err("tdmo must increase monotonically");
	if (!getparint("verbose",&verbose)) verbose=0;
	
	/* determine number of cdps to process */
	nx = cdpmax-cdpmin+1;
	
	/* allocate and zero common-offset gather p(t,x) */
	nxpad = 0.5*ABS(offset/dx);
	nxfft = npfaro(nx+nxpad,2*(nx+nxpad));
	p = ealloc2float(nt,nxfft+2);
	for (ix=0; ix<nxfft; ++ix)
		for (it=0; it<nt; ++it)
			p[ix][it] = 0.0;

	/* allocate and zero offset mix accumulator q(t,x) */
	q = ealloc2float(nt,nx);
	for (ix=0; ix<nx; ++ix)
		for (it=0; it<nt; ++it)
			q[ix][it] = 0.0;
		
	/* open temporary file for headers */
	/*hfp = tmpfile();*/
	hfp = etempfile(NULL);
	
	/* initialize */
	oldoffset = offset;
	gottrace = 1;
	done = 0;
	ntrace = 0;
	noff = 0;

	/* loop over traces */
	do {
		/* if got a trace, determine offset */
		if (gottrace) offset = tr.offset;
		
		/* if an offset is complete */
		if ((gottrace && offset!=oldoffset) || !gottrace) {
		
			/* do dmo for old common-offset gather */
			dmooff(oldoffset,nxfft,dx,nt,dt,ft,vmin,
				ntdmo,tdmo,vdmo,p);
			
			/* add dmo-corrected traces to mix */
			for (ix=0; ix<nx; ++ix)
				for (it=0; it<nt; ++it)
					q[ix][it] += p[ix][it];
			
			/* count offsets in mix */
			noff++;
							
			/* free space for old common-offset gather */
			free2float(p);
			
			/* if beginning a new offset */
			if (offset!=oldoffset) {
				
				/* allocate space for new offset */
				nxpad = 0.5*ABS(offset/dx);
				nxfft = npfaro(nx+nxpad,2*(nx+nxpad));
				p = ealloc2float(nt,nxfft+2);
				for (ix=0; ix<nxfft; ++ix)
					for (it=0; it<nt; ++it)
						p[ix][it] = 0.0;
			}
		}
		
		/* if a mix of offsets is complete */
		if (noff==noffmix || !gottrace) {
			
			/* rewind trace header file */
			fseek(hfp,0L,SEEK_SET);
			
			/* loop over all output traces */
			for (itrace=0; itrace<ntrace; ++itrace) {
			
				/* read trace header and determine cdp index */
				efread(&tro,HDRBYTES,1,hfp);
				
				/* get dmo-corrected data */
				bcopy(q[tro.cdp-cdpmin],tro.data,
					nt*sizeof(float));
				
				/* write output trace */
				puttr(&tro);
			}
			
			/* report */
			if (verbose) 
				fprintf(stderr,"\tCompleted mix of "
					"%d offsets with %d traces\n",
					noff,ntrace);
			
			/* if no more traces, break */
			if (!gottrace) break;
			
			/* rewind trace header file */
			fseek(hfp,0L,SEEK_SET);
			
			/* reset number of offsets and traces in mix */
			noff = 0;
			ntrace = 0;
			
			/* zero offset mix accumulator */
			for (ix=0; ix<nx; ++ix)
				for (it=0; it<nt; ++it)
					q[ix][it] = 0.0;
		}
			
		/* if cdp is not within range to process, skip it */
		if (tr.cdp<cdpmin || tr.cdp>cdpmax) continue;
	
		/* save trace header and update number of traces */
		efwrite(&tr,HDRBYTES,1,hfp);
		ntrace++;

		/* remember offset */
		oldoffset = offset;

		/* get trace samples */
		bcopy(tr.data,p[tr.cdp-cdpmin],nt*sizeof(float));