sumigprepspi.c

su 的源代码库

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

/* SUMIGPREPSPI: $Vision: 1.00 $ ; $Date: 2006/11/07 22:58:42 $*/

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


/*********************** self documentation ******************************/
char *sdoc[] = {
"									",
" SUMIGPREPSPI --- The 2-D PREstack commom-shot Phase-Shift-Plus 	",
"			interpolation depth MIGration.			",
"									",
"   sumigprepspi <indata >outfile [parameters] 				", 
"									",
" Required Parameters:						   	",
"									",
" nxo=	   number of total horizontal output samples			",
" nxshot=	number of shot gathers to be migrated		   	",
" nz=	    number of depth sapmles				 	",
" dx=	    horizontal sampling interval			  	",   
" dz=	    depth sampling interval				 	",
" vfile=	 velocity profile, it must be binary format.		",
" 									",
" Optional Parameters:						   	",
" fmax=25	the peak frequency of Ricker wavelet used as source wavelet",
" f1=5,f2=10,f3=40,f4=50	 frequencies to build a Hamming window  ",
" lpad=9999,rpad=9999	    number of zero traces padded on both	",
" 			    sides of depth section to determine the	",
"			    migration aperature, the default values	",
"			    are using the full aperature.		",
" verbose=0		silent, =1 additional runtime information	",
"									",  
" Notes:								",
" The input velocity file \'vfile\' consists of C-style binary floats.	",  
" The structure of this file is vfile[iz][ix]. Note that this means that",
" the x-direction is the fastest direction instead of z-direction! Such a",
" structure is more convenient for the downward continuation type	",
" migration algorithm than using z as fastest dimension as in other SU  ", 
" programs.						   		",
"									",
" Because most of the tools in the SU package (such as  unif2, unisam2, ", 
" and makevel) produce output with the structure vfile[ix][iz], you will",
" need to transpose the velocity files created by these programs. You may",
" use the SU program \'transp\' in SU to transpose such files into the  ",
" required vfile[iz][ix] structure.					",
"									",
" Also, sx must be monotonically increasing throughout the dataset, and ",
" and gx must be monotonically increasing within a shot. You may resort ",
" your data with \'susort\', accordingly.				",
"									",
" The scalco header field is honored so this field must be set correctly.",
" See selfdocs of \'susort\', \'suchw\'. Also:   sukeyword scalco	",
"									",
NULL};

/*
 * Credits: CWP, Baoniu Han, [e-mail], April 19th, 1998
 *	  Modified: Chris Stolk, 11 Dec 2005, - changed data input
 *		    to remove erroneous time delay.
 *	  Modified: CWP, John Stockwell 26 Sept 2006 - replaced Han's
 *	  "goto-loop" in two places with "do { }while loops".
 *	  Fixed it so that sx, gx, and scalco are honored.
 *
 *
 * Trace header fields accessed: ns, dt, delrt, d2
 * Trace header fields modified: ns, dt, delrt
 */

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

/* Prototypes for subroutines used internally */
float *ricker(float Freq,float dt,int *Npoint);
void get_sx_gx(float *sx, float *gx);


segy tr;


