sumigps.c

su 的源代码库

	complex *pt,*pw,*qn,*qt;
	static int tabbed=0;
	static float fntab,*ctab,*stab,opi2=1.0/(PI*2.0);

	/* if not already built, build cosine/sine tables */
	if (!tabbed) {
		int itab,ntab=1025;
		float angle,dangle=2.0*PI/(ntab-1);
		ctab = alloc1float(ntab);
		stab = alloc1float(ntab);
		for (itab=0,angle=0.0; itab<ntab; ++itab,angle+=dangle) {
			ctab[itab] = cos(angle);
			stab[itab] = sin(angle);
		}
		tabbed = 1;
		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 vtau[],
	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
vtau		array[ntau] of velocities 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,
		vel,dvel,angle,cosa,frac,
		*p,*taumax,*taumin,*tturn,**t,**a;
	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);
		
	/* loop over slopes p */
	for (jp=0; jp<np; ++jp) {
	
		/* slope p */
		p[jp] = pj = 2.0/vtau[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;
		vel = 0.5*vtau[0];
		dvel = 0.5*(vtau[1]-vtau[0]);
		angle = asin(MIN(1.0,pj*vel));
		
		/* 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);
			
			/* 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;
			}
			
			/* update angle */
			angle += pj*dvel;
			
			/* update velocity and first difference */
			if (itau<ntau-1) {
				frac = (taui-(itau*dtau))/dtau;
				dvel = 0.5*(vtau[itau+1]-vtau[itau]);
				vel = 0.5*((1.0-frac)*vtau[itau] +
					frac*vtau[itau+1]);
			} else {
				dvel = 0.5*(vtau[ntau-1]-vtau[ntau-2]);
				vel = 0.5*vtau[ntau-1] +
					(taui-(ntau-1)*dtau)*dvel/dtau;
			}
		}
		
		/* 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;
}

/* for graceful interrupt termination */
static void closefiles(void)
{
	efclose(headerfp);
	efclose(tracefp);
	eremove(headerfile);
	eremove(tracefile);
	exit(EXIT_FAILURE);
}