kdmig3d.c

用于石油地震资料数字处理

                        } else {
			      cosine=1.0;

                              rs=sqrt((xo-gd->sx)*(xo-gd->sx)+
                                      (yo-gd->sy)*(yo-gd->sy)+
                                      zo*zo);

                              rg=sqrt((xo-gd->gx)*(xo-gd->gx)+
                                      (yo-gd->gy)*(yo-gd->gy)+
                                      zo*zo);

			      rtotal=rs+rg;

                              #ifdef DEBUG
                              if (fabs((rs+rg)*1000.0-ttotal)/((rs+rg)*1000.0)>0.05 
				    && izo==0)
                              fprintf(gd->jpfp,
				    "at: ixo=%d,iyo=%d,izo=%d: est. ttot=%f, err=%f\n",
                              ixo,iyo,izo,
                              ttotal,((rs+rg)*1000.0-ttotal)/((rs+rg)*1000.0));
                              #endif
                        }

			/*cos(theta) y_3 (r_s+rg)(rs^2+rg^2)/(v^2 rs^2 rg^2)*/
                        amp=cosine*zo*rtotal*(1.0/rs/rs+1.0/rg/rg);

                        for (im=0;im<gd->multit;im++) {
                              if (tbuf[im][iys][ixo][izo][iyts0][ixts0]>=InfByte2     ||
		                  tbuf[im][iys][ixo][izo][iyts0+1][ixts0]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iyts0][ixts0+1]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iyts0+1][ixts0+1]>=InfByte2 ||
		                  tbuf[im][iys][ixo][izo][iytg0][ixtg0]>=InfByte2     ||
		                  tbuf[im][iys][ixo][izo][iytg0+1][ixtg0]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iytg0][ixtg0+1]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iytg0+1][ixtg0+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iyts1][ixts1]>=InfByte2   ||
		                  tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iyts1][ixts1+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iytg1][ixtg1]>=InfByte2   ||
		                  tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iytg1][ixtg1+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1+1]>=InfByte2) {
                                    fprintf(gd->jpfp," skipped ixo=%d ",ixo);
                                    continue;
                              }

                              /*traveltime table at [iys][ixo][izo][nyt][nxt]*/
                              tiys0=ys0b0*(xs0a0*tbuf[im][iys][ixo][izo][iyts0][ixts0]+
                                           xs0a *tbuf[im][iys][ixo][izo][iyts0][ixts0+1])+
			            ys0b *(xs0a0*tbuf[im][iys][ixo][izo][iyts0+1][ixts0]+
                                           xs0a *tbuf[im][iys][ixo][izo][iyts0+1][ixts0+1]);

                              tiyg0=yg0b0*(xg0a0*tbuf[im][iys][ixo][izo][iytg0][ixtg0]+
                                           xg0a *tbuf[im][iys][ixo][izo][iytg0][ixtg0+1])+
			            yg0b *(xg0a0*tbuf[im][iys][ixo][izo][iytg0+1][ixtg0]+
                                           xg0a *tbuf[im][iys][ixo][izo][iytg0+1][ixtg0+1]);

                              /*traveltime table at [iys+1][ixo][izo][nyt][nxt]*/
                              tiys1=ys1b0*(xs1a0*tbuf[im][iys+1][ixo][izo][iyts1][ixts1]+
                                           xs1a *tbuf[im][iys+1][ixo][izo][iyts1][ixts1+1])+
			            ys1b *(xs1a0*tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1]+
                                       xs1a *tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1+1]);

                              tiyg1=yg1b0*(xg1a0*tbuf[im][iys+1][ixo][izo][iytg1][ixtg1]+
                                           xg1a *tbuf[im][iys+1][ixo][izo][iytg1][ixtg1+1])+
			            yg1b *(xg1a0*tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1]+
                                        xg1a *tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1+1]);
                              ts=(ybeta0*tiys0+ybeta*tiys1)/(float)InfByte2*InfTime;
                              tg=(ybeta0*tiyg0+ybeta*tiyg1)/(float)InfByte2*InfTime;
                        
                              ts=MAX(ts,0.001);
                              tg=MAX(tg,0.001);

                              ttotal=ts+tg;

                              if (ttotal>=InfTime ||
                                    ttotal<=0.0) continue;

	                      tintr=(ttotal-gd->ft)/gd->dt;
	                      itintr=(int)tintr;

	                      if (itintr-ntri>0 && itintr+ntri+1 < gd->nt){
	       	                    wtintr=tintr-itintr;
		                    wtintr0=1.0-wtintr;

 		                    mig[iyo][ixo][izo]+=amp*(wtintr0*
                           (-trf[itintr-ntri]+2.0*trf[itintr]-trf[itintr+ntri])
		                          +wtintr *
                           (-trf[itintr-ntri+1]+2.0*trf[itintr+1]-trf[itintr+ntri+1]));
			      }
                        }
                  }
	    }
      }
}