int
main (int argc, char **argv)
{
	int L=10, Bz=0;	/*BIN number for determination of reference velocity.*/
	float c[101];
	float *V,P[100];
	float Sz=0;
	float Y[101];
	double a,a1,a2,theta,theta1,theta2;


	int nt;			/* number of time samples		*/
	int nz;			/* number of migrated depth samples	*/
	int nx;			/* number of horizontal samples	 */
	
	int nxshot;		/* number of shots to be migrated	*/
	int nxshot_orig;	/* first value of nxshot		*/
	int iz,iw,ix,it,ik;	/* loop counters 			*/
	int igx;		/* integerized gx value			*/
	int ntfft,nxfft;	/* fft size				*/
	int nw,truenw,nk;	/* number of wave number and frequency 	*/	

	int dip=45;		/* dip angle				*/

	float sx,gx;		/* x source and geophone location	*/
	float gxmin=0.0,gxmax=0.0;/* x source and geophone location	*/
        float min_sx_gx;        /* min(sx,gx)                           */

	float oldgx;		/* old gx position			*/
	float oldgxmin;		/* old gx position			*/
	float oldgxmax;		/* old gx position			*/
	float oldsx=0.0;	/* old sx position			*/
	int oldigx=0;		/* old value of integerized gx value	*/
	int oldisx=0;		/* old value of integerized sx value	*/

	int isx=0,nxo;		/* index for source and geophone	*/
	int ix1,ix2,ix3,ixshot,il,ir;/* dummy index			*/
	int lpad,rpad;	/* padding on both sides of the migrated section*/

	float *wl=NULL,*wtmp=NULL;
	float fmax;
	float f1,f2,f3,f4;
	int nf1,nf2,nf3,nf4;
	int ntw;

	float dt=0.004,dz;	/* time and depth sampling intervals 	*/
	float dw,dk;		/* wavenumber and frequency sampling interval */
	float fw,fk;		/* first wavenumber and frequency	*/
	float w,k;		/* wavenumber and frequency		*/
	float dx;		/* spatial sampling interval		*/
	float **p=NULL;
	float **cresult=NULL;	/* input, output data			*/

	float vmin,v1,v2,vmax;
	double kz1;
	double phase1;
	float **v=NULL,**vp=NULL;
	complex cshift1;
	complex *wlsp=NULL;
	complex **cp=NULL;
	complex **cp1=NULL;
	complex **cq=NULL;
	complex **cq1=NULL;
	complex ***cq2=NULL;
	complex ***cq3=NULL;	/* complex input,output			*/

	char *vfile="";		/* name of file containing velocities	*/
	FILE *vfp=NULL;		/* ... its file pointer			*/

	int verbose;		/* verbose flag				*/

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(1);

	/* get required parameters */
	MUSTGETPARINT("nz",&nz);
	MUSTGETPARFLOAT("dz",&dz);
	MUSTGETPARSTRING("vfile", &vfile);
	MUSTGETPARINT("nxo",&nxo);
	MUSTGETPARINT("nxshot",&nxshot);

	/* get optional parameters */
	if (!getparfloat("fmax",&fmax)) fmax = 25.0;  
	if (!getparfloat("f1",&f1)) f1 = 10.0;
	if (!getparfloat("f2",&f2)) f2 = 20.0;
	if (!getparfloat("f3",&f3)) f3 = 40.0;
	if (!getparfloat("f4",&f4)) f4 = 50.0;

	if (!getparint("lpad",&lpad)) lpad=9999;
	if (!getparint("rpad",&rpad)) rpad=9999;
	if (!getparint("dip",&dip)) dip=65;

	if (!getparint("verbose",&verbose))     verbose = 0;

	/* allocate space */
	cresult = alloc2float(nz,nxo);
	vp=alloc2float(nxo,nz);

	/* load velocity file */
	vfp=efopen(vfile,"r");
	efread(vp[0],FSIZE,nz*nxo,vfp);
	efclose(vfp);

	/* zero out cresult array */
	memset((void *) cresult[0], 0, nxo*nz*FSIZE);

	/* save value of nxshot */
	nxshot_orig=nxshot;


		/* get info from first trace */
	if (!gettr(&tr))  err("can't get first trace");
	nt = tr.ns;
	get_sx_gx(&sx,&gx);
	min_sx_gx = MIN(sx,gx);
	sx = sx - min_sx_gx;
	gx = gx - min_sx_gx;


	/* let user give dt and/or dx from command line */
	if (!getparfloat("dt", &dt)) {
		if (tr.dt) { /* is dt field set? */
			dt = ((double) 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");
		}
	}

	do {    /* begin loop over shots */

 
		/* determine frequency sampling interval */
		ntfft = npfar(nt);
		nw = ntfft/2+1;
		dw = 2.0*PI/(ntfft*dt);

		/* compute the index of the frequency to be migrated */
		fw=2.0*PI*f1;
		nf1=fw/dw+0.5;
		 
		fw=2.0*PI*f2;
		nf2=fw/dw+0.5;

		fw=2.0*PI*f3;
		nf3=fw/dw+0.5;

		fw=2.0*PI*f4;
		nf4=fw/dw+0.5;  

		/* the number of frequency to migrated */
		truenw=nf4-nf1+1;
		fw=0.0+nf1*dw;
		if (verbose)
			warn("nf1=%d nf2=%d nf3=%d nf4=%d nw=%d",nf1,nf2,nf3,nf4,truenw);

		/* allocate space */
		wl=alloc1float(ntfft);
		wlsp=alloc1complex(nw);

		/* generate the Ricker wavelet */
		wtmp=ricker(fmax,dt,&ntw);

		/* zero out wl[] array */
		memset((void *) wl, 0, ntfft*FSIZE);
	
		/* CHANGE BY CHRIS STOLK, Dec. 11, 2005 */
		/* The next two lines are the old code, */
		/* it is erroneous because the peak of  */
		/* the wavelet occurs at positive time  */
		/* instead of time zero. 		*/
		/*
		for(it=0;it<ntw;it++)
	  		wl[it]=wtmp[it];
		*/
		/* New code: we put in the wavelet in a centered fashion */ 
		for(it=0;it<ntw;it++) 
	  		wl[(it-ntw/2+ntfft) % ntfft]=wtmp[it];
		/* End of new code */

		free1float(wtmp);

		/* fourier transform in time */
		pfarc(-1,ntfft,wl,wlsp);

		/* allocate space */
		p = alloc2float(ntfft,nxo);
		cq = alloc2complex(nw,nxo);

		/* zero out p[][] array */
		memset((void *) p[0], 0, ntfft*nxo*FSIZE);

		/* initialize a number of items before looping over traces */
		nx = 0;
		igx=0;
		oldsx=sx;
		oldgx=gx;
		oldgxmax=gxmax;
		oldgxmin=gxmin;
		do { /* begin looping over traces within a shot gather */

			igx = gx/dx;
			if (igx==oldigx) 
			   warn("repeated igx!!! check dx or scalco value!!!");
			oldigx = igx;

			memcpy( (void *) p[igx], (const void *) tr.data,nt*FSIZE);

			if(gxmin>gx)gxmin=gx;
			if(gxmax<gx)gxmax=gx;

			if(verbose)
				warn(" inside loop:  min_sx_gx %f isx %d igx %d gx %f sx %f",min_sx_gx,isx,igx,gx,sx);
			/* get sx and gx */
			get_sx_gx(&sx,&gx);
			sx = (sx - min_sx_gx);
			gx = (gx - min_sx_gx);

			/* sx, gx must increase monotonically */
			if (!(oldsx <= sx) )
			 err("sx field must be monotonically increasing!");
			if (!(oldgx <= gx) )
			 err("gx field must be monotonically increasing!");

			++nx;
		} while(gettr(&tr) && sx==oldsx);
			
		isx=oldsx/dx;
		if (isx==oldisx) 
			warn("repeated isx!!! check dx or scalco value!!!");
		oldisx=isx;
		ixshot=isx;
		if(verbose) {
			warn("sx %f, gx %f , gxmin %f  gxmax %f nx %d",sx,gx,gxmin,gxmax, nx);
			warn("isx %d igx %d ixshot %d" ,isx,igx,ixshot);
		}



		/* transform the shot gather from time to frequency domain */
		pfa2rc(1,1,ntfft,nxo,p[0],cq[0]);


		/* compute the most left and right index for the migrated */
		/* section*/
		ix1=oldsx/dx;
		ix2=oldgxmin/dx;
		ix3=oldgxmax/dx;

		if(ix1>=ix3)ix3=ix1;
		if(ix1<=ix2)ix2=ix1;
		il=ix2;
		ir=ix3;

		ix2-=lpad;
		ix3+=rpad;
		if(ix2<0)ix2=0;
		if(ix3>nxo-1)ix3=nxo-1;

		/* the total traces to be migrated */
		nx=ix3-ix2+1;
		nw=truenw;

		/* determine wavenumber sampling (for complex to complex FFT) */
		nxfft = npfa(nx);
		nk = nxfft;
		dk = 2.0*PI/(nxfft*dx);
		fk = -PI/dx;

		/* allocate space for velocity profile within the aperature */
		v=alloc2float(nx,nz);   
	
		for(iz=0;iz<nz;iz++)
			for(ix=0;ix<nx;ix++)
				v[iz][ix]=vp[iz][ix+ix2];