triseis.c

su 的源代码库

			/* adjusted time */
			ta = t+delt+0.5*cpoqr*nn-ft;

			/* skip if outside time window */
			if (ta<0.0 || ta>(nt-1)*dt) continue;

			/* compute cos and sin increments */
			cosd = cos(dw*ta);
			sind = sin(dw*ta);

			/* compute exp increment */
                        if (natten==0)
				/* no attenuation */
				expd = exp(-0.5*dw*cpoqi*nn);
			else
	                        /* include attenuation term */
				expd = exp(-0.5*dw*(cpoqi*nn+atten));

			/* initialize cos, sin, and exp */
			cosw = 1.0;
			sinw = 0.0;
			expw = 1.0;

			/* if no attenuation or noncausal attenuation */
			if (natten==0 || natten==1) {

				/* loop over frequencies */
				for (iw=0; iw<nw; ++iw) {

					/* accumulate beam */
					csyn[ig][iw].r += expw
						*(campr*cosw-campi*sinw);
					csyn[ig][iw].i += expw
						*(campr*sinw+campi*cosw);

					/* update cos, sin, and exp */
					tempd = cosw*cosd-sinw*sind;
					sinw = cosw*sind+sinw*cosd;
					cosw = tempd;
					expw *= expd;
				
					/* if amplitude too small, break */
					if (expw<0.01) break;

				}

			/* else, causal attenuation */
			} else {

				/* loop over frequencies */
				for (iw=0,w=0.0; iw<nw; ++iw,w+=dw) {

					/* compute cos and sin */
					tempd = w*(ta-lnvec[iw]);
					cosw = cos(tempd);
					sinw = sin(tempd);					

					/* accumulate beam */
					csyn[ig][iw].r += expw
						*(campr*cosw-campi*sinw);
					csyn[ig][iw].i += expw
						*(campr*sinw+campi*cosw);

					/* update exp */
					expw *= expd;

					/* if amplitude too small, break */
					if (expw<0.01) break;
				}
			}
		}
	}

	/* make complex ricker wavelet */
	for (iw=0,w=0.0; iw<nw; ++iw,w+=dw) {
		cwave[iw] = cricker(w,wpeak,delay);
		cwave[iw] = crmul(cwave[iw],sqrt(w));
	}

	/* apply wavelet to synthetics and inverse Fourier transform */
	for (ig=0; ig<ng; ++ig) {
		for (iw=0; iw<nw; ++iw) {
			csyn[ig][iw] = cmul(csyn[ig][iw],cwave[iw]);
		}
		pfacr(-1,ntfft,csyn[ig],syn);
		fwrite(syn,sizeof(float),nt,stdout);
	}
	
	/* free workspace */
	free1float(syn);
	free1float(lnvec);
	free1complex(cwave);
	free2complex(csyn);
}

static complex cricker (float w, float wpeak, float delay)
/*****************************************************************************
Compute Fourier transform of Ricker wavelet - complex function of frequency
******************************************************************************
Input:
w		frequency at which to evaluate transform
wpeak		peak (dominant) frequency in radians
delay		time shift (used for an approximately causal wavelet)
******************************************************************************
Notes:
The amplitude of the Ricker wavelet at a frequency of 2.5*wpeak is 
approximately 4 percent of that at the dominant frequency wpeak.
The Ricker wavelet effectively begins at time t = -2*PI/wpeak.  Therefore,
for practical purposes, a causal wavelet may be obtained by a time delay
of 2*PI/wpeak.
The Ricker wavelet has the shape of the second derivative of a Gaussian.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 02/28/91
******************************************************************************/
{
	return crmul(cexp(cmplx(-pow(w/wpeak,2.0),delay*w)),
		4.0*w*w*sqrt(PI)/pow(wpeak,3.0));
}

static int insideModel (Model *m, float x, float z)
/* return 1 if (x,z) inside model; 0 otherwise */
{
	return (m->ymin<=x && x<=m->ymax && m->xmin<=z && z<=m->xmax);
}

