sumigtopo2d.c

su 的源代码库

			t[ix][iz] = a1*t1[ix][iz]+a2*t2[ix][iz];
}

/* compute  reference traveltime and slowness   */
      void timeb(int nr,int nz,float dr,float dz,float fz,float z0,float a,
        float v0,float **t,float **p,float **sig,float **ang)
{
        int  ir,iz;
        float r,z,v,rc,oa,temp,rou,zc;


        if( a==0.0) {
                for(ir=0,r=0;ir<nr;++ir,r+=dr)
                        for(iz=0,z=fz-z0;iz<nz;++iz,z+=dz){
                                rou = sqrt(r*r+z*z);
                                if(rou<dz) rou = dz;
                                t[ir][iz] = rou/v0;
                                p[ir][iz] = r/(rou*v0);
                                sig[ir][iz] = v0*rou;
                                ang[ir][iz] = asin(r/rou);
                        }
        } else {
                oa = 1.0/a;     zc = v0*oa;
                for(ir=0,r=0;ir<nr;++ir,r+=dr)
                        for(iz=0,z=fz+zc-z0;iz<nz;++iz,z+=dz){
                                rou = sqrt(r*r+z*z);
                                v = v0+a*(z-zc+z0);
                                if(ir==0){
                                        t[ir][iz] = log(v/v0)*oa;
                                        p[ir][iz] = 0.0;
                                        ang[ir][iz] = 0.0;
                                        sig[ir][iz] = 0.5*(z+0.1*dz-zc+z0)
						*(v0+v);
                                } else {
                                        rc = (r*r+z*z-zc*zc)/(2.0*r);
                                        rou = sqrt(zc*zc+rc*rc);
                                        t[ir][iz] = log((v*(rou+rc))
                                                /(v0*(rou+rc-r)))*oa;
                                        p[ir][iz] = sqrt(rou*rou-rc*rc)
                                                /(rou*v0);
                                        temp = v0*p[ir][iz];
                                        if(temp>1.0) temp = 1.0;
                                        ang[ir][iz] = asin(temp);
                                        sig[ir][iz] = a*rou*r;
                                }
                        }
        }
}

void filt(float *trace,int nt,float dt,float fmax,int ls,int m,float *trf);

  void mig2d(float *trace,int nt,float ft,float dt,
	float sx,float gx,float **mig,float aperx,
  	int nx,float fx,float dx,float nz,float fz,float dz,
	int ls,int mtmax,float dxm,float fmax,float angmax,
        float **tbs,float **tbg,float **pbs,
        float **pbg,float **angbs,float **angbg,int nr,float **tsum,
        int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,float *szif,
        float nangls,float nanglg,float **sigbs,float **sigbg)
