sumiggbzo.c

su 的源代码库

		gka[ik] = gk[ik+nk-lwrap];
	for (ik=lwrap+nk; ik<lwrap+nk+lwrap; ++ik)
		gka[ik] = gk[ik-lwrap-nk];
	
	/* phase shift to account for non-centered x-axis */
	xshift = 0.5*(nx-1)*dx;
	for (ik=0,k=fka; ik<nka; ++ik,k+=dk) {
		phase = k*xshift;
		c = cos(phase);
		s = sin(phase);
		temp = gka[ik].r*c-gka[ik].i*s;
		gka[ik].i = gka[ik].r*s+gka[ik].i*c;
		gka[ik].r = temp;
	}
	
	/* compute k values at which to interpolate g(k) */
	for (ip=0,p=fp; ip<np; ++ip,p+=dp) {
		kp[ip] = w*p;
		
		/* if outside Nyquist bounds, do not interpolate */
		if (kp[ip]<fk && kp[ip]<ek)
			kp[ip] = fk-1000.0*dk;
		else if (kp[ip]>fk && kp[ip]>ek)
			kp[ip] = ek+1000.0*dk;
	}
		
	/* interpolate g(k) to obtain h(p) */
	ints8c(nka,dk,fka,gka,czero,czero,np,kp,hp);
	
	/* phase shift to account for non-centered x-axis and non-zero fx */
	xshift = -fx-0.5*(nx-1)*dx;
	for (ip=0; ip<np; ++ip) {
		phase = kp[ip]*xshift;
		c = cos(phase);
		s = sin(phase);
		h[ip].r = hp[ip].r*c-hp[ip].i*s;
		h[ip].i = hp[ip].r*s+hp[ip].i*c;
	}

	/* free workspace */
	free1complex(gka);
	free1complex(hp);
	free1float(kp);
}

static BeamData* beamData (float wmin, int nt, float dt, float ft, float *f)
/*****************************************************************************
Compute filtered complex beam data as a function of real and imaginary time.
******************************************************************************
Input:
wmin		minimum frequency (in radians per unit time)
nt		number of time samples
dt		time sampling interval
ft		first time sample
f		array[nt] containing data to be filtered

Returned:	pointer to beam data
*****************************************************************************/
{
	int ntpad,ntfft,nw,iwnyq,ntrfft,ntr,nti,nwr,it,itr,iti,iw;
	float dw,fw,dtr,ftr,dti,fti,w,ti,scale,*fa;
	complex *ca,*cb,*cfi,**cf;
	BeamData *bd;

	/* pad to avoid wraparound in Hilbert transform */
	ntpad = 25;

	/* fft sampling */
	ntfft = npfaro(nt+ntpad,2*(nt+ntpad));
	nw = ntfft/2+1;
	dw = 2.0*PI/(ntfft*dt);
	fw = 0.0;
	iwnyq = nw-1;

	/* real time sampling (oversample for future linear interpolation) */
	ntrfft = nwr = npfao(NOVERSAMPLE*ntfft,NOVERSAMPLE*ntfft+ntfft);
	dtr = dt*ntfft/ntrfft;
	ftr = ft;
	ntr = 1+(nt+ntpad-1)*dt/dtr;

	/* imaginary time sampling (exponential decay filters) */
	nti = NFILTER;
	dti = EXPMIN/(wmin*(nti-1));
	fti = 0.0;

	/* allocate space for filtered data */
	cf = alloc2complex(ntr,nti);

	/* allocate workspace */
	fa = alloc1float(ntfft);
	ca = alloc1complex(nw);
	cb = alloc1complex(ntrfft);

	/* pad data with zeros */
	for (it=0; it<nt; ++it)
		fa[it] = f[it];
	for (it=nt; it<ntfft; ++it)
		fa[it] = 0.0;

	/* Fourier transform and scale to make complex analytic signal */
	pfarc(1,ntfft,fa,ca);
	for (iw=1; iw<iwnyq; ++iw) {
		ca[iw].r *= 2.0;
		ca[iw].i *= 2.0;
	}

	/* loop over imaginary time */
	for (iti=0,ti=fti; iti<nti; ++iti,ti+=dti) {

		/* apply exponential decay filter */
		for (iw=0,w=fw; iw<nw; ++iw,w+=dw) {
			scale = exp(w*ti);
			cb[iw].r = ca[iw].r*scale;
			cb[iw].i = ca[iw].i*scale;
		}

		/* pad with zeros */
		for (iw=nw; iw<nwr; ++iw)
			cb[iw].r = cb[iw].i = 0.0;

		/* inverse Fourier transform and scale */
		pfacc(-1,ntrfft,cb);
		cfi = cf[iti];
		scale = 1.0/ntfft;
		for (itr=0; itr<ntr; ++itr) {
			cfi[itr].r = scale*cb[itr].r;
			cfi[itr].i = scale*cb[itr].i;
		}
	}

	/* free workspace */
	free1float(fa);
	free1complex(ca);
	free1complex(cb);

	/* return beam data */
	bd = (BeamData*)alloc1(1,sizeof(BeamData));
	bd->ntr = ntr;
	bd->dtr = dtr;
	bd->ftr = ftr;
	bd->nti = nti;
	bd->dti = dti;
	bd->fti = fti;
	bd->cf = cf;
	return bd;
}

