module killwave_ns
!modified from lapfact
!Gboyega Ayeni
!2009/17/01
        use killwave_params
        use npolydiv_mod
        use nhelicon_mod
        use wilson_mod
        use nbound_mod

        implicit none

        logical                         :: lscale
        real, allocatable               :: scale(:)
        real                            :: ds, os
        integer                         :: ns

contains 
!######################################
        subroutine  ns_init(dip_in,eps_in)
                integer         :: ifilt
                real            :: dip_in, eps_in
                real            :: u(2),v(2)
                
                dip=dip_in
                eps=eps_in    
                ptch=1
                n_nfilt=nfilt*3
                call nallocate(n_filt,n_nfilt,ptch)
                do ifilt=1,nfilt
                        n_filt%hlx(1)%lag(         ifilt) =        ifilt                  ! early lags (first row)
                        n_filt%hlx(1)%lag( nfilt  +ifilt) =   n1 + ifilt - nfilt          ! late lags  (second row)
                        n_filt%hlx(1)%lag( nfilt*2+ifilt) =   n1 + ifilt                  ! late late lags (second row)
                end do
                n_filt%hlx(1)%flt=0.
                call allocatehelix(lapfilt, 4)                                            ! alloc. filter
                lapfilt0 = 1. + eps                                                       ! zero lag coeff.
                lapfilt%lag = (/ 1, n1-1, n1, n1+1/)                                      ! lags of one side of the autocorellation

                if(verb > 5) write(0,*)'working on angle',degrad(1,dip)
                u(1)= cos(dip); u(2)=-sin(dip); 

                m=(u(1)-u(2))/2.; 
                p=(u(1)+u(2))/2.;
                lapfilt%flt = 0.;
                lapfilt%flt = (/2*m*p,-p*p,-2*m*p,-m*m/) 
                call wilson_init( 10 * n1 )
                call wilson_factor( 100, lapfilt0, lapfilt, n_filt0, n_filt%hlx(1))
                call wilson_close()
        end subroutine
!######################################
        subroutine ns_filt_lop(adj,add,xx,yy)
                real            :: xx(:),yy(:)
                logical         :: add,adj
                if (destroy)then
                        if (verb>5) write(0,*)'I am destroying angle: ', dip
                        call nhelicon_init(n_filt) 
                        stat= filt_nhconest(.false.,.false.,xx,yy)
                        call nhelicon_close()
                else
                        if (verb>5) write(0,*)'I am keeping only angle: ', dip
                        call npolydiv_init(n1*n2,n_filt)
                        stat = filt_npolydiv(.false.,.false.,xx,yy)
                        call npolydiv_close()
                endif
        end subroutine
!######################################
        function filt_npolydiv( adj, add, xx, yy) result(stat)
                integer                   :: stat,stat1,stat2 
                logical,intent(in)        :: adj,add 
                real,dimension(:)         :: xx,yy 
                real, dimension(size(xx)) :: tmp
                tmp=0.

                if ( adj) then
                        stat= npolydiv_lop(.false.,.false., yy,tmp)     ! tmp = grad y
                        stat= npolydiv_lop(.true. ,.false., xx,tmp)     ! x = x + grad' tmp
                else
                        stat= npolydiv_lop(.false.,.false., xx,tmp)     ! tmp = grad x
                        stat= npolydiv_lop(.true. ,.false., yy,tmp)     ! y = y + grad' tmp
                endif
                stat=stat*stat2
        end function 
!######################################
        function filt_nhconest( adj, add, xx, yy) result(stat)
                integer                   :: stat,stat1,stat2 
                logical,intent(in)        :: adj,add 
                real,dimension(:)         :: xx,yy 
                real, dimension(size(xx)) :: tmp
                tmp=0.
                if ( adj) then
                        stat= nhelicon_lop(.false.,.false., yy,tmp)     ! tmp = grad y
                        stat= nhelicon_lop(.true. ,.false., xx,tmp)     ! x = x + grad' tmp
                else
                        stat= nhelicon_lop(.false.,.false., xx,tmp)     ! tmp = grad x
                        stat= nhelicon_lop(.true. ,.false., yy,tmp)     ! y = y + grad' tmp
                endif
                stat=stat*stat2
        end function 
end module
