# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"

module cross_wilson
! Wilson's factorization of cross-correlation
  use helicon
  use polydiv
  implicit none
  integer, parameter, private :: n = 1024
  real, dimension( 2*n-1), private :: cross, ab, cc
  real, dimension( n), private :: b, c
  contains
  subroutine cross_factor( niter, s0, ss, a0, aa, b0, bb)
    integer, intent( in) :: niter
    ! Newton iterations
    real, intent( in) :: s0
    ! in: cross-correlat.
    type( filter), intent( in) :: ss !
    real, intent( in out) :: a0, b0
    ! out: min. phase
    type( filter), intent( out) :: aa, bb !
    real :: eps
    integer :: i, stat
    cross = 0.
    cross( n) = s0
    b( 1) = 1. ! initialize
    cross( n+ss%lag) = ss%flt
    ! cross-correlation
    do i = 1, niter
      call polydiv_init( 2*n-1, aa)
      ! divide polynoms
      stat = polydiv_lop( .false., .false., cross, ab) ! ab = S/A
! ab = ab/(a0*b0) # scale
      call polydiv_init( 2*n-1, bb)
      ! divide polynoms
      stat = polydiv_lop( .true., .false., cc, ab)
      ! cc = S/(AB')
      b( 2:n) = cc( n+1:2*n-1)/cc( n)
      ! b = + side
! b( 1) = 0.5*(cc( n) + 1.) # of (1 + cc)
      eps = maxval(abs(b(2:n)))
      call helicon_init( aa)
      ! mutliply polynoms
      stat = helicon_lop( .false., .false., b, c) ! c = A b
! c = c*a0; a0 = c( 1); aa%flt = c(1+aa%lag)/a0 # scale
      aa%flt = c(1+aa%lag)
      b( 2:n) = cc( n-1:1:-1)/cc( n)
      ! b = - side
      eps = max(maxval(abs(b(2:n))),eps)
      write( 0,*) i, eps ! "L1 norm"
      if ( eps < epsilon(a0)) then
        exit ! convergence
      end if
      call helicon_init( bb)
      ! mutliply polynoms
      stat = helicon_lop( .false., .false., b, c) ! c = B b
! c = c*b0; b0 = c( 1); bb%flt = c(1+bb%lag)/b0 # scale
      bb%flt = c(1+bb%lag)
    end do
    call polydiv_close()
  end subroutine
end module