static void setCell (Cells *cells, int jx, int jz)
/*****************************************************************************
Set a cell by computing its Gaussian beam complex time and amplitude.  
If the amplitude is non-zero, set neighboring cells recursively.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell to set
jz		z index of the cell to set
******************************************************************************
Notes:
To reduce the amount of memory required for recursion, the actual 
computation of complex time and amplitude is performed by the cellTimeAmp() 
function, so that no local variables are required in this function, except
for the input arguments themselves.
*****************************************************************************/
{
	/* if cell is out of bounds, return */
	if (jx<0 || jx>=cells->mx || jz<0 || jz>=cells->mz) return;

	/* if cell is live, return */
	if (cells->cell[jx][jz].live==cells->live) return;

	/* make cell live */
	cells->cell[jx][jz].live = cells->live;

	/* compute complex time and amplitude.  If amplitude is
	 * big enough, recursively set neighboring cells. */
	if (cellTimeAmp(cells,jx,jz)) {
		setCell(cells,jx+1,jz);
		setCell(cells,jx-1,jz);
		setCell(cells,jx,jz+1);
		setCell(cells,jx,jz-1);
	}
}

static int cellTimeAmp (Cells *cells, int jx, int jz)
/*****************************************************************************
Compute complex and time and amplitude for a cell.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell to set
jz		z index of the cell to set

Returned:	1 if Gaussian amplitude is significant, 0 otherwise
*****************************************************************************/
{
	int lx=cells->lx,lz=cells->lz;
	float dx=cells->dx,fx=cells->fx,dz=cells->dz,fz=cells->fz,
		wmin=cells->wmin,lmin=cells->lmin;
	Ray *ray=cells->ray;
	Cell **cell=cells->cell;
	int irs,kmah;
	float tmax,xc,zc,t0,x0,z0,x,z,tr,ti,ar,ai,
		v,q1,p1,q2,p2,e,es,scale,phase,vvmr,vvmi,c,s,ov,
		px,pz,pzpz,pxpx,pxpz,dvdx,dvdz,dvds,dvdn,
		wxxr,wxzr,wzzr,wxxi,wxzi,wzzi;
	RayStep *rs;

	/* maximum time */
	tmax = ray->rs[ray->nrs-1].t;

	/* cell coordinates */
	xc = fx+jx*lx*dx;
	zc = fz+jz*lz*dz;

	/* ray step nearest to cell */
	irs = nearestRayStep(ray,xc,zc);
	rs = &(ray->rs[irs]);
		
	/* ray time and coordinates */
	t0 = rs->t;
	x0 = rs->x;
	z0 = rs->z;
	
	/* real and imaginary parts of v*v*m = v*v*p/q */
	v = rs->v;
	q1 = rs->q1;
	p1 = rs->p1;
	q2 = rs->q2;
	p2 = rs->p2;
	e = wmin*lmin*lmin;
	es = e*e;
	scale = v*v/(q2*q2+es*q1*q1);
	vvmr = (p2*q2+es*p1*q1)*scale;
	vvmi = -e*scale;
	
	/* components of slowness vector, px and pz */
	c = rs->c;
	s = rs->s;
	ov = 1.0/v;
	px = s*ov;
	pz = c*ov;
	pzpz = pz*pz;
	pxpx = px*px;
	pxpz = px*pz;
	
	/* velocity derivatives along tangent and normal */
	dvdx = rs->dvdx;
	dvdz = rs->dvdz;
	dvds = s*dvdx+c*dvdz;
	dvdn = c*dvdx-s*dvdz;
	
	/* real part of W matrix */
	wxxr = pzpz*vvmr-2.0*pxpz*dvdn-pxpx*dvds;
	wxzr = (pxpx-pzpz)*dvdn-pxpz*(vvmr+dvds);
	wzzr = pxpx*vvmr+2.0*pxpz*dvdn-pzpz*dvds;
	
	/* imaginary part of W matrix */
	wxxi = pzpz*vvmi;
	wxzi = -pxpz*vvmi;
	wzzi = pxpx*vvmi;
	
	/* vector from ray to cell */
	x = xc-x0;
	z = zc-z0;
	
	/* real and imaginary parts of complex time */
	tr = t0+(px+0.5*wxxr*x)*x+(pz+wxzr*x+0.5*wzzr*z)*z;
	ti = 0.5*wxxi*x*x+(wxzi*x+0.5*wzzi*z)*z;

	/* real and imaginary parts of complex amplitude */
	kmah = rs->kmah;
	scale = pow(es/(q2*q2+es*q1*q1),0.25);
	phase = 0.5*(atan2(q2,q1*e)+2.0*PI*((kmah+1)/2));
	ar = scale*cos(phase);
	ai = scale*sin(phase);

	/* set cell parameters */
	cell[jx][jz].tr = tr;
	cell[jx][jz].ti = ti;
	cell[jx][jz].ar = ar;
	cell[jx][jz].ai = ai;

	/* return 1 if Gaussian amplitude is significant, 0 otherwise */
	return (wmin*ti>EXPMIN && tr<=tmax)?1:0;
}

