sumigpresp.c

su 的源代码库

			oldigx = igx;

                        if(gxmin>gx)gxmin=gx;
                        if(gxmax<gx)gxmax=gx;

                        if(verbose)
                                warn(" inside loop:  min_sx_gx %f isx %d igx %d gx %f sx %f",min_sx_gx,isx,igx,gx,sx);
                        /* get sx and gx */
                        get_sx_gx(&sx,&gx);
                        sx = (sx - min_sx_gx);
                        gx = (gx - min_sx_gx);

                        /* sx, gx must increase monotonically */
                        if (!(oldsx <= sx) )
                         err("sx field must be monotonically increasing!");
                        if (!(oldgx <= gx) )
                         err("gx field must be monotonically increasing!");

                        ++nx;
                } while(gettr(&tr) && sx==oldsx);


                isx=oldsx/dx;
		if (isx==oldisx) 
			warn("repeated isx!!! check dx or scalco value!!!");
		oldisx=isx;
                ixshot=isx;
                if(verbose) {
                        warn("sx %f, gx %f , gxmin %f  gxmax %f nx %d",sx,gx,gxmin,gxmax, nx);
                        warn("isx %d igx %d ixshot %d" ,isx,igx,ixshot);
                }


		/* transform the shot gather from time to frequency domain */
		pfa2rc(1,1,ntfft,nxo,p[0],cq[0]);

                /* compute the most left and right index for the migrated */
                /* section */
                ix1=oldsx/dx;
                ix2=oldgxmin/dx;
                ix3=oldgxmax/dx;

                if(ix1>=ix3)ix3=ix1;
                if(ix1<=ix2)ix2=ix1;

                ix2-=lpad;
                ix3+=rpad;
                if(ix2<0)ix2=0;
                if(ix3>nxo-1)ix3=nxo-1;

                /* the total traces to be migrated */
                nx=ix3-ix2+1;
                nw=truenw;

		/* determine wavenumber sampling (for complex to complex FFT) */
		nxfft = npfa(nx);
		nk = nxfft;
		dk = 2.0*PI/(nxfft*dx);
		fk = -PI/dx;
	

		/* allocate space for velocity profile within the aperature */
		v=alloc2float(nx,nz);   
		for(iz=0;iz<nz;iz++)
			for(ix=0;ix<nx;ix++)
				v[iz][ix]=vp[iz][ix+ix2];
		

		/* allocate space */
		cp = alloc2complex(nx,nw);
		cp1 = alloc2complex(nx,nw);

                /* transpose the frequency domain data from     */
                /* data[ix][iw] to data[iw][ix] and apply a     */
                /* Hamming at the same time                     */
		for (ix=0; ix<nx; ix++) {
			for (iw=0; iw<nw; iw++){
				float tmpp=0.0,tmppp=0.0;

				if(iw>=(nf1-nf1)&&iw<=(nf2-nf1)){
					tmpp=PI/(nf2-nf1);
					tmppp=tmpp*(iw-nf1)-PI;
					tmpp=0.54+0.46*cos(tmppp);
					cp[iw][ix]=crmul(cq[ix+ix2][iw+nf1],tmpp);
				} else {
					if(iw>=(nf3-nf1)&&iw<=(nf4-nf1)){
						tmpp=PI/(nf4-nf3);
						tmppp=tmpp*(iw-nf3);
						tmpp=0.54+0.46*cos(tmppp);
						cp[iw][ix]=crmul(cq[ix+ix2][iw+nf1],tmpp);
					} else {
						cp[iw][ix]=cq[ix+ix2][iw+nf1];}
				}
				cp1[iw][ix]=cmplx(0.0,0.0);
			}
		}

		for(iw=0;iw<nw;iw++){
			cp1[iw][ixshot-ix2]=wlsp[iw+nf1];
		}

                if(verbose) {
                        warn("ixshot %d ix %d ix1 %d ix2 %d ix3 %d",ixshot,ix,ix1,ix2,ix3);
                        warn("oldsx %f ",oldsx);
                }

			
		free2float(p);
		free2complex(cq);
		free1float(wl);
		free1complex(wlsp);

		/* allocating space */
		cq=alloc2complex(nxfft,nw);
		cq1=alloc2complex(nxfft,nw);


		/* loops over depth */
		for(iz=0;iz<nz;++iz){

			/* the imaging condition */
			for(ix=0;ix<nx;ix++){
				for(iw=0,w=fw;iw<nw;w+=dw,iw++){   
					complex tmp;
					float ratio=10.0;
		
					if(fabs(ix+ix2-ixshot)*dx<ratio*iz*dz)
						tmp=cmul(cp[iw][ix],cp1[iw][ix]);
					else 
						tmp=cmplx(0.0,0.0);  

					cresult[ix+ix2][iz]+=tmp.r/ntfft;
				}
			}

			/* get the minimum velocity */
			vmin=0;
			for(ix=il-ix2;ix<=ir-ix2;ix++){
				vmin+=1.0/v[iz][ix]/(ir-il+1);
			}
			vmin=1.0/vmin;
		
			/* compute the shifted wavefield */
			for (ik=0;ik<nx;++ik) {
				for (iw=0; iw<nw; ++iw) {
					cq[iw][ik] = ik%2 ? cneg(cp[iw][ik]) : cp[iw][ik];
					cq1[iw][ik] = ik%2 ? cneg(cp1[iw][ik]) : cp1[iw][ik];
				}
			}
		 
			/* zero out cq[][] cq1[][] */
			for (ik=nx; ik<nk; ++ik) {
				for (iw=0; iw<nw; ++iw) {
					cq[iw][ik] = cmplx(0.0,0.0);
					cq1[iw][ik] = cmplx(0.0,0.0);
				}
			}

			/* FFT to W-K domain */
			pfa2cc(-1,1,nk,nw,cq[0]);
			pfa2cc(-1,1,nk,nw,cq1[0]);
	
			v1=vmin;
			for(ik=0,k=fk;ik<nk;++ik,k+=dk) {
				for(iw=0,w=fw;iw<nw;++iw,w+=dw){
					if(w==0.0)w=1.0e-10/dt; 
					kz1=1.0-pow(v1*k/w,2.0);
					if(kz1>0.15){
						phase1 = -w*sqrt(kz1)*dz/v1;
						cshift1 = cmplx(cos(phase1), sin(phase1));
						cq[iw][ik] = cmul(cq[iw][ik],cshift1);
						cq1[iw][ik] = cmul(cq1[iw][ik],cshift1);
					} else {
						cq[iw][ik] = cq1[iw][ik] = cmplx(0.0,0.0);
					}
				}
			}
	
			pfa2cc(1,1,nk,nw,cq[0]);
			pfa2cc(1,1,nk,nw,cq1[0]);

			for(ix=0;ix<nx;++ix) {
				for(iw=0,w=fw;iw<nw;w+=dw,++iw){
					float a=0.015,g=1.0;
					int I=10;
				
					if(ix<=I)g=exp(-a*(I-ix)*(I-ix));
					if(ix>=nx-I)g=exp(-a*(-nx+I+ix)*(-nx+I+ix));
				 
				
					cq[iw][ix] = crmul( cq[iw][ix],1.0/nxfft);
					cq[iw][ix] =ix%2 ? cneg(cq[iw][ix]) : cq[iw][ix];
					kz2=(1.0/v1-1.0/v[iz][ix])*w*dz;
					cshift2=cmplx(cos(kz2),sin(kz2));
					cp[iw][ix]=cmul(cq[iw][ix],cshift2);
		
					cq1[iw][ix] = crmul( cq1[iw][ix],1.0/nxfft);
					cq1[iw][ix] =ix%2 ? cneg(cq1[iw][ix]) : cq1[iw][ix];
					cp1[iw][ix]=cmul(cq1[iw][ix],cshift2);
		 
				}
			}
		
				

		}

		free2complex(cp);
		free2complex(cp1);
		free2complex(cq);
		free2complex(cq1);
		free2float(v);

		--nxshot;

	} while(nxshot);



        /* restore header fields and write output */
        for(ix=0; ix<nxo; ix++){
                tr.ns = nz;
                tr.d1 = dz;
                tr.d2 = dx;
                tr.offset = 0;
                tr.cdp = tr.tracl = ix;
                memcpy( (void *) tr.data, (const void *) cresult[ix],nz*FSIZE);
                puttr(&tr);
        }


	return(CWP_Exit());	

}



