! Copyright 1999 Stanford University
!All rights reserved
!
!Author: Sergey Fomel and Robert Clapp
module npef  
! Estimate non-stationary PEF
  use nhconest
  use npolydiv
  use polydiv
  use cgstep_mod
  use solver_mod
  use dottest
  use hconest
  use wilson
  use autocorr
  use helix
  implicit none
  contains
  subroutine find_pef( dd, aa, rr, niter, eps, nh,verb)
    integer,          intent( in)    :: niter, nh      
    ! number of iterations
    real,             intent( in)    :: eps	       ! epsilon
    type( nfilter)                   :: aa             
    ! estimated PEF output.
    type( nfilter),   intent( in)    :: rr             
    ! roughening filter.
    type(filter)                     :: dummy
    real,    dimension(:), pointer   :: dd             ! input data
    real,    dimension(:,:), pointer   :: rmov             
    ! input data
    integer                          :: ip, ih, np, nr 
    ! filter lengths
    real, dimension (:), allocatable :: flt	       
    ! np*na filter coefs
    real, dimension(:), allocatable  :: temp1,temp2
    real :: dot1(2),dot2(2),scale,tempr
    logical :: verb
    integer :: i1,ierr
    np = size( aa%hlx)			! data   length
    nr = np*nh
    allocate( flt( nr))
    flt=0.
    call nhconest_init( dd, aa, nh)
    call npolydiv_init( nr, rr)
    dd=dd-sum(dd)/size(dd)
    call solver_prec( nhconest_lop, cgstep, x= flt, dat= -dd, niter=&
      & niter,prec= npolydiv_lop, nprec= nr, eps= eps,verb=verb)
    call   cgstep_close()
    call npolydiv_close()
    call nhconest_close()
    do ip = 1, np 
      do ih = 1, size( aa%hlx(ip)%lag)
        aa%hlx( ip)%flt( ih) = flt( (ip-1)*nh + ih)
      end do
    end do 
    deallocate(flt)
    return
  end subroutine 
end module 
