normray.c

su 的源代码库

	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;
	}

	/* increment number of reflections */
	++nref;

	if (s2-pxr*pxr>=0) {
		/* pre-critical reflection */
        	if (infofp!=NULL)
			fprintf(infofp,
				"Reflected pre-critically at %g degrees.\n",
				acos(-pzr/sqrt(s))*180/PI);
		temp1 = -pzr*dens2/dens1/sqrt(s);
		temp2 = sqrt(s2/s-(pxr*pxr/s));
		coeff = (temp1-temp2)/(temp1+temp2);
		ampli *= coeff;  
	} else {
		/* post-critical reflection */
		coeff = 1;
		ampliphase -= 2*atan(-sqrt(pxr*pxr-s2)*dens1/(dens2*pzr));
		if (infofp!=NULL)
			fprintf(infofp,
				"Reflected post-critically at %g degrees.\n",
				acos(-pzr/sqrt(s))*180/PI);
	}

	if (infofp!=NULL) {
		fprintf(infofp,"  --incident side: v1=%g  dens1=%g\n",
			1/sqrt(s),dens1);
	   	fprintf(infofp,"  --opposite side: v2=%g  dens2=%g\n",
			1/sqrt(s2),dens2);
                fprintf(infofp,"  --reflection coeff: %g\n",coeff);
                fprintf(infofp,"  --total ampl coeff: %g\n",ampli);
                fprintf(infofp,"  --total phase shift: %g\n",ampliphase);
	}

	/* 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->ampli = ampli;
	rs->ampliphase = ampliphase;
	rs->atten = atten;

	return rs;
}

static void evaluateDynamic (float sigma, 
	float px, float pz, float dsdx, float dsdz,
	float *q1, float *p1, float *q2, float *p2, int *kmah)
/* evaluate dynamic ray parameters */
{
	double s0,s1,s2,ss,scale,a,b,c,d,e,
		q1i,p1i,q2i,p2i,q1o,p1o,q2o,p2o;

	/* get input dynamic ray parameters */
	q1i = *q1;  p1i = *p1;  q2i = *q2;  p2i = *p2;

	/* trivial case:  constant sloth or ray parallel to sloth gradient */
	if (pz*dsdx==px*dsdz) {
		q1o = q1i+p1i*sigma;
		q2o = q2i+p2i*sigma;
		p2o = p2i;
		p1o = p1i;

	/* else, general case */
	} else {

		/* constants */
		s0 = px*px+pz*pz;
		s1 = px*dsdx+pz*dsdz;
		s2 = 0.25*(dsdx*dsdx+dsdz*dsdz);
		ss = (s2*sigma+s1)*sigma+s0;
		scale = 1.0/sqrt(s0*ss);
		a = s0+0.5*s1*sigma;
		b = (0.25*s1*s1/s0-s2)*sigma;
		c = s0+s1*sigma;
		d = 0.5*s1+2.0*b;
		b *= sigma;
		e = (0.5*s1+s2*sigma)/ss;

		/* update q1 and q2 */
		q1o = scale*(a*(q1i+p1i*sigma)+b*q1i);
		q2o = scale*(a*(q2i+p2i*sigma)+b*q2i);

		/* update p1 and p2 */
		p1o = scale*(c*p1i+d*q1i)-e*q1o;
		p2o = scale*(c*p2i+d*q2i)-e*q2o;
	}

	/* update kmah index */
	if (q2i*q2o>=0.0 && p2i*p2o<0.0 && q2i*p2i<0.0)
		*kmah += 2;
	else if (q2o==0.0 || q2i*q2o<0.0
			|| (q2i==0.0 && p2i*p2o<0.0 && q2o*p2o>0.0))
		*kmah += 1;

	/* return updated dynamic ray parameters */
	*q1 = q1o;  *p1 = p1o;  *q2 = q2o;  *p2 = p2o;
}

static void evaluateDynamic2 (float sigma, 
	float px, float pz, float dsdx, float dsdz,
	float *q1, float *p1, float *q2, float *p2,
	float s0, float s1, float s2)
