sumigpsti.c

su 的源代码库

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

/* SUMIGPSTI: $Revision: 1.2 $ ; $Date: 2003/06/09 16:17:07 $		*/

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

/*********************** self documentation ******************************/
char *sdoc[] = {
"									",
" SUMIGPSTI - MIGration by Phase Shift for TI media with turning rays	",
"									",
" sumigpsti <stdin >stdout [optional parms]				",
"									",
" Required Parameters:							",
" 	None								",
"									",
" Optional Parameters:							",
" dt=from header(dt) or .004	time sampling interval			",
" dx=from header(d2) or 1.0	distance between sucessive cdp's	",
" ffil=0,0,0.5/dt,0.5/dt  trapezoidal window of frequencies to migrate	",
" tmig=0.0	times corresponding to interval velocities in vmig	",
" vnmig=1500.0	interval NMO velocities corresponding to times in tmig	",
" vmig=1500.0	interval velocities corresponding to times in tmig	",
" etamig=0.0	interval eta values corresponding to times in tmig	",
" vnfile=	binary (non-ascii) file containing NMO velocities vn(t)	",
" vfile=	binary (non-ascii) file containing velocities v(t)	",
" etafile=	binary (non-ascii) file containing eta values eta(t)	",
" nxpad=0	number of cdps to pad with zeros before FFT		",
" ltaper=0	length of linear taper for left and right edges		", 
" verbose=0	=1 for diagnostic print					",
"									",
" Notes:								",
" Input traces must be sorted by either increasing or decreasing cdp.	",
"									",
" The tmig, vnmig, vmig and etamig arrays specify an interval values	",
" function of time. Linear interpolation and constant extrapolation is	",
" used to determine interval velocities at times not specified.  Values	",
" specified in tmig must increase monotonically.			",
" Alternatively, interval velocities may be stored in a binary file	",
" containing one velocity for every time sample.  If vnfile is specified,",
" then the tmig and vnmig arrays are ignored.				",
"									",
" The time of first sample is assumed to be zero, regardless of the value",
" of the trace header field delrt.					",
"									",
" Trace header fields accessed:  ns and dt				",
"									",
NULL};

/**************** end self doc *******************************************/

/* Credits:
 *	CWP: Dave Hale (originally call supsmig.c)
 *		modified to TI media by Tariq Alkhalifah
 */

typedef struct TATableStruct {
	int np;
	float *p;
	float *taumax;
	float *taumin;
	float *tturn;
	float **t;
	float **a;
} TATable;

void mig1k (TATable *table, float k, float ffil[4], int ntflag,
	int nt, float dt, complex *cp,
	int ntau, float dtau, complex *cq);
TATable *tableta (int np, int ntau, float dtau, float a3333[],
	float a1111[], float a1133[], float a1313[],
	int nt, float dt, int verbose);

segy tr;

