susynlvfti.c

su 的源代码库

				cr = 1.0;
			}

			/* if either cosine is negative, skip diffractor */
			if (ci<0.0 || cr<0.0) continue;

			/* two-way time and amplitude */
			time = ts+tg;
			if (er) {
				amp = sqrt(vg*vd/qg);
			} else {
				if (ls)
					amp = sqrt((vs*vd*vd*vg)/(qs*qg));
				else
					amp = sqrt((vs*vd*vd*vg)/
						(qs*qg*(qs+qg)));
			}
			amp *= (ci+cr)*ar*ds;
				
			/* add sinc wavelet to trace */
			addsinc(time,amp,nt,dt,ft,trace);
		}
	}
	
	/* allocate workspace */
	temp = ealloc1float(nt);
	
	/* apply half-derivative filter to trace */
	conv(nhd,-lhd,hd,nt,0,trace,nt,0,temp);

	/* convolve wavelet with trace */
	conv(w->lw,w->iw,w->wv,nt,0,temp,nt,0,trace);
	
	/* free workspace */
	free1float(temp);
}


static void raylvt (float v00, float dvdx, float dvdz,
	float x0, float z0, float x, float z,
	float a, float f, float l, int ntries,
	float nitmax, float epst, float epsx, float angxs,
	float *c, float *s, float *t, float *q)
