program nfilt_it
use sep 
use sep3d_struct_mod
use npolydiv
use nhelicon
use io_nfilt
implicit none
logical :: adj,add,prec
real,allocatable    :: model(:),data(:)
type(sep3d) :: mod,dat
character(len=128) :: label
type(nfilter) :: nfilt
integer :: n

call sep_init(SOURCE)
call init_sep3d("data",dat,"SCRATCH")
call init_sep3d("model",mod,"SCRATCH")
call auxclose("data")
call auxclose("model")
call auxinout("model")
call auxinout("data")
call auxin("filter")

n=product(mod%n(1:sep3d_ndims(mod)))
if (n/=product(dat%n(1:sep3d_ndims(dat)))) &
  call seperr("expecting model and data to be the same size")

if(.not. read_nfilt("filter",nfilt)) &
  call seperr("trouble reading in filter")

if(size(nfilt%pch)/=n) call seperr("filter and data size not equiv")

call from_param("add",add,.false.)
call from_param("adj",adj,.false.)
call from_param("prec",prec,.false.)

allocate(model(n),data(n))
model=0
data=0
if(adj .or. add) call sreed("data",data,size(data)*4)
if(.not. adj .or. add)call sreed("model",model,size(model)*4)
 
if(.not. prec) then
  call nhelicon_init(nfilt)
  n=nhelicon_lop(adj,add,model,data)
else
  call npolydiv_init(n,nfilt)
  n=npolydiv_lop(adj,add,model,data)
end if

if (.not. adj) then 
  call srite("data",data,size(data)*4)
else
  call srite("model",model,size(model)*4)
end if

  









end program
