sudatumfd.c

su 的源代码库

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

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

/*********************** self documentation ******************************/
char *sdoc[]={
"									",
" SUDATUMFD - 2D zero-offset Finite Difference acoustic wave-equation	",
"		 DATUMing    						",
"									",
" sudatumfd <stdin > stdout [optional parameters]			",
"									",
" Required parameters:						   	",
"									",
" nt=	   number of time samples on each trace	       			",
" nx=	   number of receivers per shot gather				",
" nsx=	  number of shot gathers				    	",
" nz=	   number of downward continuation depth steps			",
" dz=	   depth sampling interval (in meters)				",
" mx=	   number of horizontal samples in the velocity model		",
" mz=	   number of vertical samples in the velocity model		",
" vfile1=       velocity file used for thin-lens term	    		",
" vfile2=       velocity file used for diffraction term			",
" dx=           horizontal sampling interval (in meters)                ",
"									",
" Optional parameters:						   	",
"									",
" dt=.004       time sampling interval (in seconds)			",
" buff=5	number of zero traces added to each side of each   	",
"	     shot gather as a pad			       		",
" tap_len=5     taper length (in number of traces)			",
" x_0=0.0       x coordinate of leftmost position in velocity model     ",
"									",
" Notes:								",
" The algorithm is a 45-degree implicit-finite-difference scheme based  ",
" on the one-way wave equation.  It works on poststack (zero-offset)    ",
" data only.  The two velocity files, vfile1 and vfile2, are binary     ",
" files containing floats with the format v[ix][iz].  There are two     ",
" potentially different velocity files for the thin-lens and            ",
" diffraction terms to allow for the use of a zero-velocity layer       ",
" which allows for datuming from an irregular surface.                  ",
"                                                                       ",
" Source and receiver locations must be set in the header values in     ",
" order for the datuming to work properly.  The leftmost position of    ",
" of the velocity models given in vfile1 and vfile2 must also be given. ",
" 									",
" 									",
NULL};

/* 
 * Author:  Chris Robinson, 10/16/00, CWP, Colorado School of Mines
 *
 *
 * References:
 *  Beasley, C., and Lynn, W., 1992, The zero-velocity layer: migration
 *    from irregular surfaces: Geophysics, 57, 1435-1443.
 *
 *  Claerbout, J. F., 1985, Imaging the earth's interior:  Blackwell
 *    Scientific Publications.
 *
 */

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

void tlens(complex ***ifr,complex ***ifr1,int nsx,int nx_pad,int nw,
	   float dw,float dz,float v,float dx,float **vel,
	   int zvel,int **gx,int buff,float x_0);

void diff(complex ***ifr1,complex ***ifr2,int nsx,int nx_pad,int nw,
	   float dw,float dz,float v,float dx,float **vel,
	   int zvel,int **gx,int buff,float x_0);

void padsub(complex ***ifr,complex ***ifr_pad,int nsx,int nx,
	    int nw,int nx_pad,int buff); 

void unpadsub(complex ***ifr,complex ***ifr_pad,int nsx,int nx,
	      int nw,int nx_pad,int buff);

void taper(complex ***ifr_pad,complex ***ifr_pt,int nsx,int nx_pad,
	   int nw,float dw,int tap_len,int buff,int nx);

segy tr;

