# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
program mute_fill
  use sep
  use io_nfilt
  use mask1
  use npolydiv
  use sep3d_struct_mod
  use solver_mod
  use cgstep_mod
  implicit none
  type(sep3d) :: in,out,pck
  real, allocatable :: pts(:,:)
  logical, pointer ::known(:)
  type(nfilter) :: filt
  real, allocatable :: vin(:),vout(:),cube(:,:,:,:)
  integer :: n123,ndim,i,i1,i2,i3,i4,ih,niter
  real :: bias,eps
  call from_param("niter",niter,20)
  call from_param("eps",eps,.001)
  call init_sep3d("in",in,"INPUT")
  call init_sep3d("mask",pck,"INPUT")
  call init_sep3d(in,out,"OUTPUT")
  call sep3d_write_description("out",out)
  allocate(pts(pck%n(1),pck%n(2)))
  ndim=sep3d_ndims(in)
  n123=product(in%n(1:ndim))/in%n(2)
  ndim=ndim-1
  allocate(vin(n123),vout(n123),known(n123))
  call sreed("mask",pts,size(pts)*4)
  bias=sum(pts(2,:))/pck%n(2)
  pts(2,:)=pts(2,:)-bias
  vin=0
  known=.false.
  do i=1,size(pts,2)
    i1=(pts(1,i)-in%o(1))/in%d(1)+1.5
    i2=(pts(3,i)-in%o(3))/in%d(3)+1.5
    if (ndim.eq.3) then
      i3=(pts(4,i)-in%o(4))/in%d(4)+1.5
      ih=i1+(i2-1)*in%n(1)+(i3-1)*in%n(1)*in%n(3)
    else
      ih=i1+(i2-1)*in%n(1)
    end if
    vin(ih)=pts(2,i)
    known(ih)=.true.
  end do
  call mask1_init(known)
  if (.not. read_nfilt("filt",filt)) then
    call seperr("trouble reading filter")
  end if
  call npolydiv_init(size(vin),filt)
  call solver_prec(mask1_lop,cgstep,npolydiv_lop,size(vin),vout,vin,niter&
    &,eps)
  vout=vout+bias
  if (ndim.eq.3) then
    allocate(cube(in%n(1),in%n(2),in%n(3),in%n(4)))
  else
    allocate(cube(in%n(1),in%n(2),in%n(3),1))
  end if
  cube=1
  ih=0
  do i3=1,size(cube,4)
    do i2=1,size(cube,3)
      do i1=1,size(cube,1)
        ih= ih+1
        i=(vout(ih)-in%o(2))/in%d(2)+1.5
        cube(i1,i+1:,i2,i3)=0.
      end do
    end do
  end do
  call srite("out",cube,size(cube)*4)
end program
