sufdmod2_pml.c

su 的源代码库

		}
	}


        if (pml_thick != 0)
		pml_init (nx, nz, dx, dz, dt, dvv, od, verbose);


	/* loop  ver time steps */
	for (it=0,t=0.0; it<nt; ++it,t+=dt) {
	
		/* if verbose, print time step */
		if (verbose>1) warn("it=%d  t=%g",it,t);
	
		/* update source function */
		if (ns==1)
			ptsrc(xs[0],zs[0],nx,dx,fx,nz,dz,fz,dt,t,
			      fmax,fpeak,tdelay,s);
		else
			exsrc(ns,xs,zs,nx,dx,fx,nz,dz,fz,dt,t,fmax,s);
		
		/* do one time step */
		tstep2(nx,dx,nz,dz,dt,dvv,od,s,pm,p,pp,abs);
		
		/* write waves */
		if (it%mt==0) {

			/* set selected trace header fields for all traces */
			cubetr.sx = xs[0];
			cubetr.sdepth = zs[0];
			cubetr.trid = 30 ;
			cubetr.ns = nz ;
			cubetr.d1 = dz ;
			cubetr.d2 = dx ;

			/* account for delay in source starting time */
			cubetr.delrt = - 1000.0 * tdelay;

			tracl = 0 ;

			/* set selected trace header fields trace by trace */
			for (ix=0 ; ix < nx ; ++ix) {
				++tracl;
				++tracr;

				cubetr.offset = ix * dx - xs[0];
				cubetr.tracl = (int) tracl;
				cubetr.tracr = (int) tracr;

				for (iz=0 ; iz < nz ; ++iz) {	
					cubetr.data[iz] = pp[ix][iz];
				}
				/* output traces of data cube */
				fputtr(stdout, &cubetr);
			}
		}

		/* if requested, save horizontal line of seismograms */
		if (hs!=NULL) {
			for (ix=0; ix<nx; ++ix)
				hs[ix][it] = pp[ix][hs1];
		}

		/* if requested, save vertical line of seismograms */
		if (vs!=NULL) {
			for (iz=0; iz<nz; ++iz)
				vs[iz][it] = pp[vs2][iz];
		}

		/* if requested, save seismograms at source locations */
		if (ss!=NULL) {
			for (is=0; is<ns; ++is)
				ss[is][it] = pp[ixs[is]][izs[is]];
		}

		/* roll time slice pointers */
		ptemp = pm;
		pm = p;
		p = pp;
		pp = ptemp;
	}

	/* if requested, write horizontal line of seismograms */
	if (hs!=NULL) {

		horiztr.sx = xs[0];
		horiztr.sdepth = zs[0];
		horiztr.trid = 1;
		horiztr.ns = nt ;
		horiztr.dt = 1000000 * dt ;
		horiztr.d2 = dx ;

		/* account for delay in source starting time */
		horiztr.delrt = -1000.0 * tdelay ; 

		tracl = tracr = 0;
		for (ix=0 ; ix < nx ; ++ix){
			++tracl;
			++tracr;

			/* offset from first source location */
			horiztr.offset = ix * dx - xs[0];

			horiztr.tracl = (int) tracl;
			horiztr.tracr = (int) tracr;

			for (it = 0 ; it < nt ; ++it){
				horiztr.data[it] = hs[ix][it];
			}
			
			fputtr(hseisfp , &horiztr);
		}

			
		fclose(hseisfp);
	}

	/* if requested, write vertical line of seismograms */
	if (vs!=NULL) {

		verttr.trid = 1;
		verttr.ns = nt ;
		verttr.sx = xs[0];
		verttr.sdepth = zs[0];
		verttr.dt = 1000000 * dt ;
		verttr.d2 = dx ;
		/* account for delay source starting time */
		verttr.delrt = -1000.0 * tdelay ;

		tracl = tracr = 0;
		for (iz=0 ; iz < nz ; ++iz){
			++tracl;
			++tracr;

			/* vertical line implies offset in z */
			verttr.offset = iz * dz - zs[0];

			verttr.tracl = (int) tracl;
			verttr.tracr = (int) tracr;

			for (it = 0 ; it < nt ; ++it){
				verttr.data[it] = vs[iz][it];
			}
			
			fputtr(vseisfp , &verttr);
		}

		fclose(vseisfp);
	}

	/* if requested, write seismogram at source position */
	if (ss!=NULL) {

		srctr.trid = 1;
		srctr.ns = nt ;
		srctr.dt = 1000000 * dt ;
		srctr.d2 = dx ;
		srctr.delrt = -1000.0 * tdelay ;

		tracl = tracr = 0;
		for (is=0 ; is < ns ; ++is){
			++tracl;
			++tracr;

			srctr.sx = xs[is];
			srctr.sdepth = zs[is];
			srctr.tracl = (int) tracl;
			srctr.tracr = (int) tracr;

			for (it = 0 ; it < nt ; ++it){
				srctr.data[it] = ss[is][it];
			}
			
			fputtr(sseisfp , &srctr);
		}

		fclose(sseisfp);
	}

	
	/* free space before returning */
	free2float(s);
	free2float(dvv);
	free2float(pm);
	free2float(p);
	free2float(pp);
	
	if (od!=NULL) free2float(od);
	if (hs!=NULL) free2float(hs);
	if (vs!=NULL) free2float(vs);
	if (ss!=NULL) free2float(ss);
	
	return(CWP_Exit());
}