static void accCell (Cells *cells, int jx, int jz)
/*****************************************************************************
Accumulate the contribution of a Gaussian beam in a cell and its
neighbors, recursively.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell to fill
jz		z index of the cell to fill
******************************************************************************
Notes:
To reduce the amount of memory required for recursion, the actual
accumulation is performed by function cellBeam(), so that no local
variables are required in this function, except for the input
arguments themselves.
*****************************************************************************/
{
	/* if cell is out of bounds, return */
	if (jx<0 || jx>=cells->mx-1 || jz<0 || jz>=cells->mz-1) return;

	/* if cell is dead, return */
	if (cells->cell[jx][jz].dead==cells->dead) return;

	/* make cell dead */
	cells->cell[jx][jz].dead = cells->dead;

	/* if upper-left corner of cell is not live, return */
	if (cells->cell[jx][jz].live!=cells->live) return;

	/* if remaining three corners of cell are live */
	if (cells->cell[jx+1][jz].live==cells->live &&
		cells->cell[jx][jz+1].live==cells->live &&
		cells->cell[jx+1][jz+1].live==cells->live) {

		/* accumulate beam in cell */
		cellBeam(cells,jx,jz);
	}

	/* recursively accumulate in neighboring cells */
	accCell(cells,jx+1,jz);
	accCell(cells,jx-1,jz);
	accCell(cells,jx,jz+1);
	accCell(cells,jx,jz-1);
}

