#<
#nfilter_it
#
#Usage:
#nfilter_it.x  <in.H >out.H 
#
#Input Parameters:
#
#
#Output Parameters:
#
#
#
#Description:
#
#
#
#>
#%end of self-documentaiton
#-------------------------------------------------
#
#Author: Robert Clapp, ESMB 463, 7230253
#
#Date Created:Thu Feb 17 22:54:20 PST 2000
#
#Purpose: 
#

program nfilter_it{ 
	use sep
	use npef
	use build_filter_mod
	use createnhelixmod
	use mask1
  integer ::  n(2),j(10),center(2),na(2),ncoef,n123,niter,niter_prec
	real ::  o(2),d(2),hwidth,amp,eps
	character(len=128) :: label
  real, pointer :: data(:),covar(:)
	logical, pointer :: mask(:)
	type(nfilter) :: rough,filt

	write(0,*) "i die w"
	#read in parameters
	call sep_get_data_axis_par("in",1,n(1),o(1),d(1),label)
	write(0,*) "i die 2"
	call sep_get_data_axis_par("in",2,n(2),o(2),d(2),label)
	write(0,*) "i die 2a"
	call from_param("j",j,j)
	write(0,*) "i die 21"
	call from_param("na",na,(/5,2/))
	call from_param("center",center,(/ceiling(na(1)/2.),1/))
	write(0,*) "i die 2"
	call from_param("hwidth",hwidth,4.)
	write(0,*) "i die 2"
	call from_param("amp",amp,.999)
	write(0,*) "i die 2"
	call from_param("eps",eps,0.)
	write(0,*) "i die 2"
	call from_param("niter",niter,30)
	call from_param("niter_prec",niter_prec,30)

	
	
	
	#build filters
	 call init_build_filters(n,j,center,na,o,d,hwidth,&
       amp,.false.,.false.,d(1)*5,3.)

  call build_2d_filter(filt,rough,ncoef)

	n123=product(n)
	write(0,*) "i die shwoehere",n123,ncoef
	write(0,*) "rough",size(rough%hlx),rough%hlx(1)%flt
	rough%pch=1
	write(0,*) "check",rough%hlx(1)%lag
	deallocate(rough%hlx(1)%flt,rough%hlx(1)%lag)
	allocate(rough%hlx(1)%flt(2),rough%hlx(1)%lag(2))
	rough%hlx(1)%flt=(/-.5,-.5/)
	rough%hlx(1)%lag=(/ncoef,ncoef*n(1)/j(1)/)
	

	
	#find the npef
	allocate(covar(n123))
	call sreed("covar",covar,4*n123)
	write(0,*) "i die shwoehere",minval(covar),maxval(covar),j
  call find_pef(covar,filt,rough,niter,eps,ncoef,.true.)

	#read in the missing data
	allocate(data(n123))
	call sreed("in",data,n123*4)


	#apply the npef
  call npolydiv_init(n123,filt)
	allocate(mask(n123))
	mask =  (data != 0.)
  call mask1_init(mask)
  call solver_prec(mask1_lop,cgstep,niter=niter_prec,x=data, dat=data,prec=npolydiv_lop,nprec=size(mask),eps=0.,verb=.true.)


	#write it out
	call srite("out",data,size(data)*4)



}