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()
                integer	::i
                allocate (ptch(n1*n2*n3),n_nfilt(nang))
        
                call from_param('scale',lscale,.false.)
                call from_param('destroy',destroy,.false.)
        
                if (lscale)then !scale according to angle
                        call from_param("os",os,1.); call from_param("ds",ds,0.); call from_param("ns",ns,nang)
                        allocate (scale(ns))
                        do i=1,ns
                                scale(i)=os+(i-1)*ds
                        enddo
                endif
        
                ptch=0.
        
                do i=1,n1*n2*n3
                        dip = atan(dips(i))
                        idip = 1.5+(dip-oang)/dang
                        ptch(i)=idip
                        !if(ptch(i)>nang.or.(ptch(i)<0.))write(0,*)dip,idip,ptch(i)
                end do
                where (ptch>float(nang))ptch=nang
                where (ptch<1.         )ptch=1
                n_nfilt=nfilt*3
                call nallocate(n_filt,n_nfilt,ptch)
        end subroutine
        !######################################
        subroutine ns_lop()
                real, dimension(2)      :: u,v
                do iang=1,nang

                        do ifilt=1,nfilt
                                n_filt%hlx(iang)%lag(         ifilt) = ifilt			! early lags (first row)
                                n_filt%hlx(iang)%lag( nfilt  +ifilt) = n1 + ifilt - nfilt 	! late lags  (second row)
                                n_filt%hlx(iang)%lag( nfilt*2+ifilt) = n1 + ifilt 		! late late lags (second row)
                        end do
                        !call nboundn( iang, (/n1,n2/), (/10, 1/), n_filt)                       !deal with boundary issues here|| Use tentn??
                end do
                call allocatehelix(lapfilt, 4)             ! alloc. filter
                eps=0.015
                lapfilt0 = 1. + eps                        ! zero lag coeff.
                lapfilt%lag = (/ 1, n1-1, n1, n1+1/)       ! lag(1)= 1; lag(2)=n1  # one side only
        


                do iang=1,nang
                        dip=oang+(iang-1)*dang
                        write(0,*)'working on angle',iang,dip,degrad(1,dip)
                        !dip=pi/120.
                        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(iang))
                        call wilson_close()
                        n_filt%hlx(iang)%flt=n_filt%hlx(iang)%flt!/n_filt0
                end do	
        
                if (lscale.and.nang==ns)then    !scale according to the angle!!!!!! useful for zeroing specific angles
                        do iang=1,nang
                                n_filt%hlx(iang)%flt=n_filt%hlx(iang)%flt*scale(iang)
                        end do
                endif
        
                if (destroy)then
                        if (verb>5) write(0,*)'I am destroying all events!!!!!!!'
                        call nhelicon_init(n_filt)
                        stat= filt_nhconest(.false.,.false.,infile,outfile)
                        call nhelicon_close()
                else
                        if (verb>5) write(0,*)'I am precondition all events!!!!!!!'
                        call npolydiv_init(n1*n2*n3,n_filt)
                        stat = filt_npolydiv(.false.,.false.,infile,outfile)
                        call npolydiv_close()
                endif
                if (verb>5) write(0,*)'I am writing the output'
                call srite("outfile", outfile, 4*n1*n2*n3)
        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
