sudatumfd.c

su 的源代码库

   {
      for(ix=0;ix<nx;ix++)
      {
	 tr.ns = nt;
	 tr.dt = dt*1000000;
	 tr.gx = gx[ix][isx];
	 tr.sx = sx[ix][isx];
	 tr.cdp = cdp[ix][isx];
	 tr.offset = offset[ix][isx];
	 memcpy( (void *) tr.data, (void *) outtime[ix][isx], nt*4);
	 puttr(&tr);
      }
   }
   free3float(outtime);
   /*****************************************************************/

   free2int(gx);
   free2int(sx);
   free2int(cdp);
   free2int(offset);

   return(CWP_Exit());
}




/**************************thin lens subroutine*****************************/
void tlens(complex ***val_ifr_pad,complex ***val_ifr1,int val_nsx,int val_nx_pad,
	   int val_nw,float val_dw,float val_dz,float val_v,float val_dx,
	   float **val_vel,int val_zlev,
	   int **val_gx,int val_buff,float val_x_0)
{
  int i,j,k,a;
  complex mult;
  float ig;
  a = 0;
  for(i=0;i<val_nsx;i++)
  {
    for(j=0;j<val_nx_pad;j++)
    {
      if(j<val_buff)
      {
	a = val_gx[0][i];
      }
      if(j>val_nx_pad-1-val_buff)
      {
	a = val_gx[val_nx_pad-1-2*val_buff][i];
      }
      if(j>=val_buff && j<=val_nx_pad-1-val_buff)
      {
	a = val_gx[j-val_buff][i];
      }
      val_v = 0.5*val_vel[(int)((a-val_x_0)/val_dx)][val_zlev];
      ig = (-1.0/val_v)*val_dz;
      for(k=0;k<val_nw;k++)
      {
	 mult = cmplx(cos(k*val_dw*ig),sin(k*val_dw*ig));
	 val_ifr1[j][i][k] = cmul(val_ifr_pad[j][i][k],mult);
      }
    }
  }
}
/*************************************************************************/



/*************************diffraction term subroutine********************/
void diff(complex ***val_ifr1,complex ***val_ifr2,int val_nsx,int val_nx_pad,
	  int val_nw,float val_dw,float val_dz,float val_v,float val_dx,
	  float **val_vel,int val_zlev,
	  int **val_gx,int val_buff,float val_x_0)
{
  int i,j,k,m,n,a;
  float w;
  complex *e,*f,*d,*ohm,*alpha,*beta,*gamma;
  complex num1,num2,num3,num4,num5;
  complex num6,num7,num8,num9,num10,num11;
  float sigma,tau;

  sigma = 0.5;
  tau = 0.25;

  e = alloc1complex(val_nx_pad);
  f = alloc1complex(val_nx_pad);
  d = alloc1complex(val_nx_pad);
  ohm = alloc1complex(val_nx_pad);
  alpha = alloc1complex(val_nx_pad);
  beta = alloc1complex(val_nx_pad);
  gamma = alloc1complex(val_nx_pad);

  a = 0;
  for(i=0;i<val_nsx;i++)
  {
    for(k=0;k<val_nw;k++)
    {
      w = k*val_dw;
      if(w<5.0) w=5.0;
      for(n=0;n<val_nx_pad;n++)
      {
	if(n<val_buff)
	{
	  a = val_gx[0][i];
	}
	if(n>val_nx_pad-1-val_buff)
	{
	  a = val_gx[val_nx_pad-1-2*val_buff][i];
	}
	if(n>=val_buff && n<=val_nx_pad-1-val_buff)
	{
	  a = val_gx[n-val_buff][i];
	}
	val_v = 0.5*val_vel[(int)((a-val_x_0)/val_dx)][val_zlev];
	alpha[n] = cmplx(0.0,-1.0*val_v*sigma*val_dz/(val_dx*val_dx*2.0*w));
	beta[n] = cmul(alpha[n],cadd(cmplx(1.0,0.0),
		  cmplx(0.0,-2.0*val_v*tau/(w*sigma*val_dz))));
	gamma[n] = cmul(alpha[n],cadd(cmplx(1.0,0.0),
		  cmplx(0.0,2.0*val_v*tau/(w*sigma*val_dz))));
      }
      /**********************Calculate d[n]'s**********************/
      num6 = cmul(cadd(cmplx(1.0,0.0),cmul(cmplx(-2.0,0.0),gamma[0]))
		  ,val_ifr1[0][i][k]);
      num7 = cmul(gamma[0],val_ifr1[1][i][k]);
      d[0] = cadd(num6,num7);
      ohm[0] = cmul(cmplx(-1.0,0.0),cdiv(d[0],beta[0])); 
      for(n=1;n<(val_nx_pad-1);n++)
      {
	num3 = cmul(gamma[n],val_ifr1[n-1][i][k]);
	num4 = cmul(cadd(cmplx(1.0,0.0),cmul(cmplx(-2.0,0.0),gamma[n]))
		    ,val_ifr1[n][i][k]);
	num5 = cmul(gamma[n],val_ifr1[n+1][i][k]);

	d[n] = cadd(num3,cadd(num4,num5));
	ohm[n] = cmul(cmplx(-1.0,0.0),cdiv(d[n],beta[n]));
      }
      num8 = cmul(gamma[val_nx_pad-1],val_ifr1[val_nx_pad-2][i][k]);
      num9 = cmul(cadd(cmplx(1.0,0.0),cmul(cmplx(-2.0,0.0),gamma[val_nx_pad-1]))
		  ,val_ifr1[val_nx_pad-1][i][k]);
      d[val_nx_pad-1] = cadd(num8,num9);
      ohm[val_nx_pad-1] = cmul(cmplx(-1.0,0.0),cdiv(d[0],beta[0]));
      /***********************************************************/ 
      /*******************Calculate e's and f's*******************/
      e[1] = cadd(cdiv(cmplx(1.0,0.0),beta[1]),cmplx(2.0,0.0));
      f[1] = ohm[1];
      for(m=2;m<val_nx_pad-1;m++)
      {
	 num1 = cadd(cdiv(cmplx(1.0,0.0),beta[m]),cmplx(2.0,0.0));
	 num2 = cdiv(cmplx(-1.0,0.0),e[m-1]);
	 e[m] = cadd(num1,num2);
	 f[m] = cadd(ohm[m],cdiv(f[m-1],e[m-1]));
      }
      e[0] = cmplx(0.0,0.0);
      e[val_nx_pad-1] = cmplx(0.0,0.0);
      /**********************************************************/
      /*******Downward continue one freq. on one shot gather*****/
      val_ifr2[val_nx_pad-1][i][k] = cmplx(0.0,0.0);
      num10 = cdiv(f[val_nx_pad-2],e[val_nx_pad-2]);
      val_ifr2[val_nx_pad-2][i][k] = cmul(cmplx(-1.0,0.0),num10);
      for(j=val_nx_pad-3;j>0;j--)
      {
	 num11 = cadd(val_ifr2[j+1][i][k],cmul(cmplx(-1.0,0.0),f[j]));
	 val_ifr2[j][i][k] = cdiv(num11,e[j]);
      }
      val_ifr2[0][i][k] = cmplx(0.0,0.0);
      /**********************************************************/
    }/*end freq loop*/
  }/*end source loop*/

  free1complex(e);
  free1complex(f);
  free1complex(d);
  free1complex(ohm);
  free1complex(alpha);
  free1complex(beta);
  free1complex(gamma);
}
/***********************************************************************/





