sufctanismod.c

su 的源代码库

*	 equation: Geophysics, vol. 39, p. 834-842.
*
* The waveform here differs from the one in Alford, Kelly, and Boore in
* that there is a time shift and an arbitrary amplitude scaling factor of 60.
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void tforce_akb(int nt, float *tforce, float dt, float fpeak)
{
  int	it;
  float	t1; 
  float 	t0=1.8/fpeak;

  for (it=0; it<nt; it++) {
    t1=it*dt;
    tforce[it] = -60.0*(t1-t0)*exp(-2.0*fpeak*fpeak
				   *(t1-t0)*(t1-t0));
  }
}

/************************************************************************
* tforce_spike -- Compute the	time response of a source function as
*	a spike.	
*************************************************************************
*************************************************************************
* Input: 
*	int nt		number of time step
*	float dt 	time step
*	float fpeak	(not used) peak frequency
*	
*************************************************************************
* Output: float *tforce		source array
*
*************************************************************************
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void tforce_spike(int nt, float *tforce, float dt, float fpeak)
{
  int	it;
  float	t1;

  for (it=0; it<nt; it++) {
    t1=it*dt;
    tforce[it] = 0.0;
  }

  tforce[1]=1.0;
}


/************************************************************************
* tforce_unit -- Compute the	time response of a source function as
*	a constant unit shift.	
*************************************************************************
*************************************************************************
* Input: 
*	int nt		number of time step
*	float dt 	time step
*	float fpeak	(not used) peak frequency 
*	
*************************************************************************
* Output: float *tforce		source array
*
*************************************************************************
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void tforce_unit(int nt, float *tforce, float dt, float fpeak)
{
  int	it;
  float	t1;

  for (it=0; it<nt; it++) {
    t1=it*dt;
    tforce[it] = 1.0;
  }

}


/************************************************************************
* locate_source --	specify the source location	
*			user can change to give many source locations
*************************************************************************
*************************************************************************
* Input: 
*	int nx		number of grids in x direction 
*	int nz		number of grids in z direction
*	int sx		single source x location
*	int sz		single source z location
*	int source 	1 for single source
*			2 many sources 
*	
*************************************************************************
* Output: float **xzsource	source positions array
*
*************************************************************************
* Note: 
*	User may modify this function to re-arrange the source locations
*	for their purpose
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	locate_source(int nx, int nz, int sx, int sz, 
		      float **xzsource, int source)
{
  int	ix, iz;

  for (ix=0; ix<nx; ix++)
    for (iz=0; iz<nz; iz++)
      { 
	xzsource[ix][iz]=0.0;
      }

  if( source==1 ) 
    xzsource[sx][sz]=1.0;

  if( source==2 ) {
    for(ix=0; ix<nx/14; ix++)
      xzsource[ix][nz/3]=1.0;
    for(ix=nx/14; ix<nx*17/70; ix++)
      xzsource[ix][(int)(nz/3+tan(PI/4.)*(ix-nx/14.))]=1.0;
    for(ix=17*nx/70; ix<80*nx/80; ix++)
      xzsource[ix][22*nz/30]=1.0;
    /*		for(ix=nx*63/80; ix<nx*15/16; ix++)
		xzsource[ix][(int)(nz*22/30-tan(PI/4.)
		*(ix-nx*63/80.))]=1.0;
		for(ix=15*nx/16; ix<nx; ix++)
		xzsource[ix][nz/3]=1.0;
    */	}

  if( source==3 ) {
    for(ix=0; ix<nx/14; ix++)
      xzsource[ix][nz/3]=1.0;
    for(ix=nx/14; ix<nx*34/70; ix++)
      xzsource[ix][(int)(nz/3+tan(PI/18.)*(ix-nx/14.))]=1.0;
    for(ix=17*nx/70; ix<80*nx/80; ix++)
      xzsource[ix][22*nz/30]=1.0;
  }
}



/************************************************************************
* moving_bc -- Computes the boundary containing the wavefield.
*		Finite difference computations are performed only inside
*		this boundary.
*************************************************************************
* Input: 
*	int it		time index
*	int nx		number of x values
*	int nz		number of z values
*	int sx		source x location
* 	int sz 		source z location
*	float dx	x increment
*	float dz	z increment
*	int impulse	1 if single source  
*	int movebc 	1 for moving boundary, 0 no moving boundary
*	float *t 	time array
*	float vmax	maximum velocity
*	
*************************************************************************
* Output: 	int *mbx1	right boundary
*		int *mbx2	left boundary
*		int *mbz1	top boundary
*		int *mbz2	bottom boundary
* These delimit the boundary within which the wavefield is computed.
*************************************************************************
* Notes: the moving boundary is used to improve the efficiency of the
*	modeling or migration, by constraining the area over which
*	the finite-differencing is computed.
* Caveat: this improves computational speed for forward modeling of
*		i.e. generating snapshots, this does not save much time
*		for the migration.
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	moving_bc (int it, int nx, int nz, int sx, int sz, 
		   float dx, float dz, int impulse, int movebc, 
		   float *t, float vmax, int *mbx1, int *mbz1, int *mbx2, int *mbz2)
{
  float	temp, dzdx=dz/dx;

  vmax *= 1.5;

  temp=vmax*t[it]/dz;

  if (movebc) {
    if (impulse) {
      if (0 > sz - temp - 0.04*nz) {
	/*endring foer: mbz1=0*/
	*mbz1=0;
      } else {
	*mbz1 = sz - temp - 0.04*nz;
      }

      if (nz > sz+temp+0.04*nz) {
	*mbz2=sz+temp+0.04*nz;
      } else {
	*mbz2=nz;
      }

      if (0 > sx - temp*dzdx - 0.04*nx) {
	*mbx1=0;
      } else {
	*mbx1 = sx - temp*dzdx - 0.04*nx;
      } 

      if (nx > sx + temp*dzdx + 0.04*nx) {
	*mbx2=sx + temp*dzdx + 0.04*nx; 
      } else {
	*mbx2 = nx;
      }

    } else {
      /*endring mbz1=0*/
      *mbz1=0;
      if (nz > temp+0.04*nz) {
	*mbz2=temp+0.04*nz;
      } else {
	*mbz2=nz;
      }

      *mbx1=0;
      *mbx2=nx;
    }

  } else {
    *mbx1=0;
    *mbx2=nx;
    *mbz1=0;
    *mbz2=nz;
  }

}