static void cellBeam (Cells *cells, int jx, int jz)
/*****************************************************************************
Accumulate Gaussian beam for one cell.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell in which to accumulate beam
jz		z index of the cell in which to accumulate beam
*****************************************************************************/
{
	int lx=cells->lx,lz=cells->lz,nx=cells->nx,nz=cells->nz;
	float **g=cells->g;
	Cell **cell=cells->cell;
	BeamData *bd=cells->bd;
	int ntr=bd->ntr,nti=bd->nti;
	float dtr=bd->dtr,ftr=bd->ftr,dti=bd->dti,fti=bd->fti;
	complex **cf=bd->cf;
	int kxlo,kxhi,kzlo,kzhi,kx,kz,itr,iti;
	float odtr,odti,t00r,t01r,t10r,t11r,t00i,t01i,t10i,t11i,
		a00r,a01r,a10r,a11r,a00i,a01i,a10i,a11i,
		tx0r,tx1r,tx0i,tx1i,ax0r,ax1r,ax0i,ax1i,
		txzr,txzi,axzr,axzi,
		dtx0r,dtx0i,dtx1r,dtx1i,dax0r,dax0i,dax1r,dax1i,
		dtxzr,dtxzi,daxzr,daxzi,xdelta,zdelta,
		trn,tin,trfrac,mtrfrac,tifrac,mtifrac,
		cf0r,cf0i,cf1r,cf1i,cfr,cfi;
	complex *cf0,*cf1;

	/* inverse of time sampling intervals */
	odtr = 1.0/dtr;
	odti = 1.0/dti;

	/* complex time and amplitude for each corner */
	t00r = cell[jx][jz].tr;
	t01r = cell[jx][jz+1].tr;
	t10r = cell[jx+1][jz].tr;
	t11r = cell[jx+1][jz+1].tr;
	t00i = cell[jx][jz].ti;
	t01i = cell[jx][jz+1].ti;
	t10i = cell[jx+1][jz].ti;
	t11i = cell[jx+1][jz+1].ti;
	a00r = cell[jx][jz].ar;
	a01r = cell[jx][jz+1].ar;
	a10r = cell[jx+1][jz].ar;
	a11r = cell[jx+1][jz+1].ar;
	a00i = cell[jx][jz].ai;
	a01i = cell[jx][jz+1].ai;
	a10i = cell[jx+1][jz].ai;
	a11i = cell[jx+1][jz+1].ai;

	/* x and z samples for cell */
	kxlo = jx*lx;
	kxhi = kxlo+lx;
	if (kxhi>nx) kxhi = nx;
	kzlo = jz*lz;
	kzhi = kzlo+lz;
	if (kzhi>nz) kzhi = nz;

	/* fractional increments for linear interpolation */
	xdelta = 1.0/lx;
	zdelta = 1.0/lz;

	/* increments for times and amplitudes at top and bottom of cell */
	dtx0r = (t10r-t00r)*xdelta;
	dtx1r = (t11r-t01r)*xdelta;
	dtx0i = (t10i-t00i)*xdelta;
	dtx1i = (t11i-t01i)*xdelta;
	dax0r = (a10r-a00r)*xdelta;
	dax1r = (a11r-a01r)*xdelta;
	dax0i = (a10i-a00i)*xdelta;
	dax1i = (a11i-a01i)*xdelta;

	/* times and amplitudes at top-left and bottom-left of cell */
	tx0r = t00r;
	tx1r = t01r;
	tx0i = t00i;
	tx1i = t01i;
	ax0r = a00r;
	ax1r = a01r;
	ax0i = a00i;
	ax1i = a01i;

	/* loop over x samples */
	for (kx=kxlo; kx<kxhi; ++kx) {

		/* increments for time and amplitude */
		dtxzr = (tx1r-tx0r)*zdelta;
		dtxzi = (tx1i-tx0i)*zdelta;
		daxzr = (ax1r-ax0r)*zdelta;
		daxzi = (ax1i-ax0i)*zdelta;

		/* time and amplitude at top of cell */
		txzr = tx0r;
		txzi = tx0i;
		axzr = ax0r;
		axzi = ax0i;

		/* loop over z samples */
		for (kz=kzlo; kz<kzhi; ++kz) {

			/* index of imaginary time */
			iti = tin = (txzi-fti)*odti;
			if (iti<0 || iti>=nti-1) continue;

			/* pointers to left and right imaginary time samples */
			cf0 = cf[iti];
			cf1 = cf[iti+1];

			/* imaginary time linear interpolation coefficients */
			tifrac = tin-iti;
			mtifrac = 1.0-tifrac;

			/* index of real time */
			itr = trn = (txzr-ftr)*odtr;
			if (itr<0 || itr>=ntr-1) continue;

			/* real time linear interpolation coefficients */
			trfrac = trn-itr;
			mtrfrac = 1.0-trfrac;

			/* real and imaginary parts of complex beam data */
			cf0r = mtrfrac*cf0[itr].r+trfrac*cf0[itr+1].r;
			cf1r = mtrfrac*cf1[itr].r+trfrac*cf1[itr+1].r;
			cfr = mtifrac*cf0r+tifrac*cf1r;
			cf0i = mtrfrac*cf0[itr].i+trfrac*cf0[itr+1].i;
			cf1i = mtrfrac*cf1[itr].i+trfrac*cf1[itr+1].i;
			cfi = mtifrac*cf0i+tifrac*cf1i;

			/* accumulate beam */
			g[kx][kz] += axzr*cfr-axzi*cfi;

			/* increment time and amplitude */
			txzr += dtxzr;
			txzi += dtxzi;
			axzr += daxzr;
			axzi += daxzi;
		}

		/* increment times and amplitudes at top and bottom of cell */
		tx0r += dtx0r;
		tx1r += dtx1r;
		tx0i += dtx0i;
		tx1i += dtx1i;
		ax0r += dax0r;
		ax1r += dax1r;
		ax0i += dax0i;
		ax1i += dax1i;
	}
}