sumigpreffd.c

su 的源代码库

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

		
			/* get the minimum velocity */
			vmin=v[iz][0];
			for(ix=0;ix<nx;ix++){
				if(v[iz][ix]<vmin)vmin=v[iz][ix];
			}
		
			/* compute time-invariant 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 parts of the cq[][] and cq1[][] arrays */
			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;

			/* apply phase shift */
			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);
					}
				}
			}
	
			/* fourier transform */
			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);
		 
				}
			}
				
			/* apply fdmig algorithm */
			fdmig( cp, nx, nw,v[iz],fw,dw,dz,dx,dt,v1,dip);
			fdmig( cp1,nx, nw,v[iz],fw,dw,dz,dx,dt,v1,dip);

		}

		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 fdmig( complex **cp, int nx, int nw, float *v,float fw,float
	dw,float dz,float dx,float dt,float vc,int dip)
{
	int iw,ix;
	float *p,*s1,*s2,w,coefa,coefb,v1,vn,trick=0.1;
	complex cp2,cp3,cpnm1,cpnm2;
	complex a1,a2,b1,b2;
	complex endl,endr;
	complex *data,*d,*a,*b,*c;

	p=alloc1float(nx);
	s1=alloc1float(nx);
	s2=alloc1float(nx);

	data=alloc1complex(nx);
	d=alloc1complex(nx);
	a=alloc1complex(nx);
	b=alloc1complex(nx);
	c=alloc1complex(nx);

	for(ix=0;ix<nx;ix++){
		p[ix]=vc/v[ix];
		p[ix]=(p[ix]*p[ix]+p[ix]+1.0);
	}

	
	if(dip!=65){
		coefa=0.5;coefb=0.25;
	} else {
		coefa=0.4784689;
		coefb=0.37607656;
	}

	v1=v[0];
	vn=v[nx-1];

	for(iw=0,w=fw;iw<nw;iw++,w+=dw){
		if(fabs(w)<=1.0e-10)w=1.0e-10/dt; 

		for(ix=0;ix<nx;ix++){
			s1[ix]=(v[ix]*v[ix])*p[ix]*coefb/(dx*dx*w*w)+trick;
			s2[ix]=-(1-vc/v[ix])*v[ix]*dz*coefa/(w*dx*dx)*0.5;
		}

		for(ix=0;ix<nx;ix++){
			data[ix]=cp[iw][ix];
		}

		cp2=data[1];
		cp3=data[2];
		cpnm1=data[nx-2];
		cpnm2=data[nx-3];
		a1=crmul(cmul(cp2,conjg(cp3)),2.0);
		b1=cadd(cmul(cp2,conjg(cp2)),cmul(cp3,conjg(cp3)));

		if(b1.r==0.0 && b1.i==0.0)
			a1=cexp(cmplx(0.0,-w*dx*0.5/v1));
		else
			a1=cdiv(a1,b1);

		if(a1.i>0.0)
			a1=cexp(cmplx(0.0,-w*dx*0.5/v1));

		a2=crmul(cmul(cpnm1,conjg(cpnm2)),2.0);
		b2=cadd(cmul(cpnm1,conjg(cpnm1)),cmul(cpnm2,conjg(cpnm2)));

		if(b2.r==0.0 && b2.i==0.0)
			a2=cexp(cmplx(0.0,-w*dx*0.5/vn));
		else
			a2=cdiv(a2,b2);

		if(a2.i>0.0)
			a2=cexp(cmplx(0.0,-w*dx*0.5/vn));

		for(ix=0;ix<nx;ix++){
			a[ix]=cmplx(s1[ix],s2[ix]);
			b[ix]=cmplx(1.0-2.0*s1[ix],-2.0*s2[ix]);
		}

		for(ix=1;ix<nx-1;ix++){

			d[ix]=cadd(cadd(cmul(data[ix+1],a[ix+1]),
					cmul(data[ix-1],a[ix-1])),
					cmul(data[ix],b[ix]));
		}

		d[0]=cadd(cmul(cadd(b[0],cmul(a[0],a1)),
				data[0]),cmul(data[1],a[1]));

		d[nx-1]=cadd(cmul(cadd(b[nx-1],
			cmul(a[nx-1],a2)),data[nx-1]),
			cmul(data[nx-2],a[nx-2]));

		for(ix=0;ix<nx;ix++){
			data[ix]=cmplx(s1[ix],-s2[ix]);
			b[ix]=cmplx(1.0-2.0*s1[ix],2.0*s2[ix]);
		}
		endl=cadd(b[0],cmul(data[0],a1));
		endr=cadd(b[nx-1],cmul(data[nx-1],a2));

		
		for(ix=1;ix<nx-1;ix++){
			a[ix]=data[ix+1];
			c[ix]=data[ix-1];
		}
		a[0]=data[1];
		c[nx-1]=data[nx-2];
			
		retris(data,a,c,b,endl,endr,nx,d);

		for(ix=0;ix<nx;ix++){
			cp[iw][ix]=data[ix];
		}

	}


	free1complex(data);
	free1float(p);
	free1complex(d);
	free1complex(b);
	free1complex(c);
	free1complex(a);
	free1float(s1);
	free1float(s2);
		
	return;
}
		 

void retris(complex *data,complex *a,complex *c, complex *b,
		complex endl,complex endr, int nx, complex *d)
{
		 
	int ix;
	complex *e,den;
	complex *f;

	e=alloc1complex(nx);
	f=alloc1complex(nx);
	e[0]=cdiv(cneg(a[0]),endl);
	f[0]=cdiv(d[0],endl);

	for(ix=1;ix<nx-1;++ix){
		den=cadd(b[ix],cmul(c[ix],e[ix-1]));
		e[ix]=cdiv(cneg(a[ix]),den);
		f[ix]=cdiv(csub(d[ix],cmul(f[ix-1],c[ix])),den);
	}
		 

	data[nx-1]=cdiv(csub(d[nx-1],cmul(f[nx-2],c[nx-2])),
			cadd(endr,cmul(c[nx-2],e[nx-2])));
		
	for(ix=nx-2;ix>-1;--ix)
	data[ix]=cadd(cmul(data[ix+1],e[ix]),f[ix]);

	free1complex(e);
	free1complex(f);
	return;  
}

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