module filt_fact_mod  
        use  wilson_mod

        implicit none
        type( filter),   private  :: aa,aa1 
        contains
!######################################
        subroutine filt_fact_init (n1,n2,n3,na,eps,filt0,filt,hfilt )
                integer, intent (in) :: n1, n2, n3, na
                real,    intent (in) :: eps,filt0
		real,    intent (in) ::     filt(:)
	        type( filter)	     :: hfilt
	        aa = filt_fact2( eps, n1, na, filt0, filt)
		hfilt=aa
                call helicon_init( aa)
                call polydiv_init(n1*n2*n3, aa)
        end subroutine
!######################################
        function filt_fact_lop ( adj, add, pp, qq) result(stat)
        integer            :: stat 
        logical,intent(in) :: adj,add 
        real,dimension(:)  :: pp,qq 
                call adjnull (adj,add,pp,qq )
                call filt_fact_lop2(adj,add,pp,qq )
                stat=0
        end function 
!######################################
        subroutine filt_fact_lop2(adj,add,pp,qq)
        logical,intent(in) :: adj,add 
	        real, dimension (:)  :: pp 
        	real, dimension (:)  :: qq 
	        logical              :: opt
        	integer stat1
                	stat1 = helicon_lop( adj, .false., pp, qq,opt)
        end subroutine

!######################################
        function filt_conv_lop3( 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= helicon_lop(.false.,.false., yy,tmp)     ! tmp = grad y
                        stat= helicon_lop(.true. ,.false., xx,tmp)     ! x = x + grad' tmp
                else
                        stat= helicon_lop(.false.,.false., xx,tmp)     ! tmp = grad x
                        stat= helicon_lop(.true. ,.false., yy,tmp)     ! y = y + grad' tmp
                endif
                stat=stat*stat2
        end function 
!######################################
        function filt_decon_lop3( 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= polydiv_lop(.false.,.false., yy,tmp)     ! tmp = grad y
                        stat= polydiv_lop(.true. ,.false., xx,tmp)     ! x = x + grad' tmp
                else
                        stat= polydiv_lop(.false.,.false., xx,tmp)     ! tmp = grad x
                        stat= polydiv_lop(.true. ,.false., yy,tmp)     ! y = y + grad' tmp
                endif
                stat=stat*stat2
        end function 
!######################################
        subroutine  helderiv_close()
                call deallocatehelix( aa)
        end subroutine 
!######################################
        function filt_fact2(eps, n1, na, filt0, filt, verb)   result (aa)
                type( filter)                       :: aa, lap
                real,dimension(4),      intent( in) :: filt
                real,                   intent( in) :: eps, filt0
                integer,                intent( in) :: n1, na
                integer                             :: i
                real                                :: a0, lap0
                logical                             :: verb 
                optional                            :: verb    
                        call allocatehelix( lap, 4)             ! alloc. filter
                        lap0 = filt0 + eps                      ! zero lag coeff.
                        lap%lag = (/ 1, n1-1, n1, n1+1/)        ! lag(1)= 1; lag(2)=n1  # one side only
                        lap%flt = filt				! put filter on helix
                        call allocatehelix( aa, 3*na)           ! laplacian derivative
                        aa%flt = 0.                             ! probably done already in allocation.
                        do i = 1, na 
                                aa%lag(      i) =      i        ! early lags (first row)
                                aa%lag(   na+i) = n1 + i - na   ! late lags (second row)
                                aa%lag( 2*na+i) = n1 + i        ! late late lags (second row)
                        end do 
                        call wilson_init( 10 * n1 )
                        call wilson_factor( 100, lap0, lap, a0, aa,.false.)  
                        call wilson_close()
                        call deallocatehelix( lap)
                        if(present(verb))then; if (verb)then
                                write(0,*)'+++++++++++++++filtr'
                                write(0,*)aa%lag
                                write(0,*)aa%flt
                                write(0,*)'+++++++++++++++'     
                        endif; endif
  	end function 


end module 