void exsrc (int ns, float *xs, float *zs,
	int nx, float dx, float fx,
	int nz, float dz, float fz,
	float dt, float t, float fmax, float **s)
/*****************************************************************************
exsrc - update source pressure function for an extended source
******************************************************************************
Input:
ns		number of x,z coordinates for extended source
xs		array[ns] of x coordinates of extended source
zs		array[ns] of z coordinates of extended source
nx		number of x samples
dx		x sampling interval
fx		first x sample
nz		number of z samples
dz		z sampling interval
fz		first z sample
dt		time step (ignored)
t		time at which to compute source function
fmax		maximum frequency

Output:
s		array[nx][nz] of source pressure at time t+dt
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 03/01/90
******************************************************************************/
{
	int ix,iz,ixv,izv,is;
	float sigma,tbias,ascale,tscale,ts,xn,zn,
		v,xv,zv,dxdv,dzdv,xvn,zvn,amp,dv,dist,distprev;
	static float *vs,(*xsd)[4],(*zsd)[4];
	static int made=0;
	
	/* if not already made, make spline coefficients */
	if (!made) {
		vs = alloc1float(ns);
		xsd = (float(*)[4])alloc1float(ns*4);
		zsd = (float(*)[4])alloc1float(ns*4);
		for (is=0; is<ns; ++is)
			vs[is] = is;
		cmonot(ns,vs,xs,xsd);
		cmonot(ns,vs,zs,zsd);
		made = 1;
	}
	
	/* zero source array */
	for (ix=0; ix<nx; ++ix)
		for (iz=0; iz<nz; ++iz)
			s[ix][iz] = 0.0 *dt ;
	
	/* compute time-dependent part of source function */
	sigma = 0.25/fmax;
	tbias = 3.0*sigma;
	ascale = -exp(0.5)/sigma;
	tscale = 0.5/(sigma*sigma);
	if (t>2.0*tbias) return;
	ts = ascale*(t-tbias)*exp(-tscale*(t-tbias)*(t-tbias));
	
	/* loop over extended source locations */
	for (v=vs[0],distprev=0.0,dv=1.0; dv!=0.0; distprev=dist,v+=dv) {
		
		/* determine x(v), z(v), dx/dv, and dz/dv along source */
		intcub(0,ns,vs,xsd,1,&v,&xv);
		intcub(0,ns,vs,zsd,1,&v,&zv);
		intcub(1,ns,vs,xsd,1,&v,&dxdv);
		intcub(1,ns,vs,zsd,1,&v,&dzdv);
		
		/* determine increment along extended source */
		if (dxdv==0.0)
			dv = dz/ABS(dzdv);
		else if (dzdv==0.0)
			dv = dx/ABS(dxdv);
		else
			dv = MIN(dz/ABS(dzdv),dx/ABS(dxdv));
		if (v+dv>vs[ns-1]) dv = vs[ns-1]-v;
		dist = dv*sqrt(dzdv*dzdv+dxdv*dxdv)/sqrt(dx*dx+dz*dz);
		
		/* determine source amplitude */
		amp = (dist+distprev)/2.0;
		
		/* let source contribute within limited distance */
		xvn = (xv-fx)/dx;
		zvn = (zv-fz)/dz;
		ixv = NINT(xvn); 
		izv = NINT(zvn);
		for (ix=MAX(0,ixv-3); ix<=MIN(nx-1,ixv+3); ++ix) {
			for (iz=MAX(0,izv-3); iz<=MIN(nz-1,izv+3); ++iz) {
				xn = ix-xvn;
				zn = iz-zvn;
				s[ix][iz] += ts*amp*exp(-xn*xn-zn*zn);
			}
		}
	}
}

/* prototype of subroutine used internally */
static float ricker (float t, float fpeak);

void ptsrc (float xs, float zs,
	int nx, float dx, float fx,
	int nz, float dz, float fz,
	float dt, float t, float fmax, float fpeak, float tdelay, float **s)
