kdmig3d.c

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

			      tro.tracr=1+ixo;
			      tro.tracl=1+ixoffset+nxoffset*ixo;
			      tro.cdp=gd->fxo+ixo*gd->dxo;
			      tro.cdpt=1+ixoffset;

			      /*write out*/
			      fputtr(outfp,&tro);
			}
		  }
	    }
      }

      fprintf(gd->jpfp,"\n");
      fprintf(gd->jpfp,"Output done\n");
      fflush(gd->jpfp);

      efclose(gd->jpfp);
      efclose(outfp);

      free6ushort(tbuf);
      free5uchar(cbuf);
      free5uchar(rbuf);
      free6float(ttab);
      free5float(ctab);
      free5float(rtab);
      free5float(mig);
      return(CWP_Exit());
}

void mig3d( struct GD *gd, float *trace, float *trf, float ***mig,
      unsigned short ******tbuf, unsigned char  *****cbuf,
      unsigned char  *****rbuf)
/******************************************************************
mig3d - migrate a single trace
*******************************************************************
Function prototype:
void mig3d( struct GD *gd, float *trace, float *trf, float ***mig,
      unsigned short ******tbuf, unsigned char  *****cbuf,
      unsigned char  *****rbuf);
*******************************************************************
Inputs:
struct *gd   grid information
float *trace trace data
float *trf   integrated trace
unsigned short *******tbuf traveltime buffer
unsigned char  *****cbuf  cosine buffer
unsigned char  *****rbuf  path buffer
Outputs:
float ***mig migrated section

*******************************************************************
Notes:

For formulations:
		    |    ps  + pg    |
	h(y,xi)=det | d(ps+pg)/d xi1 |
		    | d(ps+pg)/d xi2 |

	       = 2 cos(theta)^2 z (rs+rg)(rs^2+rg^2)/(c^3 rs^3 rg^3);

	a(y,xi)=A(x,xs)A(x,xg)
	       =1/(16 pi^2 rs rg);

	|\Delta_y phi|= 2cos(theta)/c(y);

	here y is (y_1,y_2,y_3) meaning the outputing points.

	The correct correction factor is:

	| h |/a /|\Delta_y phi| =
	16 pi^2 cos(theta) y_3 (r_s+rg)(rs^2+rg^2)/(v^2 rs^2 rg^2); 

*******************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
*******************************************************************/

