Operator definition for constant-velocity 3-D prestack migration

  Following is an example of implementing the transformation function for 3-D prestack migration with constant velocity:

function mig3d_const (want_info,inv_map,inv_amp,inp_geom,out_geom, &
      first_samp, last_samp, max_freq, map, amp, inp_cent_coord) result (info)
type(genkir3d_info_type) :: info
logical, intent(in)      :: want_info,inv_map,inv_amp
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,max_freq
real, pointer,  optional,  dimension (:,:,:)  :: inp_cent_coord

logical in_bound
integer :: n_t_mod,n_t_dat,i_samp_mod,n_tr_mod,n_tr_dat,i_tr_dat
integer :: n_dim,i_dim
real ::  t_mod,t_mod_sq,t_dat,inv_d_t_dat,inv_d_t_mod
real ::  t_0_mod,t_0_dat,d_t_mod,dpos,dneg
real :: dm(3)

if(.not. want_info) then
!{

n_dim=out_geom%irr%n_dim
n_tr_mod=out_geom%irr%n_tr
n_t_mod=out_geom%irr%n_t
t_0_mod=out_geom%irr%t_0 
d_t_mod=out_geom%irr%d_t
inv_d_t_mod=1./out_geom%irr%d_t

n_tr_dat=inp_geom%irr%n_tr
n_t_dat=inp_geom%irr%n_t
t_0_dat=inp_geom%irr%t_0
inv_d_t_dat=1./inp_geom%irr%d_t

max_freq=1.

if(.not. inv_map) then
!{
  do i_tr_dat=1,n_tr_dat
  !{

dpos=0.
    dneg=0.
    do i_dim=1,n_dim
      dm(i_dim)=(out_geom%irr%cmp(1,i_dim) - &
                 inp_geom%irr%cmp(i_tr_dat,i_dim))
      dpos = dpos + (dm(i_dim) + &
                     .5*inp_geom%irr%off(i_tr_dat,i_dim))**2
      dneg = dneg + (dm(i_dim) - &
                     .5*inp_geom%irr%off(i_tr_dat,i_dim))**2
    end do
    dpos=dpos*slow_sq
    dneg=dneg*slow_sq

in_bound = .true.
    first_samp(i_tr_dat,1)=1
    last_samp(i_tr_dat,1)=n_t_mod
    i_samp_mod = first_samp(i_tr_dat,1)
    do while ((i_samp_mod < n_t_mod ) .and. in_bound)
      t_mod=.5*(t_0_mod + (i_samp_mod-1)*d_t_mod)
      t_mod_sq=t_mod**2
      t_dat=sqrt(t_mod_sq + dpos) + sqrt(t_mod_sq + dneg) 
      map(i_samp_mod,i_tr_dat,1)=(t_dat-t_0_dat)*inv_d_t_dat + 1

if(inv_amp) then
        if(t_mod_sq > 0.) then
          amp(i_samp_mod,i_tr_dat,1)=t_dat/t_mod
        else
          amp(i_samp_mod,i_tr_dat,1)=1.
        end if
      else
        if(t_dat > 0.) then
          amp(i_samp_mod,i_tr_dat,1)=t_mod/t_dat
        else
          amp(i_samp_mod,i_tr_dat,1)=1.
        end if
      end if

if(map(i_samp_mod,i_tr_dat,1) > n_t_dat) then
        in_bound=.false.
        last_samp(i_tr_dat,1)=i_samp_mod-1
      end if
      i_samp_mod=i_samp_mod + 1
    end do
  !}
  end do
!}
else
!{
  call seperr('mig3d_const: cannot implement inv_map with mig3d')
!}
end if

info=info_default
  info%status=0
!}
else
!{
 info=info_default
!}
end if

return
end function mig3d_const

.