/* evaluate dynamic ray parameters at each point of the ray */
{
	double ss,scale,a,b,c,d,e,
		q1i,p1i,q2i,p2i,q1o,p1o,q2o,p2o;
	
	/* get input dynamic ray parameters */
	q1i = *q1;  p1i = *p1;  q2i = *q2;  p2i = *p2;
		
	/* trivial case:  constant sloth or ray parallel to sloth gradient */
	if (pz*dsdx==px*dsdz) {
		q1o = q1i+p1i*sigma;
		q2o = q2i+p2i*sigma;
		p2o = p2i;
		p1o = p1i;

	/* else, general case */
	} else {

		/* constants */
		ss = (s2*sigma+s1)*sigma+s0;
		scale = 1.0/sqrt(s0*ss);
		a = s0+0.5*s1*sigma;
		b = (0.25*s1*s1/s0-s2)*sigma;
		c = s0+s1*sigma;
		d = 0.5*s1+2.0*b;
		b *= sigma;
		e = (0.5*s1+s2*sigma)/ss;

		/* update q1 and q2 */
		q1o = scale*(a*(q1i+p1i*sigma)+b*q1i);
		q2o = scale*(a*(q2i+p2i*sigma)+b*q2i);

		/* update p1 and p2 */
		p1o = scale*(c*p1i+d*q1i)-e*q1o;
		p2o = scale*(c*p2i+d*q2i)-e*q2o;
	}
	
	/* return updated dynamic ray parameters */
	*q1 = q1o;  *p1 = p1o;  *q2 = q2o;  *p2 = p2o;
}		

int checkIfSourceOnEdge (Face *tris, float zs, float xs)
{
	float x1,x2,x3,z1,z2,z3,m12,m13,m23,b12,b13,b23,eps;
	EdgeUse *eu;

	eu = tris->eu;
	eps = tris->m->eps;

	/* get vertices */
	x1 = eu->vu->v->y;
	z1 = eu->vu->v->x;
	x2 = eu->euCW->vu->v->y;
	z2 = eu->euCW->vu->v->x;
	x3 = eu->euCCW->vu->v->y;
	z3 = eu->euCCW->vu->v->x;

	/* source is sitting on vertex */
	if ((xs-x1)*(xs-x1)+(zs-z1)*(zs-z1)<eps
		|| (xs-x2)*(xs-x2)+(zs-z2)*(zs-z2)<eps
		|| (xs-x3)*(xs-x3)+(zs-z3)*(zs-z3)<eps) return 2;

	/* check special cases and compute slope of edges */
	if (ABS(x1-x2)<0.1*eps) {
		if (ABS(xs-x1)<0.1*eps) {
			return 1;
		} else {
			m12 = 0.01*FLT_MAX;
		}
	} else {
		m12 = (z1-z2)/(x1-x2);
	}
	if (ABS(x1-x3)<0.1*eps) {
		if (ABS(xs-x1)<0.1*eps) {
			return 1;
		} else {
			m13 = 0.01*FLT_MAX;
		}
	} else {
		m13 = (z1-z3)/(x1-x3);
	}
	if (ABS(x2-x3)<0.1*eps) {
		if (ABS(xs-x2)<0.1*eps) {
			return 1;
		} else {
			m23 = 0.01*FLT_MAX;
		}
	} else {
		m23 = (z2-z3)/(x2-x3);
	}
	b12 = z1-m12*x1;
	b13 = z1-m13*x1;
	b23 = z2-m23*x2;

	/* source on edge? */
	if (ABS(zs-m12*xs-b12)<0.1*eps
		|| ABS(zs-m13*xs-b13)<0.1*eps
		|| ABS(zs-m23*xs-b23)<0.1*eps) return 1;

	return 0;
}

void writeFresnel (Model *m, FILE *rayfp, FILE *infofp, int nxz, int nray,
	RayStep *rs[], RayEnd re[], int krecord, float freq, int prim,
	FILE *fresnelfp)