int main(int argc, char **argv)
{
   FILE *vfp1,*vfp2;
   int nt,nz,nx,nx_pad,nsx;
   int mx,mz;
   int it,ix,isx,j;
   int **gx,**sx,**cdp,**offset;
   int buff,tap_len;

   float dt,dx,dz,x_0;
   float ***intime,***outtime;
   complex ***ifr,***ofr,***ofr1,***ifr3;
   complex ***ifr1,***ifr2;
   float v,**vel1,**vel2;

   int ntfft,nw,iw;
   float dw,*storage1,*p1,*p2,pi=PI;
   complex *cp1,*cp2;

   char *vfile1="";
   char *vfile2="";

   initargs(argc,argv);
   requestdoc(1);

   if (!getparint("nt",&nt)) err("must specify nt");
   if (!getparint("nx",&nx)) err("must specify nx");
   if (!getparint("nz",&nz)) err("must specify nz");
   if (!getparint("nsx",&nsx)) err("must specify nsx");
   if (!getparfloat("v",&v)) v=1000.0;
   if (!getparfloat("dz",&dz)) err("must specify dz");
   if (!getparfloat("x_0",&x_0)) err("must specify x_0");
   if (!getparint("buff",&buff)) err("must specify buff");
   if (!getparint("tap_len",&tap_len)) err("must specify tap_len");
   if (!getparint("mx",&mx)) err("must specify mx");
   if (!getparint("mz",&mz)) err("must specify mz");
   if (!getparstring("vfile1",&vfile1)) err("must specify vfile1");
   if (!getparstring("vfile2",&vfile2)) err("must specify vfile2");

   vfp1 = fopen(vfile1,"r");   
   vel1 = alloc2float(mz,mx);
   fread(vel1[0],sizeof(float),mz*mx,vfp1);
   fclose(vfp1);
   vfp2 = fopen(vfile2,"r");
   vel2 = alloc2float(mz,mx);
   fread(vel2[0],sizeof(float),mz*mx,vfp2);
   fclose(vfp2);

   nx_pad = nx+2*buff;
  
   intime = alloc3float(nt,nsx,nx);

   /*********************Initialize Array*******************/
   for(ix=0;ix<nx;ix++)
   for(it=0;it<nt;it++)
   for(isx=0;isx<nsx;isx++)
   {
      intime[ix][isx][it] = 0.0;
   }
   /*********************************************************/

   gx = alloc2int(nsx,nx);
   sx = alloc2int(nsx,nx); 
   cdp = alloc2int(nsx,nx);
   offset = alloc2int(nsx,nx);

   /*********************Read Traces into Array******************/
   ix=0; isx=0;
   while(gettr(&tr)&&ix<nx)
   {
      nt = (int) tr.ns;
      if(!getparfloat("dt",&dt))
      {
	 if(tr.dt)
	 {
	    dt = ((double) tr.dt)/1000000.0;
	 }
	 else
	 {
	    dt = 0.004;
	    warn("tr.dt not set, assuming dt=0.004");
	 }
      }/*end if*/
      if(!getparfloat("dx",&dx))
      {
	 if(tr.d2)
	 {
	    dx = tr.d2;
	 }
	 else
	 {
	    dx=1.0;
	    warn("tr.d2 not set, assuming d2=1.0");   
	 }
      }/*end if*/
      gx[ix][isx] = tr.gx;
      sx[ix][isx] = tr.sx;
      cdp[ix][isx] = tr.cdp;
      offset[ix][isx] = tr.offset;
      memcpy( (void *) intime[ix][isx], (void *) tr.data, nt*4);
      ix++;
      if(ix==nx) {ix=0; isx++;}
   }/*end while*/
   /**********************End Read Traces*****************************/


   /*******************Data Processing*******************************/
   /*****To set outtime=intime without looping use the following*****/
   /*memcpy((void *) &outtime[0][0][0],(void *) &intime[0][0][0],nx*nsx*nt*4);*/

   /**************Temporal Fourier Transform Data*******************/ 
   ntfft = npfar(2*nt);
   nw=ntfft/2+1;
   dw=2.0*pi/((ntfft-1)*dt);

   ifr = alloc3complex(nw,nsx,nx);
   ofr = alloc3complex(nw,nsx,nx_pad);

   for(isx=0;isx<nsx;isx++)
   {
      for(ix=0;ix<nx;ix++)
      {
	 storage1 = alloc1float(nt);
	 for(it=0;it<nt;it++) storage1[it] = intime[ix][isx][it];
	 p1 = alloc1float(ntfft);
	 cp1 = alloc1complex(nw);
	 memcpy(p1,storage1,nt*4);
	 for(it=nt;it<ntfft;it++) p1[it]=0.0;
	 pfarc(1,ntfft,p1,cp1);
	 for(iw=0;iw<nw;iw++)
	 {
	    ifr[ix][isx][iw] = cp1[iw];
	 }
	 free1float(storage1);
	 free1float(p1);
	 free1complex(cp1);
      }
   }
   free3float(intime);
   /*********End Temporal Fourier Transform**************************/

   ifr1 = alloc3complex(nw,nsx,nx_pad);
   ifr2 = alloc3complex(nw,nsx,nx_pad);

   padsub(ifr,ifr1,nsx,nx,nw,nx_pad,buff);
   free3complex(ifr);

   for(j=0;j<nz;j++)
   {
      /*warn("j=%d\n",j);*/
      taper(ifr1,ifr2,nsx,nx_pad,nw,dw,tap_len,buff,nx);
      tlens(ifr2,ifr1,nsx,nx_pad,nw,dw,dz,v,dx,vel1,
	    j,gx,buff,x_0);
      diff(ifr1,ifr2,nsx,nx_pad,nw,dw,dz,v,dx,vel2,
	   j,gx,buff,x_0);
      memcpy((void *) &ifr1[0][0][0],(void *) &ifr2[0][0][0],nx_pad*nsx*nw*8);
   }
   free3complex(ifr1);

   ifr3 = alloc3complex(nw,nsx,nx);
   unpadsub(ifr2,ifr3,nsx,nx,nw,nx_pad,buff);
   free3complex(ifr2);

   ofr1 = alloc3complex(nw,nsx,nx);
   memcpy((void *) &ofr1[0][0][0],(void *) &ifr3[0][0][0],nx*nsx*nw*8);
   free3complex(ifr3);

   /*********Temporal Inverse Fourier Transform Data*****************/
   outtime = alloc3float(nt,nsx,nx);
   for(isx=0;isx<nsx;isx++)
   {
      for(ix=0;ix<nx;ix++)
      {
	 p2 = alloc1float(ntfft);
	 cp2 = alloc1complex(nw);
	 for(iw=0;iw<nw;iw++) 
	 {
	    cp2[iw] = ofr1[ix][isx][iw];
	 }
	 pfacr(-1,ntfft,cp2,p2);
	 for(it=0;it<nt;it++) outtime[ix][isx][it] = p2[it]/(1.0*ntfft); 
	 free1float(p2);
	 free1complex(cp2);
      }
   }
   free3complex(ofr);
   /*********End Temporal Inverse Fourier Transform********************/



   /**********************End Data Processing*************************/


   /******************Output Data into Trace Format******************/
   for(isx=0;isx<nsx;isx++)