/************************************************************************
* moving_fctbc -- Computes the boundary containing the area where do
*	       the FCT correction.
*	       FCT correction are performed only inside
*	       this boundary.
*************************************************************************
* Input: 
* 	  int mbx1	right boundary containing the wavefield
*	  int mbx2	left boundary containing the wavefield
*	  int mbz1	top boundary containing the wavefield
*	  int mbz2	bottom boundary containing the wavefield
* 	  int nxcc1	right boundary for doing FCT
*	  int nxcc2	left boundary for doing FCT
*	  int nzcc1	top boundary for doing FCT
*	  int nzcc2	bottom boundary for doing FCT
*	
*************************************************************************
* Output: int *fctxbeg	right boundary
*	  int *fctxend	left boundary
*	  int *fctzbeg	top boundary
*	  int *fctzend	bottom boundary
* These delimit the boundary within which the FCT is applied.
*************************************************************************
* Notes: the moving boundary is used to improve the efficiency of the
*       modeling or migration, by constraining the area over which
*       the finite-differencing is computed.
* Caveat: this improves computational speed for forward modeling of
*	  i.e. generating snapshot 
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	moving_fctbc (int mbx1, int mbz1, 
		      int mbx2, int mbz2, 
		      int nxcc1, int nzcc1, 
		      int nxcc2, int nzcc2, 
		      int *fctxbeg, int *fctzbeg, 
		      int *fctxend, int *fctzend)
{
  /* FCT correction is localized to the area bounded by */
  /* (fctxend - fctxbeg) by (fctzend-fctzbeg) */
  /* contain the FCT correction boundary by the model boundary */
  if (mbx1 > nxcc1) {  /* left boundary */
    *fctxbeg=mbx1;
  } else {
    *fctxbeg=nxcc1;
  }
  if (mbx2 < nxcc2) {  /* right boundary */
    *fctxend=mbx2;
  } else {
    *fctxend=nxcc2;
  }
  if (mbz1 > nzcc1) {  /* top boundary */
    *fctzbeg=mbz1;
  } else {
    *fctzbeg=nzcc1;
  }
  if (mbz2 < nzcc2) {  /* bottom boundary */
    *fctzend=mbz2;
  } else {
    *fctzend=nzcc2;
  }
}




/************************************************************************
* strain2_x --  use 2nd-order FD to compute the derivative in x-direction
*		(the strain tensor e11 from input wavefield). 
*		e11 has the size (nx-1)*nz
*************************************************************************
*************************************************************************
* Input: 
* 	  int mbx1	right boundary containing the wavefield
*	  int mbx2	left boundary containing the wavefield
*	  int mbz1	top boundary containing the wavefield
*	  int mbz2	bottom boundary containing the wavefield
*	  float dx	spatial step in x-direction 	
*	  float **a	x-component wavefield;  size nx*nz
*	
*************************************************************************
* Output: float **da	derivative field (tensor e11)
*
*************************************************************************
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	strain2_x(int mbx1, int mbx2, int mbz1, int mbz2, 
		  float dx, float **a, float **da)
{
  int 	i, k;
  int	rdx = 1.0/dx;

  for (i=mbx1; i < mbx2-1; i++)
    for (k=mbz1; k < mbz2; k++)
      {
	da[i][k]=rdx*(a[i+1][k]-a[i][k]);
      }

}

/************************************************************************
* strain2_z --  use 2nd-order FD to compute the derivative in z-direction
*		(the strain tensor e33 from input wavefield). 
*		e33 has the size (nx-1)*nz
*************************************************************************
*************************************************************************
* Input: 
* 	  int mbx1	right boundary containing the wavefield
*	  int mbx2	left boundary containing the wavefield
*	  int mbz1	top boundary containing the wavefield
*	  int mbz2	bottom boundary containing the wavefield
*	  float dx	spatial step in x-direction 	
*	  float **a	z-component wavefield;  size (nx-1)*(nz-1)
*	
*************************************************************************
* Output: float **da	derivative field (tensor e33)
*
*************************************************************************
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	strain2_z(int mbx1, int mbx2, int mbz1, int mbz2, 
		  float dz, float **a, float **da)
{
  int 	i, k;
  int	rdz = 1.0/dz;

  for (i=mbx1; i < mbx2-1; i++)
    for (k=mbz1+1; k < mbz2-1; k++)
      {
	da[i][k]=rdz*(a[i][k]-a[i][k-1]);
      }

}

/************************************************************************
* strain2_xy --  use 2nd-order FD to compute the derivative in x-direction
*		(the strain tensor e12 from input wavefield). 
*		e12 has the size nx*(nz-1)
*************************************************************************
*************************************************************************
* Input: 
* 	  int mbx1	right boundary containing the wavefield
*	  int mbx2	left boundary containing the wavefield
*	  int mbz1	top boundary containing the wavefield
*	  int mbz2	bottom boundary containing the wavefield
*	  float dx	spatial step in x-direction 	
*	  float **ay	y-component wavefield;  size (nx-1)*(nz-1)
*	
*************************************************************************
* Output: float **da	derivative field (tensor e12)
*
*****************************************************************