/*******************Padding subroutine********************************/
void padsub(complex ***val_ifr, complex ***val_ifr_pad,int val_nsx,
	    int val_nx,int val_nw,int val_nx_pad,int val_buff)
{
  int i,j,k;

  for(i=0;i<val_nsx;i++)
  {
    for(j=0;j<val_nx_pad;j++)
    {
      for(k=0;k<val_nw;k++)
      {
	if(j<val_nx+val_buff && j>-1+val_buff )
	  val_ifr_pad[j][i][k] = val_ifr[j-val_buff][i][k];
	else 
	{
	  val_ifr_pad[j][i][k] = cmplx(0.0,0.0); 
	}
      }
    }
  }
}
/*********************************************************************/



/*********************Unpadding routine*******************************/
void unpadsub(complex ***val_ifr, complex ***val_ifr_upad,int val_nsx,
	      int val_nx,int val_nw,int val_nx_pad,int val_buff)
{
  int i,j,k;

  for(i=0;i<val_nsx;i++)
  {
    for(j=0;j<val_nx;j++)
    {
      for(k=0;k<val_nw;k++)
      {
	val_ifr_upad[j][i][k] = val_ifr[j+val_buff][i][k];
      }
    }
  }
}
/********************************************************************/



/****************Tapering subroutine**********************************/
void taper(complex ***val_ifr_pad,complex ***val_ifr_pt,int val_nsx,
	   int val_nx_pad,int val_nw,float val_dw,int val_tap_len,
	   int val_buff,int val_nx)
{
  int i,j,k,m;
  complex *tap_func;
  float pi=3.1415927;

  tap_func = alloc1complex(val_nx_pad);

  if(val_tap_len==0)
  {
    for(m=0;m<val_nx_pad;m++) tap_func[m]=cmplx(1.0,0.0);
  }
  if(val_tap_len>0)
  {
    for(m=0;m<val_nx_pad;m++)
    {
      if(m<val_buff) tap_func[m]=cmplx(0.0,0.0);
      if(m>=val_buff && m<(val_buff+val_tap_len-1))
      {
	tap_func[m]=cmplx(sin(2.0*pi*(m+1-val_buff)/(4.0*val_tap_len)),0.0);
      }
      if(m>=(val_buff+val_tap_len-1) && m<=(val_buff+val_nx-1-(val_tap_len-1)))
      {
	tap_func[m]=cmplx(1.0,0.0);
      }
      if(m>(val_buff+val_nx-1-(val_tap_len-1)) && m<=(val_buff+val_nx-1))
      {
	tap_func[m]=cmplx(cos(2.0*pi*(m-(val_buff+val_nx-1-(val_tap_len-1)))/
			  (4.0*val_tap_len)),0.0);
      }
      if(m>(val_buff+val_nx-1)) tap_func[m]=cmplx(0.0,0.0);
    }
  }

  for(i=0;i<val_nsx;i++)
  {
    for(j=0;j<val_nx_pad;j++)
    {
      for(k=0;k<val_nw;k++)
      {
	 val_ifr_pt[j][i][k] = cmul(tap_func[j],val_ifr_pad[j][i][k]);
      }
    }
  }

  free1complex(tap_func);
}
/********************************************************************/