/* for each ray, write x,z pairs and FresnelVolume */
{
	int iray,nrs,irs,nxzw,ns,is,icount;
	float sigmamax=0.0,sigma1,sigma2,x,z,px,pz,dsdx,dsdz,dummy,
		ds,xs,zs,dsigma,mfa,mfb,scale;
        float q1end=0.0,q2end=0.0,q1,q2,p1,p2,s0,s1,s2,q1i,q2i,p1i,p2i;
        float xfresnel1,xfresnel2,zfresnel1,zfresnel2;
        float *xa,*za,*r;
	RayStep *rsi;
	FaceAttributes *fa;

	/* allocate space for auxiliary arrays */
	xa = ealloc1float(nxz);
	za = ealloc1float(nxz);
	r = ealloc1float(nxz);

	/* initialize counter */
	icount = 0;

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

		if (krecord==INT_MAX || (krecord==re[iray].kend
			&& (prim==INT_MAX || prim==re[iray].nref))) { 

			/* count rays */
			++icount;

			/* count ray steps while determining maximum sigma */
			/* and the spreading values at the ray ends */
			for (nrs=0,rsi=rs[iray]; rsi!=NULL;
				++nrs,rsi=rsi->rsNext) {
				sigmamax = rsi->sigma;
				q1end = rsi->q1;
				q2end = rsi->q2;
			}

			/* 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;

			        /* constants for nontrivial case */
				s0 = px*px+pz*pz;
				s1 = px*dsdx+pz*dsdz;
				s2 = 0.25*(dsdx*dsdx+dsdz*dsdz);

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

				/* q1, q2, p1, p2 at this step */
				q1i = rsi->q1; 
				q2i = rsi->q2;
				p1i = rsi->p1; 
				p2i = rsi->p2;

				/* 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);
                        	        xa[nxzw] = xs;
                                	za[nxzw] = zs;
					if (efwrite(&zs,sizeof(float),1,rayfp)
						!=1) err("error writing "
						"ray file!\n");
					if (efwrite(&xs,sizeof(float),1,rayfp)
						!=1) err("error writing "
						"ray file!\n");

					q1 = q1i;
					q2 = q2i;
					p1 = p1i;
					p2 = p2i;

					/* compute and write p1,p2,q1,q2 */
                	                evaluateDynamic2(dsigma,px,pz,
						dsdx,dsdz,&q1,&p1,&q2,&p2,
						s0,s1,s2);

					/* compute m(Of,A) */
					if (q2==0)
						mfa = FLT_MAX;
					else
						mfa = p2/q2;

					/* compute m(Of,B) */
					mfb = -q1*q2end+q2*q1end;
					if (mfb==0)
						mfb = FLT_MAX;
					else
						mfb = (p2*q1end-p1*q2end)/mfb;

					/* compute fresnel radius */
					dummy = ABS(1.0/((mfa-mfb)*freq));
					r[nxzw] = sqrt(dummy);

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

			/* if necessary, repeat last (x,z) and r */
			while (nxzw<nxz) {
				xa[nxzw] = xs;
				za[nxzw] = zs;
				r[nxzw] = 0;
				if (efwrite(&zs,sizeof(float),1,rayfp)!=1)
					err("error writing ray file!\n");
				if (efwrite(&xs,sizeof(float),1,rayfp)!=1)
					err("error writing ray file!\n");
				++nxzw;
			}

			/* compute fresnel points */
			if (infofp!=NULL)
				fprintf(infofp,"Fresnel points:\n");
			for (is=0; is<nxz; ++is) {
				if (is==0 || is==nxz-1) {
					xfresnel1 = xfresnel2 = xa[is];
					zfresnel1 = zfresnel2 = za[is];
				} else {
					scale = sqrt((za[is-1]-za[is+1])
						*(za[is-1]-za[is+1])
						+(xa[is+1]-xa[is-1])
						*(xa[is+1]-xa[is-1]));
					if (scale==0 || is==0 || is==nxz-1) {
						xfresnel1 = xfresnel2 = xa[is];
						zfresnel1 = zfresnel2 = za[is];
					} else {
						dummy = 1.0/scale
							*(za[is-1]-za[is+1])
							*r[is];
						xfresnel1 = xa[is]+dummy;
						xfresnel2 = xa[is]-dummy;
						dummy = 1.0/scale
							*(xa[is+1]-xa[is-1])
							*r[is];
						zfresnel1 = za[is]+dummy;
						zfresnel2 = za[is]-dummy;
					}
				}
				if (infofp!=NULL)
					fprintf(infofp,
						"x1=%g  z1=%g  "
						"x2=%g  z2=%g  r=%g\n",
						xfresnel1,zfresnel1,
						xfresnel2,zfresnel2,r[is]);
				if (efwrite(&zfresnel1,sizeof(float),1,
					fresnelfp)!=1)
					err("error writing Fresnel file!\n");
				if (efwrite(&xfresnel1,sizeof(float),1,
					fresnelfp)!=1)
					err("error writing Fresnel file!\n");
				if (efwrite(&zfresnel2,sizeof(float),1,
					fresnelfp)!=1)
					err("error writing Fresnel file!\n");
				if (efwrite(&xfresnel2,sizeof(float),1,
					fresnelfp)!=1)
					err("error writing Fresnel file!\n");

			} /* end loop over fresnel points */

		} else {
			/* do not write rays and fresnel volumes; */
			/* manipulate rayend information */
			re[iray].kend = -1;
		}

        } /* close loop over rays */

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


	/* free space */
        free1float(xa);
        free1float(za);
        free1float(r);
}