sumigpsti.c

su 的源代码库

		fntab = (float)ntab;
	}
	
	/* parts of time/amplitude table */
	np = table->np;
	p = table->p;
	taumax = table->taumax;
	taumin = table->taumin;
	tturn = table->tturn;
	t = table->t;
	a = table->a;
	
	/* frequency sampling */
	ntfft = npfao(2*nt,4*nt);
	nw = ntfft;
	dw = 2.0*PI/(ntfft*dt);
	
	/* allocate workspace */
	pt = pw = alloc1complex(ntfft);
	qn = alloc1complex(ntau);
	qt = alloc1complex(ntau);
	
	/* nyquist frequency */
	wnyq = PI/dt;
	
	/* index of smallest frequency >= nyquist */
	iwnyq = (ntfft%2)?ntfft/2:ntfft/2+1;
	
	/* determine frequency filter sample indices */
	iwfil0 = MAX(0,MIN(iwnyq,NINT(2.0*PI*ffil[0]/dw)));
	iwfil1 = MAX(0,MIN(iwnyq,NINT(2.0*PI*ffil[1]/dw)));
	iwfil2 = MAX(0,MIN(iwnyq,NINT(2.0*PI*ffil[2]/dw)));
	iwfil3 = MAX(0,MIN(iwnyq,NINT(2.0*PI*ffil[3]/dw)));
	
	/* pad p(t) with zeros */
	for (it=0; it<nt; ++it)
		pt[it] = cp[it];
	for (it=nt; it<ntfft; ++it)
		pt[it].r = pt[it].i = 0.0;
	
	/* Fourier transform p(t) to p(w) (include scaling) */
	pfacc(1,ntfft,pt);
	for (iw=0; iw<nw; ++iw) {
		pw[iw].r *= 1.0/ntfft;
		pw[iw].i *= 1.0/ntfft;
	}
	
	/* initially zero normal and turned images */
	for (itau=0; itau<ntau; ++itau)
		qn[itau].r = qn[itau].i = qt[itau].r = qt[itau].i = 0.0;
	
	/* minimum and maximum frequency indices */
	wmin = ABS(k)/p[np-1];
	iwmin = MAX(MAX(1,iwfil0),(int)(wmin/dw));
	if (p[0]<=ABS(k)/wnyq)
		wmax = wnyq;
	else
		wmax = ABS(k)/p[0];
	iwmax = MIN(MIN(iwnyq-1,iwfil3),(int)(wmax/dw));
	
	/* loop over frequencies */
	for (iw=iwmin,w=iwmin*dw; iw<=iwmax; ++iw,w+=dw) {
			
		/* if slope not within range of table, continue */
		kow = ABS(k)/w;
		if (kow>p[np-1]) continue;
		
		/* amplitude of frequency filter */
		if (iwfil0<=iw && iw<iwfil1)
			ampf = (float)(iw-iwfil0)/(float)(iwfil1-iwfil0);
		else if (iwfil2<iw && iw<=iwfil3)
			ampf = (float)(iwfil3-iw)/(float)(iwfil3-iwfil2);
		else
			ampf = 1.0;
		
		/* weights for interpolation in table */
		xindex(np,p,kow,&ip1);
		ip1 = MIN(ip1,np-2);
		ip2 = ip1+1;
		p1 = p[ip1];
		p2 = p[ip2];
		a1 = (p2*p2-kow*kow)/(p2*p2-p1*p1);
		a2 = (kow*kow-p1*p1)/(p2*p2-p1*p1);
		wa1 = w*a1;
		wa2 = w*a2;
		
		/* maximum tau index for normal image */
		taumaxi = a1*taumax[ip1]+a2*taumax[ip2];
		itaumax = MIN(ntau-1,NINT(taumaxi/dtau));
		
		/* minimum tau index for turned image */
		taumini = a1*taumin[ip1]+a2*taumin[ip2];
		itaumin = MAX(0,NINT(taumini/dtau));
		
		/* phase at turning point */
		phst = 2.0*(wa1*tturn[ip1]+wa2*tturn[ip2]);
		
		/* filtered sum and differences for positive and negative w */
		pwsr = ampf*(pw[iw].r+pw[nw-iw].r);
		pwsi = ampf*(pw[iw].i+pw[nw-iw].i);
		pwdr = ampf*(pw[iw].r-pw[nw-iw].r);
		pwdi = ampf*(pw[iw].i-pw[nw-iw].i);
		
		/* accumulate normal image */
		if (ntflag&1) {
		for (itau=0; itau<itaumax; ++itau) {
			ampi = a1*a[ip1][itau]+a2*a[ip2][itau];
			phsi = wa1*t[ip1][itau]+wa2*t[ip2][itau];
			phsn = phsi*opi2;
			itab = fntab*(phsn-(int)phsn);
			qn[itau].r = qn[itau].r +
				ampi*(pwsr*ctab[itab]+pwdi*stab[itab]);
			qn[itau].i = qn[itau].i +
				ampi*(pwsi*ctab[itab]-pwdr*stab[itab]);
		}
		}
		
		/* accumulate turned image */
		if (ntflag&2) {
		for (itau=itaumin+1; itau<itaumax; ++itau) {
			ampi = a1*a[ip1][itau]+a2*a[ip2][itau];
			phsi = phst-(wa1*t[ip1][itau]+wa2*t[ip2][itau]);
			phsn = phsi*opi2;
			itab = fntab*(phsn-(int)phsn);
			qt[itau].r = qt[itau].r +
				ampi*(pwsr*stab[itab]-pwdi*ctab[itab]);
			qt[itau].i = qt[itau].i +
				ampi*(pwsi*stab[itab]+pwdr*ctab[itab]);
		}
		}
	}
	
	/* sum normal and turned images */
	for (itau=0; itau<ntau; ++itau) {
		cq[itau].r = qn[itau].r+qt[itau].r;
		cq[itau].i = qn[itau].i+qt[itau].i;
	}
	
	/* free workspace */
	free1complex(pt);
	free1complex(qn);
	free1complex(qt);
}	

