module polydiv1  
! Polynomial division (recursive filtering)
  use adj_mod
  implicit none
  real ,private,  dimension (:), pointer  :: aa 
  contains
  subroutine polydiv1_init ( aa_in )
    real ,   dimension (:), pointer :: aa_in 
    aa => aa_in
  end subroutine
  function polydiv1_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 polydiv1_lop2(adj,add,xx,yy )
    stat=0
  end function 
  subroutine polydiv1_lop2(adj,add,xx,yy)
    logical,intent(in) :: adj,add 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    integer  ia, ix, iy
    real     tt
    if ( adj) then
      do ix= size(xx), 1, -1  
        tt = yy( ix)
        do ia = 1, min( size(aa), size (xx) - ix)
          iy = ix + ia
          tt =			tt -  aa( ia) * xx( iy)
        end do 
        xx( ix) = xx( ix) + tt
      end do
    else
      do iy= 1, size(xx)
        tt = xx( iy)
        do ia = 1, min( size(aa), iy-1)
          ix = iy - ia
          tt =			tt -  aa( ia) * yy( ix)
        end do 
        yy( iy) =  yy( iy) + tt
      end do
    end if
  end subroutine 
  subroutine polydiv1_close()
  end subroutine 
end module 
