normray.c

su 的源代码库

				/* have been written */
				if (nxzw>=nxz) break;

				/* ray parameters at this step */
				x = rsi->x;
				z = rsi->z;
				px = rsi->px;
				pz = rsi->pz;
				fa = rsi->f->fa;
				dsdx = fa->dsdx;
				dsdz = fa->dsdz;

				/* sigma at this step and next step */
				sigma1 = rsi->sigma;
				sigma2 = (irs<nrs-1)?rsi->rsNext->sigma:sigma1;

				/* number of sigma between this */
				/* step and next */
				ns = 1+(sigma2-sigma1)*(nxz-1)/(2.0*sigmamax);

				/* increment in sigma */
				ds = (sigma2-sigma1)/ns;

				/* loop over sigma */
				for (is=0,dsigma=0.0; is<ns; ++is,dsigma+=ds) {

					/* compute and write x and z */
					xs = x+dsigma*(px+dsigma*0.25*dsdx);
					zs = z+dsigma*(pz+dsigma*0.25*dsdz);
					if (efwrite(&zs,sizeof(float),1,rayfp)
						!=1)
						err("error writing rays!\n");
					if (efwrite(&xs,sizeof(float),1,rayfp)
						!=1)
						err("error writing rays!\n");


					/* if sufficient number written, */
					/* break */
					if (++nxzw>=nxz) break;
				}
			}

			/* if necessary, repeat last (x,z) */
			while (nxzw<nxz) {
				if (efwrite(&zs,sizeof(float),1,rayfp)!=1)
					err("error writing rays!\n");
				if (efwrite(&xs,sizeof(float),1,rayfp)!=1)
					err("error writing rays!\n");
				++nxzw;
			}

		} else {

			/* manipulate rayend information */
			re[iray].kend = -1;

		}

	    }

	} /* end loop over rays */

	if (infofp!=NULL)
		fprintf(infofp,"\nWrote %d rays to ray file.\n",icount);

	/* write number of rays to parameter file */
}

void writeWaves (FILE *wavefp, int nt, int nray, RayStep *rs[])
/* for each ray, write nt z(t),x(t) pairs uniformly sampled in time t */
{
	int iray,it;
	float tmax,dt,t1,t2,x1,z1,px1,pz1,sigma1,sigma2,s00,dsdx,dsdz,
		ax,bx,az,bz,at,bt,ct,ti,sigmah,sigmam,tm,
		dsigma,sigma,t,x,z,xstart,zstart;
	RayStep *rsi,*rs1,*rs2;
	FaceAttributes *fa;

	/* determine maximum time for all rays */
	for (iray=0,tmax=0.0; iray<nray; ++iray)
		for (rsi=rs[iray]; rsi!=NULL; rsi=rsi->rsNext)
			tmax = MAX(tmax,rsi->t);

	/* determine time sampling interval */
	dt = tmax/((nt>1)?nt-1:1);
	fprintf(stderr,"wavefront time increment = %g\n",dt);

	/* loop over rays */
	for (iray=0; iray<nray; ++iray) {

		/* initialize */
		rs1 = rs[iray];
		rs2 = rs1->rsNext;
		t1 = rs1->t;
		t2 = rs2->t;
		x1 = rs1->x;
		z1 = rs1->z;
		px1 = rs1->px;
		pz1 = rs1->pz;
		sigma1 = rs1->sigma;
		sigma2 = rs2->sigma;
		fa = rs1->f->fa;
		s00 = fa->s00;
		dsdx = fa->dsdx;
		dsdz = fa->dsdz;
		ax = px1;
		bx = 0.25*dsdx;
		az = pz1;
		bz = 0.25*dsdz;
		at = s00+dsdx*x1+dsdz*z1;
		bt = 0.5*(dsdx*px1+dsdz*pz1);
		ct = 0.0833333*(dsdx*dsdx+dsdz*dsdz);
		sigma = sigma1;
		t = t1;

		/* remember start x,z */
		xstart = x1;
		zstart = z1;

		/* loop over times */
		for (it=0,ti=0.0; it<nt; ++it,ti+=dt) {

			/* if necessary, go to next ray step */
			if (t2<ti) {
				do {
					rs1 = rs2;
					rs2 = rs1->rsNext;
					if (rs2==NULL) break;
				} while (rs2->t<ti);
				if (rs2==NULL) break;
				t1 = rs1->t;
				t2 = rs2->t;
				x1 = rs1->x;
				z1 = rs1->z;
				px1 = rs1->px;
				pz1 = rs1->pz;
				sigma1 = rs1->sigma;
				sigma2 = rs2->sigma;
				fa = rs1->f->fa;
				s00 = fa->s00;
				dsdx = fa->dsdx;
				dsdz = fa->dsdz;
				ax = px1;
				bx = 0.25*dsdx;
				az = pz1;
				bz = 0.25*dsdz;
				at = s00+dsdx*x1+dsdz*z1;
				bt = 0.5*(dsdx*px1+dsdz*pz1);
				ct = 0.0833333*(dsdx*dsdx+dsdz*dsdz);
				sigma = sigma1;
				t = t1;
			}

			/* determine dsigma for this time by bisection */
			sigmah = sigma2;
			while (t+0.01*dt<ti) {
				sigmam = 0.5*(sigma+sigmah);
				dsigma = sigmam-sigma1;
				tm = t1+dsigma*(at+dsigma*(bt+dsigma*ct));
				if (tm<ti) {
					t = tm;
					sigma = sigmam;
				} else {
					sigmah = sigmam;
				}
			}
			dsigma = sigma-sigma1;

			/* compute x and z */
			x = x1+dsigma*(ax+dsigma*bx);
			z = z1+dsigma*(az+dsigma*bz);

			/* write x and z */
			if (efwrite(&z,sizeof(float),1,wavefp)!=1)
				err("error writing waves!\n");
			if (efwrite(&x,sizeof(float),1,wavefp)!=1)
				err("error writing waves!\n");
		}

		/* finish writing x and z */
		while (it<nt) {
			if (efwrite(&zstart,sizeof(float),1,wavefp)!=1)
				err("error writing waves!\n");
			if (efwrite(&xstart,sizeof(float),1,wavefp)!=1)
				err("error writing waves!\n");
			++it;
		}
	}
}

