gbbeam.c

该程序主要用于三角网

				"Ray %d is out of the aperature\n",ire);
			continue;
		}

		/* ray end parameters */
		sigma = 1/sqrt(re[ire].sigma);
		x = re[ire].x;
		z = re[ire].z;
		t = re[ire].t;
		pz = re[ire].pz;
		q1 = re[ire].q1;
		p1 = re[ire].p1;
		q2 = re[ire].q2;
		p2 = re[ire].p2;
		kmah = re[ire].kmah;
		atten = re[ire].atten;
		ampli = re[ire].ampli;
		ampliphase = re[ire].ampliphase;
		dangle = re[ire].dangle;
		dsdx = re[ire].dsdxe;
		dsdz = re[ire].dsdze;
		v0 = 1.0/sqrt(re[ire].sb);
		a0 = re[ire].ab;

		/* compute ray dependant constants */
		v2 = v*v;
		v4 = v2*v2;
		dvds = -0.5*v4*(dsdx*px+dsdz*pz);
		dvdn = -0.5*v4*(dsdx*pz-dsdz*px);

		/* useful ray dependant information */
                if (infofp!=NULL) {
			fprintf(infofp,
				"\n***Ray emerges at (x=%g,z=%g) "
				"on interface %d.\n",x,z,kend); 		
			fprintf(infofp,"   Emergence angle=%g    time=%g\n",
				asin(v*px)*180.0/PI,t); 		
			fprintf(infofp,"   v=%g    dsdx=%g    dsdz=%g\n",
				v,dsdx,dsdz); 		
			fprintf(infofp,"   dvds=%g    dvdn=%g\n",dvds,dvdn); 
                	fprintf(infofp,"   ampli=%g    ampliphase=%g\n",
				ampli,ampliphase);
		}

		if (dangle==0.0) dangle = 1.0;
		scale = sigma*dangle*sqrt(v/(v0*2*PI));

		/* angularly dependent amplitude */
		intlin(nang,ang,amp,amp[0],amp[nang-1],1,&a0,&ampa);

		/* complex beam epsilon */
		if (bw!=0.0) {
			eps1 = 0.0;
			eps2 = -0.5*wpeak*bw*bw;
		} else {
			temp = c * ABS(p2/q2) + eps;
			eps1 = -q2/q1 * (p2*p1/(q2*q1) + temp*temp);
			eps1 = eps1/(p1*p1/(q1*q1) + temp*temp);
			eps2 = -temp/(p1*p1+q1*q1*temp*temp);		
		}
		ceps = cmplx(eps1,eps2);

		/* complex p, q, and p/q */
		cp = crmul(ceps,p1);  cp.r += p2;
		cq = crmul(ceps,q1);  cq.r += q2;
		cpoq = cdiv(cp,cq);
		cpoqr = cpoq.r;
		cpoqi = cpoq.i;

		/* if cpoqi negative, skip ray */
		/* (cpoqi<0 leads to exponential growth!) */
		if (cpoqi<0.0) {
			if (infofp!=NULL)
			 fprintf(infofp,"cpoqi negative, skipping ray!");
			continue;
		}

		/* complex q, adjusted for reflections */
		cqr = (nref%2) ? cneg(cq) : cq;
		
		/* frequency-independent amplitude factor */
		phase = -PI*((kmah+1)/2);
		camp = csqrt(cdiv(ceps,cqr));

                if (infofp!=NULL) {
			temp = sqrt(2*((eps1*q1+q2)*(eps1*q1+q2)
				+ eps2*eps2*q1*q1)/(-wpeak*eps2));
			fprintf(infofp,"   Gaussian beam information:\n"); 
			fprintf(infofp,"   Re(M)=%g \t Im(M)=%g\n",
				cpoqr,cpoqi); 		
			fprintf(infofp,"   const beam phase=%g degrees\n",
				180.0/PI*atan(camp.i/(camp.r+eps)));
			fprintf(infofp,"   epsilon1=%g \t epsilon2=%g\n",
				eps1,eps2); 	
			fprintf(infofp,"   ray-end beam width=%g\n",temp); 		
		}

		/* complex reflection/transmission coefficient */
		if (reftrans) {
			phase += ampliphase;
			camp = crmul(camp,ampli);
		}

		/* phase correction for 2.5-D spreading */
		phase += PI/4.0;

                camp = crmul(camp,scale*ampa);
		camp = cmul(camp,cmplx(cos(phase),sin(phase)));
		campr = camp.r;
		campi = camp.i;

		/* precompute causal attenuation */
		if (natten==2) {
			for (iw=0; iw<nw; ++iw)
				attenvec[iw] = atten*lnvec[iw];
		}

		/* loop over receiver locations */
		for (ig=0; ig<ng; ++ig) {

			/* delta x and delta z */
			dx = xg[ig]-x;
			dz = zg[ig]-z;
			disq = dx*dx+dz*dz;

			/* delta t and delta s */
			n = v*(pz*dx-px*dz);
			nn = n*n;

			/* do not extrapolate too far */
			if (lscale!=FLT_MAX && disq>sqlscale) continue;

			/* count beams that influence this receiver */
			++count[ig];

			dels = v*(px*dx+pz*dz)*(1.0-dvdn/v*n);
			delt = dels/v-0.5*dvds*dels*dels/v2;

			/* adjusted time */
			ta = t+delt+0.5*cpoqr*nn-ft;

			/* uncomment for extrapolation information */
/*			if (infofp!=NULL)
				fprintf(infofp,
					"   n=%g   ds=%g   dt=%g   tt=%g\n",
					n,dels,delt,ta);
*/

			/* skip if outside time window */
			if (ta<0.0 || ta>(nt-1)*dt) continue;

			/* compute cos and sin increments */
			cosd = cos(dw*ta);
			sind = sin(dw*ta);
                      
			/* compute exp increment */
			if (natten==0)
				/* no attenuation */
				expd = exp(-0.5*dw*cpoqi*nn);
                        else
	                        /* include attenuation term */
				expd = exp(-0.5*dw*(cpoqi*nn+atten));

			/* initialize cos, sin, and exp */
			cosw = 1.0;
			sinw = 0.0;
			expw = 1.0;

			/* if no attenuation or noncausal attenuation */
			if (natten==0 || natten==1) {

				/* loop over frequencies */
				for (iw=0; iw<nw; ++iw) {

					/* accumulate beam */
					csyn[ig][iw].r += expw
						*(campr*cosw-campi*sinw);
					csyn[ig][iw].i += expw
						*(campr*sinw+campi*cosw);

					/* update cos, sin, and exp */
					tempd = cosw*cosd-sinw*sind;
					sinw = cosw*sind+sinw*cosd;
					cosw = tempd;
					expw *= expd;

					/* if amplitude too small, break */
					if (expw<0.01) break;
				}

			/* else, causal attenuation */
			} else {

				/* loop over frequencies */
				for (iw=0,w=0.0; iw<nw; ++iw,w+=dw) {

					/* compute cos and sin */
					tempd = w*(ta-attenvec[iw]);
					cosw = cos(tempd);					
					sinw = sin(tempd);					

					/* accumulate beam */
					csyn[ig][iw].r += expw
						*(campr*cosw-campi*sinw);
					csyn[ig][iw].i += expw
						*(campr*sinw+campi*cosw);

					/* update exp */
					expw *= expd;

					/* if amplitude too small, break */
					if (expw<0.01) break;
				}
			}
		}
	}

	/* make complex ricker wavelet */
	for (iw=0,w=0.0; iw<nw; ++iw,w+=dw) {
		cwave[iw] = cricker(w,wpeak,delay);
		cwave[iw] = crmul(cwave[iw],sqrt(w));
	}

	/* apply wavelet to synthetics and inverse Fourier transform */
	for (ig=0; ig<ng; ++ig) {
		for (iw=0; iw<nw; ++iw)
			csyn[ig][iw] = cmul(csyn[ig][iw],cwave[iw]);
		pfacr(-1,ntfft,csyn[ig],syn);
		if (efwrite(syn,sizeof(float),nt,stdout)!=nt)
			err("error writing output traces!\n");
	}

	/* receiver contibution check */
        if (infofp!=NULL) {
		fprintf(infofp,"\nBeam contribution to each receiver "
				"*************************\n");
		fprintf(infofp,"Receiver No\t x_coord. \t z_coord. \t "
				"No. of beams\n");
	        for (ig=0; ig<ng; ++ig) 
			fprintf(infofp,"%d \t\t %g \t\t %g \t\t %d\n",
				ig,xg[ig],zg[ig],count[ig]);
        }

	/* free workspace */
        free1float(syn);
        free1float(lnvec);
        free1float(attenvec);
        free1complex(cwave);
	free2complex(csyn);
}

static complex cricker (float w, float wpeak, float delay)
/*****************************************************************************
Compute Fourier transform of Ricker wavelet - complex function of frequency
******************************************************************************
Input:
w		frequency at which to evaluate transform
wpeak		peak (dominant) frequency in radians
delay		time shift (used for an approximately causal wavelet)
******************************************************************************
Notes:
The amplitude of the Ricker wavelet at a frequency of 2.5*wpeak is 
approximately 4 percent of that at the dominant frequency wpeak.
The Ricker wavelet effectively begins at time t = -2*PI/wpeak.  Therefore,
for practical purposes, a causal wavelet may be obtained by a time delay
of 2*PI/wpeak.
The Ricker wavelet has the shape of the second derivative of a Gaussian.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 02/28/91
******************************************************************************/
{
	return crmul(cexp(cmplx(-pow(w/wpeak,2.0),delay*w)),
		4.0*w*w*sqrt(PI)/pow(wpeak,3.0));
}