Operator definition for 3-D DMO

  Following is an example of implementing the operator domain definition function for 3-D DMO:

function dmo3d_domain (inv_map,inv_amp,inp_geom,out_geom, &
                       inp_cent_coord,out_coord) result (stat)
logical                 :: inv_map,inv_amp
integer                 :: stat
type (seis3d_geom_type), intent(in) :: inp_geom
type (seis3d_geom_type), intent(inout) :: out_geom

real, pointer,  dimension (:,:,:)  :: inp_cent_coord
real, pointer,  dimension (:,:,:)  :: out_coord
! no cross-line extension for operator
nh_z2=0

sin_inp=sin(inp_geom%irr%azim(1))
cos_inp=cos(inp_geom%irr%azim(1))
h_inp=.5*inp_geom%irr%aoff(1)
h_cos_inp=h_inp*cos_inp
h_sin_inp=h_inp*sin_inp
! Operator aperture limit according to min velocity
t_bar=(2.*h_inp*max_slow)/inp_geom%irr%d_t
! Estimate sampling rate along z1 and corresponding number of samples
d_z1 = sqrt((cos_inp*out_geom%reg%d_cmpx)**2  + &
            (sin_inp*out_geom%reg%d_cmpy)**2) / h_inp
nh_z1=(1./d_z1)-1

n_z1=1+2*nh_z1
n_z2=1+2*nh_z2
n_z1z2=n_z1*n_z2
n_dim=inp_geom%irr%n_dim
! allocate array for operator domain
allocate(inp_cent_coord(-nh_z1:nh_z1,-nh_z2:nh_z2,2), stat=ierr)
if(ierr .ne. 0) call seperr('point_domain: Troubles in allocating cordinate arrays')
! allocate arrays for output geometry info n
allocate( out_geom%irr%cmp(n_z1z2,inp_geom%irr%n_dim), stat=ierr)
if(ierr .ne. 0) call seperr('point_domain: Troubles in allocating output geometry arrays')
! loop over operator domain
do i_z2=-nh_z2,nh_z2
  do i_z1=-nh_z1,nh_z1
    z1 = i_z1*d_z1
! computer input_centered coordinates
    inp_cent_coord(i_z1,i_z2,1) = z1
    inp_cent_coord(i_z1,i_z2,2) = 0.
    i_z1z2=(i_z1+nh_z1)+(i_z2+nh_z2)*n_z1+1
! compute corresonding output space coordinates
    out_geom%irr%cmp(i_z1z2,1) = inp_geom%irr%cmp(1,1) + z1 * h_cos_inp
    out_geom%irr%cmp(i_z1z2,2) = inp_geom%irr%cmp(1,2) + z1 * h_sin_inp
  end do
end do

stat=0

return
end function dmo3d_domain

Following is an example of implementing the transformation function for 3-D DMO:

function dmo3d_map_amp (inv_map,inv_amp,inp_geom,out_geom, &
      first_samp, last_samp, map, amp, inp_cent_coord) result (stat)
logical                 :: inv_map,inv_amp
integer                 :: stat
type (seis3d_geom_type), intent(in) :: inp_geom
type (seis3d_geom_type), intent(in) :: out_geom
integer, dimension(:,:), intent(out) :: first_samp,last_samp
real,  dimension (:,:,:), intent(out)  :: map,amp
real, pointer,  optional,  dimension (:,:,:)  :: inp_cent_coord

logical in_bound
integer :: i_samp_out,i_tr_out,i_z1
real aoff_sq,t_out,t_inp,inv_d_t_inp,rat,z1,ratsq
real inv_d_t_out,t0_out
real ::  epsi,max_t_inp,t0_inp

if(.not. present(inp_cent_coord)) then
  call seperr('dmo3d_map_amp: inp_cent_coord should be present')
end if

epsi=epsilon(z1)

if(.not. inv_map) then

inv_d_t_inp = 1./inp_geom%irr%d_t
t0_inp = inv_d_t_inp * inp_geom%irr%t_0
! loop over operator domain
  do i_z1=lbound(inp_cent_coord,dim=1),ubound(inp_cent_coord,dim=1)

i_tr_out=i_z1-lbound(inp_cent_coord,dim=1)+1
! compute stretching factor
    z1=inp_cent_coord(i_z1,0,1)
    ratsq=1./(1. - z1*z1)
    rat=sqrt(ratsq)
! set first output sample
    first_samp(1,i_tr_out)=ceiling(t0_inp/rat)
! set input time limits according to aperture limitation
    if(abs(z1) > epsi) then
      max_t_inp=min(inp_geom%irr%n_t,(t_bar/abs(z1))-t0_inp)
    else
      max_t_inp=inp_geom%irr%n_t
    end if
! set values before looping
    last_samp(1,i_tr_out)=out_geom%irr%n_t
    i_samp_out=0
    in_bound = .true.
 !loop over output samples
    do while ((i_samp_out < out_geom%irr%n_t) .and. in_bound)
      i_samp_out=i_samp_out+1
      t_out=out_geom%irr%t_0 + (i_samp_out-1)*out_geom%irr%d_t
! compute time map
      t_inp=t_out*rat
      map(i_samp_out,1,i_tr_out)=(t_inp - inp_geom%irr%t_0)*inv_d_t_inp +1
! compute amplitude
      if(inv_amp) then
        amp(i_samp_out,1,i_tr_out)=1.
      else
        amp(i_samp_out,1,i_tr_out)=(2.*ratsq-1.)
      end if
! check to be within input trace limits
      if(map(i_samp_out,1,i_tr_out) > max_t_inp) then
        in_bound=.false.
        last_samp(1,i_tr_out)=i_samp_out-1
      end if
    end do
  end do
else

inv_d_t_inp = 1./inp_geom%irr%d_t
inv_d_t_out = 1./out_geom%irr%d_t
t0_out = inv_d_t_out * out_geom%irr%t_0
! loop over operator domain
  do i_z1=lbound(inp_cent_coord,dim=1),ubound(inp_cent_coord,dim=1)

i_tr_out=i_z1-lbound(inp_cent_coord,dim=1)+1
! compute stretching factor
    z1=inp_cent_coord(i_z1,0,1)
    ratsq=(1. - z1*z1)
    rat=sqrt(ratsq)
! set output output sample according to aperture limitation
    if(abs(z1) > epsi) then
      last_samp(1,i_tr_out)=min(out_geom%irr%n_t, (t_bar/abs(z1))-t0_out)
    else
      last_samp(1,i_tr_out)=out_geom%irr%n_t
    end if
! set values before looping
    first_samp(1,i_tr_out)=1
    i_samp_out=last_samp(1,i_tr_out)
    in_bound = .true.
 !loop over output samples
    do while ((i_samp_out > 1 ) .and. in_bound)
      i_samp_out=i_samp_out - 1
      t_out=out_geom%irr%t_0 + (i_samp_out-1)*out_geom%irr%d_t
! compute time map
      t_inp=t_out*rat
      map(i_samp_out,1,i_tr_out)=(t_inp - inp_geom%irr%t_0)*inv_d_t_inp +1
! compute amplitude
      if(inv_amp) then
        amp(i_samp_out,1,i_tr_out)=1.
      else
        amp(i_samp_out,1,i_tr_out)=(2.*ratsq-1.)
      end if
! check to be within input trace limits
      if(map(i_samp_out,1,i_tr_out) < 1.) then
        in_bound=.false.
        first_samp(1,i_tr_out)=i_samp_out + 1
      end if
    end do
  end do
end if

stat=0

return
end function dmo3d_map_amp