/*****************************************************************************
Migrate one trace 
******************************************************************************
Input:
*trace		one seismic trace 
nt		number of time samples in seismic trace
ft		first time sample of seismic trace
dt		time sampleing interval in seismic trace
sx,gx		lateral coordinates of source and geophone 
aperx		lateral aperature in migration
nx,fx,dx,nz,fz,dz  dimension parameters of migration region
ls		=1 for line source; =0 for point source
mtmax		number of time samples in triangle filter
dxm		midpoint sampling interval
fmax		frequency-highcut for input trace	 
angmax		migration angle aperature from vertical 	 
(tb,pb,sigb,angb)s,g	reference traveltime, lateral slowness, cosine of 
		incident angle, and emergent angle
nr		number of lateral samples in reference quantities
tsum		sum of residual traveltimes from shot and receiver
nxt,fxt,dxt,nzt,fzt,dzt		dimension parameters of traveltime table
szif            array[] of the z component of the recording surface (ignored)
nangls,nanglg   Angles that the normal unit vector of the current surface
                forms with the vertical at source and geophone locations.
Output:
mig		migrated section
*****************************************************************************/
{
	int nxf,nxe,nxtf,nxte,ix,iz,iz0,izt0,nzp,jrs,jrg,jz,jt,mt,jx;
	float xm,x,dis,rxz,ar,srs,srg,srs0,srg0,sigp,z0,rdz,ampd,res0,
	      angs,angg,ax,ax0,pmin,odt=1.0/dt,pd,az,sz,sz0,sigmas,sigmag,
	      at,td,res,temp,sigs,sigg;
	float **tmt,**ampt,**ampti,*tm,*amp,*ampi,*tzt,*trf,*zpt;

	tmt = alloc2float(nzt,nxt);
	ampt = alloc2float(nzt,nxt);
	ampti = alloc2float(nzt,nxt);
	tm = szif; /* dummy to keep compiler happy */
	tm = alloc1float(nzt);
	tzt = alloc1float(nzt);
	amp = alloc1float(nzt);
	ampi = alloc1float(nzt);
	zpt = alloc1float(nxt);
	trf = alloc1float(nt+2*mtmax);

	z0 = (fz-fzt)/dzt;
	rdz = dz/dzt;
	pmin = 1.0/(2.0*dxm*fmax);
	
	filt(trace,nt,dt,fmax,ls,mtmax,trf);

	xm = 0.5*(sx+gx);
	rxz = (angmax==90)?0.0:1.0/tan(angmax*PI/180.);
	nxtf = (xm-aperx-fxt)/dxt;
	if(nxtf<0) nxtf = 0;
	nxte = (xm+aperx-fxt)/dxt+1;
	if(nxte>=nxt) nxte = nxt-1;

	/* compute amplitudes and filter length	*/
	for(ix=nxtf; ix<=nxte; ++ix){
		x = fxt+ix*dxt;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = ((dis-dxt)*rxz-fzt)/dzt-1;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;

		ar = (sx>=x)?(sx-x)/dx:(x-sx)/dx;
		jrs = (int)ar;
		if(jrs>nr-2) jrs = nr-2;
		srs = ar-jrs;
		srs0 = 1.0-srs;
		ar = (gx>=x)?(gx-x)/dx:(x-gx)/dx;
		jrg = (int)ar;
		if(jrg>nr-2) jrg = nr-2;
		srg = ar-jrg;
		srg0 = 1.0-srg;
		sigp = ((sx-x)*(gx-x)>0)?1.0:-1.0;
                sigs = ((sx-x)<0)?1.0:-1.0;
                sigg = ((gx-x)<0)?1.0:-1.0;
		zpt[ix] = fzt+(nzt-1)*dzt;

		for(iz=izt0; iz<nzt; ++iz){
			angs = srs0*angbs[jrs][iz]+srs*angbs[jrs+1][iz];
                        angs = ABS(sigs*angs - nangls);
			angg = srg0*angbg[jrg][iz]+srg*angbg[jrg+1][iz];
                        angg = ABS(sigg*angg - nanglg);
			sigmas = srs0*sigbs[jrs][iz]+srs*sigbs[jrs+1][iz]; 
			sigmag = srg0*sigbg[jrg][iz]+srg*sigbg[jrg+1][iz]; 
                        ampd = cos(angs)*sigmag + cos(angg)*sigmas;
			ampd /= (cos(nangls)*cos(nanglg)*sqrt(sigmas*sigmag));

        /* Filter of 90 degrees to the operator (Peels,1988)

           Peels, G. L., 1988, True amplitude wave field extrapolation
           with applications in  seismic shot record redatuming, PhD
           thesis, Delft University of Technology.              */

                        if(ABS(angs)>=(PI/2) || ABS(angg)>=(PI/2)) {
                                ampd=0.0;
                        }

			if(ampd<0.0) ampd = -ampd;
			ampt[ix][iz] = ampd;

			pd = srs0*pbs[jrs][iz]+srs*pbs[jrs+1][iz]+sigp 
			     *(srg0*pbg[jrg][iz]+srg*pbg[jrg+1][iz]);
			if(pd<0.0) pd = -pd;
			temp = pd*dxm*odt;
			if(temp<1) temp = 1.0;
			if(temp>mtmax) temp = mtmax;
			ampti[ix][iz] = ampd/(temp*temp);
			tmt[ix][iz] = temp;
			if(pd<pmin && zpt[ix]>fzt+(nzt-1.1)*dzt) 
				zpt[ix] = fzt+iz*dzt;

		}
	}

	nxf = (xm-aperx-fx)/dx+0.5;
	if(nxf<0) nxf = 0;
	nxe = (xm+aperx-fx)/dx+0.5;
	if(nxe>=nx) nxe = nx-1;
	
	/* interpolate amplitudes and filter length along lateral	*/
	for(ix=nxf; ix<=nxe; ++ix){
		x = fx+ix*dx;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = (dis*rxz-fzt)/dzt;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;
		iz0 = (dis*rxz-fz)/dz;
		if(iz0<0) iz0 = 0;
		if(iz0>=nz) iz0 = nz-1;

		ax = (x-fxt)/dxt;
		jx = (int)ax;
		ax = ax-jx;
		if(ax<=0.01) ax = 0.;
		if(ax>=0.99) ax = 1.0;
		ax0 = 1.0-ax;
		if(jx>nxte-1) jx = nxte-1;
		if(jx<nxtf) jx = nxtf;

		ar = (sx>=x)?(sx-x)/dx:(x-sx)/dx;
		jrs = (int)ar;
		if(jrs>nr-2) jrs = nr-2;
		srs = ar-jrs;
		srs0 = 1.0-srs;
		ar = (gx>=x)?(gx-x)/dx:(x-gx)/dx;
		jrg = (int)ar;
		if(jrg>nr-2) jrg = nr-2;
		srg = ar-jrg;
		srg0 = 1.0-srg;

		for(iz=izt0; iz<nzt; ++iz){
		    tzt[iz] = ax0*tsum[jx][iz]+ax*tsum[jx+1][iz]
				+srs0*tbs[jrs][iz]+srs*tbs[jrs+1][iz]
				+srg0*tbg[jrg][iz]+srg*tbg[jrg+1][iz];

		    amp[iz] = ax0*ampt[jx][iz]+ax*ampt[jx+1][iz];
		    ampi[iz] = ax0*ampti[jx][iz]+ax*ampti[jx+1][iz];
		    tm[iz] = ax0*tmt[jx][iz]+ax*tmt[jx+1][iz];

		}

		nzp = (ax0*zpt[jx]+ax*zpt[jx+1]-fz)/dz+1.5;
		if(nzp<iz0) nzp = iz0;
		if(nzp>nz) nzp = nz;

		/* interpolate along depth if operater aliasing 	*/
		for(iz=iz0; iz<nzp; ++iz) {
			az = z0+iz*rdz;
			jz = (int)az;
			if(jz>=nzt-1) jz = nzt-2;
			sz = az-jz;
			sz0 = 1.0-sz;
			td = sz0*tzt[jz]+sz*tzt[jz+1];
			at = (td-ft)*odt+mtmax;
			jt = (int)at;
			if(jt > mtmax && jt < nt+mtmax-1){
			    ampd = sz0*ampi[jz]+sz*ampi[jz+1];
			    mt = (int)(0.5+sz0*tm[jz]+sz*tm[jz+1]);
			    res = at-jt;
			    res0 = 1.0-res;
 			    temp = (res0*(-trf[jt-mt]+2.0*trf[jt]-trf[jt+mt]) 
				+res*(-trf[jt-mt+1]+2.0*trf[jt+1]
				-trf[jt+mt+1]))*ampd;				
			    mig[ix][iz] += temp;

			}
		}

		/* interpolate along depth if not operater aliasing 	*/
		for(iz=nzp; iz<nz; ++iz) {
			az = z0+iz*rdz;
			jz = (int)az;
			if(jz>=nzt-1) jz = nzt-2;
			sz = az-jz;
			sz0 = 1.0-sz;
			td = sz0*tzt[jz]+sz*tzt[jz+1];
			at = (td-ft)*odt;
			jt = (int)at;
			if(jt > 0 && jt < nt-1){
			    ampd = sz0*amp[jz]+sz*amp[jz+1];
			    res = at-jt;
			    res0 = 1.0-res;
 			    temp = (res0*trace[jt]+res*trace[jt+1])*ampd; 
			    mig[ix][iz] += temp;
			}
		}

	}

	free2float(ampt);
	free2float(ampti);
	free2float(tmt);
	free1float(amp);
	free1float(ampi);
	free1float(zpt);
	free1float(tm);
	free1float(tzt);
	free1float(trf);
}