int
main(int argc, char **argv)
{
	int nt,nx,nxfft,nxpad,ix,it,nk,ik,nfil,ntflag,
		ltaper,ntmig,nvmig,nvnmig,netamig,itmig,verbose,np;
	float dt,dx,dk,taper,t,k,fftscl,ffil[4],
		*a1111,*a3333,*a1313,*a1133,
		*tmig,*vmig,*vnmig,*etamig,*vt,*vnt,*etat,**gtx;
	float vnt2,vt2,delta;
	complex **gtk;
	char *vfile="",*vnfile="",*etafile="";
	FILE *hfp,*tfp;
	TATable *table;

	/* 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;

	/* let user give dt and/or dx from command line */
	if (!getparfloat("dt", &dt)) {
		if (tr.dt) { /* is dt field set? */
			dt = (float) tr.dt / 1000000.0;
		} else { /* dt not set, assume 4 ms */
			dt = 0.004;
			warn("tr.dt not set, assuming dt=0.004");
		}
	}
	if (!getparfloat("dx",&dx)) {
		if (tr.d2) { /* is d2 field set? */
			dx = tr.d2;
		} else {
			dx = 1.0;
			warn("tr.d2 not set, assuming d2=1.0");
		}
	}


	/* get parameters */
	if (!(nfil=getparfloat("ffil",ffil))) {
		ffil[0] = ffil[1] = 0.0;
		ffil[2] = ffil[3] = 0.5/dt;
	} else if (nfil!=4) {
		err("if ffil specified, exactly 4 values must be provided");
	}
	if (!getparint("nxpad",&nxpad)) nxpad=0;
	if (!getparint("ltaper",&ltaper)) ltaper=0;
	if (!getparint("verbose",&verbose)) verbose=0;
	if (!getparint("np",&np)) np=200;
	if (!getparint("ntflag",&ntflag)) ntflag=1;

	/* determine NMO velocity function vnm(t) */
	vnt = ealloc1float(nt);
        if (!getparstring("vnfile",&vnfile)) {
                ntmig = countparval("tmig");
                if (ntmig==0) ntmig = 1;
                tmig = ealloc1float(ntmig);
                if (!getparfloat("tmig",tmig)) tmig[0] = 0.0;
                nvnmig = countparval("vnmig");
                if (nvnmig==0) nvnmig = 1;
                if (nvnmig!=ntmig) err("number of tmig and vnmig must be equal");
                vnmig = ealloc1float(nvnmig);
                if (!getparfloat("vnmig",vnmig)) vnmig[0] = 1500.0;
                for (itmig=1; itmig<ntmig; ++itmig)
                        if (tmig[itmig]<=tmig[itmig-1])
                                err("tmig must increase monotonically");
                for (it=0,t=0.0; it<nt; ++it,t+=dt)
                        intlin(ntmig,tmig,vnmig,vnmig[0],vnmig[ntmig-1],
                                1,&t,&vnt[it]);
        } else {
                if (fread(vnt,sizeof(float),nt,fopen(vnfile,"r"))!=nt)
                        err("cannot read %d velocities from file %s",nt,vnfile);
        }

	/* determine function eta(t) */
	etat = ealloc1float(nt);
        if (!getparstring("etafile",&etafile)) {
                ntmig = countparval("tmig");
                if (ntmig==0) ntmig = 1;
                tmig = ealloc1float(ntmig);
                if (!getparfloat("tmig",tmig)) tmig[0] = 0.0;
                netamig = countparval("etamig");
                if (netamig==0) netamig = 1;
                if (netamig!=ntmig) err("number of tmig and etamig must be equal");
                etamig = ealloc1float(netamig);
                if (!getparfloat("etamig",etamig)) etamig[0] = 0.0;
                for (itmig=1; itmig<ntmig; ++itmig)
                        if (tmig[itmig]<=tmig[itmig-1])
                                err("tmig must increase monotonically");
                for (it=0,t=0.0; it<nt; ++it,t+=dt)
                        intlin(ntmig,tmig,etamig,etamig[0],etamig[ntmig-1],
                                1,&t,&etat[it]);
        } else {
                if (fread(etat,sizeof(float),nt,fopen(etafile,"r"))!=nt)
                        err("cannot read %d velocities from file %s",nt,etafile);
        }
	
	/* determine velocity function v(t) */
	vt = ealloc1float(nt);
	if (!getparstring("vfile",&vfile)) {
		ntmig = countparval("tmig");
		if (ntmig==0) ntmig = 1;
		tmig = ealloc1float(ntmig);
		if (!getparfloat("tmig",tmig)) tmig[0] = 0.0;
		nvmig = countparval("vmig");
		if (nvmig!=ntmig) {
			for(it=0; it<nt; ++it)
				vt[it] = vnt[it];
		} else {
			vmig = ealloc1float(nvmig);
			if (!getparfloat("vmig",vmig)) vmig[0] = vnt[0];
			for (itmig=1; itmig<ntmig; ++itmig)
				if (tmig[itmig]<=tmig[itmig-1])
					err("tmig must increase monotonically");
			for (it=0,t=0.0; it<nt; ++it,t+=dt)
				intlin(ntmig,tmig,vmig,vmig[0],vmig[ntmig-1],
					1,&t,&vt[it]);
		}
	} else {
		if (fread(vt,sizeof(float),nt,fopen(vfile,"r"))!=nt)
			err("cannot read %d velocities from file %s",nt,vfile);
	}

	/* allocate space for velocity arrays */
	a1111 = ealloc1float(nt);
	a3333 = ealloc1float(nt);
	a1313 = ealloc1float(nt);
	a1133 = ealloc1float(nt);

	/*convert eta and NMO velocity values to elastic coeffecients*/
	for(it=0; it<nt; ++it){
		vnt2   =vnt[it]*vnt[it];
		vt2    =vt[it]*vt[it];
		a1111[it]=vnt2*(1+2*etat[it]);
		a3333[it]=vt2;
		a1313[it]=0.25*vt2;
		delta= 0.5*(vnt2/vt2-1);
		a1133[it]=sqrt(2*delta*a3333[it]*
			(a3333[it]-a1313[it])+(a3333[it]-
			a1313[it])*(a3333[it]-a1313[it]))-a1313[it];
	}
	
	/* copy traces and headers to temporary files */
	tfp = tmpfile();
	hfp = tmpfile();
	nx = 0;
	do {
		nx++;
		fwrite(&tr,HDRBYTES,1,hfp);
		fwrite(tr.data,sizeof(float),nt,tfp);
	} while(gettr(&tr));
	efseeko(hfp,(off_t) 0,SEEK_SET);
	efseeko(tfp,(off_t) 0,SEEK_SET);
	if (verbose) fprintf(stderr,"\t%d traces input\n",nx);
	
	/* determine wavenumber sampling */
	nxfft = npfaro(nx+nxpad,2*(nx+nxpad));
	nk = nxfft/2+1;
	dk = 2.0*PI/(nxfft*dx);

	/* allocate space for Fourier transform */
	gtk = ealloc2complex(nt,nk);
	gtx = ealloc1(nxfft,sizeof(float*));
	for (ix=0; ix<nxfft; ++ix)
		gtx[ix] = (float*)gtk[0]+ix*nt;

	/* read and apply fft scaling to traces and pad with zeros */
	fftscl = 1.0/nxfft;
	for (ix=0; ix<nx; ++ix) {
		efread(gtx[ix],sizeof(float),nt,tfp);
		for (it=0; it<nt; ++it)
			gtx[ix][it] *= fftscl;
		if (ix<ltaper) {
			taper = (float)(ix+1)/(float)(ltaper+1);
			for (it=0; it<nt; ++it)
				gtx[ix][it] *= taper;
		} else if (ix>=nx-ltaper) {
			taper = (float)(nx-ix)/(float)(ltaper+1);
			for (it=0; it<nt; ++it)
				gtx[ix][it] *= taper;
		}
	}
	for (ix=nx; ix<nxfft; ++ix)
		for (it=0; it<nt; ++it)
			gtx[ix][it] = 0.0;
	
	/* Fourier transform g(t,x) to g(t,k) */
	pfa2rc(-1,2,nt,nxfft,gtx[0],gtk[0]);
	if (verbose) fprintf(stderr,"\tFourier transform done\n");
	
	/* build time/amplitude table */
	table = tableta(np,nt,dt,a3333,a1111,a1133,a1313,nt,dt,verbose);
	if (verbose) fprintf(stderr,"\tTime/amplitude table built\n");

	/*free space*/
	free1float(vt);
	free1float(vnt);
	free1float(etat);
	free1float(a1111);
	free1float(a3333);
	free1float(a1133);
	free1float(a1313);
	
	/* loop over wavenumbers */
	for (ik=0,k=0.0; ik<nk; ++ik,k+=dk) {
	
		/* report */
		if (verbose && ik%(nk/10>0?nk/10:1)==0)
			fprintf(stderr,"\t%d of %d wavenumbers done\n",
				ik,nk);
		
		/* migrate */
		mig1k(table,k,ffil,ntflag,nt,dt,gtk[ik],nt,dt,gtk[ik]);
	}
	
	/* Fourier transform g(t,k) to g(t,x) */
	pfa2cr(1,2,nt,nxfft,gtk[0],gtx[0]);
	if (verbose) fprintf(stderr,"\tinverse Fourier transform done\n");
	
	/* output migrated traces with headers */
	for (ix=0; ix<nx; ++ix) {
		efread(&tr,HDRBYTES,1,hfp);
		memcpy((void *) tr.data,
				(const void *) gtx[ix], nt*sizeof(float));
		puttr(&tr);
	}

	return(CWP_Exit());
}

