! Copyright 1999 Stanford University
!All rights reserved
!
!Author: Sergey Fomel and Robert Clapp
module npolydiv  
! Helix polynomial division
  use adj_mod
  use nhelix
  implicit none
  integer, private :: nd 
  type( nfilter), private :: aa 
  real ,private,  dimension (:), allocatable  :: tt 
  integer ,private,  dimension(:), pointer  :: lag 
  real ,private,     dimension(:), pointer  :: flt 
  contains
  subroutine npolydiv_init ( nd_in,aa_in )
    integer    :: nd_in 
    type( nfilter)    :: aa_in 
    nd = nd_in 
    aa = aa_in 
    if (.not. allocated(tt)) then
      allocate(tt ( nd)) 
    end if
  end subroutine
  function npolydiv_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 npolydiv_lop2(adj,xx,yy )
    stat=0
  end function 
  subroutine npolydiv_lop2(adj,xx,yy)
    logical,intent(in) :: adj 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    integer                        :: ia, ix, iy, ip
    tt = 0.
!    allocate(flt(size(aa%hlx(1)%lag)))
!    allocate(lag(size(aa%hlx(1)%flt)))
    if ( adj) then
      tt = yy
      do iy= nd, 1, -1  
!		if(associated(aa%mis)){if(aa%mis(iy)){ tt(iy)=0.; next}}
        if (associated(aa%mis)) then
          if (aa%mis(iy)) then
            tt(iy)=0.
            cycle
          end if
        end if
        ip = aa%pch( iy)
        lag => aa%hlx( ip)%lag
        flt => aa%hlx( ip)%flt
        do ia = 1, size( lag)
          ix = iy - lag( ia)
          if ( ix < 1) then
            cycle
          end if
!			if(ix <0 || ix >size(xx)) write(0,*) "1x error",ix,size(xx)
          tt( ix) =			tt( ix) -  flt( ia) * tt( iy)
        end do
      end do 
      xx =	xx + tt
    else
      tt = xx
      do iy= 1, nd  
!		if(associated(aa%mis)){if(aa%mis(iy)){ tt(iy)=0.; next}}
        if (associated(aa%mis)) then
          if (aa%mis(iy)) then
            tt(iy)=0.
            cycle
          end if
        end if
        ip = aa%pch( iy)
        lag => aa%hlx( ip)%lag
        flt => aa%hlx( ip)%flt
        do ia = 1, size( lag)
          ix = iy - lag( ia)
          if ( ix < 1) then
            cycle
          end if
!			if(ix <0 || ix >size(xx)) write(0,*) "2x error",ix,size(xx),size(yy),lag(ia)
          tt( iy) =			tt( iy) -  flt( ia) * tt( ix)
        end do
      end do 
      yy =	yy + tt
    end if 
!    deallocate(lag,flt)
  end subroutine 
  subroutine npolydiv_close()
    deallocate( tt )
  end subroutine 
end module 