/*****************************************************************************
Trace ray between two points, for linear velocity v(x,z) = v00+dvdx*x+dvdz*z
in a transversely isotropic medium
Reference:
Alkhalifah, T., 1993, Efficient synthetic seismograms for transversely isotropic media with constant velocity gradient:CWP project review. 
******************************************************************************
Input:
v00	     velocity at (x=0,z=0)
dvdx	    derivative dv/dx of velocity with respect to x
dvdz	    derivative dv/dz of velocity with respect to z
x0		x coordinate at beginning of ray
z0		z coordinate at beginning of ray
x		x coordinate at end of ray
z		z coordinate at end of ray
ntries	  number of iterations in Snell's law search
epsx	    lateral offset tolerance
a,f,l	   ratios of velocities for the FAI transversely isotropic medium
angxs		angle of symmetry axis of anisotriopy

Output:
c		cosine of propagation angle at end of ray
s		sine of propagation angle at end of ray
t		time along raypath
q		integral of velocity along raypath
*****************************************************************************/
{

	double px0,v0,aa,alfa,beta=0.0,vcp=0.0,oa;
	double v,c1,c2,f1,fx,cr,sr;
	double r,or,vcpi=0.0,vcps=0.0,p11,so,co=0.0,ci=0.0,cz2,
			dso,co1,so1=0.0,ss=0.0,cs=0.0;
	double si2,co2,sc,sx=0.0,cx=0.0,so3,diff,so4,diff1=0.0,sz2;
	double w1,w2,w3,px2,pz2,betaa,betab,alpha,pz0,betat,xx,dvcp,bbb,soo;
	double umax,time=0.0,fctn1,fctn2,fctnu,temp,mx,sm=0.0,vg,gang=0.0;
	double rt,timenew=0.0,tnm,u,sum,del,si=0.0,vrat,zz,uv,cw=0.0,sw=0.0,
			uv0,px,pz;
	int j,count,it=0,nit,sig,zneg,err,tfturn,sign;

	/* if (x,z) same as (x0,z0) */
	if ((x==x0) && (z==z0)) {
		*c = 1.0;
		*s = 0.0;
		*t = 0.0;
		*q = 0.0;
		return;
	}

	/*defining variables variables*/
	c1 = 1 + l;
	c2 = a + l;
	w1	= a-l;
	w2     = 1-l;
	w3     = (f+l)*(f+l);

	/*if necessary, rotate coordinates until the symmetry axis is in z direction*/
	if (angxs != 0.) {
		sx    = sin(angxs);
		cx    = cos(angxs);
		co2  = (cx)*(cx);
		si2  = (sx)*(sx);
		alfa = c2*si2+c1*co2;
		beta = sqrt(pow(w1*si2-w2*co2,2)+4*si2*co2*w3);
		vrat  = sqrt(.5*(alfa+beta));
		temp  = cx*x0 - sx*z0;
		z0    = cx*z0 + sx*x0;
		x0    = temp;
		temp  = cx*x - sx*z;
		z     = cx*z + sx*x;
		x     = temp;
		temp  = cx*dvdx - sx*dvdz;
		dvdz  = (cx*dvdz + sx*dvdx)/vrat;
		dvdx  = temp/vrat;
		v00   = v00/vrat;
	}

	 /* if constant velocity */
	if ((dvdx==0.0) && (dvdz==0.0)) {
		x  -= x0;
		z  -= z0;
		r   = sqrt(x*x+z*z);
		or  = 1.0/r;
		*s  = x*or;
		so  = *s;
		dso = .1;
		so4 =0;
		
		/*secant search for takeoff ray angle*/
		for(j=0; j<ntries; ++j) {
			if(so>1)
				so = .5*(1+so1);
			if(so<-1)
				so = .5*(-1+so1);
			so4 = so;
			co   = sqrt(1-so*so);
			co2  = (co)*(co);
			si2  = (so)*(so);
			sc   = (so)*(co);
			alfa = c2*si2+c1*co2;
			beta = sqrt(pow(w1*si2-w2*co2,2)+4*si2*co2*w3);
			vcp  = sqrt(.5*(alfa+beta));
			dvcp = (.5*(a+l)*sc-.5*(1+l)*sc+.25*(2*((a-l)*si2-(1-l)*co2)*(
				(a-l)*sc+(1-l)*sc)+4*sc*co2*(l+f)*(l+f)-4*
				si2*sc*(l+f)*(l+f))/beta)/vcp;
			gang = atan(dvcp/vcp);
			diff = z*tan(gang+asin(so))-x;
			if(j ==0) {
				if(ABS(diff) < epsx){
					err=0;
					break;
				}
				so1 = so;
				so  = so + dso;
				diff1  = diff;
			} else {
				count += 1;
				so3 = so;
				if(ABS(diff)<epsx ) {
					err  = 0;
					break;
				}
				so  = (diff1*so3-diff*so1)/(diff1-diff);
				diff1  = diff;
				so1 = so3;
			}
			/*fprintf(stderr,"diff=%f so=%f count=%d/n",diff,so,count);*/
			}
		*s = sin(asin(*s)+gang);
		*c = sqrt(1-(*s)*(*s));
		vg   = vcp/cos(gang);
		*t = r/(v00*vg);
		*q = r*v00*vg;

		/* if coordinates were rotated, unrotate cosine and sine */
		if (angxs != 0.) {
			temp  = cx*(*s) + sx*(*c);
			*c    = cx*(*c) - sx*(*s);
			*s    = temp;
		}
		return;
	}

	/*define distances before axes rotation*/
	xx = x-x0;
	zz = z-z0;

	/*if not rotated*/
	cr  = 1.0;
	sr  = 0.0;

	/* if necessary, rotate coordinates so that v(x,z) = v0+a*(z-z0) */
	if (dvdx!=0.0 || dvdz<0) {
		aa  = sqrt(dvdx*dvdx+dvdz*dvdz);
		oa  = 1.0/aa;
		cr  = dvdz*oa;
		sr  = dvdx*oa;
		temp = cr*x0-sr*z0;
		z0 = cr*z0+sr*x0;
		x0 = temp;
		temp  = cr*x-sr*z;
		z  = cr*z+sr*x;
		x  = temp;
	} else {
		aa = dvdz;
	}

	/*velocities and distances after possible rotation of coordinates*/
	v0  = v00 + aa*z0;
	v   = v00 + aa*z;
	x  -= x0;
	z  -= z0;
	z0  = v0/aa;
	
	/* defining regions depending on wether we needed the equation reversed */
	if(((x<0.) && (aa>0.)) || ((x>0.) && (aa<0.))) 
		sig=-1;
	else
		sig=1;
	zneg=1;

	

	/* if raypath is parallel to velocity gradient */
	if (x*x<0.00001*z*z) {
		so = sr;
		dso= .1;
		so4=0;

		/*secant search for takeoff ray angle*/
		for(j=0; j<ntries; ++j) {
			if(ABS(so)>1){
				so = .5*(sig*sm+so4);
			}
			so4 = so;
			co   = sqrt(1-so*so);
			co2  = (co)*(co);
			si2  = (so)*(so);
			sc   = (so)*(co);
			alfa = c2*si2+c1*co2;
			beta = sqrt(pow(w1*si2-w2*co2,2)+4*si2*co2*w3);
			vcp  = sqrt(.5*(alfa+beta));
			dvcp = (.5*(a+l)*sc-.5*(1+l)*sc+.25*(2*((a-l)*si2-(1-l)*co2)*(
				(a-l)*sc+(1-l)*sc)+4*sc*co2*(l+f)*(l+f)-4*
				si2*sc*(l+f)*(l+f))/beta)/vcp;
			gang = atan(dvcp/vcp);
			diff = zz*tan(gang+asin(so))-xx;
			if(j ==0) {
				if(ABS(diff) < epsx){
					err=0;
					break;
				}
				so1 = so;
				so  = so + dso;
				diff1  = diff;
			} else {
				count += 1;
				so3 = so;
				if(ABS(diff)<epsx ) {
					err  = 0;
					break;
				}
				so  = (diff1*so3-diff*so1)/(diff1-diff);
				diff1  = diff;
				so1 = so3;
			}
			}
		/* if ray is going down */
		if (z>0.0) {
			*s = sin(asin(sr)-gang);
			*c = sqrt(1-(*s)*(*s));
			vg   = vcp/cos(gang);
			*t = log(v/v0)/(aa*vg);
			*q = 0.5*z*(v+v0)*vg;

		/* else, if ray is going up */
		} else {
			*s =  sin(asin(sr)-gang);
			*c = -sqrt(1-(*s)*(*s));
			vg   = vcp/cos(gang);
			*t = -log(v/v0)/(aa*vg);
			*q = -0.5*z*(v+v0)*vg;
		}

	/* else raypath is circular with finite radius */
	} else {


		/*calculating the center for a circular raypath for isotropic*/
		z0   = v0/aa;
		z   += z0;
		x0   = -(x*x+z*z-z0*z0)/(2.0*x);
		rt    = sqrt(x0*x0+z0*z0);
		f1    = 1/(aa*aa*z0*z0/(rt*rt*v0*v0));
		so    = sig*z0/rt;
		soo   = so;

		/*intial trial parameters*/
		count = 1;
		err  = 1;
		dso   = .1;
		so4   = 0.;
		sm    = 1;

		/*for search effeciency*/
		if(dvdz*cx<0) {
			zneg = -1;
		}
	
		/* looping over ntries to find the spatial root through Secant method*/
		for(j=0; j<ntries; ++j) {
			if(ABS(so)>sm){
				so = .5*(sig*sm+so4);
			}
			if(sig*so<0)
				so = .5*(so4+0.);
			so4  = so;
			co   = zneg*sqrt(1-so*so);
			sw   = cr*so+sr*co;
			cw   = cr*co-sr*so;
			sz2  = sw*sw;
			cz2  = cw*cw;
			alfa = c2*sz2+c1*cz2;
			bbb  = w1*sz2-w2*cz2;
			beta = sqrt(bbb*bbb+4*sz2*cz2*w3);
			vcpi = sqrt(.5*(alfa+beta));
			if((z0-x*so/co)==0) {
				ss=sig;
			} else {
				p11  = z*so/(z0*co-x*so);
				ss   = sig*sqrt(p11*p11/(1+p11*p11));
			}
			/*if beyond the ray turning point*/
			if((z0-x*so/co)<0) tfturn=-1;
			else tfturn=1;
			cs   = zneg*tfturn*sqrt(1-ss*ss);
			si   = cr*ss+sr*cs;
			ci   = cr*cs-sr*ss;
			sz2  = si*si;
			cz2  = ci*ci;
			alfa = c2*sz2+c1*cz2;
			bbb  = w1*sz2-w2*cz2;
			beta = sqrt(bbb*bbb+4*sz2*cz2*w3);
			vcps = sqrt(.5*(alfa+beta));
			diff = z/ss - z0*(vcpi/vcps)/so;
			if(j ==0) {
				if(ABS(diff) < epsx){
					err=0;
					break;
				}
				co1 = co;
				so1 = so;
				so  = so + dso*tfturn*sig;
				diff1  = diff;
			} else {
				count += 1;
				so3 = so;
				if(ABS(diff)<epsx ) {
					err  = 0;
					break;
				}
				if(j==(ntries-1)){
					zneg *= -1;
					j     = -1;
					so    = soo;
					if(err == 3){
						err = 1;
						break;
					}
					err = 3;
				}else{
					so  = (diff1*so3-diff*so1)/(diff1-diff);
					diff1  = diff;
					so1 = so3;
				}
			}
			} 
		/*ignoring unfound rays*/
		if(err){ 
			*t=0;
			*q=10000000;
			*s=0;
			*c=1;
		}else{
			fx  = -z*cs/ss;
			f1  = vcps*vcps*(fx*fx+z*z);

			/* signed radius of raypath; if ray going clockwise, r > 0 */
			if (x<0.0){
				r = sqrt(fx*fx+z*z);
			}else {
				r = -sqrt(fx*fx+z*z);
			}

			/*initial ray parameters for traveltime calculations*/
			uv0  = 1/(vcpi*v0);
			uv   = 1/(vcps*v);
			px0  = sw*uv0;
			pz0  = cw*uv0;
			px   = si*uv;
			pz   = ci*uv;

			/*calculating the time*/
			if(dvdx==0) umax = (pz0-pz)/dvdz;
			else umax = (px0-px)/dvdx;

	    		for(nit=0; nit<nitmax; nit++)   {
				if(nit==0)   {
					u     =  0.;
						px2   = px0*px0;
					pz2   = pz0*pz0;
					alpha = c2*px2 + c1*pz2;
					betaa = (w1*px2 - w2*pz2)*(w1*px2 - w2*pz2);
					betab = 4*w3*px2*pz2;
					betat = sqrt(betaa+betab);
					fctn1 = 1/sqrt(.5*(alpha+betat));
					/*fctn1 = v0;*/
					u     = umax;
					px2   = (px0-dvdx*u)*(px0-dvdx*u);
					pz2   = (pz0-dvdz*u)*(pz0-dvdz*u);
					alpha = c2*px2 + c1*pz2;
					betaa = (w1*px2 - w2*pz2)*(w1*px2 - w2*pz2);
					betab = 4*w3*px2*pz2;
					betat = sqrt(betaa+betab);
					fctn2 = 1/sqrt(.5*(alpha+betat));
					/*fctn2 = v;*/

		    		timenew = 0.5*umax*(fctn2+fctn1);
		    		it	= 1;
				}
				else   {
		    			tnm = 1./it;
		    			del = umax*tnm;
		    			u   = 0.5*del;

		    			sum = 0.;
		    			for(j=0; j<it; j++)   {
						px2   = (px0-dvdx*u)*(px0-dvdx*u);
						pz2   = (pz0-dvdz*u)*(pz0-dvdz*u);
						alpha = c2*px2 + c1*pz2;
						betaa = (w1*px2 - w2*pz2)*(w1*px2 - w2*pz2);
						betab = 4*w3*px2*pz2;
						betat = sqrt(betaa+betab);
						fctnu = 1/sqrt(.5*(alpha+betat));
						sum  += fctnu;
						u    += del;
		    			}
		    			timenew = 0.5*(timenew+umax*sum*tnm);
		    			it	= 2*it;
				}
				if(ABS(timenew-time)<epst)   {
		    			time = ABS(timenew);
		    			break;
				}
				time = timenew;
	    		}
			*t = ABS(time);

			/*calculating the geometrical spreading*/
			sign = umax/ABS(umax);
			*s   = si*sign;
			*c   = ci*sign;
			mx   = sqrt(f1)*((cw-ci)*cr+(sw-si)*sr);
			co2  = ci*ci;
			si2  = si*si;
			sc   = si*ci;
			dvcp = (.5*c2*sc-.5*c1*sc+.25*(2*(w1*si2-w2*co2)*(
				w1*sc+w2*sc)+4*sc*co2*w3-4*
				si2*sc*w3)/beta)/vcps; 
			*q   = ABS(aa*cos(atan(dvcp/vcps))*r*mx);
		}
	}

	/* if coordinates were rotated, unrotate cosine and sine */
	if (angxs != 0.) {
		temp  = cx*(*s) + sx*(*c);
		*c    = cx*(*c) - sx*(*s);
		*s    = temp;
	}
}