module tomo_reg_prec_mod
use sep
use nhelicon2
use npolydiv2
use nhelicon
use tomo_eps_mod
use npolydiv
use helicon
use chain_mod
use sep3d_struct_mod
use io_nfilt
implicit none
logical, private, save :: prec,dual,new
real, private,save,pointer :: space(:),space2(:)
type(nfilter) :: filt1_prec,filt2_prec
type(nfilter) :: filt1_reg,filt2_reg,filta,filtb
type(filter) :: filt1
integer, private, save :: nmod



contains


logical function load_reg_prec(slow_s,reg0) 
type(sep3d)   :: slow_s,reg_s
real, pointer :: reg0(:),junk(:)
integer, external :: auxin,sreed
integer :: ierr
integer :: i
integer :: ns
real :: sm

load_reg_prec=.false.

nmod=product(slow_s%n(1:sep3d_ndims(slow_s)))

  allocate(reg0(nmod))
if(-1==auxin("reg0")) then
  reg0=0.
else

  call init_sep3d("reg0",reg_s,"INPUT")
  if(product(reg_s%n) /= size(reg0)) then
    write(0,*) "size of reg0 must be the same as the size of the slowness space"
    return
  end if
    
  if(4*size(reg0)/=sreed("reg0",reg0,size(reg0)*4)) then
    write(0,*) "trouble reading in reg0"
    return
  end if
    call srite("bb.H",reg0,size(reg0)*4)

  reg0=-reg0
  ns=0;
  sm=0
  do i=1,size(reg0)
    if(abs(reg0(i)) <.3.and. i < 200) then
      sm=sm+reg0(i)
      ns=ns+1
    end if
  end do
  sm=sm/ns
 ! sm=.06
   do i=1,size(reg0)
    if(abs(reg0(i)) <.3 .and. i < 200) then
    reg0(i)=reg0(i)-sm
   ! write(0,*) "SUBTA",reg0(i),sm
    else
      reg0(i)=0
    end if
  end do
  !write(0,*) "CHECK SM",sm,maxval(reg0),minval(reg0)
end if
    call srite("bb.H",reg0,size(reg0)*4)

call from_param("new",new,.true.)
call from_param("prec",prec,.true.)
if(0==auxin("filt_reg")) then
  write(0,*) "filt not supplied"
  return
end if


if(.not. read_nfilt("filt_reg",filt1_reg)) then
  write(0,*) "trouble reading nfilter"
  return
end if
dual=.false.
if(-1/=auxin("filt2_reg")) then
  dual=.true.
  if(.not. read_nfilt("filt2_reg",filt2_reg)) then
    write(0,*) "trouble reading nfilter"
    return
  end if
 ! allocate(space(nmod))

end if
  if(new)then
    if(.not. read_nfilt("filt1",filtb)) then
    write(0,*) "trouble reading nfilter"
    return
  end if 
  if(.not. read_nfilt("filt2",filta)) then
    write(0,*) "trouble reading nfilter"
    return
  end if
end if
  allocate(space(nmod),space2(nmod))
if(prec) then
  if(.not. read_nfilt("filt_prec",filt1_prec)) then
    write(0,*) "trouble reading nfilter"
    return
  end if
  dual=.false.
  if(-1/=auxin("filt2_prec")) then
    dual=.true.
    if(.not. read_nfilt("filt2_prec",filt2_prec)) then
      write(0,*) "trouble reading nfilter"
      return
    end if
 !   allocate(space(nmod))
  end if
end if 

call allocatehelix(filt1,3)

ierr=slow_s%n(1)
filt1%lag=(/ierr-1,ierr,ierr/)
filt1%flt=(/-.33,-.33,-.33/)
call helicon_init(filt1)
if(.not. new) then
  call nhelicon_init(filt1_reg)

  if(dual) call nhelicon2_init(filt2_reg)

  if(prec) then
    call npolydiv_init(nmod,filt1_prec)
    if(dual) call npolydiv2_init(nmod,filt2_prec)
  end if
else
  call npolydiv_init(nmod,filta)
   call nhelicon_init(filtb)
end if

!allocate(junk(size(reg0)))
!junk=reg0
!ierr=eps_op(.false.,.false.,junk,reg0)
!deallocate(junk)
!reg0=0.
load_reg_prec=.true.
end function

integer function prec_op(adj,add,model,data)
logical, intent(in) :: adj,add
real :: model(:),data(:)
prec_op=0
if(new) then
  !call chain(npolydiv_lop,nhelicon_lop,adj,add,model,data,space)
 ! call adjnull(adj,add,model,data)
if(1==1) then
 prec_op=npolydiv_lop(adj,add,model,data)
endif

!write(0,*) "BB",adj
 if(1==2) then
 if(.not. adj) then
   data=model+data
  else
      model=data+model
  end if
  end if
else
if (prec) then
  if(dual) then
    call chain(npolydiv_lop,npolydiv2_lop,adj,add,model,data,space)
  else
    prec_op=npolydiv_lop(adj,add,model,data)
  end if
else 
  if(.not. adj) then
    data=model+data
  else
     model=data+model
  end if
end if
end if
end function

  
integer function reg_range()
reg_range=nmod
end function

integer function reg_domain()
reg_domain=nmod
end function

integer function prec_range()
prec_range=nmod
end function

integer function prec_domain()
prec_domain=nmod
end function


integer function reg_op(adj,add,model,data)
logical, intent(in) :: adj,add
real :: model(:),data(:)
integer :: ierr

reg_op=0;
if(new) then
!write(0,*) "in reg",adj
 ! call chain(nhelicon_lop,prec_op,adj,add,model,data,space2)
ierr=nhelicon_lop(adj,add,model,data)
else
if(dual) then
  call chain(nhelicon_lop,nhelicon2_lop,adj,add,model,data,space)
else
  ierr=nhelicon_lop(adj,add,model,data)
 ! ierr=helicon_lop(adj,add,model,data)
end if
end if
end function


end module
