# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!factor_filt_mod
!
!Usage:
!call factor_filt_mod (arg1, arg2 arg3)
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Mon Oct 29 14:33:42 PST 2001
!
!Purpose:
!
!
module factor_filt_mod
  use patch_bob
  use mkwallwt
  use nhelicon_patch
  use npolydiv_patch
  use sep3d_struct_mod
  use sep
  use autocorr
  use filt_lagrange_mod
  use wilson_bob
  implicit none
  integer, dimension( :), pointer,save,private :: npatch, nwall,&
    & nwind
  real, dimension( :), pointer,save,private :: windwt
  contains
  logical function factor_filt_init(filt,data,dip12,dip13)
    type(sep3d) :: data
    type(nfilter) :: filt
    real :: dip12(:),dip13(:)
    real, pointer::wt(:)
    optional :: dip13
    logical :: is2d
    integer :: ndim
    real, pointer :: wind_wt(:)
    real :: odip12,odip13,ddip12,ddip13,mx,dip_13,dip_12
    integer, pointer :: n(:),nd(:)
    integer, pointer :: np(:)
    integer :: ndip12,ndip13,i12,i13,ncoef13,ncoef12,ic,i
    type(filter) :: f12,f13,a12,a13,tot,final
    logical :: valid
    real :: z12,z13,ztot,eps,s0,a0
    integer,allocatable :: lag(:)
    factor_filt_init=.false.
!check for stupid errors
    if (present(dip13)) then
      if (sep3d_ndims(data).ne.3) then
        write(0,*) "data is not 3-d but 2 dips supplied"
        return
      end if
      is2d=.false.
      ndim=3
      call from_param("ndip13",ndip13,100)
      call from_param("min_dip13",odip13,minval(dip13))
      call from_param("max_dip13",mx,maxval(dip13))
      ddip13=(mx-odip13)/(ndip13-1)
    else
      if (sep3d_ndims(data).ne.2) then
        write(0,*) "data is not 2-d but 1 dip supplied"
        return
      end if
      is2d=.true.
      ndim=2
      ndip13=1
    end if
    allocate(lag(product(data%n(1:ndim-1)+7)))
    do i=1,product(data%n(1:ndim-1)+7)
      lag(i)=i
    end do
    call wilson_bob_init( size(lag),1)
    call from_param("ncoef12",ncoef12,5)
    call from_param("ncoef13",ncoef13,5)
    call from_param("ndip12",ndip12,100)
    call from_param("min_dip12",odip12,minval(dip12))
    call from_param("max_dip12",mx,maxval(dip12))
    ddip12=(mx-odip12)/(ndip12-1)
!grab patch size
    allocate(n(ndim),np(ndim),nd(ndim))
    n=25
    call from_param("nwind_patch",n,n)
    np=ceiling(real(data%n(1:ndim))/(n*1.2))
    call from_param("npatch",np,np)
    nd=data%n(1:ndim)
!
!FIX THIS FIX THIS
!DO SOME BOUNDING
!
    allocate(wind_wt(product(n)))
    wind_wt=1.
    call patch_init_op(np,nd,n,wind_wt)
!time to construct the filters
    allocate(filt%hlx(ndip13*ndip12),filt%pch(product(data%n(1:ndim)))&
      &)
    ic=0
    do i13=1,ndip13
      dip_13=odip13+ddip13*(i13-1)
      if (.not. is2d) then
        if (.not. lagrange_it(n(1)*n(2),f13,dip_13,ncoef12)) then
          call erexit("trouble initializing lagrange ")
        end if
        a13=autocorr1(f13,z13,1.)
        call deallocatehelix(f13)
      end if
      do i12=1,ndip12
        dip_12=odip12+ddip12*(i12-1)
        ic= ic+1
        if (.not. lagrange_it(n(1),f12,dip_12,ncoef12)) then
          call erexit("trouble initializing lagrange ")
        end if
        if (is2d) then
          tot=autocorr1(f12,ztot,1.)
