# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
module autocorr
  use helix
  use compress
  implicit none
  contains
  function autocorr1( aa, s0, a0) result ( ss)
    type( filter) :: aa, ss
    real, intent (in), optional :: a0
    real, intent (out) :: s0
    integer :: i, j, k, n, na
    if (present (a0)) then
      s0 = a0*a0 + dot_product (aa%flt,aa%flt)
    else
      s0 = 1. + dot_product (aa%flt,aa%flt)
    end if
    na = size (aa%lag)
    call allocatehelix (ss, na*(na+1)/2)
    k = na
    if (present (a0)) then
      ss%flt (:na) = aa%flt * a0
    else
      ss%flt (:na) = aa%flt
    end if
    ss%lag (:na) = aa%lag
    do i = 1, na
      do j = i+1, na
        k = k + 1
        ss%flt (k) = aa%flt (j) * aa%flt (i)
        ss%lag (k) = aa%lag (j) - aa%lag (i)
        do n = 1, k-1
          if (ss%lag (n) .eq. ss%lag (k) .and. ss%flt (n) /= 0.) then
            ss%flt (n) = ss%flt (n) + ss%flt (k)
            ss%flt (k) = 0.
          end if
        end do
      end do
    end do
    call resize (ss)
  end function
end module