# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module steer_filter_new_mod
  use nhelix
  use sep
  use sep3d_struct_mod
  implicit none
  type(sep3d), private , save :: mod
  integer, private, save :: n1,n2
  real, private, pointer, save :: dips(:),func(:)
  integer, private, pointer, save :: lags(:)
  logical, private, pointer, save :: have(:)
  logical, private,save :: causal
  real, private, save :: dip_range
  integer, private, save :: n_integrate,iover,idown
  real, private, save :: ddip,pow
  interface interp_filter
    module procedure interp2d_filter
  end interface
  interface create_steer_integrate_op
    module procedure create_steer_integrate_op_2d
  end interface
  contains
  logical function init_filt_interp(nsearch,diprange,nintegrate,d_dip&
    &,model,axis,caus)
    logical :: caus
    integer :: i1,i2,n2s,axis
    real :: diprange,d_dip
    integer ::nintegrate,nsearch
    type(sep3d) :: model
    init_filt_interp=.false.
    mod=model
    n1=mod%n(1)
    if (axis.eq.12) then
      iover=mod%n(1)
      idown=1
    else if (axis.eq.13) then
      iover=mod%n(1)*mod%n(2)
      idown=1
    else
      write(0,*) "axis",axis,"don't understand axis specification"
      return
    end if
    call from_param("sin_pow",pow,2.)
    causal=caus
    ddip=d_dip
    n_integrate=nintegrate
    dip_range=diprange
    dip_range=dip_range*atan(1.)/45.
    write(0,*) "CHECK DUIP",dip_range
    if (causal) then
      init_filt_interp=init_causal_interp(nsearch)
    else
      init_filt_interp=init_noncausal_interp()
    end if
  end function
  logical function init_noncausal_interp()
    integer :: i1,i2
    allocate(lags(24),func(24),dips(24),have(24))
    lags( 1:8 )=(/-iover,-idown-iover,-idown,iover-idown,iover,iover&
      &+idown,idown,idown-iover/)
    dips( 1:8 )=(/-540.,-495.,-450.,-405.,-360.,-315.,-270.,-225./)
    lags( 9:16)=(/-iover,-idown-iover,-idown,iover-idown,iover,iover&
      &+idown,idown,idown-iover/)
    dips( 9:16)=(/-180.,-135.,-90.,-45.,0.,45.,90.,135./)
    lags(17:24)=(/-iover,-idown-iover,-idown,iover-idown,iover,iover&
      &+idown,idown,idown-iover/)
    dips(17:24)=(/180.,225.,270.,315.,360.,405.,450.,495./)
    dips=dips*atan(1.)/45.
    init_noncausal_interp=.true.
  end function
  logical function init_causal_interp(nsearch)
    integer :: i1,i2,nsearch
    allocate(lags(nsearch*2+2))
    allocate(func(nsearch*2+2))
    allocate(dips(nsearch*2+2))
    allocate(have(nsearch*2+2))
!DO THE DIPS BOTH CASES
    do i2=1,nsearch*2+1
      dips(i2)=atan(real(i2-nsearch-1))
    end do
    dips(nsearch*2+2)=atan(1.)*2.
    do i1=1,nsearch*2+1
      lags(i1)=iover+(i1-nsearch-1)*idown
    end do
    lags(nsearch*2+2)=idown
    init_causal_interp=.true.
  end function
  logical function clean_filt_interp()
    deallocate(lags,func,dips,have)
    clean_filt_interp=.true.
  end function
  logical function assign_filt_database_2d(filt,mydip,ndip)
    type(nfilter) :: filt
    real :: mydip(:,:,:)
    integer :: idip,iaz,i,i1,i2,i3,ndip
    real :: odip,oaz
    assign_filt_database_2d=.false.
    odip=minval(mydip)
    allocate(filt%pch(size(mydip)))
    i=0
    do i3=1,size(mydip,3)
      do i2=1,size(mydip,2)
        do i1=1,size(mydip,1)
          i= i+1
          idip=min(ndip,max(1,nint( (mydip(i1,i2,i3) - odip)/ddip) +1)&
            &)
          filt%pch(i)=idip
        end do
      end do
    end do
    assign_filt_database_2d=.true.
  end function
  logical function create_steer_integrate_op_2d(filt,p1)
    type(nfilter) :: filt
    integer :: ndip
    real :: p1(:,:,:)
    real, pointer :: mydip(:,:,:)
    create_steer_integrate_op_2d=.false.
    allocate(mydip(size(p1,1),size(p1,2),size(p1,3)))
    mydip=atan(p1)*45./atan(1.)
    if (.not. create_filt_database_2d(filt,mydip,ndip)) then
      write(0,*) "trouble creating filt database"
      return
    end if
    if (.not. assign_filt_database_2d(filt,mydip,ndip)) then
      write(0,*) "trouble assigning filters"
      return
    end if
    deallocate(mydip)
    create_steer_integrate_op_2d=.true.
  end function
  logical function create_filt_database_2d(filt,mydip,ndip)
    type(nfilter) :: filt
    real :: mydip(:,:,:)
    integer :: ndip
    real :: odip,d
    integer :: idip,i
    create_filt_database_2d=.false.
    odip=minval(mydip)
    ndip=ceiling((maxval(mydip)-odip)/ddip)+1
    allocate(filt%hlx(ndip))
    nullify(filt%mis)
    do idip=1,ndip
      d=odip+(idip-1)*ddip
      if (.not. interp2d_filter(d,filt%hlx(idip))) then
        write(0,*) "trouble creating filter", idip,d
        return
      end if
    end do