void get_bufs( struct GD *gd, float ******ttab, float *****ctab,
      float *****rtab, unsigned short ******tbuf, unsigned char *****cbuf,
      unsigned char  *****rbuf)
/**********************************************************************
get_bufs - Traveltime buffer. This will dramatically speed up the code.
If you could even give a larger buffer (use nyo instead of
nys), you can further speed up the code.
***********************************************************************
Function prototype:
void get_bufs( struct GD *gd, float ******ttab, float *****ctab,
      float *****rtab, unsigned short ******tbuf, unsigned char *****cbuf,
      unsigned char  *****rbuf);
***********************************************************************
Inputs:
struct *gd   grid information
float ******ttab multivalued traveltime table

outputs:
float *****ctab           cosine buffer
float *****rtab           ray path buffer
unsigned short ******tbuf multivalued traveltime buffer
***********************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
***********************************************************************/
{
      int iyt;	        /*index for nyt*/
      int ixt;	        /*index for nxt*/
      int iys;          /*index for nys*/
      int ixs;          /*index for nxs*/
      int izs;	        /*index for nzs*/
      int izo;          /*index for nzo*/
      int ixo;          /*index for nxo*/      
      int im;

      float xalpha,xalpha0;
      float zgamma,zgamma0;	

      float xo,xs;
      float zo,zs;
      float t,c,r;

      for (ixo=0;ixo<gd->nxo;ixo++) {
            xo=gd->fxo+ixo*gd->dxo;
    	    xs=(xo-gd->fxs)/gd->dxs;
            ixs=(int)xs;
            ixs=MIN(ixs,gd->nxs-2);
            ixs=MAX(ixs,0);
            xalpha=xs-ixs;
            xalpha0=1.0-xalpha;

            for (izo=0;izo<gd->nzo;izo++) {
                  zo=gd->fzo+izo*gd->dzo;
    	          zs=(zo-gd->fzs)/gd->dzs;
                  izs=(int)zs;
	          izs=MIN(izs,gd->nzs-2);
	          izs=MAX(izs,0);
                  zgamma=zs-izs;
                  zgamma0=1.0-zgamma;

                  for (iys=0;iys<gd->nys;iys++) {
		        for (iyt=0;iyt<gd->nyt;iyt++) {
			      for (ixt=0;ixt<gd->nxt;ixt++) {
                                    if (ttab[0][iys][ixs  ][izs][iyt][ixt]>=InfTime ||
					ttab[0][iys][ixs+1][izs][iyt][ixt]>=InfTime ||
					ttab[0][iys][ixs  ][izs+1][iyt][ixt]>=InfTime ||
					ttab[0][iys][ixs+1][izs+1][iyt][ixt]>=InfTime) {
				          tbuf[0][iys][ixo][izo][iyt][ixt]=InfByte2;
                                          cbuf[iys][ixo][izo][iyt][ixt]=0;
					  rbuf[iys][ixo][izo][iyt][ixt]=InfByte1;
				    } else {
                                          t=zgamma0*(
                                              xalpha0*ttab[0][iys][ixs  ][izs][iyt][ixt]+
                                              xalpha *ttab[0][iys][ixs+1][izs][iyt][ixt])+
                                            zgamma *(
					      xalpha0*ttab[0][iys][ixs  ][izs+1][iyt][ixt]+
                                              xalpha *ttab[0][iys][ixs+1][izs+1][iyt][ixt]);
                                          tbuf[0][iys][ixo][izo][iyt][ixt]=
					        (unsigned short)((t/InfTime)*InfByte2);

					  if (gd->crfp==NULL) continue;
					  c=zgamma0*(
                                              xalpha0*ctab[iys][ixs  ][izs][iyt][ixt]+
                                              xalpha *ctab[iys][ixs+1][izs][iyt][ixt])+
                                            zgamma *(
					      xalpha0*ctab[iys][ixs  ][izs+1][iyt][ixt]+
                                              xalpha *ctab[iys][ixs+1][izs+1][iyt][ixt]);
					  c=MIN(c,1.0);
					  cbuf[iys][ixo][izo][iyt][ixt]=
					        (unsigned char)(c*InfByte1);

					  r=zgamma0*(
                                              xalpha0*rtab[iys][ixs  ][izs][iyt][ixt]+
                                              xalpha *rtab[iys][ixs+1][izs][iyt][ixt])+
                                            zgamma *(
					      xalpha0*rtab[iys][ixs  ][izs+1][iyt][ixt]+
                                              xalpha *rtab[iys][ixs+1][izs+1][iyt][ixt]);
					  r=MIN(r,InfDistance);
					  rbuf[iys][ixo][izo][iyt][ixt]=
					      (unsigned char)((r/InfDistance)*InfByte1);
				    }
				    for (im=1;im<gd->multit;im++) {
                                        if (ttab[im][iys][ixs  ][izs  ][iyt][ixt]>=InfTime ||
				            ttab[im][iys][ixs+1][izs  ][iyt][ixt]>=InfTime ||
				            ttab[im][iys][ixs  ][izs+1][iyt][ixt]>=InfTime ||
				            ttab[im][iys][ixs+1][izs+1][iyt][ixt]>=InfTime) 
				              tbuf[im][iys][ixo][izo][iyt][ixt]=InfByte2;
				        else {
                                            t=zgamma0*(
                                            xalpha0*ttab[im][iys][ixs  ][izs][iyt][ixt]+
                                            xalpha *ttab[im][iys][ixs+1][izs][iyt][ixt])+
                                              zgamma *(
			                    xalpha0*ttab[im][iys][ixs  ][izs+1][iyt][ixt]+
                                            xalpha *ttab[im][iys][ixs+1][izs+1][iyt][ixt]);
                                            tbuf[im][iys][ixo][izo][iyt][ixt]=
			      	            (unsigned short)((t/InfTime)*InfByte2);
					  }
				    }
			      }
                        }
                  }
	    }
      }
}

