sumigprepspi.c

su 的源代码库

				vmax=v[0][0];vmin=v[0][0];
	
		for(iz=0;iz<nz;++iz) {
			for(ix=0;ix<nx;++ix) {
	 			if(v[iz][ix]>=vmax) vmax=v[iz][ix];
	 			if(v[iz][ix]<=vmin) vmin=v[iz][ix];
			}
		}


		/* 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);

		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;
				}
			}

			for (ik=0; ik<nx; ++ik) {
				for (iw=0,w=fw; iw<nw;w+=dw, ++iw){
					cp[iw][ik]=cmul(cp[iw][ik],
							cexp(cmplx(0.0,-w*dz/v[iz][ik])));
					cp1[iw][ik]=cmul(cp1[iw][ik],
							cexp(cmplx(0.0,-w*dz/v[iz][ik])));
					cq[iw][ik] = ik%2 ? cneg(cp[iw][ik]) : cp[iw][ik];
					cq1[iw][ik] = ik%2 ? cneg(cp1[iw][ik]) : cp1[iw][ik];
				}
			}
	
			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]);
	
			/* The second time phase shift */
			v1=vmin;
			v2=vmax;
	
			if((v2-v1)/v1<0.01){
		
				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.0){
							phase1 = -w*sqrt(kz1)*dz/v1+w*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 {
							phase1 = -w*sqrt(-kz1)*dz/v1;
							cshift1=cexp(cmplx(phase1,w*dz/v1));
							cq[iw][ik] = cmul(cq[iw][ik],cshift1);
							cq1[iw][ik] = cmul(cq1[iw][ik],cshift1);
						}
					}		
				}
				pfa2cc(1,1,nk,nw,cq[0]);
				pfa2cc(1,1,nk,nw,cq1[0]);
	
				for(ix=0;ix<nx;++ix){
					for(iw=0;iw<nw;++iw){
						cq[iw][ix] = crmul( cq[iw][ix], 1.0/nxfft);
						cp[iw][ix] =ix%2 ? cneg(cq[iw][ix]) : cq[iw][ix];
						cq1[iw][ix] = crmul( cq1[iw][ix], 1.0/nxfft);
						cp1[iw][ix] =ix%2 ? cneg(cq1[iw][ix]) : cq1[iw][ix];
		 
					}
				}
			} else {
		
				for(ik=0;ik<=L;ik++)
					c[ik]=vmin+ik*1.0*(vmax-vmin)/((float) L);
			 
					memset((void *) p, 0, L*FSIZE);


		for(ix=0;ix<nx;ix++){
			for(ik=0;ik<L;ik++){
				if(((v[iz][ix]>=c[ik])
					&&(v[iz][ix]<c[ik+1]))
					||((ik==L-1)
					&&(v[iz][ix]==vmax))){
						P[ik]+=1.0/nx; break;
				}
			}
		}
		Sz=0.0;
		for(ik=0;ik<L;ik++) {
			if(P[ik]!=0.00) Sz=Sz-P[ik]*log(P[ik]);
		}

		Bz=exp(Sz)+0.5;
		Y[0]=0.0;
		Y[L]=1.0;
 
		for(ik=1;ik<L;ik++) {
			Y[ik]=0.0;
			for(ix=0;ix<ik;ix++)
				for(ix=0;ix<ik;ix++)
					Y[ik]=Y[ik]+P[ix];
		}
		 
		V=alloc1float(Bz+1);
		V[0]=vmin;
	
		for(ix=1;ix<=Bz;ix++){
			for(ik=0;ik<L;ik++){
				if((ix*1.0/Bz>Y[ik])&&(ix*1.0/Bz<=Y[ik+1])){	
					V[ix]=c[ik]+(ix*1.0/Bz-Y[ik])*(c[ik+1]-c[ik])/(Y[ik+1]-Y[ik]);
				break;
				}
			}
	
		}

	 
		 
		V[Bz]=V[Bz]*1.005;  
		cq2=ealloc3complex(nk,nw,Bz+1);
		cq3=ealloc3complex(nk,nw,Bz+1);

		for(ix=0;ix<Bz+1;ix++){
			for(iw=0,w=fw;iw<nw;++iw,w+=dw)
				for(ik=0,k=fk;ik<nk;++ik,k+=dk){

					if(w==0.0)w=1.0e-10/dt;

					kz1=1.0-pow(V[ix]*k/w,2.0);
					if(kz1>=0.00){
						phase1 =-w*sqrt(kz1)*dz/V[ix]+w*dz/V[ix];
						cshift1 = cexp(cmplx(0.0,phase1));
						cq2[ix][iw][ik] = cmul(cq[iw][ik],cshift1);
						cq3[ix][iw][ik] = cmul(cq1[iw][ik],cshift1);
	 
					} else {
						phase1 =-w*sqrt(-kz1)*dz/V[ix];
						cshift1 =cexp(cmplx(phase1,w*dz/V[ix]));
						cq2[ix][iw][ik] = cmul(cq[iw][ik],cshift1);
						cq3[ix][iw][ik] = cmul(cq1[iw][ik],cshift1);
		
					}

				}

		
				pfa2cc(1,1,nk,nw,cq2[ix][0]);
				pfa2cc(1,1,nk,nw,cq3[ix][0]);
	
				for(ik=0;ik<nx;++ik)
					for(iw=0,w=fw;iw<nw;w+=dw,++iw){
						float a=0.015,g=1.0;
							int I=10;
		
							if(ik<=I)g=exp(-a*(I-ik)*(I-ik));
							if(ik>=nx-I)g=exp(-a*(-nx+I+ik)*(-nx+I+ik));
							g=1.0;

							cq2[ix][iw][ik] = crmul( cq2[ix][iw][ik], g*1.0/nxfft);
							cq2[ix][iw][ik] =ik%2 ? cneg(cq2[ix][iw][ik]) : cq2[ix][iw][ik];
							cq3[ix][iw][ik] = crmul( cq3[ix][iw][ik], g*1.0/nxfft);
							cq3[ix][iw][ik] =ik%2 ? cneg(cq3[ix][iw][ik]) : cq3[ix][iw][ik];
	 
						}
					}

					for(ix=0;ix<nx;++ix) {
						for(ik=0;ik<Bz;++ik){
	
							if(((v[iz][ix]>=V[ik])&&(v[iz][ix]<V[ik+1]))) {


								v1=V[ik];v2=V[ik+1];
		
						for(iw=0,w=fw;iw<nw;w+=dw,++iw){

							a1=cq2[ik][iw][ix].r;
							a2=cq2[ik+1][iw][ix].r;
							theta1=cq2[ik][iw][ix].i ;theta2=cq2[ik+1][iw][ix].i;  
		
							a= a1*(v2-v[iz][ix])/(v2-v1)+a2*(v[iz][ix]-v1)/(v2-v1);
							theta=theta1*(v2-v[iz][ix])/(v2-v1)+theta2*(v[iz][ix]-v1)/(v2-v1);
							cp[iw][ix] =cmplx(a,theta);
	 
		
							a1=cq3[ik][iw][ix].r;a2=cq3[ik+1][iw][ix].r;
							theta1=cq3[ik][iw][ix].i ;
							theta2=cq3[ik+1][iw][ix].i;
	
							a= a1*(v2-v[iz][ix])/(v2-v1)+a2*(v[iz][ix]-v1)/(v2-v1);
							theta=theta1*(v2-v[iz][ix])/(v2-v1)+theta2*(v[iz][ix]-v1)/(v2-v1);
							cp1[iw][ix] =cmplx(a,theta);
		
		
								}

							break;
		
							}	
						}
					}

				free3complex(cq2);
				free3complex(cq3);
				free1float(V);
		
			}
		
				

		}

		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;
}