module model_covar_mod
use nbound_bob_mod
use sep3d_struct_mod
use lapfac
use adj_mod
use build_steering2d_filters
use basic_util_mod
use build_steering3d_filters
use steer_filter_new_mod
use sep
implicit none
integer, private, save :: ns(3)
real, private, save :: ds(3),os(3)
character(len=128),private,save  :: method
type(sep3d),private,save ::  mod




contains 

logical function init_model_covar(struct,meth)
  type(sep3d) :: struct
  integer          :: n_s(3)
  real             :: o_s(3),d_s(3)
  character(len=*) :: meth
  init_model_covar=.false.
  mod=struct
  call set_n_o_d_3d(mod,ns,os,ds)
  method=meth
  init_model_covar=.true.
end function

logical function build_model_covar(filt1,filt2,nfilt)
  type(nfilter) :: filt1,filt2
  integer       :: nfilt
  build_model_covar=.false.

  if(method(1:6)=="lap_1d") then
write(0,*) "AA 1"
    if(.not. calc_lap_1pass(filt1)) then
      call seperr("trouble calculating 1-D laplacian")
      return
    end if
    nfilt=1
  else if(method(1:6)=="lap_2d") then
write(0,*) "AA 2"
    if(.not. calc_lap_2pass(filt1,filt2)) then
      call seperr("trouble calculating 2-D laplacian")
      return
    end if
    nfilt=2
  else if(method(1:5)=="steer") then
write(0,*) "AA 3"
    if(.not. calc_steer(filt1,filt2,nfilt)) then
      call seperr("trouble calculating steer filter")
      return
    end if
  else 
    write(0,*) "invalid model covariance description",trim(method)
    return
  end if
!  call nbound_bob(ns,filt1)
  if(nfilt==2) call nbound_bob(ns,filt2)
  build_model_covar=.true.
end function


logical function calc_lap_1pass(myfilt)
  type(nfilter) :: myfilt
  type(filter)  :: filt,filt_out
  type(sep3d)   :: tomo_s
  logical       :: is_2d
  integer       :: i,n123,npts_lap
  calc_lap_1pass=.false.
                                                                                         
write(0,*) "AAA x"
  is_2d=.false.; 
  if(ns(3)==1) is_2d=.true.
  allocate(myfilt%hlx(1))
  if(is_2d) then
write(0,*) "AAA 1"
    call from_param("npts_lap",npts_lap,6)
write(0,*) "AAA a",npts_lap,tomo_s%n(1)
    myfilt%hlx(1)=lapfac2(.0001,tomo_s%n(1),npts_lap)
write(0,*) "AAA b"
    myfilt%hlx(1)%flt=-myfilt%hlx(1)%flt/sum(myfilt%hlx(1)%flt)
    n123=ns(1)*ns(2)