TATable *tableta (int np, int ntau, float dtau, float a3333[],
	float a1111[], float a1133[], float a1313[],
	int nt, float dt, int verbose)
/*****************************************************************************
tabulate time shifts and amplitudes for use in phase-shift migration
******************************************************************************
Input:
np		number of slopes p
ntau		number of vertical two-way times tau
dtau		vertical time sampling interval
a1111		array[ntau] of a1111 elastic coeffecient as a function of tau
a3333		array[ntau] of a3333 elastic coeffecient as a function of tau
a1133		array[ntau] of a1133 elastic coeffecient as a function of tau
a1313		array[ntau] of a1313 elastic coeffecient as a function of tau
nt		number of time samples
dt		time sampling interval
verbose	 non-zero to print diagnostic info on stderr

Returned: 	pointer to the TATable
******************************************************************************/
{
	int jp,ktau,itau,jtau,ltau,it;
	float pj,taui,taul,tauj,ti,tl,ai,al,
		angle,cosa,frac,
		*p,*taumax,*taumin,*tturn,**t,**a;
	float f1,f2,f3,f4,f5,f6,px2,pz2,vp,eps,
		alpha,beta,gamma,det,rad,signbeta,q;
	float a1111t,a3333t,a1313t,a1133t;
	TATable *table;
	
	/* allocate table */
	table = alloc1(1,sizeof(TATable));
	table->np = np;
	table->p = p = alloc1float(np);
	table->taumax = taumax = alloc1float(np);
	table->taumin = taumin =alloc1float(np);
	table->tturn = tturn = alloc1float(np);
	table->t = t = alloc2float(ntau,np);
	table->a = a = alloc2float(ntau,np);
	
	for (it=0; it<nt; ++it){
		a1111[it] = .25*a1111[it];
		a3333[it] = .25*a3333[it];
		a1313[it] = .25*a1313[it];
		a1133[it] = .25*a1133[it];
	}
		
	/* loop over slopes p */
	for (jp=0; jp<np; ++jp) {
	
		/* slope p */
		p[jp] = pj = 2.0/sqrt(a1111[0])*sqrt((float)(jp)/(float)(np-1));
			
		/* time and amplitude at tau = 0 */
		t[jp][0] = 0.0;
		a[jp][0] = 1.0;
		
		/* tau index */
		ktau = 0;
		
		/* initial ray tracing parameters */
		taui = 0.0;
		ti = 0.0;
		ai = 1.0;

		a1111t = a1111[0];
		a3333t = a3333[0];
		a1313t = a1313[0];
		a1133t = a1133[0];

		f1 = a1111t+a1313t;
		f2 = a3333t+a1313t;
		f3 = a1111t-a1313t;
		f4 = a3333t-a1313t;
		f5 = 2.*(a1133t+a1313t)*(a1133t+a1313t);
		f6 = 2;

		eps   = .0001;
		px2   = pj*pj;
	    	alpha = f2*f2-f4*f4;
	    	beta  = 2*((f1*f2+f3*f4-f5)*px2-f2*f6);
	    	gamma = f6*f6-(2.*f1*f6-(f1*f1-f3*f3)*px2)*px2;
	    	det   = beta*beta-4.*alpha*gamma;

		if (det<0) continue;
		rad = sqrt(det);
		if(ABS(beta)>eps)   signbeta = ABS(beta)/beta;
	    	else                signbeta = 1.;
		q    = -.5*(beta+signbeta*rad);
		pz2  = gamma/q;
		if(pz2<0) continue;
		vp   = 1/sqrt(px2+pz2);

		angle = asin(MIN(1.0,pj*vp));
		
		/* loop over times t */
		for (it=1; it<nt; ++it) {
			
			/* remember last tau, t, and a */
			taul = taui;
			tl = ti;
			al = ai;
			
			/* update cosine of propagation angle */
			cosa = cos(angle)*sqrt(a3333t)/vp;
			
			/* update tau, t, and a */
			taui += dt*cosa;
			ti += dt*cosa*cosa;
			ai = 1.0;
			
			/* if ray emerges at surface, break */
			if (taui<0.0) break;
			
			/* update turning time and max,min tau */
			if (taui>=taul) {
				tturn[jp] = ti;
				taumax[jp] = taui;
				taumin[jp] = taui;
			} else {
				taumin[jp] = taui;
			}
			
			/* compute tau sample indices */
			itau = (int)(taui/dtau);
			ltau = (int)(taul/dtau);
			
			/* loop over tau samples crossed by ray */
			for (jtau=ltau+1; jtau<=MIN(itau,ntau-1); ++jtau) {
				
				/* tau of sample crossed */
				tauj = jtau*dtau;
				
				/* time and amp via linear interpolation */
				frac = (tauj-taul)/(taui-taul);
				ktau++;
				t[jp][ktau] = (1.0-frac)*tl+frac*ti;
				a[jp][ktau] = (1.0-frac)*al+frac*ai;
			}

			if (itau<ntau-1) {
				frac = (taui-(itau*dtau))/dtau;
				a1111t = (1.0-frac)*a1111[itau] +
					frac*a1111[itau+1];
				a3333t = (1.0-frac)*a3333[itau] +
					frac*a3333[itau+1];
				a1313t = (1.0-frac)*a1313[itau] +
					frac*a1313[itau+1];
				a1133t = (1.0-frac)*a1133[itau] +
					frac*a1133[itau+1];
			} else {
				a1111t = a1111[ntau-1] +
					(taui-(ntau-1)*dtau)*
					(a1111[itau]-a1111[itau-1])/dtau;
				a3333t = a3333[ntau-1] +
					(taui-(ntau-1)*dtau)*
					(a3333[itau]-a3333[itau-1])/dtau;
				a1313t = a1313[ntau-1] +
					(taui-(ntau-1)*dtau)*
					(a1313[itau]-a1313[itau-1])/dtau;
				a1133t = a1133[ntau-1] +
					(taui-(ntau-1)*dtau)*
					(a1133[itau]-a1133[itau-1])/dtau;
			}
			
		

			f1 = a1111t+a1313t;
			f2 = a3333t+a1313t;
			f3 = a1111t-a1313t;
			f4 = a3333t-a1313t;
			f5 = 2.*(a1133t+a1313t)*(a1133t+a1313t);
			f6 = 2;

			eps   = .0001;
	    		alpha = f2*f2-f4*f4;
	    		beta  = 2*((f1*f2+f3*f4-f5)*px2-f2*f6);
	    		gamma = f6*f6-(2.*f1*f6-(f1*f1-f3*f3)*px2)*px2;
	    		det   = beta*beta-4.*alpha*gamma;

			if (det<0) continue;
			rad = sqrt(det);
			if(ABS(beta)>eps)   signbeta = ABS(beta)/beta;
	    		else                signbeta = 1.;
			q    = -.5*(beta+signbeta*rad);
			pz2  = gamma/q;
			if(pz2<0) continue;
			vp   = 1/sqrt(px2+pz2);
			angle = asin(MIN(1.0,pj*vp));

		}

		
		/* extrapolate times and amplitudes for interpolation */
		for (jtau=ktau+1; jtau<ntau; ++jtau) {
			t[jp][jtau] = t[jp][ktau];
			a[jp][jtau] = a[jp][ktau];
		}
		
		/* print turning time and max,min tau */
		/*
		if (verbose)
			fprintf(stderr,"p=%g tturn=%g taumax=%g taumin=%g\n",
				p[jp],tturn[jp],taumax[jp],taumin[jp]);
		*/
	}
	return table;
}