static RayStep* traceRayInTri (RayStep *rs)
/* Trace ray in triangle.  Return pointer to 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 */
	rs->rsNext = (RayStep*)ealloc1(1,sizeof(RayStep));
	rs = rs->rsNext;
	rs->sigma = sigma;
	rs->x = x;
	rs->z = z;
	rs->px = px;
	rs->pz = pz;
	rs->t = t;
	rs->q1 = q1;
	rs->p1 = p1;
	rs->q2 = q2;
	rs->p2 = p2;
	rs->kmah = kmah;
	rs->nref = nref;
	rs->eu = eu;
	rs->f = f;
	rs->rsNext = NULL;
	rs->atten = atten;
	rs->ampli = ampli;
	rs->ampliphase = ampliphase;

	return rs;
}

static RayStep* traceRayAcrossEdge (RayStep *rs, FILE *infofp)
/* Trace ray across edge. */
/* If successful, return pointer to new RayStep. */
/* If unsuccessful, return NULL. */
/* 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,temp1,temp2,
		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,taper;
	float ampli,atten,dens1,dens2,ampliphase,coeff;
	EdgeUse *eu,*eum;
	EdgeUseAttributes *eua,*euma;
	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;
	atten = rs->atten;
	ampli = rs->ampli;
	ampliphase = rs->ampliphase;

	/* check for boundary */
	if (eu->euMate->f==NULL) {
		if (infofp!=NULL)
			fprintf(infofp,"Transmitted outside model.  "
					"Ray stopped.\n");
		return NULL;
	}

	/* determine sloth on this side of edge */
	fa = f->fa;
	s00 = fa->s00;
	ds1dx = fa->dsdx;
	ds1dz = fa->dsdz;
	s1 = s00+ds1dx*x+ds1dz*z;

	/* determine density on this side */
	dens1 = fa->dens;

	/* determine sloth on other side of edge */
	eum = eu->euMate;
	f = eum->f;
	fa = f->fa;
	s00 = fa->s00;
	ds2dx = fa->dsdx;
	ds2dz = fa->dsdz;
	s2 = s00+ds2dx*x+ds2dz*z;

	/* determine density on other side of the edge */
	dens2 = fa->dens;

	/* edge vector */
	dx = eum->vu->v->y-eu->vu->v->y;
	dz = eum->vu->v->x-eu->vu->v->x;

	/* fractional distance along edge */
	frac = (ABS(dx)>ABS(dz)) ? (x-eu->vu->v->y)/dx : (z-eu->vu->v->x)/dz;

	/* here is speeding up potential */

	/* linearly interpolate unit vector g tangent to edge */
	eua = eu->eua;
	euma = eum->eua;
	if (eua!=NULL && euma!=NULL) {
		gx = frac*euma->tx-(1.0-frac)*eua->tx;
		gz = frac*euma->tz-(1.0-frac)*eua->tz;
	} else {
		gx = -dx;
		gz = -dz;
	}
	scale = 1.0/sqrt(gx*gx+gz*gz);
	gx *= scale;
	gz *= scale;

	/* unit vector h normal to edge */
	hx = -gz;
	h_z = gx;

	/* remember ray parameters on this side */
	px1 = px;
	pz1 = pz;

	/* rotated ray parameters on this side */
	px1r = px*h_z-pz*hx;
	pz1r = px*hx+pz*h_z;

	/* rotated ray parameters on other side */
	px2r = px1r;
	pz2rs = s2-px2r*px2r;

	/* post-critical incidence*/
	if (pz2rs<=0.0) {
		if (infofp!=NULL)
			fprintf(infofp,"Post-critical transmission.  "
					"Ray stopped.\n");
		return NULL;
	}

	pz2r = sqrt(pz2rs);

	/* ray parameters on other side */
	px = px2r*h_z+pz2r*hx;
	pz = pz2r*h_z-px2r*hx;

	/* curvature term */
	c1ov1 = pz1r;
	c2ov2 = pz2r;
	oc2s = s2/pz2rs;
	if (eua!=NULL && euma!=NULL) {
		g = frac*euma->c-(1.0-frac)*eua->c;
		cterm = g*oc2s*(c1ov1-c2ov2);
	} else {
		cterm = 0.0;
	}

	/* update dynamic ray parameters */
	scale = (pz2r*sqrt(s1))/(pz1r*sqrt(s2));
	q1 = q1*scale;
	p1 = p1/scale+cterm*q1;
	q2 = q2*scale;
	p2 = p2/scale+cterm*q2;

	/* velocity derivatives tangent and normal to ray */
	fac1 = -0.5/(s1*s1);
	fac2 = -0.5/(s2*s2);
	dv1dl = fac1*(px1*ds1dx+pz1*ds1dz);
	dv2dl = fac2*(px*ds2dx+pz*ds2dz);
	dv1dm = fac1*(pz1*ds1dx-px1*ds1dz);
	dv2dm = fac2*(pz*ds2dx-px*ds2dz);

	/* inhomogeneity term */
	iterm = -px1r*oc2s
		*(2.0*(dv1dm*c1ov1-dv2dm*c2ov2)+px1r*(dv1dl-dv2dl));

	/* update dynamic ray parameters */
	p1 += iterm*q1;
	p2 += iterm*q2;

	/* transmission effects on amplitudes */
	if (ABS(s1-s2)>0.001 || dens1!=dens2) {
		if (infofp!=NULL)
			fprintf(infofp,"Transmitted WITH "
					"influence on amplitude.\n");
		if (dens1==FLT_MAX)
			dens1 = dens2 = 1.0;
		temp1 = dens2/dens1;
		temp2 = pz2r/pz1r;
		coeff = 1.0+(temp1-temp2)/(temp1+temp2);

		/* check if too close to critical */
		taper = (px1r*px1r)/s2;
		taper = (taper>TAPER) ? cos((taper-TAPER)*5.0*PI) : 1.0;

		/* complete amplitude coeff */
		ampli *= coeff*sqrt(pz2r/pz1r)*taper;

		if (infofp!=NULL) {
	   		fprintf(infofp,
				"  --incident side: v1=%g  dens1=%g  "
				"incidence angle=%g\n",
				1/sqrt(s1),dens1,acos(pz1r/sqrt(s1))*180/PI);
			fprintf(infofp,
				"  --opposite side: v2=%g  dens2=%g  "
				"refracted angle=%g\n",
				1/sqrt(s2),dens2,acos(pz2r/sqrt(s2))*180/PI);
			fprintf(infofp,
				"  --transmission coeff: %g\n",coeff);
			fprintf(infofp,
				"  --critical transm. taper: %g\n",taper);
			fprintf(infofp,
				"  --total ampli. coeff: %g\n",ampli);
		}

		/* grazing incidence*/
		if (pz1r*pz1r<=0.008*s1) {
			if (infofp!=NULL)
				fprintf(infofp,"Grazing incidence.  "
						"Ray stopped.\n");
			return NULL;
		}

	} else if (infofp!=NULL) {
		fprintf(infofp,"Transmitted WITHOUT influencing amplitude.\n");
	}

	/* return new raystep */
	rs->rsNext = (RayStep*)ealloc1(1,sizeof(RayStep));
	rs = rs->rsNext;
	rs->sigma = sigma;
	rs->x = x;
	rs->z = z;
	rs->px = px;
	rs->pz = pz;
	rs->t = t;
	rs->q1 = q1;