write(0,*) "AAA 2"
  else
    !
    !  WARNING WARNING WARNING
    !
    !  This is probably wrong
    !
    call from_param("npts_lap",npts_lap,ns(1)*ns(2)*2)
    call allocatehelix(filt,3)
    filt%lag=(/1,ns(1),ns(1)*ns(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
  calc_lap_1pass=.true.
end function
                                                                                         
logical function calc_lap_2pass(filt1,filt2)
  type(nfilter) :: filt1,filt2
  type(filter)  :: filt,filt_out
  type(sep3d)   :: tomo_s
  logical       :: is_2d
  integer       :: i,n123,npts_lap
                                                                                         
  calc_lap_2pass=.false.
                                                                                         
  if(ns(3)==1) then
    write(0,*) "calling 2 pass laplacian doesn't make sense for 2-D case"
    return
  end if
                                                                                         
  allocate(filt1%hlx(1))
  allocate(filt2%hlx(2))
  call allocatehelix(filt,2)
  filt%lag=(/1,ns(1)/)
  filt%flt=-1.
  call from_param("npts_lap",npts_lap,6)
  call allocatehelix(filt1%hlx(1),npts_lap*2)
  call allocatehelix(filt2%hlx(1),npts_lap*2)
  do i=1,npts_lap
    filt2%hlx(1)%lag(i)=i
    filt2%hlx(1)%lag(i+npts_lap)=ns(1)*ns(2)+i-npts_lap
    filt1%hlx(1)%lag(i)=i
    filt1%hlx(1)%lag(i+npts_lap)=ns(1)+i-npts_lap
  end do
  if(.not. factor_filt(filt,4.00001,filt1)) then
    write(0,*) "trouble factoring filter"
    return
  end if
  n123=ns(1)*ns(2)*ns(3)
  allocate(filt2%pch(n123)); filt2%pch=1
  allocate(filt1%pch(n123)); filt1%pch=1
  do i=1,size(filt1%hlx(1)%lag)
    filt2%hlx(1)%flt(i)=filt%flt(i)
    filt1%hlx(1)%flt(i)=filt%flt(i)
  end do
  call deallocatehelix(filt)
  calc_lap_2pass=.true.
end function
                                                                                         

logical function calc_steer(filt1,filt2,nfilt)
  type(nfilter) :: filt1,filt2
  type(filter)  :: filt,filt_out
  type(sep3d)   :: tomo_s,ref_s,dip_s,dip2_s
  real, pointer :: ref3d(:,:,:),amp(:,:,:),p(:,:,:),p2(:,:,:)
  complex, pointer :: ref2d(:,:)
  logical       :: is_2d
  integer       :: i,n123,nref_use,i1,nfilt
  logical       :: steer_integrate,causal
  real          :: dip_range,azimuth_range
  real          :: min_steer_amp,max_steer_amp,hwidth_steer
  integer       :: nwind_steer_smooth,integrate_npts,integrate_noff
  integer       :: noff_steer,ntri_steer
  real          :: daz,ddip
                                                                                         
  calc_steer=.false.
                                                                                         
  call init_sep3d("dip_xz",dip_s,"INPUT")
  is_2d=.false.
  if(sep3d_ndims(dip_s)==2) is_2d=.true.
                                                                                         
                                                                                         
  nfilt=1
  if(.not. is_2d) call init_sep3d("dip_yz",dip2_s,"INPUT")
                                                                                         

  call from_param("steer_daz",daz,.5)
  call from_param("steer_ddip",ddip,.5)
  call from_param("steer_noff",noff_steer,8)
  call from_param("nwind_steer_smooth" ,nwind_steer_smooth,2)
  call from_param("ntri_steer" ,ntri_steer,2)
  call from_param("integrate_npts",integrate_npts,50)
  call from_param("integrate_noff",integrate_noff,8)
  call from_param("hwidth_steer" ,hwidth_steer,2.5)
  call from_param("min_steer_amp" ,min_steer_amp,.1)
  call from_param("max_steer_amp" ,max_steer_amp,.9)
  call from_param("azimuth_range",azimuth_range,30.)
  call from_param("dip_range",dip_range,30.)
  call from_param("causal",causal,.true.)
  call from_param("steer_integrate",steer_integrate,.true.)
  if(steer_integrate) then

    if(.not. init_filt_interp(noff_steer,dip_range,integrate_npts,ddip,mod,12,causal)) then
      write(0,*) "trouble initializing integration filters"
      return
    end if
  end if
                                                                                         
  if(is_2d) then
    call default_steering_filters_2d()
    call init_steering_filters_2d(n_in=dip_s%n,o_in=dip_s%o,d_in=dip_s%d, &
             nwind_in=nwind_steer_smooth,hwidth_in=hwidth_steer, ns_in=ntri_steer,&
       min_amp_in=min_steer_amp,max_amp_in=max_steer_amp)

    allocate(amp(dip_s%n(1),dip_s%n(2),1))
    call init_sep3d("dip_xz",dip_s,"INPUT")
    allocate(p(dip_s%n(1),dip_s%n(2),1))
    allocate(p2(dip_s%n(1),dip_s%n(2),1))
    if(.not. sreed_field("dip_xz",p)) then
      write(0,*) "TROUBLE READING IN DIP_XZ"
      return
    end if
    amp=max_steer_amp

    if(.not.steer_integrate) then
      call create_2d_steer_op(p(:,:,1)*dip_s%d(2)/dip_s%d(1),amp(:,:,1),filt1)
    else
      if(.not. create_steer_integrate_op_2d(filt1,p*dip_s%d(2)/dip_s%d(1))) then
        write(0,*) "trouble creating 2-D steer op"
        return
      end if
    end if
    deallocate(p,p2,amp)
    nullify(filt1%mis)
  else
    call default_steering_filters_3d()
    call init_steering_filters_3d(n_in=dip_s%n,o_in=dip_s%o,d_in=dip_s%d,&
             nwind_in=nwind_steer_smooth,hwidth_in=hwidth_steer, ns_in=ntri_steer,&
       min_amp_in=min_steer_amp,max_amp_in=max_steer_amp)
    allocate(amp(dip_s%n(1),dip_s%n(2),dip_s%n(3)))
    amp=max_steer_amp
    call init_sep3d("dip_xz",dip_s,"INPUT")
    allocate(p(dip_s%n(1),dip_s%n(2),dip_s%n(3)))
    allocate(p2(dip_s%n(1),dip_s%n(2),dip_s%n(3)))
    if(.not. sreed_field("dip_xz",p)) then
      write(0,*) "TROUBLE READING IN DIP_XZ"
      return
    end if
    if(.not. sreed_field("dip_yz",p2)) then
      write(0,*) "TROUBLE READING IN DIP_YZ"
      return
    end if
    if(.not.steer_integrate) then
      call create_3d_steer_2filt_op(p/dip_s%d(2)*dip_s%d(1),p2&
           /dip_s%d(3)*dip_s%d(1),amp,filt1,filt2)
    else 
      call seperr("Put it back in")
    end if
    if(.not.steer_integrate) nfilt=2
    deallocate(p,p2,amp)
  end if
  calc_steer=.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



end module
