normray.c

该程序主要用于三角网

						break;
					}
				}
			}

			if (infofp!=NULL)
				fprintf(infofp,"^^^\n"
					"Interaction with interface %d at "
					"(x=%g,z=%g).\n",
					temp,rslast->x,rslast->z);

			/* if ray at stopping edge, stop */
			if (stop) {
				if (infofp!=NULL)
					fprintf(infofp,"Hit stopping edge.  "
							"Ray stopped.\n");
				break;
			}

			/* else if ray at reflecting edge, reflect */
			else if (refl)
				rs = reflectRayFromEdge(rs,infofp);

			/* else trace ray across edge */
			else
				rs = traceRayAcrossEdge(rs,infofp);
			
		} while (rs!=NULL);
		
		/* reinitialize RayChecks for new ray */
		for (i=0; i<kk; ++i) {
			rc[i].hitseq = rc[i].hitseq-rc[i].nhits;
			rc[i].nhits = 0;
		}

		/* 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].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 = temp;
		re[iangle].ampli = rslast->ampli;
		re[iangle].ampliphase = rslast->ampliphase;
		re[iangle].atten = rslast->atten;
		re[iangle].dangle = dangle;

		/* optional output of rayend information */
		if (infofp!=NULL) {
			fprintf(infofp,
				"---------------------------------------"
				"------------\n\n"
				"The following information is stored at "
				"the rayend:\n\n");
			fprintf(infofp,
				"Ray stops at (%g,%g) at interface %d\n",
				re[iangle].x,re[iangle].z,re[iangle].kend);
			fprintf(infofp,
				"Takeoff angle=%g   Emergence angle=%g\n",
				re[iangle].ab*180./PI,-asin(re[iangle].px
				*1./sqrt(re[iangle].se))*180./PI);
			fprintf(infofp,"Takeoff velocity=%g   "
				"Emergence velocity=%g\n",
				1./sqrt(sb),1./sqrt(re[iangle].se));
			fprintf(infofp,"Slowness components px=%g   pz=%g\n",
				re[iangle].px,re[iangle].pz);
			fprintf(infofp,"Traveltime=%g   Sigma=%g\n",
				re[iangle].t,re[iangle].sigma);
			fprintf(infofp,"kmah index=%d   "
				"Number of reflections=%d\n",
				re[iangle].kmah,re[iangle].nref);
			fprintf(infofp,"Ray propagator q1=%g   q2=%g\n",
				re[iangle].q1,re[iangle].q2);
			fprintf(infofp,"Ray propagator p1=%g   p2=%g\n",
				re[iangle].p1,re[iangle].p2);
			fprintf(infofp,"Attenuation factor=%g   "
				"Angle increment=%g\n",
				re[iangle].atten,dangle*180./PI);
			fprintf(infofp,"Refl/Transm Amplitude=%g   Phase=%g\n",
				re[iangle].ampli,re[iangle].ampliphase); 
			fprintf(infofp,
				"---------------------------------------"
				"------------\n\n");
		}
	}
}

void writeRays (FILE *rayfp, int nxz, int nray, RayStep *rs[],
	RayEnd re[], int krecord, int prim, FILE *infofp)
/* for each ray, write x,z pairs */
{
	int iray,nrs,irs,nxzw,ns,is,icount;
	float sigmamax=0.0,sigma1,sigma2,x,z,px,pz,dsdx,dsdz,
		ds,xs,zs,dsigma;
	RayStep *rsi;
	FaceAttributes *fa;

	/* initialize counter */
	icount = 0;

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

		if(caustic==0)
		   re[iray].kmah=1;

		if(nonsurface==1)
		   re[iray].z=0;

		if(re[iray].kmah!=0 && re[iray].z<0.000001){
		if (krecord==INT_MAX || (krecord==re[iray].kend
			&& (prim==INT_MAX || prim==re[iray].nref))) {

			/* count rays */
			++icount;
			++raycount;

			/* count ray steps while determining maximum sigma */
			for (nrs=0,rsi=rs[iray]; rsi!=NULL;
				++nrs,rsi=rsi->rsNext)
				sigmamax = rsi->sigma;

			/* initialize number of (x,z) written for this ray */
			nxzw = 0;

			/* loop over ray steps */
			for (irs=0,rsi=rs[iray]; irs<nrs;
				++irs,rsi=rsi->rsNext) {

				/* if sufficient number of (x,z) */
				/* 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 */