!i=ndip ;
!filt%hlx(i)%lag=0 ; filt%hlx(i)%flt=0.
!call allocatehelix(filt%hlx(i),1)
    create_filt_database_2d=.true.
  end function
  logical function interp2d_filter(dip,filt)
    real :: dip,dip2
    type(filter) :: filt
    interp2d_filter=.false.
    dip2=dip*atan(1.)/45.
    if (.not. calc_2d_func(dip2)) then
      write(0,*) "trouble calculating 2d filter"
      return
    end if
!if(.not. calc_2d_lagrange(dip2)){
! write(0,*) "trouble calculating 2d filter"
! return
!}
    if (.not. make_filter(filt)) then
      write(0,*) "trouble making filter"
      return
    end if
    interp2d_filter=.true.
  end function
  logical function calc_2d_func(dip)
    real :: dip
    real :: odip,myd,delta_dip,ddip
    real :: val,f
    integer :: i,i1
    calc_2d_func=.false.
    odip=dip-dip_range/2.
    ddip=dip_range/(n_integrate-1)
    delta_dip=3.1415926536/dip_range
    func=0.
    have=.false.
    do i=1,n_integrate
      myd=(i-1)*ddip+odip
      i1=1
!FIND THE LARGEST DIPS ARRAY VALUE LESS THEN THE CURRENT DIP
      if (myd >= dips(size(dips,1))) then
        i1=size(dips,1)
      else
        do while (dips(i1) < myd)
          i1= i1+1
        end do
        i1= i1-1
      end if
      if (i1 > 0 .and. i1 < size(dips,1)) then
        val=sin(delta_dip*(myd-odip))**pow
        f=(myd-dips(i1))/(dips(i1+1)-dips(i1))
        func(i1)= func(i1)+(1.-f)*val
        func(i1+1) = func(i1+1) +f *val
        have(i1:i1+1)=.true.
      end if
    end do
    func=-func/sum(func)
    calc_2d_func=.true.
  end function
  logical function make_filter(filt)
    type(filter) :: filt
    integer :: i1,i2,nc
    make_filter=.false.
    nc=0
    do i1=1,size(have)
      if (abs(func(i1))>.0001) then
        nc= nc+1
      end if
    end do
    call allocatehelix(filt,nc)
    nc=0
    do i1=1,size(have)
      if (abs(func(i1)) > .0001) then
        nc= nc+1
        filt%flt(nc)=func(i1)
        filt%lag(nc)=lags(i1)
      end if
    end do
    filt%flt=-1.*filt%flt/sum(filt%flt)
    make_filter=.true.
  end function
  logical function lagrange_2d_tab(x_out,y_out,x_in,y_in,tab)
    real :: x_out,y_out,x_in(:),y_in(:),tab(:,:)
    integer :: i1,i2
    real :: num,den
    real :: tabx(size(x_in)),taby(size(y_in))
    lagrange_2d_tab=.false.
    if (size(y_in).ne.size(tab,2) .or. size(x_in).ne.size(tab,1)) then
      write(0,*) "expecting size(x_in==tab(1,:)) and size(y_in(tab(:,1&
        &)))"
      return
    end if
    tabx=0
    taby=0
    if (.not. lagrange_1d_tab(x_out,x_in,tabx)) then
      write(0,*) "trouble interpolating x"
      return
    end if
    if (.not. lagrange_1d_tab(y_out,y_in,taby)) then
      write(0,*) "trouble interpolating y"
      return
    end if
    do i2=1,size(tab,2)
      do i1=1,size(tab,1)
        tab(i1,i2)=tabx(i1)*taby(i2)
      end do
    end do
    lagrange_2d_tab=.true.
  end function
  logical function lagrange_1d_tab(x_out,x_in,tab)
    real :: x_out,x_in(:),tab(:)
    integer :: i1,i2
    real :: num,den
! lagrange_1d_tab=.false.
! tab=1
!
! do i2=1,size(x_in){
! do i1=1,size(x_in){
! if(i1!=i2) tab(i2)=tab(i2)*(x_out-x_in(i1))/(x_in(i2)-x_in(i1))
! }
! }
    lagrange_1d_tab=.true.
  end function
end module
