! Copyright 1999 Stanford University
!All rights reserved
!
!Author: Sergey Fomel and Robert Clapp
module npolydiv2
! 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 npolydiv2_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 npolydiv2_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 npolydiv2_lop2(adj,xx,yy )
    stat=0
  end function 
  subroutine npolydiv2_lop2(adj,xx,yy)
    logical,intent(in) :: adj 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    real :: f
    integer                        :: ia, ix, iy, ip,ilag
    tt = 0.
    if ( adj) then
     if(associated(aa%mis)) then
       where(.not. aa%mis) tt = yy
     else 
       tt=yy
     end if 
      do iy= nd, 1, -1  
        ip = aa%pch( iy)
        flt => aa%hlx( ip)%flt
        do ia = 1, size( flt)
          f=flt(ia)
          ix=iy-aa%hlx(ip)%lag(ia)
          if(associated(aa%mis)) then
              if(aa%mis(iy)) f=0
          end if
          if ( ix < 1) cycle
          tt( ix) =			tt( ix) -  f * tt( iy)
        end do
      end do 
      xx =	xx + tt

    else
     if(associated(aa%mis)) then
       where(.not. aa%mis) tt = xx
     else
       tt=xx
     end if 
      do iy= 1, nd  
        ip = aa%pch( iy)
        flt => aa%hlx( ip)%flt
        do ia = 1, size( flt)
          ix=iy-aa%hlx(ip)%lag(ia)
          f=flt(ia)
          if(associated(aa%mis)) then
              if(aa%mis(iy)) f=0
          end if
!              if(f/=0 .and. ia==1) write(0,*) "IN APPY",ia,iy
          if ( ix < 1) cycle
          tt( iy) =			tt( iy) -  f * tt( ix)
        end do
      end do 
      yy =	yy + tt
    end if 
  end subroutine 
  subroutine npolydiv2_close()
    deallocate( tt )
  end subroutine 
end module 