! write(0,*) "1HECK ",dip_12,f12%flt,f12%lag
! write(0,*) "2HECK ",dip_12,tot%flt,tot%lag
        else
          a12=autocorr1(f12,z12,1.)
          call allocatehelix(tot,size(a13%lag)+size(a12%lag))
          do i=1,size(a13%lag)
            tot%flt(i)=a13%flt(i)
            tot%lag(i)=a13%lag(i)
          end do
          do i=1,size(a12%lag)
            tot%flt(i+size(a12%lag))=a12%flt(i)
            tot%lag(i+size(a12%lag))=a12%lag(i)
          end do
          ztot=z13+z12
          call deallocatehelix(a12)
        end if
        call deallocatehelix(f12)
        valid=.false.
        eps=.002
        do while (.not. valid)
          s0=ztot+eps
          call allocatehelix(filt%hlx(ic),size(lag))
          filt%hlx(ic)%lag=lag
          if ( wilson_bob_factor (15, s0, tot, a0,filt%hlx(ic),.true.)&
            &) then
            valid=.true.
          end if
          eps=eps*2.
        end do
        call resize (filt%hlx(ic), .001)
        filt%hlx(ic)%flt = 0.
        if ( wilson_bob_factor (10, s0, tot, a0, filt%hlx(ic))) then
          valid=.true.
        end if
! if(.not. lagrange_it(n(1),filt%hlx(ic),dip_12,ncoef12))
! call erexit("trouble initializing lagrange ")
! write(0,*) "CHECK ",dip_12,filt%hlx(ic)%flt,filt%hlx(ic)%lag
        call deallocatehelix(tot)
      end do
      if (.not. is2d) then
        call deallocatehelix(a13)
      end if
    end do
    do i=1,size(filt%pch)
      i12=nint((dip12(i)-odip12)/ddip12+1.)
      if (is2d) then
        i13=1
      else
        i13=nint((dip13(i)-odip13)/ddip13+1.)
      end if
      filt%pch(i)=(i13-1)*ndip12+i12
    end do
    deallocate(lag)
    factor_filt_init=.true.
  end function
!patching init op
  subroutine patch_init_op(n_patch,n_wall,n_wind,wind_wt)
    integer, dimension( :), pointer :: n_patch, n_wall, n_wind
    real, dimension(:), pointer :: wind_wt
    allocate(npatch(size(n_patch)))
    allocate(nwall(size(n_wall)))
    allocate(nwind(size(n_wind)))
    allocate(windwt(size(wind_wt)))
    npatch=n_patch
    nwall=n_wall
    nwind=n_wind
    windwt=wind_wt
  end subroutine
  subroutine patch_op_close()
    deallocate(npatch,nwind,nwall,windwt)
  end subroutine
!Helicon op
  integer function patch_nhelicon_op(adj,add,model,data)
    logical :: adj,add
    real :: model(:),data(:)
    patch_nhelicon_op=patch_me_op(nhelicon_patch_op,adj,add,model,data&
      &)
  end function
!polydiv op
  integer function patch_npolydiv_op(adj,add,model,data)
    logical :: adj,add
    real :: model(:),data(:)
    patch_npolydiv_op=patch_me_op(npolydiv_patch_op,adj,add,model,data&
      &)
  end function
!filter construction
!filter division
!patching op
  integer function patch_me_op(oper,adj,add, modl, data)
    interface
      integer function oper(adj, add, modl, data,index)
        logical, intent (in) :: adj, add
        real, dimension( :) :: modl, data
        integer,dimension(:) :: index
      end function
    end interface
    logical :: adj,add
    real, dimension( :) :: modl,data
    real, dimension( size( modl)) :: wallwt,junk
    real, dimension( size( windwt)) :: winmodl, windata
    integer, dimension( size( windwt)) :: index
    integer :: i, stat1, stat2
    call patch_bob_init( npatch, nwall, nwind)
    call adjnull(adj,add,modl,data)
    call wallwtn( npatch, nwall, nwind, windwt, wallwt)
! call seperr("");
    if (adj) then
      junk=data
      data = data * wallwt
    else
      junk=data
      data=0
    end if
    do i = 1, product( npatch)
      call patch_index_extract(index)
      if (adj) then
        stat1 = patch_bob_lop( .false., .false., data, windata)
        windata=windata*windwt
      else
        stat1 = patch_bob_lop( .false., .false., modl, winmodl)
      end if
      if (adj) then
        winmodl=0.
      else
        windata=0.
      end if
      stat2 = oper( adj,add,winmodl, windata,index)
      if (adj) then
        stat1 = patch_bob_lop( .true. , .true. , modl, winmodl)
      else
        stat1 = patch_bob_lop( .true. , .true. , data, windwt * windata&
          &)
      end if
      call patch_bob_close ()
    end do
    if (adj) then
      data=junk
    else
      data = junk+ data * wallwt
    end if
    patch_me_op=0
  end function
end module
