sudmotx.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* SUDMOTX: $Revision: 1.14 $ ; $Date: 2003/06/09 16:17:07 $		*/

#include "su.h"
#include "segy.h"
#include "header.h"
#include <signal.h>

/*********************** self documentation ******************************/
char *sdoc[] = {
"									",
" SUDMOTX - DMO via T-X domain (Kirchhoff) method for common-offset gathers",
"									",
" sudmotx <stdin >stdout cdpmin= cdpmax= dxcdp= noffmix= [optional parms]",
"									",
" Required Parameters:							",
" cdpmin                  minimum cdp (integer number) for which to apply DMO",
" cdpmax                  maximum cdp (integer number) for which to apply DMO",
" dxcdp                   distance between successive cdps		",
" noffmix                 number of offsets to mix (see notes)		",
"									",
" Optional Parameters:							",
" offmax=3000.0           maximum offset				",
" tmute=2.0               mute time at maximum offset offmax		",
" vrms=1500.0             RMS velocity at mute time tmute		",
" verbose=0               =1 for diagnostic print			",
" tmpdir=	if non-empty, use the value as a directory path	prefix	",
"		for storing temporary files; else if the CWP_TMPDIR	",
"		environment variable is set use	its value for the path;	",
"		else use tmpfile()					",
"									",
"									",
" Notes:								",
" Input traces should be sorted into common-offset gathers.  One common-",
" offset gather ends and another begins when the offset field of the trace",
" headers changes.							",
" 									",
" The cdp field of the input trace headers must be the cdp bin NUMBER, NOT",
" the cdp location expressed in units of meters or feet.		",
"									",
" The number of offsets to mix (noffmix) should typically equal the ratio of",
" the shotpoint spacing to the cdp spacing.  This choice ensures that every",
" cdp will be represented in each offset mix.  Traces in each mix will	",
" contribute through DMO to other traces in adjacent cdps within that mix.",
"									",
" The defaults for offmax and vrms are appropriate only for metric units.",
" If distances are measured in feet, then these parameters should be	",
" specified explicitly.							",
"									",
" offmax, tmute, and vrms need not be specified precisely.		",
" If these values are unknown, then one should overestimate offmax	",
" and underestimate tmute and vrms.					",
"									",
" No muting is actually performed.  The tmute parameter is used only to	",
" determine parameters required to perform DMO.				",
NULL};

/* Credits:
 *	CWP: Dave Hale
 *
 * Technical Reference:
 *      A non-aliased integral method for dip-moveout
 *      Dave Hale
 *      submitted to Geophysics, June, 1990
 *
 * Trace header fields accessed:  ns, dt, delrt, offset, cdp.
 */
/**************** end self doc *******************************************/

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

static void closefiles(void);

/* Globals (so can trap signal) defining temporary disk files */
char headerfile[BUFSIZ];/* filename for the file of headers	*/
FILE *headerfp;		/* fp for header storage file		*/

segy tr,tro;

int
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=0.0;/* source-receiver offset of current trace */
	float oldoffset;/* offset of previous trace */
	int cdp=0;	/* 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 */
	char *tmpdir;	/* directory path for tmp files	*/
	cwp_Bool istmpdir=cwp_false;/* true for user given path */

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

	/* get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = ((double) tr.dt)/1000000.0;
	ft = 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;

	/* Look for user-supplied tmpdir */
	if (!getparstring("tmpdir",&tmpdir) &&
	    !(tmpdir = getenv("CWP_TMPDIR"))) tmpdir="";
	if (!STREQ(tmpdir, "") && access(tmpdir, WRITE_OK))
		err("you can't write in %s (or it doesn't exist)", tmpdir);
	
	/* 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 */
	if (STREQ(tmpdir,"")) {
		headerfp = etmpfile();
		if (verbose) warn("using tmpfile() call");
	} else { /* user-supplied tmpdir */
		char directory[BUFSIZ];
		strcpy(directory, tmpdir);
		strcpy(headerfile, temporary_filename(directory));
		/* Trap signals so can remove temp files */
		signal(SIGINT,  (void (*) (int)) closefiles);
		signal(SIGQUIT, (void (*) (int)) closefiles);
		signal(SIGHUP,  (void (*) (int)) closefiles);
		signal(SIGTERM, (void (*) (int)) closefiles);
		headerfp = efopen(headerfile, "w+");
      		istmpdir=cwp_true;		
		if (verbose)
			warn("putting temporary header file in %s", directory);
	}
	
	/* 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 */
 			memcpy( (void *) p,
				  (const void *) tr.data, 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);
				memcpy( (void *) q[icdp],
					(const void *) temp, nt*sizeof(float));
			}
			
			/* rewind trace header file */
			erewind(headerfp);
			
			/* loop over all output traces */
			for (itrace=0; itrace<ntrace; ++itrace) {
			
				/* read trace header and determine cdp index */
				efread(&tro,HDRBYTES,1,headerfp);
				icdp = tro.cdp-cdpmin;
				
				/* get dmo-corrected data */
				memcpy((void *) tro.data,
				      (const void *) q[icdp],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 */
			erewind(headerfp);
			
			/* 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 within range of cdps to process */