triseis.c

该程序主要用于三角网

	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) return 0;

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

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

	/* if sloth function not same on both sides of edge */
	if (s1!=s2 || ds1dx!=ds2dx || ds1dz!=ds2dz || dens1!=dens2) {

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

		/* 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 */
		if (pz2rs<=0.0) return 0;

		/* grazing incidence */
		if (pz1r*pz1r <= 0.008*s1) return 0;

		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 (s1!=s2 || dens1!=dens2) {
			taper = (px1r*px1r)/s2;

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

			coeff1 = dens2/dens1; 
			coeff2 = pz2r/pz1r; 
			coeff = 1.0+(coeff1-coeff2)/(coeff1+coeff2);

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

	/* 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 = eum;
	rsnew->f = f;
        rsnew->ampli = ampli;
        rsnew->ampliphase = ampliphase;
        rsnew->atten = atten;

	return 1;
}

static void reflectRayFromEdge (RayStep *rs, RayStep *rsnew)
/* Reflect ray from edge and return new RayStep. */
{
	int kmah,nref;
	float sigma,x,z,px,pz,t,q1,p1,q2,p2,
		s00,dsdx,dsdz,s,
		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;
	atten = rs->atten;
	ampli = rs->ampli;
	ampliphase = rs->ampliphase;

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

        /* pre-critical */
	if (s2-pxr*pxr>=0) {
		temp1 = -pzr*dens2/dens1/sqrt(s);
		temp2 = sqrt(s2/s-(pxr*pxr/s));
		ampli *= (temp1-temp2)/(temp1+temp2);  
	} else {
		ampliphase -= 2*atan(-sqrt(pxr*pxr-s2)*dens1/(dens2*pzr));
	}

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

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 via modified midpoint method */
{
	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;
}

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