void filt(float *trace,int nt,float dt,float fmax,int ls,int m,float *trf)
/*********************************************************************
  Low-pass filter, integration and phase shift for input data	 
  input: 
    trace(nt)	single seismic trace
    fmax	high cut frequency
    ls		ls=1, line source; ls=0, point source
  output:
    trace(nt) 	filtered and phase-shifted seismic trace 
    tracei(nt) 	filtered, integrated and phase-shifted seismic trace 
**********************************************************************/
{
	static int nfft=0, itaper, nw, nwf;
	static float *taper, *amp, *ampi, dw;
	int  it, iw, itemp;
	float temp, ftaper, const2, *rt;
	complex *ct;

	fmax *= 2.0*PI;
	ftaper = 0.1*fmax;
	const2 = 0.5*sqrt(2.0);

	if(nfft==0) {
        	/* Set up FFT parameters */
        	nfft = npfaro(nt+2*m, 2*(nt+2*m));
        	if (nfft >= SU_NFLTS || nfft >= 720720)
                	err("Padded nt=%d -- too big", nfft);

        	nw = nfft/2 + 1;
		dw = 2.0*PI/(nfft*dt);

		itaper = 0.5+ftaper/dw;
		taper = ealloc1float(2*itaper+1);
		for(iw=-itaper; iw<=itaper; ++iw){
			temp = (float)iw/(1.0+itaper); 
			taper[iw+itaper] = (1-temp)*(1-temp)*(temp+2)/4;
		}

		nwf = 0.5+fmax/dw;
		if(nwf>nw-itaper-1) nwf = nw-itaper-1;
		amp = ealloc1float(nwf+itaper+1);
		ampi = ealloc1float(nwf+itaper+1);
		amp[0] = ampi[0] = 0.;
		for(iw=1; iw<=nwf+itaper; ++iw){
			amp[iw] = sqrt(dw*iw)/nfft;
			ampi[iw] = 0.5/(1-cos(iw*dw*dt));
		}
	}

        /* Allocate fft arrays */
        rt   = ealloc1float(nfft);
        ct   = ealloc1complex(nw);

        memcpy(rt, trace, nt*FSIZE);
        memset((void *) (rt + nt), 0, (nfft-nt)*FSIZE); 
        pfarc(1, nfft, rt, ct);

	for(iw=nwf-itaper;iw<=nwf+itaper;++iw){
		itemp = iw-(nwf-itaper);
		ct[iw].r = taper[itemp]*ct[iw].r; 
		ct[iw].i = taper[itemp]*ct[iw].i; 
	}
	for(iw=nwf+itaper+1;iw<nw;++iw){
		ct[iw].r = 0.; 
		ct[iw].i = 0.; 
	}

       	if(!ls){
		for(iw=0; iw<=nwf+itaper; ++iw){
			/* phase shifts PI/4 + PI/4 (Half derivative)  */
			temp = -ct[iw].i*amp[iw];
			ct[iw].i = ct[iw].r*amp[iw];
			ct[iw].r = temp;
		    }
	} else {
		for(iw=0; iw<=nwf+itaper; ++iw){
                        /* phase shifts PI/4 -Half derivative   */
                        temp = (ct[iw].r-ct[iw].i)*amp[iw]*const2;
                        ct[iw].i = (ct[iw].r+ct[iw].i)*amp[iw]*const2;
                        ct[iw].r = temp;
		}
	}              
        pfacr(-1, nfft, ct, rt);
		
        /* Load traces back in */
	for (it=0; it<nt; ++it) trace[it] = rt[it];

        /* Integrate traces   */
	for(iw=0; iw<=nwf+itaper; ++iw){
		ct[iw].i = ct[iw].i*ampi[iw];
		ct[iw].r = ct[iw].r*ampi[iw];
	}
        pfacr(-1, nfft, ct, rt);
        for (it=0; it<m; ++it)  trf[it] = rt[nfft-m+it];
        for (it=0; it<nt+m; ++it)  trf[it+m] = rt[it];

	free1float(rt);
	free1complex(ct);
}