float * ricker(float Freq,float dt,int *Npoint)
{
	int i; /* they are the dummy counter*/
	float Bpar,t,u,*Amp;
	int Np1,N;
	
	if(Freq==0.0)Freq=30.0;
	if(dt==0.0)dt=0.004;
	Bpar=sqrt(6.0)/(PI*Freq);
	N=ceil(1.35*Bpar/dt);
	Np1=N;
	*Npoint=2*N+1;
	 
	Amp=alloc1float(*Npoint);
	
	Amp[Np1]=1.0;
  
	for(i=1;i<=N;i++) {
		t=dt*(float)i;
		u=2.0*sqrt(6.0)*t/Bpar;
		Amp[Np1+i]=Amp[Np1-i]=0.5*(2.0-u*u)*exp(-u*u/4.0);
	}

	return Amp;

}

void get_sx_gx(float *sx, float *gx)
{ 
/*****************************************************************************
get_sx_gx - get sx and gx from headrs
*****************************************************************************/

	float sy;		/* source coordinates */
	float gy;		/* geophone coordinates */

		if (tr.scalco) { /* if tr.scalco is set, apply value */
			if (tr.scalco>0) {
				*sx = tr.sx*tr.scalco;
				*gx = tr.gx*tr.scalco;
				sy = tr.sy*tr.scalco;
				gy = tr.gy*tr.scalco;
			} else { /* if tr.scalco is negative divide */
				*sx = tr.sx/ABS(tr.scalco);
				*gx = tr.gx/ABS(tr.scalco);
				sy = tr.sy/ABS(tr.scalco);
				gy = tr.gy/ABS(tr.scalco);
			}

		} else {
				*sx = tr.sx;
				*gx = tr.gx;
				sy = tr.sy;
				gy = tr.gy;
		}

		
		/* use pythagorean theorem to remap radial direction */
		/* to x-direction */
		*sx = SGN(*sx-sy)*sqrt((*sx)*(*sx) + sy*sy);
		*gx = SGN(*gx-gy)*sqrt((*gx)*(*gx) + gy*gy);

	return;
}