sustolt.c

seismic software,very useful

/* SUSTOLT: $Revision: 1.2 $ ; $Date: 92/10/23 15:35:08 $		*/

/*----------------------------------------------------------------------
 * 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 ******************************/
char *sdoc[] = {
"									",
" SUSTOLT - Stolt migration for stacked data or common-offset gathers	",
"									",
" sustolt <stdin >stdout cdpmin= cdpmax= dxcdp= noffmix= [...]		",
"									",
" 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 adjacent cdp bins (m)	",
"									",
" Optional Parameters:							",
" noffmix=1               number of offsets to mix (for unstacked data only)",
" tmig=0.0                times corresponding to rms velocities in vmig (s)",
" vmig=1500.0             rms velocities corresponding to times in tmig (m/s)",
" smig=1.0                stretch factor (0.6 typical if vrms increasing)",
" vscale=1.0              scale factor to apply to velocities		",
" fmax=Nyquist            maximum frequency in input data (Hz)		",
" lstaper=0               length of side tapers (# of traces)		",
" lbtaper=0               length of bottom taper (# of samples)		",
" verbose=0               =1 for diagnostic print			",
" 									",
" Notes:								",
" If unstacked traces are input, they 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 be specified for	",
" unstacked data only.  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 migration to other traces in adjacent cdps within	",
" that mix.								",
"									",
" The tmig and vmig arrays specify a velocity function of time that is	",
" used to implement Stolt's stretch for depth-variable velocity.  The	",
" stretch factor smig is often referred to as the \"W\" factor.		",
" The times in tmig must be monotonically increasing.			",
"									",
" Trace header fields accessed:  nt, dt, delrt, offset, cdp.		",
NULL};
/**************** end self doc *******************************************/

/* Credits:
 *	CWP: Dave Hale
 */

static void makev (int nmig, float *tmig, float *vmig, float vscale,
	int nt, float dt, float ft, float **v, float *vmin, float *vmax);
static void makeut (float vstolt, float fmax, float *vt,
	int nt, float dt, float **ut, int *nu, float *du, float **tu);
static void makeu (float vstolt, float *v, int nt, float dt, float *u);
static void stolt1k (float k, float v, float s, float fmax, int nt, float dt,
	complex *p, complex *q);
static void taper (int lxtaper, int lbtaper, 
	int nx, int ix, int nt, float *trace);

segy tr,tro;

