sudmotx.c

seismic software,very useful

/* SUDMOTX: $Revision: 1.6 $ ; $Date: 90/11/04 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"
"SUDMOTX - DMO via T-X domain (Kirchhoff) method for common-offset gathers\n"
"\n"
"sudmotx <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"
"offmax=3000.0           maximum offset\n"
"tmute=2.0               mute time at maximum offset offmax\n"
"vrms=1500.0             RMS velocity at mute time tmute\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 defaults for offmax and vrms are appropriate only for metric units.\n"
"If distances are measured in feet, then these parameters should be\n"
"specified explicitly.\n"
"\n"
"offmax, tmute, and vrms need not be specified precisely.\n"
"If these values are unknown, then one should overestimate offmax\n"
"and underestimate tmute and vrms.\n"
"\n"
"No muting is actually performed.  The tmute parameter is used only to\n"
"determine parameters required to perform DMO.\n"
"\n"
"Trace header fields accessed:  nt, dt, delrt, offset, cdp.\n"
"\n";
/**************** end self doc *******************************************/

/* Credits:
 *	CWP: Dave
 *
 * Technical Reference:
 *      A non-aliased integral method for dip-moveout
 *      Dave Hale
 *      submitted to Geophysics, June, 1990
 */

void maketa (float dx, float dt, float offmax, float tmute, float vrms,
	int nsmax, int *nsp, float *ts, float *as);
void makeds (int ns, float *ts, float *as, int lds, int ifds, float *ds);
void dmotx (int ns, float *ts, float *as, float offset, float x, float dx,
	int itmute, int nt, float dt, float ft, float *p, float *q);

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 dxcdp;	/* cdp sampling interval */
	int noffmix;	/* number of offsets to mix */
	float offmax;	/* maximum offset */
	float tmute;	/* mute time at far offset */
	float vrms;	/* rms velocity at mute time */
	int nsmax;	/* maximum number of time shifts per trace in DMO */
	int ns;		/* actual number of time shifts per trace in DMO */
	float *p;	/* input trace */
	float **q;	/* output DMO-corrected traces */
	float *temp;	/* temporary array */
	float *ts;	/* table of time shifts for DMO */
	float *as;	/* table of amplitudes for DMO */
	float offset;	/* source-receiver offset of current trace */
	float oldoffset;/* offset of previous trace */
	int cdp;	/* cdp number of current trace */
	int ncdp;	/* number of cdps */
	int icdp;	/* cdp index */
	int jcdp;	/* cdp index */
	int jcdplo;	/* lower bound for jcdp */
	int jcdphi;	/* upper bound for jcdp */
	int ntrace;	/* number of traces processed in current mix */
	int itrace;	/* trace index */
	int noff;	/* number of offsets processed in current mix */
	int gottrace;	/* non-zero if an input trace was read */
	int done;	/* non-zero if done */
	float *ds;	/* shaping filter to complete DMO processing */
	int lds=125;	/* length of shaping filter */
	int ifds=-100;	/* time index of first sample in shaping filter */
	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;

	/* get parameters */
	if (!getparint("cdpmin",&cdpmin)) err("must specify cdpmin");
	if (!getparint("cdpmax",&cdpmax)) err("must specify cdpmax");
	if (cdpmin>cdpmax) err("cdpmin must be less than cdpmax");
	if (!getparfloat("dxcdp",&dxcdp)) err("must specify dxcdp");
	if (!getparint("noffmix",&noffmix)) err("must specify noffmix");
	if (!getparfloat("offmax",&offmax)) offmax=3000.0;
	if (!getparfloat("tmute",&tmute)) tmute=2.0;
	if (!getparfloat("vrms",&vrms)) vrms=1500.0;
	if (!getparint("nsmax",&nsmax)) nsmax=400;
	if (!getparint("verbose",&verbose)) verbose=0;
	
	/* determine number of cdps */
	ncdp = cdpmax-cdpmin+1;

	/* allocate workspace */
	q = ealloc2float(nt,ncdp);
	p = ealloc1float(nt);
	temp = ealloc1float(nt);
	ts = ealloc1float(nsmax);
	as = ealloc1float(nsmax);
	ds = ealloc1float(lds);
	
	/* tabulate time shifts and amplitudes for dmo */
	maketa(dxcdp,dt,offmax,tmute,vrms,nsmax,&ns,ts,as);
	if (verbose) 
		fprintf(stderr,"\tDMO will be performed via %d time shifts\n",
			ns);
	
	/* compute shaping filter for dmo horizontal reflection response */
	makeds(ns,ts,as,lds,ifds,ds);
	
	/* open temporary file for headers */
	/*hfp = tmpfile();*/
	hfp = etempfile(NULL);
	
	/* initialize */
	oldoffset = tr.offset;
	gottrace = 1;
	done = 0;
	ntrace = 0;
	noff = 0;
	for (icdp=0; icdp<ncdp; ++icdp)
		for (it=0; it<nt; ++it)
			q[icdp][it] = 0.0;

	/* loop over traces */
	do {
		
		/* if got a trace */
		if (gottrace) {

			/* determine offset and cdp */
			offset = tr.offset;
			cdp = tr.cdp;
		
			/* update number of offsets mixed */
			if (offset!=oldoffset) noff++;

			/* get trace samples */
			bcopy(tr.data,p,nt*sizeof(float));
		}
		
		/* if a mix of offsets is complete */
		if (noff==noffmix || !gottrace) {
			
			/* update number of offsets mixed */
			if (!gottrace) noff++; 
			
			/* apply shaping filter to complete dmo processing */
			for (icdp=0; icdp<ncdp; ++icdp) {
				conv(lds,ifds,ds,nt,0,q[icdp],nt,0,temp);
				bcopy(temp,q[icdp],nt*sizeof(float));
			}
			
			/* 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);
				icdp = tro.cdp-cdpmin;
				
				/* get dmo-corrected data */
				bcopy(q[icdp],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 */
			noff = 0;
			ntrace = 0;
			
			/* zero dmo accumulators */
			for (icdp=0; icdp<ncdp; ++icdp)
				for (it=0; it<nt; ++it)
					q[icdp][it] = 0.0;
		}
				
		/* if cdp is not within range of cdps to process, skip it */
		if (cdp<cdpmin || cdp>cdpmax) continue;
		
		/* save trace header and update number of traces */
		efwrite(&tr,HDRBYTES,1,hfp);
		ntrace++;
		
		/* determine output traces potentially modified by input */