# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!reflector_mod
!
!Usage:
!call reflector_mod (arg1, arg2 arg3)
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Mon Jan 15 18:06:56 PST 2001
!
!Purpose:
!
!
module reflector_mod
  use tomo_space_conv_mod
  use tau_z_conv_3d_mod
  implicit none
  type(sep3d) :: refs
  logical :: ref_is_3d
  complex, pointer,private,save :: ref_2d(:,:)
  real,pointer,private,save :: ref_3d(:,:,:)
  contains
  logical function return_3d_conv_refs(myreflect,convert)
    real :: myreflect(:,:,:)
    logical :: convert_me,tau
    logical, optional :: convert
    type(sep3d) :: vel_conv
    return_3d_conv_refs=.true.
    if (present(convert)) then
      convert_me=convert
    else
      convert_me=.true.
    end if
    if (.not.ref_is_3d) then
      write(0,*) "reflector is 2-D requested 3-D reflector"
      return_3d_conv_refs=.false.
      return
    end if
    myreflect=ref_3d
    if (.not. vel_tomo_structure(vel_conv)) then
      write(0,*) "trouble grabbing converted vel"
      return_3d_conv_refs=.false.
      return
    end if
    call from_param("tau",tau,.false.)
    if (tau .and. convert_me) then
      if (.not. convert_ref_z_to_tau_3d(myreflect, refs%o(1:2),refs%d&
        &(1:2),-1.)) then
        write(0,*) "trouble converting reflector from z to tau"
        return
      end if
    end if
  end function
  logical function return_2d_conv_refs(myreflect,convert)
    complex, allocatable ::tempc(:,:)
    complex :: myreflect(:,:)
    real, pointer :: velt(:,:),mapt(:,:)
    logical :: tau,convert_me
    logical, optional :: convert
    type(sep3d) :: vel_conv
    return_2d_conv_refs=.true.
    if (present(convert)) then
      convert_me=convert
    else
      convert_me=.true.
    end if
    if (ref_is_3d) then
      write(0,*) "reflector is 3-D requested 2-D reflector"
      return_2d_conv_refs=.false.
    else
      myreflect=ref_2d
      if (.not. vel_tomo_structure(vel_conv)) then
        write(0,*) "trouble grabbing converted vel"
        return_2d_conv_refs=.false.
      else
        call from_param("tau",tau,.false.)
        if (tau .and. convert_me) then
          if (.not. convert_ref_z_to_tau_2d(myreflect)) then
            write(0,*) "trouble converting reflector from z to tau"
            return
          end if
        end if
      end if
    end if
  end function
  logical function init_reflector()
    integer, external :: auxin,getch,sreed
    type(sep3d) :: vel
    real, allocatable :: tempr(:,:,:)
    complex,allocatable :: tempc(:,:)
    character(len=1024) :: temp_ch
    real :: j(8)
    logical ::verb
    if (.not. vel_tomo_structure(vel)) then
      write(0,*) "invalid velocity structure"
      return
    end if
    init_reflector=.false.
    if (-1.eq.auxin("reflectors")) then
      call auxclose("reflectors")
      call auxinout("reflectors")
      call from_param("verb",verb,.false.)
      if (verb) then
        write(0,*) "    Assuming no reflector information"
      end if
      if (sep3d_ndims(vel).eq.2) then
        call init_sep3d(refs,"OUTPUT","COMPLEX","REGULAR",ctag="reflectors"&
          &)
        call axis_allocate(refs,2)
        refs%n(:)=(/2,2/)
      else
        call init_sep3d(refs,"OUTPUT","FLOAT","REGULAR",ctag="reflectors"&
          &)
        call axis_allocate(refs,3)
        refs%n(1:3)=(/2,2,2/)
      end if
      refs%ntraces=1
      call sep3d_set_sep3d(refs)
      call sep3d_set_write_status(refs,.true.,.true.)
      call sep3d_write_description("reflectors",refs)
      j=0
      call srite("reflectors",j,32)
      call sep3d_delete("reflectors")
      call auxclose("reflectors")
    end if
    call init_sep3d("reflectors",refs,"INPUT")
    if (.not. valid_structure(refs)) then
      write(0,*) "invalid reflector init"
      return
    end if
    if (sep3d_ndims(vel).eq.2) then
      if (sep3d_ndims(refs).ne.2 .and. refs%data_format(1:3).ne."COM")&
        & then
        write(0,*) "velocity is 2d reflector file must be 2-D complex"
        write(0,*) sep3d_ndims(refs),"=refs%ndims refs%data_format="&
          &,refs%data_format
        return
      end if
      ref_is_3d=.false.
      allocate(tempc(refs%n(1),refs%n(2)))
      allocate(ref_2d(refs%n(1),refs%n(2)))
      if (8*product(refs%n).ne. sreed("reflectors",tempc,size(tempc)*8&
        &)) then
        write(0,*) "trouble reading in reflectors"
        return
      end if
      ref_2d=tempc
      deallocate(tempc)
    else if (sep3d_ndims(refs).ne.3) then
      write(0,*) "    only understand velocity of dim 2 or 3"
      return
    else
      if (sep3d_ndims(refs).ne.3 .and. refs%data_format(1:3).ne."REA")&
        & then
        write(0,*) "velocity is 3d reflector file must be 3-D REAL"
        write(0,*) sep3d_ndims(refs),"=refs%ndims refs%data_format="&
          &,refs%data_format
        return
      end if
      ref_is_3d=.true.
      allocate(tempr(refs%n(1),refs%n(2),refs%n(3)))
      allocate(ref_3d(refs%n(1),refs%n(2),refs%n(3)))
      if (4*product(refs%n).ne. sreed("reflectors",tempr,size(tempc)*4&
        &)) then
        write(0,*) "trouble reading in reflectors"
        init_reflector=.false.
      end if
      ref_3d=tempr
      deallocate(tempr)
    end if
    init_reflector=.true.
  end function
  logical function return_ref(myref)
    type(sep3d) :: myref
    integer, external :: auxin
    return_ref=.true.
    myref=refs
    return_ref=valid_structure(refs)
  end function
  logical function ref_complex_real(ref_complex,ref_real,n,o,d,default&
    &)
    integer :: n,i,iref,icrp
    integer,allocatable :: ibeg(:)
    real :: o,d,ref_real(:,:,:),f,crp,default
    complex :: ref_complex(:,:)
    ref_complex_real=.false.
    if (size(ref_complex,2).ne.size(ref_real,3)) then
      write(0,*) "nrefs of ref_complex and ref_real not equal"
      return
    end if
    allocate(ibeg(size(ref_complex,2)))
    ibeg=1
    ref_real=default
    do icrp=1,size(ref_real,1)
      crp=o+d*(icrp-1)
      ibeg=1
      do iref=1,size(ref_real,3)
        if (crp< real(ref_complex(1,iref)) .or. crp > real(ref_complex&
          &(size(ref_complex,1),iref))) then
          cycle
        end if
        do while (ibeg(iref) < size(ref_complex,1)-1 .and. real&
          &(ref_complex(ibeg(iref)+1,iref))< crp)
          ibeg(iref)= ibeg(iref)+1
        end do
        f=(crp-real(ref_complex(ibeg(iref),iref)))/ (real(ref_complex&
          &(ibeg(iref)+1,iref))-real(ref_complex(ibeg(iref),iref)))
        ref_real(icrp,1,iref)= (1.-f)*aimag(ref_complex(ibeg(iref),iref&
          &))+f*aimag(ref_complex(ibeg(iref),iref))
      end do
    end do
    deallocate(ibeg)
    ref_complex_real=.true.
  end function
  subroutine smooth_ref_3d(d,ref_in,ref_out,npts)
    real :: ref_in(:,:),ref_out(:,:),d(2)
    integer :: npts,iref,nbeg,ndist
    integer :: i1,i2,dist,val,id,nfound,ib1,ib2,ie1,ie2,ib,ia,tot_dist
    logical :: found
    nbeg=sqrt(1.*npts)/2+1
    do i2=1,size(ref_in,2)
      do i1=1,size(ref_in,1)
        ndist=nbeg
        found=.false.
        do while (.not. found .and. ndist < max(size(ref_in,1),size&
          &(ref_in,2)))
          ib2=max(1,i2-ndist)
          ib1=max(1,i1-ndist)
          ie2=min(size(ref_in,2),i2+ndist)
          ie1=min(size(ref_in,1),i1+ndist)
          nfound=0
          val=0
          tot_dist=0.
          do ib=ib2,ie2
            do ia=ib1,ie1
              if (ref_in(ia,ib) > 0.) then
                nfound= nfound+1
                if (ia.eq.i1 .and. ib.eq.i2) then
                  dist=minval(d)/4.
                else
                  dist=sqrt((d(1)*(ia-i1))**2+(d(2)*(ib-i2))**2)
                end if
                val= val+ref_in(ia,ib)/dist
                tot_dist= tot_dist+1./dist
              end if
            end do
          end do
          if (nfound > npts) then
            found=.true.
          end if
        end do
        ref_out(i1,i2)=val/tot_dist
      end do
    end do
  end subroutine
end module
