!!$
!!$=head1 NAME
!!$
!!$ Lloyd_vel - Select reference velocities 
!!$
!!$=head1 SYNOPSIS
!!$
!!$Lloyd_vel  < vel.H >out.H
!!$
!!$
!!$=over 4
!!$
!!$
!!$=item nref (integer)
!!$
!!$  Number of reference velocities
!!$
!!$
!!$=item min_region_pct (real)
!!$
!!$  [1.]  Minimum percentage of the velocity field associated with a given reference
!!$
!!$=item min_dev 
!!$
!!$  [6.]  Minimum deviation (in percent) between reference velocities
!!$
!!$=item niter_lloyd 
!!$
!!$  [30] Number of iteratons of the lloyd to run
!!$
!!$=item When doing offset gathers  ahy_n,ahy_d,ahy_o = axis 4 of D and R
!!$
!!$=item nzs must be a multiple of az__n
!!$
!!$=item nws*nmpi_threads*nomp_threads must be a multiple of aw__n
!!$
!!$=cut
!!$
!!$

program lloyd
use sep3d_struct_mod
use lloyd_mod
use sep
implicit none
type(sep3d) :: attr_s
real, allocatable :: attr(:,:),ref_attr(:,:),slice(:,:),junk(:,:),scales(:)
integer,allocatable :: iregion(:),nstart(:)
integer :: nref,i,iref,izs,ia,ib,imap,n(4),iter0,jiter,niter,ntr,it

call sep_init(SOURCE)

call from_param("nref",nref,20)
call init_lloyd();


call init_sep3d("in",attr_s,"INPUT")

!if(sep3d_ndims(attr_s)/=4) &
!  call seperr("expecting 4-d cube attr,x,y,z")


n=1
n(1:sep3d_ndims(attr_s))=attr_s%n(1:sep3d_ndims(attr_s))
write(0,*) 'ndim',n

allocate(scales(n(1)))
allocate(nstart(n(1)))
call from_param("nstart",nstart)

!write(0,*) nstart,'hi brad'
!call seperr("asdf")

scales=1
call from_param("scales",scales,scales)
if(product(nstart) > nref) &
  call seperr("product(nstart) > nref")

call to_history("n2",n(1))
call to_history("n1",nref)
call to_history("n3",n(3))
write(0,*) "Xyy",nref
call from_param("iter0",iter0,0)
call from_param("niter",niter,20
call from_param("jiter",jiter,niter)
call sep_put_data_axis_par("out",4,niter,real(iter0),real(jiter),"Iteration")



allocate(iregion(n(2)*n(3)),junk(n(1),n(2)*n(3)),slice(n(2)*n(3),n(1)), attr(nref,n(1)))

do it=1,niter
   ntr=iter0+(it-1)*jiter
   call sreed("in",junk,size(junk)*4)
   slice=transpose(junk)
   do i=1,n(1)
     slice(:,i)=slice(:,i)*scales(i)
   end do
   iref= lloyd_go(slice, attr, nref, iregion,nstart,ntr)
   do i=1,n(1)
     slice(:,i)=slice(:,i)/scales(i)
     attr(:,i)=attr(:,i)/scales(i)
   end do

     if(iref < nref) then
       attr(iref+1:,:)=-1
     end if
write(0,*) "WRITING OUT",size(attr)
     call srite("out",attr,size(attr)*4)
end do

end program
