# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module nhelicon
! Nonstationary convolution, inverse to deconvolution.
  use adj_mod
! Requires the filter be causal with an implicit "1." at the onset.
  use nhelix
  use adj_mod
  implicit none
  type( nfilter), private :: aa
  contains
  subroutine nhelicon_init ( aa_in )
    type( nfilter) :: aa_in
    aa = aa_in
  end subroutine
  function nhelicon_lop ( adj, add, xx, yy) result(stat)
    integer :: stat
    logical,intent(in) :: adj,add
    real,dimension(:) :: xx,yy
    call adjnull (adj,add,xx,yy )
    call nhelicon_lop2(adj,add,xx,yy )
    stat=0
  end function
  subroutine nhelicon_lop2(adj,add,xx,yy)
    logical,intent(in) :: adj,add
    real, dimension (:) :: xx
    real, dimension (:) :: yy
    integer :: iy, ix, ia, ip
    integer, dimension( :), pointer :: lag
    real, dimension( :), pointer :: flt
    real :: f
    if ( adj) then
! zero lag
      if (associated(aa%mis)) then
        where(.not. aa%mis)
          xx = xx + yy
end where
      else
        xx = xx + yy
      end if
    else
!this
      if (associated(aa%mis)) then
        where(.not. aa%mis)
          yy = yy + xx
end where
      else
        yy = yy + xx
      end if
    end if
    do iy = 1, size( yy)
      ip = aa%pch( iy)
      lag => aa%hlx( ip)%lag
      flt => aa%hlx( ip)%flt
      do ia = 1, size( lag)
        ix = iy - lag( ia)
        f=flt(ia)
        if ( associated( aa%mis)) then
          if (aa%mis(iy)) then
            f=0.
          end if
        end if
        if ( ix < 1) then
          cycle
        end if
        if ( ix >= size(xx)) then
          cycle
        end if
        if ( adj) then
          xx(ix) = xx(ix) + yy(iy) * f
        else
          yy(iy) = yy(iy) + xx(ix) * f
        end if
      end do
    end do
  end subroutine
  subroutine nhelicon_close()
  end subroutine
end module
