normray.c

该程序主要用于三角网

		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;
	rs->p1 = p1;
	rs->q2 = q2;
	rs->p2 = p2;
	rs->kmah = kmah;
	rs->nref = nref;
	rs->eu = eum;
	rs->f = f;
	rs->rsNext = NULL; 
	rs->ampli = ampli;
	rs->ampliphase = ampliphase;
	rs->atten = atten;

	return rs;
}

static RayStep* reflectRayFromEdge (RayStep *rs, FILE *infofp)
/* Reflect ray from edge and  return pointer to new RayStep. */
{
	int kmah,nref;
	float sigma,x,z,px,pz,t,q1,p1,q2,p2,
		s00,dsdx,dsdz,s,coeff,
		px1,pz1,pxr,pzr,
		c1ov1,c2ov2,oc2s,g,cterm,iterm,
		dv1dl,dv2dl,dv1dm,dv2dm,
		scale,gx,gz,hx,h_z,frac,dx,dz;
	float atten,ampli,ampliphase,dens1,dens2,s2,
		ds2dx,ds2dz,temp1,temp2;
	EdgeUse *eu,*eum;
	EdgeUseAttributes *eua,*euma;
	Face *f,*fn;
	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 on incident side of edge */
	fa = f->fa;
	s00 = fa->s00;
	dsdx = fa->dsdx;
	dsdz = fa->dsdz;
	s = s00+dsdx*x+dsdz*z;

	/* determine dens on incident side of edge */
	dens1 = fa->dens;
	
	/* edge vector */
	eum = eu->euMate;

	/* determine sloth on other side of edge */
	fn = eum->f;
	fa = fn->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;

	if (dens1==FLT_MAX) dens1 = 1.0;
	if (dens2==FLT_MAX) dens2 = 1.0;

	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;

	/* linearly interpolate unit vector g tangent to edge */
	eua = eu->eua;
	euma = eum->eua;
	gx = frac*euma->tx-(1.0-frac)*eua->tx;
	gz = frac*euma->tz-(1.0-frac)*eua->tz;
	scale = 1.0/sqrt(gx*gx+gz*gz);
	gx *= scale;
	gz *= scale;

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

	/* remember incident ray parameters */
	px1 = px;
	pz1 = pz;

	/* rotated incident ray parameters */
	pxr = px*h_z-pz*hx;
	pzr = px*hx+pz*h_z;

	/* rotated reflected ray parameters */
	pxr = pxr;
	pzr = -pzr;

	/* reflected ray parameters */
	px = pxr*h_z+pzr*hx;
	pz = pzr*h_z-pxr*hx;

	/* curvature term */
	c1ov1 = c2ov2 = -pzr;
	oc2s = s/(pzr*pzr);
	g = frac*euma->c-(1.0-frac)*eua->c;
	cterm = g*oc2s*(c1ov1+c2ov2);

	/* update dynamic ray parameters */
	scale = -c2ov2/c1ov1;
	q1 = q1*scale;
	p1 = p1/scale+cterm*q1;
	q2 = q2*scale;
	p2 = p2/scale+cterm*q2;

	/* if sloth not constant */
	if (dsdx!=0.0 || dsdz!=0.0) {

		/* velocity derivatives tangent and normal to ray */
		scale = -0.5/(s*s);
		dv1dl = scale*(px1*dsdx+pz1*dsdz);
		dv2dl = scale*(px*dsdx+pz*dsdz);
		dv1dm = scale*(pz1*dsdx-px1*dsdz);
		dv2dm = scale*(pz*dsdx-px*dsdz);

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

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