main(int argc, char **argv)
{
	int nt,it,cdpmin,cdpmax,nx,nxfft,nxpad,ix,noffmix,
		lstaper,lbtaper,ntmig,nvmig,itmig,nk,nu,iu,
		noff,ntrace,itrace,gottrace,done,verbose;
	float dt,ft,dx,*tmig,*vmig,vscale,vstolt,smig,fmax,
		offset,oldoffset,vmin,vmax,dk,du,
		*v,*ut,*tu,**px,**q;
	complex **pk;
	FILE *hfp;

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

	/* get information from the first header */
	if (!gettr(&tr)) err("cannot get first trace");
	nt = tr.ns;
	dt = (float)tr.dt/1000000.0;
	ft = tr.delrt/1000.0;
	if (ft!=0.0) err("cannot handle non-zero time of first sample");
	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)) noffmix = 1;
	ntmig = countparval("tmig");
	if (ntmig==0) ntmig = 1;
	tmig = ealloc1float(ntmig);
	if (!getparfloat("tmig",tmig)) tmig[0] = 0.0;
	nvmig = countparval("vmig");
	if (nvmig==0) nvmig = 1;
	if (nvmig!=ntmig) err("number of tmig and vmig must be equal");
	vmig = ealloc1float(nvmig);
	if (!getparfloat("vmig",vmig)) vmig[0] = 2000.0;
	for (itmig=1; itmig<ntmig; ++itmig)
		if (tmig[itmig]<=tmig[itmig-1])
			err("tmig must increase monotonically");
	if (!getparfloat("smig",&smig)) smig = 1.0;
	if (!getparfloat("fmax",&fmax)) fmax = 0.5/dt;
	fmax = MIN(fmax,0.5/dt);
	if (!getparfloat("vscale",&vscale)) vscale = 1.0;
	if (!getparint("lstaper",&lstaper)) lstaper=0;
	if (!getparint("lbtaper",&lbtaper)) lbtaper=0;
	if (!getparint("verbose",&verbose)) verbose=0;
	
	/* make uniformly sampled rms velocity function of time */
	makev(ntmig,tmig,vmig,vscale,nt,dt,ft,&v,&vmin,&vmax);
	
	/* Stolt migration velocity is the minimum velocity */
	vstolt = vmin;

	/* make u(t) and t(u) for Stolt stretch */
	makeut(vstolt,fmax,v,nt,dt,&ut,&nu,&du,&tu);
	free1float(v);
	
	/* determine number of cdps to process */
	nx = cdpmax-cdpmin+1;

	/* wavenumber (k) sampling */
	nxpad = 0.5*vmax*nt*dt/dx;
	nxfft = npfar(nx+nxpad);
	nk = nxfft/2+1;
	dk = 2.0*PI/(nxfft*dx);
		
	/* allocate and zero common-offset gather p(t,x) */
	pk = alloc2complex(nu,nk);
	px = alloc1(nxfft,sizeof(float*));
	px[0] = (float*)pk[0];
	for (ix=1; ix<nxfft; ++ix)
		px[ix] = px[0]+ix*nu;
	for (ix=0; ix<nxfft; ++ix)
		for (iu=0; iu<nu; ++iu)
			px[ix][iu] = 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();
	
	/* 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, migrate it */
		if ((gottrace && offset!=oldoffset) || !gottrace) {
		
			/* local variables */
			int ik,jx,iu;
			float scale=1.0/nxfft,k;
			
			/* apply side and bottom tapers */
			for (ix=0; ix<nx; ++ix)
				taper(lstaper,lbtaper,nx,ix,nt,px[ix]);
			
			/* if necessary, stretch */
			if (nu!=nt && tu!=NULL) {
				float *temp=ealloc1float(nu);
				for (ix=0; ix<nx; ++ix) {
					ints8r(nt,dt,0.0,px[ix],0.0,0.0,
						nu,tu,temp);
					for (iu=0; iu<nu; ++iu)
						px[ix][iu] = temp[iu];
				}
				free1float(temp);
			}

			/* Fourier transform p(u,x) to p(u,k) */
			pfa2rc(-1,2,nu,nxfft,px[0],pk[0]);

			/* migrate each wavenumber */
			for (ik=1,k=dk; ik<nk; ++ik,k+=dk)
				stolt1k(k,vstolt,smig,fmax,
					nu,du,pk[ik],pk[ik]);

			/* Fourier transform p(u,k) to p(u,x) and scale */
			pfa2cr(1,2,nu,nxfft,pk[0],px[0]);
			for (jx=0; jx<nx; ++jx)
				for (iu=0; iu<nu; ++iu)
					px[jx][iu] *= scale;

			/* if necessary, unstretch */
			if (nu!=nt && ut!=NULL) {
				float *temp=ealloc1float(nu);
				for (ix=0; ix<nx; ++ix) {
					ints8r(nu,du,0.0,px[ix],0.0,0.0,
						nt,ut,temp);
					for (it=0; it<nt; ++it)
						px[ix][it] = temp[it];
				}
				free1float(temp);
			}
			
			/* add migrated traces to mix */
			for (ix=0; ix<nx; ++ix)
				for (it=0; it<nt; ++it)
					q[ix][it] += px[ix][it];
			
			/* zero common-offset gather */
			for (ix=0; ix<nxfft; ++ix)
				for (iu=0; iu<nu; ++iu)
					px[ix][iu] = 0.0;
			
			/* count offsets in mix */
			noff++;
		}
		
		/* 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);
				
				/* index of cdp, zero-based */
				ix = tro.cdp-cdpmin;
				
				/* get migrated data */
				for (it=0; it<nt; ++it)
					tro.data[it] = q[ix][it];
				
				/* 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 within range to process */
		if (tr.cdp>=cdpmin && tr.cdp<=cdpmax) {
	
			/* save trace header and update number of traces */
			efwrite(&tr,HDRBYTES,1,hfp);
			ntrace++;

			/* remember offset */
			oldoffset = offset;
		
			/* index of cdp, zero-based */
			ix = tr.cdp-cdpmin;

			/* save trace */
			for (it=0; it<nt; ++it)
				px[ix][it] = tr.data[it];
		}

		/* get next trace (if there is one) */
		if (!gettr(&tr)) gottrace = 0;
		
	} while (!done);

	return EXIT_SUCCESS;
}
	
static void makev (int nmig, float *tmig, float *vmig, float vscale,