void filt( float *trace, struct GD *gd, float *trf)
/**********************************************************************
filt - Low-pass filter, integration and phase shift for input data
***********************************************************************
Function prototype:
void filt( float *trace, struct GD *gd, float *trf);
***********************************************************************
Inputs:
float *trace single seismic trace
struct *gd   grid information

outputs:
float *trf  filtered and phase-shifted seismic trace
***********************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
***********************************************************************/
{
      static int nfft=0; 
      static int itaper, nw, nwf;
      static float *taper, *amp, *ampi, dw;

      int  it, iw, itemp;
      float temp, ftaper, const2, *rt;
      complex *ct;

      gd->fmax *= 2.0*PI;
      ftaper = 0.1*gd->fmax;
      const2 = sqrt(2.0);

      /***************************************************************
      this will only be done once, nice! All the static values are set
      ****************************************************************/
      if (nfft==0) {
            /* Set up FFT parameters */
            nfft = npfaro(gd->nt+gd->ntrimax, 2 * (gd->nt+gd->ntrimax));
            if (nfft >= SU_NFLTS || nfft >= 720720)
                  err("Padded nt=%d -- too big", nfft);

            nw = nfft/2 + 1;
            dw = 2.0*PI/(nfft*gd->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.0;
            }

            nwf = 0.5+gd->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*gd->dt));
            }
      }

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

      memcpy(rt, trace, gd->nt*FSIZE);
      memset((void *) (rt + gd->nt), (int) '\0', (nfft- gd->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(!gd->ls){
            for(iw=0; iw<=nwf+itaper; ++iw){
                  /* phase shifts PI/4    */
                  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;
            }
      } else {
            for(iw=0; iw<=nwf+itaper; ++iw){
                  ct[iw].i = ct[iw].i*amp[iw];
                  ct[iw].r = ct[iw].r*amp[iw];
            }
      }
      pfacr(-1, nfft, ct, rt);

      /* Load traces back in */
      for (it=0; it<gd->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<gd->ntrimax; ++it)  trf[it] = rt[nfft-gd->ntrimax+it];
      for (it=0; it<gd->nt + gd->ntrimax; ++it)  trf[it+gd->ntrimax] = rt[it];

      free1float(rt);
      free1complex(ct);
}