# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!Thu Jan 9 16:28:30 PST 2003
!
!Purpose:
!
! Covariance filter that describes the inverse model covariance
!
!
!
module lap_filter_mod
  use param_check_mod
  use wilson
  use nhelix
  use helix
  use io_nfilt
  use tomo_space_conv_mod
  use lapfac
  implicit none
  integer, private, save :: npts_lap
  real, private, save :: scale_lap
  contains
  logical function calc_lap_2pass()
    type(nfilter) :: myfilt
    type(filter) :: filt,filt_out
    type(sep3d) :: tomo_s
    logical :: is_2d
    integer :: i,n123
    calc_lap_2pass=.false.
    if (.not. vel_tomo_structure(tomo_s)) then
      write(0,*) "could not grab converted velocity structure"
      return
    end if
    if (sep3d_ndims(tomo_s).eq.2) then
      write(0,*) "calling 2 pass laplacian doesn't make sense for 2-D&
        & case"
      return
    end if
    allocate(myfilt%hlx(1))
    call allocatehelix(filt,2)
    filt%lag=(/1,tomo_s%n(1)/)
    filt%flt=-1.
    call allocatehelix(myfilt%hlx(1),npts_lap*2)
    do i=1,npts_lap
      myfilt%hlx(1)%lag(i)=i
      myfilt%hlx(1)%lag(i+npts_lap)=tomo_s%n(1)+i-npts_lap
    end do
    if (.not. factor_filt(filt,4.00001,myfilt)) then
      write(0,*) "trouble factoring filter"
      return
    end if
    n123=product(tomo_s%n(1:3))
    allocate(myfilt%pch(n123))
    myfilt%pch=1
    if (.not. rite_nfilt("filt",myfilt,tomo_s%n,tomo_s%o,tomo_s%d))&
      & then
      write(0,*) "      ERROR: writing filter"
      return
    end if
    do i=1,npts_lap
      myfilt%hlx(1)%lag(i)=i
      myfilt%hlx(1)%lag(i+npts_lap)=tomo_s%n(1)*tomo_s%n(2)+i-npts_lap
    end do
    if (.not. rite_nfilt("filt2",myfilt,tomo_s%n,tomo_s%o,tomo_s%d))&
      & then
      write(0,*) "      ERROR: writing filter"
      return
    end if
    call deallocatehelix(filt)
    call ndeallocate(myfilt)
    calc_lap_2pass=.true.
  end function
  logical function factor_filt(filt_in,filt0,filt_out)
    type( nfilter) :: filt_out
    type( filter) :: aa, filt_in
    real :: eps
    integer:: n1, na
    integer :: i
    real :: a0, filt0
    factor_filt=.false.
    call wilson_init( 10 * maxval(filt_in%lag) )
    call wilson_factor( 20, filt0, filt_in, a0, filt_out%hlx(1))
    call wilson_close()
    factor_filt=.true.
  end function
  logical function calc_lap_1pass()
    type(nfilter) :: myfilt
    type(filter) :: filt,filt_out
    type(sep3d) :: tomo_s
    logical :: is_2d
    integer :: i,n123
    calc_lap_1pass=.false.
    if (.not. vel_tomo_structure(tomo_s)) then
      write(0,*) "      ERROR:Could not grab converted velocity&
        & structure"
      return
    end if
    is_2d=.false.
    if (sep3d_ndims(tomo_s).eq.2) then
      is_2d=.true.
    end if
    allocate(myfilt%hlx(1))
    if (is_2d) then
! call allocatehelix(myfilt%hlx(1),2)
! myfilt%hlx(1)%lag=(/1,tomo_s%n(1)/)
! myfilt%hlx(1)%flt=-.5
      myfilt%hlx(1)=lapfac2(.0001,tomo_s%n(1),npts_lap)
! call allocatehelix(filt,2)
! filt%lag=(/1,tomo_s%n(1)/)
! filt%flt=-1.
! call allocatehelix(myfilt%hlx(1),npts_lap*2)
! do i=1,npts_lap{
! myfilt%hlx(1)%lag(i)=i
! myfilt%hlx(1)%lag(i+npts_lap)=tomo_s%n(1)+i-npts_lap
! }
! if(.not. factor_filt(filt,4.00001,myfilt)){
! write(0,*) "      ERROR: factoring filter"
! return
! }
      write(0,*) "CHECK FILT",myfilt%hlx(1)%lag
      write(0,*) "CHECK FILT",myfilt%hlx(1)%flt,sum(myfilt%hlx(1)%flt)
      myfilt%hlx(1)%flt=-myfilt%hlx(1)%flt/sum(myfilt%hlx(1)%flt)
!write(0,*) "CHECK FILT",myfilt%hlx(1)%flt,sum(myfilt%hlx(1)%flt)
      n123=product(tomo_s%n(1:2))
    else
!
! WARNING WARNING WARNING
!
! This is probably wrong
!
      call allocatehelix(filt,3)
      filt%lag=(/1,tomo_s%n(1),tomo_s%n(1)*tomo_s%n(2)/)
      filt%flt=-1.
      call allocatehelix(myfilt%hlx(1),npts_lap*3)
      do i=1,npts_lap
        myfilt%hlx(1)%lag(i)=i
        myfilt%hlx(1)%lag(i+npts_lap)=tomo_s%n(1)+i-npts_lap
        myfilt%hlx(1)%lag(i+npts_lap*2)=tomo_s%n(2)*tomo_s%n(1)+i&
          &-npts_lap
      end do
      if (.not. factor_filt(filt,6.00001,myfilt)) then
        write(0,*) "      ERROR: factoring filter"
        return
      end if
      n123=product(tomo_s%n(1:3))
    end if
    allocate(myfilt%pch(n123))
    myfilt%pch=1
! myfilt%hlx(1)%flt= myfilt%hlx(1)%flt*.9
    if (.not. rite_nfilt("filt",myfilt,tomo_s%n,tomo_s%o,tomo_s%d))&
      & then
      write(0,*) "      ERROR: writing filter"
      return
    end if
    call ndeallocate(myfilt)
    calc_lap_1pass=.true.
  end function
  logical function check_pars_lap_filter(tag,struct,nfilt,norerun)
    integer :: nfilt
    type(sep3d) :: struct
    character(len=*) :: tag
    logical :: stat_ok(20),stat_rerun(20),norerun
    stat_ok=.true.
    stat_rerun=.true.
    check_pars_lap_filter=.false.
    stat_rerun(1)=check_par_file("scale_lap" ,scale_lap,tag,1.)
    stat_rerun(2)=check_par_file("npts_lap" ,npts_lap,tag,6)
    nfilt=1
    if (all(stat_ok)) then
      check_pars_lap_filter=.true.
    end if
    if (all(stat_rerun)) then
      norerun=.true.
    else
      norerun=.false.
    end if
  end function
end module