static void shootRays (Model *m, float xs, float zs,
	int nangle, float dangle, float fangle,
	int kstop, int kreflect,
	RayEnd re[])
/*****************************************************************************
Shoot rays via dynamic ray tracing for a sloth model
******************************************************************************
Input:
m		trianglulated sloth model
xs		horizontal coordinate of source
zs		vertical coordinate (depth) of source
nangle		number of ray takeoff angles
dangle		increment in ray takeoff angle
fangle		first ray takeoff angle
kstop		index of edge at which to stop ray
kreflect	index of edge at which to reflect ray

Output:
re		array[nangle] of RayEnds
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 02/28/91
modified: Andreas Rueger, Colorado School of Mines, 09/13/93
******************************************************************************/
{
	int iangle,dec;
	float angle,sb,zmax,xmax,eps;
	Tri *tris;
	FaceAttributes *fa;
	EdgeAttributes *ea;
	RayStep *rs,*rslast,*rsnext,*rs1,*rs2;

	/* determine triangle in model containing source location */
	tris = insideTriInModel(m,NULL,zs,xs);

        /* if source is on edge, move it by eps */
	dec = checkIfSourceOnEdge(tris,zs,xs);
	if (dec!=0) {
		zmax = m->xmax;
		eps = m->eps;
		if (dec==2) eps *= 5;
		zs = (zs<zmax-eps) ? zs+eps : zs-eps;

		/* if we are still on an edge, move x value */
		if (checkIfSourceOnEdge(tris,zs,xs) != 0) {
			xmax=m->ymax;
			xs = (xs<xmax-eps) ? xs+eps : xs-eps;
		}

		/* find new triangle */
		tris = insideTriInModel(m,NULL,zs,xs);
        }

	/* determine sloth at source location (beginning of ray) */
	fa = tris->fa;
	sb = fa->s00+fa->dsdx*xs+fa->dsdz*zs;

	/* allocate space for ray steps */
	rs1 = (RayStep*)malloc(sizeof(RayStep));
	rs2 = (RayStep*)malloc(sizeof(RayStep));
	
	/* loop over takeoff angles */
	for (iangle=0,angle=fangle; iangle<nangle; ++iangle,angle+=dangle) {
	
		/* initialize ray step pointers */
		rs = rs1;
		rsnext = rs2;
	
		/* initialize linked list of ray steps */
		rs->sigma = 0.0;
		rs->x = xs;
		rs->z = zs;
		rs->px = sin(angle)*sqrt(sb);
		rs->pz = cos(angle)*sqrt(sb);
		rs->q1 = 1.0;
		rs->p1 = 0.0;
		rs->q2 = 0.0;
		rs->p2 = 1.0;
		rs->t = 0.0;
		rs->kmah = 0;
		rs->nref = 0;
		rs->eu = NULL;
		rs->f = tris;
		rs->ampli = 1;
                rs->ampliphase = 0;
                rs->atten = 0;

		/* trace ray */
		do {
			/* trace ray in triangle */
			traceRayInTri(rs,rsnext);
			
			/* remember last ray step */
			rslast = rs;
			rs = rsnext;
			rsnext = rslast;
			
			/* edge attributes */
			ea = rs->eu->e->ea;
			
			/* if ray at stopping edge, stop */
			if (ea!=NULL && ea->k==kstop) {
				rslast = rs;
				break;
				
			/* else if ray at reflecting edge, reflect */
			} else if (ea!=NULL && ea->k==kreflect) {
				reflectRayFromEdge(rs,rsnext);
			
			/* else attempt to trace ray across edge */
			} else if (!traceRayAcrossEdge(rs,rsnext)) {
					rsnext = NULL;
			}
			
			/* remember last ray step */
			rslast = rs;
			rs = rsnext;
			rsnext = rslast;
			
		} while(rs!=NULL);
		
		/* save ray end parameters */
		fa = rslast->f->fa;
		re[iangle].sigma = rslast->sigma;
		re[iangle].x = rslast->x;
		re[iangle].z = rslast->z;
		re[iangle].px = rslast->px;
		re[iangle].pz = rslast->pz;
		re[iangle].t = rslast->t;
		re[iangle].q1 = rslast->q1;
		re[iangle].p1 = rslast->p1;
		re[iangle].q2 = rslast->q2;
		re[iangle].p2 = rslast->p2;
		re[iangle].kmah = rslast->kmah;
		re[iangle].nref = rslast->nref;
                re[iangle].ampli = rslast->ampli;
		if (rslast->ampliphase<-3.0 && kstop==1 && rslast->nref==0)
			fprintf(stderr,"phase=%f    ampli=%f\n",
				rslast->ampliphase,rslast->ampli);
                re[iangle].ampliphase = rslast->ampliphase;
                re[iangle].atten = rslast->atten;
		re[iangle].sb = sb;
		re[iangle].se = fa->s00+fa->dsdx*rslast->x+fa->dsdz*rslast->z;
		re[iangle].dsdxe = fa->dsdx;
		re[iangle].dsdze = fa->dsdz;
		re[iangle].ab = angle;
		re[iangle].kend = (rs==NULL) ? -1 : kstop;
		re[iangle].dangle = dangle;
	}
}