{
      int nyof;	        /*starting index for y scope: 0<=nyof<nyo*/
      int nyoe;	        /*ending index for y scope: 0<=nyoe<nyo*/
      int nxof;	        /*starting index for x scope: 0<=nxof<nx */
      int nxoe;	        /*ending index for x scope: 0<=nxoe<nx */
      int iyo;	        /*index for nyof and nyoe*/
      int ixo;	        /*index for nxof and nxoe*/
      float ym;	        /*y midpoint position*/
      float xm;	        /*x midpoint position*/
      float yo; 	/*y position */
      float xo;	        /*x position */

      float xs,ys,ybeta,ybeta0,xalpha,xalpha0;
      int ixs,iys,izo;
      float xt,yt;
      float xs0a,xs0a0,xg0a,xg0a0,xs1a,xs1a0,xg1a,xg1a0;
      float ys0b,ys0b0,yg0b,yg0b0,ys1b,ys1b0,yg1b,yg1b0;
      float tiys0,tiys1,ttotal=0.0;
      float tiyg0,tiyg1;
      float ciys0,ciys1;
      float riys0,riys1,riyg0,riyg1,rtotal;
      float cosine;
      int itintr;
      float tintr,wtintr,wtintr0;
      float amp=1.0;
      int ixts0,ixtg0,ixts1,ixtg1;
      int iyts0,iytg0,iyts1,iytg1;
      float rs,rg;
      float zo;
      int im;

      /*for the meaning of these varibles, refer to Claerbout's
      triabgle filter notation*/
      float **rhos,**rhog;
      int ntri;    /*triangle filter length*/
      float drho,ddxg,ddyg,ddxs,ddys;
      float temp,ts=0.0,tg=0.0;

      rhos=ealloc2float(gd->nxo,gd->nyo);
      rhog=ealloc2float(gd->nxo,gd->nyo);

      /**************************************************************
      Geometry:

		fz _______________________
		   |			 |
		   |fzs_______________   |
		   |   |	     |   |
		   |   |	     |   |
		   |   |   Ttimes    |   |
		   |   |	     |   |	
		   |   |_____________|   |
		   |			 |
		   |			 |
		   -----------------------
	
      and the y dimension is similar to the x dimension
      **************************************************************/

      ym=0.5*(gd->sy+gd->gy);		/*mid point in y*/
      xm=0.5*(gd->sx+gd->gx);		/*mid point in x*/

      #ifdef DEBUG
      fprintf(gd->jpfp,"xm,ym=%g %g\n",xm,ym);
      #endif

      /*******************************************************
      For all ixo from nxof to nxoe and iyo from nyof to nyoe:
      ON OUTPUT GRIDS
      ********************************************************/
      nyof=(ym-gd->yaper-gd->fyo)/gd->dyo+0.5;
      if (nyof<0) nyof=0;
      nyoe=(ym+gd->yaper-gd->fyo)/gd->dyo+0.5;
      if (nyoe>=gd->nyo) nyoe=gd->nyo-1;

      #ifdef DEBUG
      fprintf(gd->jpfp,"This trace: fyo,dyo,nyof,nyoe=%g %g %d %d\n",
      gd->fyo,gd->dyo,nyof,nyoe);
      #endif

      nxof=(xm-gd->xaper-gd->fxo)/gd->dxo+0.5;
      if (nxof<0) nxof=0;
      nxoe=(xm+gd->xaper-gd->fxo)/gd->dxo+0.5;
      if (nxoe>=gd->nxo) nxoe=gd->nxo-1;

      #ifdef DEBUG
      fprintf(gd->jpfp,"This trace: fxo,dxo,nxof,nxoe=%g %g %d %d\n",
      gd->fxo,gd->dxo,nxof,nxoe);
      #endif

      for (iyo=nyof;iyo<=nyoe;iyo++){	
	    yo=gd->fyo+iyo*gd->dyo;
    	    ys=(yo-gd->fys)/gd->dys;
    	    iys=(int)ys;
	    iys=MIN(iys,gd->nys-2);
	    iys=MAX(0,iys);
	    ybeta=ys-iys;
            ybeta0=1.0-ybeta;

            yt=(gd->sy-gd->fygd-ybeta*gd->dys)/gd->dygd;
            iyts0=(int)yt;
            iyts0=MIN(iyts0,gd->nygd-2);
            iyts0=MAX(iyts0,0);
            ys0b=yt-iyts0;
            ys0b0=1.0-ys0b;
	    iyts0-=(gd->iysf+gd->iysr*iys-gd->nyt/2);
	    iyts0=MIN(iyts0,gd->nyt-2);
            iyts0=MAX(iyts0,0);

            yt=(gd->gy-gd->fygd-ybeta*gd->dys)/gd->dygd;
            iytg0=(int)yt;
            iytg0=MIN(iytg0,gd->nygd-2);
            iytg0=MAX(iytg0,0);
            yg0b=yt-iytg0;
            yg0b0=1.0-yg0b;
	    iytg0-=(gd->iysf+gd->iysr*iys-gd->nyt/2);
	    iytg0=MIN(iytg0,gd->nyt-2);
            iytg0=MAX(iytg0,0);

            yt=(gd->sy-gd->fygd+ybeta0*gd->dys)/gd->dygd;
            iyts1=(int)yt;
            iyts1=MIN(iyts1,gd->nygd-2);
            iyts1=MAX(iyts1,0);
            ys1b=yt-iyts1;
            ys1b0=1.0-ys1b;
	    iyts1-=(gd->iysf+gd->iysr*(iys+1)-gd->nyt/2);
	    iyts1=MIN(iyts1,gd->nyt-2);
            iyts1=MAX(iyts1,0);

            yt=(gd->gy-gd->fygd+ybeta0*gd->dys)/gd->dygd;
            iytg1=(int)yt;
            iytg1=MIN(iytg1,gd->nygd-2);
            iytg1=MAX(iytg1,0);
            yg1b=yt-iytg1;
            yg1b0=1.0-yg1b;
	    iytg1-=(gd->iysf+gd->iysr*(iys+1)-gd->nyt/2);
	    iytg1=MIN(iytg1,gd->nyt-2);
            iytg1=MAX(iytg1,0);

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

                  /*to understand this following paragraph, see
                   Claerbout's notation*/

                  rhos[iyo][ixo]=sqrt((xo-gd->sx)*(xo-gd->sx)+(yo-gd->sy)*(yo-gd->sy));
                  rhog[iyo][ixo]=sqrt((xo-gd->gx)*(xo-gd->gx)+(yo-gd->gy)*(yo-gd->gy));

                  ddxs=fabs(xo-gd->sx)*gd->dxo/rhos[iyo][ixo];
                  ddys=fabs(yo-gd->sy)*gd->dyo/rhog[iyo][ixo];

                  ddxg=fabs(xo-gd->gx)*gd->dxo/rhos[iyo][ixo];
                  ddyg=fabs(yo-gd->gy)*gd->dyo/rhog[iyo][ixo];

                  drho=sqrt(0.5*(ddxs*ddxs+ddys*ddys+ddxg*ddxg+ddyg*ddyg));
                  temp=4.0/gd->v0/gd->v0*drho/gd->dt;
                  rhos[iyo][ixo]*=temp;
                  rhog[iyo][ixo]*=temp;

                  xt=(gd->sx-gd->fxgd-xalpha*gd->dxs)/gd->dxgd;
                  ixts0=(int)xt;
	      	  ixts0=MIN(ixts0,gd->nxgd-2);
                  ixts0=MAX(ixts0,0);
                  xs0a=xt-ixts0;
                  xs0a0=1.0-xs0a;
                  ixts0-=(gd->ixsf+gd->ixsr*ixs-gd->nxt/2);
		  ixts0=MIN(ixts0,gd->nxt-2);
                  ixts0=MAX(ixts0,0);

                  xt=(gd->gx-gd->fxgd-xalpha*gd->dxs)/gd->dxgd;
                  ixtg0=(int)xt;
		  ixtg0=MIN(ixtg0,gd->nxgd-2);
                  ixtg0=MAX(ixtg0,0);
                  xg0a=xt-ixtg0;
                  xg0a0=1.0-xg0a;
                  ixtg0-=(gd->ixsf+gd->ixsr*ixs-gd->nxt/2);
		  ixtg0=MIN(ixtg0,gd->nxt-2);
                  ixtg0=MAX(ixtg0,0);

                  xt=(gd->sx-gd->fxgd+xalpha0*gd->dxs)/gd->dxgd;
                  ixts1=(int)xt;
		  ixts1=MIN(ixts1,gd->nxgd-2);
                  ixts1=MAX(ixts1,0);
                  xs1a=xt-ixts1;
                  xs1a0=1.0-xs1a;
                  ixts1-=(gd->ixsf+gd->ixsr*(ixs+1)-gd->nxt/2);
		  ixts1=MIN(ixts1,gd->nxt-2);
                  ixts1=MAX(ixts1,0);

                  xt=(gd->gx-gd->fxgd+xalpha0*gd->dxs)/gd->dxgd;
                  ixtg1=(int)xt;
		  ixtg1=MIN(ixtg1,gd->nxgd-2);
                  ixtg1=MAX(ixtg1,0);
                  xg1a=xt-ixtg1;
                  xg1a0=1.0-xg1a;
                  ixtg1-=(gd->ixsf+gd->ixsr*(ixs+1)-gd->nxt/2);
		  ixtg1=MIN(ixtg1,gd->nxt-2);
                  ixtg1=MAX(ixtg1,0);

		  /********************************
		  Apply anti-aliasing operator here
		  using scaling functions
		  ********************************/
                  for(izo=0;izo<gd->nzo;izo++) {
                        ntri=MAX(1,(rhos[iyo][ixo]/ts+rhog[iyo][ixo]/tg)-1);

                        zo=gd->fzo+izo*gd->dzo;

                        if (gd->crfp!=NULL) {
                              /*cosine table at [iys][ixo][izo][nyt][nxt]*/
                              ciys0=yg0b0*(xg0a0*cbuf[iys][ixo][izo][iytg0][ixtg0]+
                                         xg0a *cbuf[iys][ixo][izo][iytg0][ixtg0+1])+
			            yg0b *(xg0a0*cbuf[iys][ixo][izo][iytg0+1][ixtg0]+
                                         xg0a *cbuf[iys][ixo][izo][iytg0+1][ixtg0+1]);

                              /*cos table at [iys+1][ixo][izo][nyt][nxt]*/
                              ciys1=yg1b0*(xg1a0*cbuf[iys+1][ixo][izo][iytg1][ixtg1]+
                                          xg1a *cbuf[iys+1][ixo][izo][iytg1][ixtg1+1])+
			             yg1b *(xg1a0*cbuf[iys+1][ixo][izo][iytg1+1][ixtg1]+
                                          xg1a *cbuf[iys+1][ixo][izo][iytg1+1][ixtg1+1]);

			      cosine=(ybeta0*ciys0+ybeta*ciys1)/(float)InfByte1;

                              /*r table at [iys][ixo][izo][nyt][nxt]*/
                              riys0=ys0b0*(xs0a0*rbuf[iys][ixo][izo][iyts0][ixts0]+
                                         xs0a *rbuf[iys][ixo][izo][iyts0][ixts0+1])+
			            ys0b *(xs0a0*rbuf[iys][ixo][izo][iyts0+1][ixts0]+
                                         xs0a *rbuf[iys][ixo][izo][iyts0+1][ixts0+1]);

                              riyg0=yg0b0*(xg0a0*rbuf[iys][ixo][izo][iytg0][ixtg0]+
                                         xg0a *rbuf[iys][ixo][izo][iytg0][ixtg0+1])+
			            yg0b *(xg0a0*rbuf[iys][ixo][izo][iytg0+1][ixtg0]+
                                         xg0a *rbuf[iys][ixo][izo][iytg0+1][ixtg0+1]);

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

                              riyg1=yg1b0*(xg1a0*rbuf[iys+1][ixo][izo][iytg1][ixtg1]+
                                          xg1a *rbuf[iys+1][ixo][izo][iytg1][ixtg1+1])+
			             yg1b *(xg1a0*rbuf[iys+1][ixo][izo][iytg1+1][ixtg1]+
                                          xg1a *rbuf[iys+1][ixo][izo][iytg1+1][ixtg1+1]);
                              
			      rs=(ybeta0*riys0+ybeta*riys1)/
				     (float)InfByte1*InfDistance;

			      rg=(ybeta0*riyg0+ybeta*riyg1)/
				     (float)InfByte1*InfDistance;

                              rtotal=rs+rg;