void mig1k (TATable *table, float k, float ffil[4], int ntflag,
	int nt, float dt, complex *cp,
	int ntau, float dtau, complex *cq)
/*****************************************************************************
phase shift migration for one wavenumber k via tabulated times and amps
******************************************************************************
Input:
table	 	pointer to the TATable
k		wavenumber
ffil[4]		four frequencies (in Hz) defining trapezoidal filter
ntflag		=1 for normal, 2 for turned, 3 for normal+turned
nt		number of time samples
dt		time sampling interval
cp		array[nt] of data to be migrated
ntau		number of vertical two-way times tau
dtau		vertical time sampling interval

Output:
cq		array[ntau] of migrated data
******************************************************************************/
{
	int ntfft,nw,iwnyq,iwmin,iwmax,iw,itau,it,np,
		ip1,ip2,iwfil0,iwfil1,iwfil2,iwfil3,
		itaumax,itaumin,itab;
	float dw,wnyq,wmin,wmax,w,ampf,p1,p2,a1,a2,wa1,wa2,
		taumaxi,taumini,kow,phst,pwsr,pwsi,pwdr,pwdi,
		ampi,phsi,phsn,
		*p,*taumax,*taumin,*tturn,**t,**a;
	complex *pt,*pw,*qn,*qt;
	static int tabbed=0;
	static float fntab,*ctab,*stab,opi2=1.0/(PI*2.0);

	/* if not already built, build cosine/sine tables */
	if (!tabbed) {
		int itab,ntab=1025;
		float angle,dangle=2.0*PI/(ntab-1);
		ctab = alloc1float(ntab);
		stab = alloc1float(ntab);
		for (itab=0,angle=0.0; itab<ntab; ++itab,angle+=dangle) {
			ctab[itab] = cos(angle);
			stab[itab] = sin(angle);
		}
		tabbed = 1;