static void traceRayInTri (RayStep *rs, RayStep *rsnew)
/* Trace ray in triangle.  Return new RayStep. */
{
	int kmah,nref;
	float x,z,px,pz,t,q1,p1,q2,p2,ddt,
		s00,dsdx,dsdz,xa,za,xb,zb,dx,dz,a,b,c,
		ds,q,sigma,dsigma,dsigma1,dsigma2,small,frac;
	float ampli,ampliphase,atten,qfac;
	EdgeUse *eu,*eut,*eusmall,*euzero=NULL;
	Face *f;
	FaceAttributes *fa;

	/* get input parameters */
	sigma =	rs->sigma;
	x = rs->x;
	z = rs->z;
	px = rs->px;
	pz = rs->pz;
	t = rs->t;
	q1 = rs->q1;
	p1 = rs->p1;
	q2 = rs->q2;
	p2 = rs->p2;
	kmah = rs->kmah;
	nref = rs->nref;
	eu = rs->eu;
	f = rs->f;
	ampli = rs->ampli;
	ampliphase = rs->ampliphase;
	atten = rs->atten;

	/* determine sloth function */
	fa = f->fa;
	s00 = fa->s00;
	dsdx = fa->dsdx;
	dsdz = fa->dsdz;

	/* determine Q-factor */
	qfac = fa->qfac;

	/* determine edge intersection with smallest positive dsigma */
	eusmall = NULL;
	small = FLT_MAX;
	eut = f->eu;
	do {
		/* edge endpoints */
		xa = eut->vu->v->y;  za = eut->vu->v->x;
		xb = eut->euCW->vu->v->y;  zb = eut->euCW->vu->v->x;

		/* coefficients b and c of equation to be solved for dsigma */
		dx = xb-xa;  dz = zb-za;
		b = dx*pz-dz*px;
		c = (eut==eu) ? 0.0 : dx*(z-za)-dz*(x-xa); 

		/* if sloth constant, solve linear equation */
		if (dsdx==0.0 && dsdz==0.0) {
			if (b!=0.0)
				dsigma1 = -c/b;
			else
				dsigma1 = FLT_MAX;
			dsigma2 = FLT_MAX;
		
		/* else, if sloth not constant, solve quadratic equation */
		} else {
			a = 0.25*(dx*dsdz-dz*dsdx);
			ds = b*b-4.0*a*c;
			if (ds<0.0) {
				dsigma1 = dsigma2 = FLT_MAX;
			} else {
				q = -0.5*((b<0.0)?b-sqrt(ds):b+sqrt(ds));
				if (a!=0.0)
					dsigma1 = q/a;
				else
					dsigma1 = FLT_MAX;
				if (q!=0.0)
					dsigma2 = c/q;
				else
					dsigma2 = FLT_MAX;
			}
		}
		
		/* remember edge with smallest positive dsigma */
		if (0.0<dsigma1 && dsigma1<small) {
			small = dsigma1;
			eusmall = eut;
		}
		if (0.0<dsigma2 && dsigma2<small) {
			small = dsigma2;
			eusmall = eut;
		}

		/* remember edge for which dsigma is zero */
		if (dsigma1==0.0) euzero = eut;
		if (dsigma2==0.0) euzero = eut;

		/* next edge use */
		eut = eut->euCW;

	} while (eut!=f->eu);

	/* ray exits at edge with smallest positive dsigma */
	if (eusmall!=NULL) {
		dsigma = small;
		eu = eusmall;

	/* but if no dsigma>0, choose the edge we are up against */
	} else {
		dsigma = 0.0;
		eu = euzero;
	}

	/* update dynamic ray parameters */
	evaluateDynamic(dsigma,px,pz,dsdx,dsdz,&q1,&p1,&q2,&p2,&kmah);	

	/* update ray parameters */
	sigma += dsigma;
	ddt = dsigma*(s00+dsdx*x+dsdz*z
		+dsigma*(0.5*(dsdx*px+dsdz*pz)
		+dsigma*(0.0833333*(dsdx*dsdx+dsdz*dsdz))));
	t += ddt;
	x += dsigma*(px+0.25*dsdx*dsigma);
	z += dsigma*(pz+0.25*dsdz*dsigma);
	px += 0.5*dsdx*dsigma;
	pz += 0.5*dsdz*dsigma;

	/* don't let ray exit too close to a vertex */
	xa = eu->vu->v->y;  dx = eu->euCW->vu->v->y-xa;
	za = eu->vu->v->x;  dz = eu->euCW->vu->v->x-za;
	frac = (ABS(dx)>ABS(dz))?(x-xa)/dx:(z-za)/dz;
	if (frac<0.0001) {
		x = xa+0.0001*dx;
		z = za+0.0001*dz;
	} else if (frac>0.9999) {
		x = xa+0.9999*dx;
		z = za+0.9999*dz;
	}

	/* compute contribution due to attenuation */
	if (qfac!=FLT_MAX)
		atten += ddt/qfac;
	else
		atten = 0;

	/* return new raystep */
	rsnew->sigma = sigma;
	rsnew->x = x;
	rsnew->z = z;
	rsnew->px = px;
	rsnew->pz = pz;
	rsnew->t = t;
	rsnew->q1 = q1;
	rsnew->p1 = p1;
	rsnew->q2 = q2;
	rsnew->p2 = p2;
	rsnew->kmah = kmah;
	rsnew->nref = nref;
	rsnew->eu = eu;
	rsnew->f = f;
	rsnew->atten = atten;
	rsnew->ampli = ampli;
	rsnew->ampliphase = ampliphase;
}

static int traceRayAcrossEdge (RayStep *rs, RayStep *rsnew)
/* Trace ray across edge. */
/* If successful, return 1, along with new RayStep. */
/* If unsuccessful, return 0. */
/* Failure to trace across an edge is because: */
/* (1) the ray was incident with angle greater than the critical angle, or */
/* (2) the ray is incident at a boundary edge */
{
	int kmah,nref;
	float sigma,x,z,px,pz,t,q1,p1,q2,p2,
		s00,s1,ds1dx,ds1dz,s2,ds2dx,ds2dz,
		px1,pz1,px1r,pz1r,px2r,pz2r,pz2rs,
		c1ov1,c2ov2,oc2s,g,cterm,iterm,
		fac1,fac2,dv1dl,dv2dl,dv1dm,dv2dm,
		scale,gx,gz,hx,h_z,frac,dx,dz;
        float ampli,atten,dens1,dens2,ampliphase,coeff;
	float coeff1, coeff2, taper;
	EdgeUse *eu,*eum;