/* 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, bhan@dix.mines.edu, 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);