# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module steering
! Bob Clapp's steering filter
  use helix
  implicit none
  integer, private :: n ! helix period
  real, private :: hw, dc ! triangle half-width, filter d.c.
  logical, private,save :: sergey
  contains
  subroutine steer_init( hwidth, amp, n1,sergey_in)
    real, intent( in) :: hwidth, amp
    integer, intent( in) :: n1
    logical, optional :: sergey_in
    if (present(sergey_in)) then
      sergey=sergey_in
    else
      sergey=.false.
    end if
    hw = hwidth
    dc = amp
    n = n1
  end subroutine
  subroutine steering2( slope, aa)
    type(filter) :: aa
    real, intent( in) :: slope
    integer :: h, nh,p0i
    real :: p0,center ,p
    integer :: i, j, nf, nc,ia,ip
    real :: f, xi
    nh = size( aa%lag)
    center=n+slope
    if (.not. sergey) then
      do h = 1, nh
        aa%lag(h) = ceiling(center- real(nh)/2.)-1+h
        aa%flt(h) = abs(hw-abs(real(aa%lag(h))-center))/hw
      end do
    else
      p=slope
      ip=floor(slope)
      p0i = (nh-1)/2
      p=slope-ip
      nc = (nh+1)*0.5
      do ia = 1, nh
        aa%lag (ia) = n+ia-p0i-1+ip
      end do
      do i = 1, nh
        f = 1.
        xi = p + nc - i
        do j = 1, nh
          if (i /= j) then
            f = f * (1. + xi / (i - j))
          end if
        end do
        aa%flt (i) = f
      end do
    end if
    aa%flt = -dc/sum( aa%flt)*aa%flt
  end subroutine
end module