/*****************************************************************************
ptsrc - update source pressure function for a point source
******************************************************************************
Input:
xs		x coordinate of point source
zs		z coordinate of point source
nx		number of x samples
dx		x sampling interval
fx		first x sample
nz		number of z samples
dz		z sampling interval
fz		first z sample
dt		time step (ignored)
t		time at which to compute source function
fmax		maximum frequency (ignored)
fpeak		peak frequency

Output:
tdelay		time delay of beginning of source function
s		array[nx][nz] of source pressure at time t+dt
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 03/01/90
******************************************************************************/
{
	int ix,iz,ixs,izs;
	float ts,xn,zn,xsn,zsn;
	
	/* zero source array */
	for (ix=0; ix<nx; ++ix)
		for (iz=0; iz<nz; ++iz)
			s[ix][iz] = 0.0 * dt*fmax;
	
	/* compute time-dependent part of source function */
	/* fpeak = 0.5*fmax;  this is now getparred */

	tdelay = 1.0/fpeak;
	if (t>2.0*tdelay) return;
	ts = ricker(t-tdelay,fpeak);
	
	/* let source contribute within limited distance */
	xsn = (xs-fx)/dx;
	zsn = (zs-fz)/dz;
	ixs = NINT(xsn);
	izs = NINT(zsn);
	for (ix=MAX(0,ixs-3); ix<=MIN(nx-1,ixs+3); ++ix) {
		for (iz=MAX(0,izs-3); iz<=MIN(nz-1,izs+3); ++iz) {
			xn = ix-xsn;
			zn = iz-zsn;
			s[ix][iz] = ts*exp(-xn*xn-zn*zn);
		}
	}
}

static float ricker (float t, float fpeak)
/*****************************************************************************
ricker - Compute Ricker wavelet as a function of time
******************************************************************************
Input:
t		time at which to evaluate Ricker wavelet
fpeak		peak (dominant) frequency of wavelet
******************************************************************************
Notes:
The amplitude of the Ricker wavelet at a frequency of 2.5*fpeak is 
approximately 4 percent of that at the dominant frequency fpeak.
The Ricker wavelet effectively begins at time t = -1.0/fpeak.  Therefore,
for practical purposes, a causal wavelet may be obtained by a time delay
of 1.0/fpeak.
The Ricker wavelet has the shape of the second derivative of a Gaussian.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 04/29/90
******************************************************************************/
{
	float x,xx;
	
	x = PI*fpeak*t;
	xx = x*x;
	/* return (-6.0+24.0*xx-8.0*xx*xx)*exp(-xx); */
	/* return PI*fpeak*(4.0*xx*x-6.0*x)*exp(-xx); */
	return exp(-xx)*(1.0-2.0*xx);
}

/* 2D finite differencing subroutine */

/* functions declared and used internally */
static void star1 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp);
static void star2 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp);
static void star3 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp);
static void star4 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp);
static void absorb (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **pm, float **p, float **pp,
	int *abs);

void tstep2 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp, int *abs)
/*****************************************************************************
One time step of FD solution (2nd order in space) to acoustic wave equation
******************************************************************************
Input:
nx		number of x samples
dx		x sampling interval
nz		number of z samples
dz		z sampling interval
dt		time step
dvv		array[nx][nz] of density*velocity^2
od		array[nx][nz] of 1/density (NULL for constant density=1.0)
s		array[nx][nz] of source pressure at time t+dt
pm		array[nx][nz] of pressure at time t-dt
p		array[nx][nz] of pressure at time t

Output:
pp		array[nx][nz] of pressure at time t+dt
******************************************************************************
Notes:
This function is optimized for special cases of constant density=1 and/or
equal spatial sampling intervals dx=dz.  The slowest case is variable
density and dx!=dz.  The fastest case is density=1.0 (od==NULL) and dx==dz.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 03/13/90
******************************************************************************/
{
	/* convolve with finite-difference star (special cases for speed) */
	if (od!=NULL && dx!=dz) {
		star1(nx,dx,nz,dz,dt,dvv,od,s,pm,p,pp);
	} else if (od!=NULL && dx==dz) {
		star2(nx,dx,nz,dz,dt,dvv,od,s,pm,p,pp);
	} else if (od==NULL && dx!=dz) {
		star3(nx,dx,nz,dz,dt,dvv,od,s,pm,p,pp);
	} else {
		star4(nx,dx,nz,dz,dt,dvv,od,s,pm,p,pp);
	}
	
	/* absorb along boundaries */
        if (pml_thick == 0) {
	   absorb(nx,dx,nz,dz,dt,dvv,od,pm,p,pp,abs);
        } else {
           pml_absorb(nx,dx,nz,dz,dt,dvv,od,pm,p,pp,abs);
	}
}

/* convolve with finite-difference star for variable density and dx!=dz */
static void star1 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp)
{