migprefd.c

用于石油地震资料数字处理

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

/* SUMIGPREFD: $Vision: 1.00 $ ; $Date: 2004/12/23 22:41:54 $       */

#include "su.h"
#include "segy.h"
#include "header.h"
#include <signal.h>
/* #include <time.h> */


/*********************** self documentation ******************************/
char *sdoc[] = {
"                                                                       ",
" SUMIGPREFD - The 2-D prestack common-shot 45-90 degree   		",
"                     finite-difference migration.       		",
"   sumigprefd <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:                                                   ",
"                                                                       ",
" dip=79       the maximum dip to migrate, it can be 45,65,79,80,87,89,90", 
"	      The computation cost is 45=65=79<80<87<89<90		",
"                                                                       ",
" 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.           ",
"                                                                       ",
" 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.                                     ",

NULL};

/*
 * Credits: CWP, Baoniu Han, [e-mail], April 19th, 1998
 *
 *
 * 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 retris(complex *data,complex *a,complex *c,complex *b,complex
                endl,complex endr, int nx, complex *d);

void fdmig( complex **cp, int nx, int nw, float *v,float fw,float
                dw,float dz,float dx,float dt,int dip);
segy tr;


/* static time_t t1,t2; */

int main (int argc, char **argv)
{
	int nt;			/* number of time samples */
	int nz;			/* number of migrated depth samples */
	int nx,nxshot,oldsx;	/* number of midpoints 	*/
	int iz,iw,ix,it;	/* loop counters 	*/
	int ntfft;		/* fft size		*/
	int nw,truenw;		/* number of wave numbers */	
	int dip=79;		/* dip angle	*/

	/* prestack goes here */
	float sx,gxmin,gxmax;	/* source and geophone location	*/
	int isx,nxo,ifx=0;	/* index for source and geophone*/
	int ix1,ix2,ix3,ixshot; /* dummy index			*/
	int lpad,rpad; /* padding on both sides of the migrated section */
	int flag=1;		/* flag control for feet or meters */

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

	float dt=0.004,dz;	/* time and depth sampling interval 	*/
	float dw;		/* frequency  sampling interval		*/
	float fw;		/* first frequency 			*/
	float w;		/* frequency				*/
	float dx;		/* spatial sampling interval		*/
	float **p,**cresult;	/* input, output data			*/
	float v1;		/* average velocity			*/
	double kz2;	
	float **v,**vp;		/* pointer for the velocity profile	*/
	complex cshift2;
	complex *wlsp,**cp,**cp1,**cq;	/*complex input,output*/
	char *vfile="";		/* name of file containing velocities */
	FILE *vfp;

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

        /* get optional parameters */
        if (!getparint("nz",&nz)) err("nz must be specified");
        if (!getparfloat("dz",&dz)) err("dz must be specified");
        if (!getparstring("vfile", &vfile)) err("vfile must be specified");
        if (!getparint("nxo",&nxo)) err("nxo must be specified");
        if (!getparint("nxshot",&nxshot)) err("nshot must be specified");
        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("flag",&flag)) flag=1;
        if (!getparint("dip",&dip)) dip=79;

	cresult = alloc2float(nz,nxo);
	vp=alloc2float(nxo,nz);

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

	for(ix=0;ix<nxo;ix++)
	for(iz=0;iz<nz;iz++)
	cresult[ix][iz]=0.0;
			

	/* get info from first trace */
loop:
/*	time(&t1); */

	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 = ((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");
		}
	}

	sx=tr.sx;
	isx=sx/dx;
	gxmin=gxmax=tr.gx;
	oldsx=sx;
 
        /* 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;
        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);

        for(it=0;it<ntfft;it++)
        wl[it]=0.0;  
        
        for(it=0;it<ntw;it++)
        wl[it]=wtmp[it];
        free1float( wtmp);

        pfarc(-1,ntfft,wl,wlsp);

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

        for (ix=0; ix<nxo; ix++)
                for (it=0; it<ntfft; it++)
                        p[ix][it] = 0.0;


        /*read in a single shot gather*/
        if (tr.gx < 0 ) {
                ix=tr.gx/dx + nxo;
        } else {
                ix=tr.gx/dx ;
        }


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

	nx = 0;

        while(gettr(&tr)){
                        int igx=0;

                        if(tr.sx!=oldsx){ efseeko(stdin,(off_t)(-240-nt*4),SEEK_CUR); break;}
        		if (tr.gx < 0 ) {
                		igx=tr.gx/dx + nxo;
        		} else {
                		igx=tr.gx/dx ;
        		}

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

                        if(gxmin>tr.gx)gxmin=tr.gx;
                        if(gxmax<tr.gx)gxmax=tr.gx;
                        nx++;
                        oldsx=tr.sx;
                        }

        warn("sx %f , gxmin %f  gxmax %f",sx,gxmin,gxmax);

        /*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=sx/dx;
        ix2=gxmin/dx;
        ix3=gxmax/dx;

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

        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;

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


        /*allocate space*/
        cp = alloc2complex(nx,nw);
        cp1 = alloc2complex(nx,nw);

        /*transpose the frequency domain data from data[ix][iw] to data[iw][ix] and
        apply a Hamming at the same time*/

        for (ix=0; ix<nx; ix++)
        for (iw=0; iw<nw; iw++){

        float tmpp=0.0,tmppp=0.0;

	if(iw>=(nf1-nf1)&&iw<=(nf2-nf1)){
	tmpp=PI/(nf2-nf1);tmppp=tmpp*(iw-nf1)-PI;tmpp=0.54+0.46*cos(tmppp);
        cp[iw][ix]=crmul(cq[ix+ix2][iw+nf1],tmpp);}
        else{
	if(iw>=(nf3-nf1)&&iw<=(nf4-nf1)){
	tmpp=PI/(nf4-nf3);tmppp=tmpp*(iw-nf3);tmpp=0.54+0.46*cos(tmppp);
        cp[iw][ix]=crmul(cq[ix+ix2][iw+nf1],tmpp);}
        else{
	cp[iw][ix]=cq[ix+ix2][iw+nf1];}
        }
	cp1[iw][ix]=cmplx(0.0,0.0);
	}

        ix=sx/dx-ifx;
	ixshot=ix;
        warn("ix %d",ix);

        for(iw=0;iw<nw;iw++){
        cp1[iw][ix-ix2]=wlsp[iw+nf1];
        }

                        
        free2float(p);
        free2complex(cq);
        free1float(wl);
        free1complex(wlsp);