module lapfac
! Factor 2-D Laplacian.
  use wilson
!  use helix
  implicit none
  contains
  function lapfac2( eps, n1, na) result (aa)
    type( filter) :: aa, lap
    real, intent( in) :: eps
    integer, intent( in) :: n1, na
    integer :: i
    real :: a0, lap0
    write(0,*) "Ab1"
    call allocatehelix( lap, 2) ! laplacian filter
    write(0,*) "0A1"
    lap0 = 4. + eps ! zero lag coeff.
    write(0,*) "1A1"
    lap%lag = (/ 1, n1 /) ! lag(1)= 1; lag(2)=n1 # one side only
    write(0,*) "2A1"
    lap%flt = -1. ! flt(1)=-1; flt(2)=-1
    write(0,*) "3A1"
    call allocatehelix( aa, 2*na) ! laplacian derivative
    write(0,*) "4A1"
    aa%flt = 0.
! probably done already in allocation.
    write(0,*) "5A1"
    do i = 1, na
      aa%lag( i ) = i ! early lags (first row)
      aa%lag( na+i) = n1 + i - na ! late lags (second row)
    end do
    write(0,*) "AA1"
    call wilson_init( 10 * n1 )
    write(0,*) "AA2"
    call wilson_factor( 20, lap0, lap, a0, aa)
    write(0,*) "AA3"
    call wilson_close()
    write(0,*) "AA4"
    call deallocatehelix( lap)
    write(0,*) "AA5"
  end function
end module
