program dipit
  use sep
  use wilson
  use helix
  use nhelicon
  use autocorr
  use npolydiv
  implicit none
  integer :: n1,n2,na,i,ia,ip,ns,idip,loc,is
  real :: a0, lap0,eps,dip,m,p,s0,os,scale,ds,s
  type(filter) :: filt,lap,ss,bb
  type(nfilter) :: nfilt,nfilt2
  real,allocatable :: in(:),out(:),sv(:),dips(:)
   
   integer ,allocatable :: pch(:),n_s(:)
  
  call sep_init();
  call from_history("n1",n1)
  call from_history("n2",n2)
  call from_param("dip",dip);
  call from_param("os",os)
  call from_param("ds",ds)
  call from_param("ns",ns)
  call to_history("n1",n1)
  call to_history("n2",n2)
  call to_history("n3",ns)
  allocate(in(n1*n2),out(n1*n2))
  allocate(dips(n1*n2),sv(n1*n2))
  call sreed("in",sv,size(sv)*4)
 ! sv=0
 ! sv(n1/2+n2*n1/2)=5.

  call sreed("dip",dips,size(dips)*4)
 dips=-dips
!dips=0
  m=cos(dip*atan(1.)/45.)-sin(dip*atan(1.)/45.)
  p=cos(dip*atan(1.)/45.)+sin(dip*atan(1.)/45.)
  m=m/2.
  p=p/2.
   na=7

   allocate(pch(n1*n2),n_s(90))
   do i=1,n1*n2
     dip=45./atan(1.)*atan(dips(i))
     idip=(dip+90)/2.+1
     pch(i)=max(1,min(90,idip))
     !call srite("out",dip,4)
    !write(0,*) "CHECK IT",i,idip,dip,dips(i)
    
   end do
   n_s=na*2
   call nallocate(nfilt,n_s,pch)
   call nallocate(nfilt2,n_s,pch)
   do i=1,90
     do ia=1,na
       nfilt%hlx(i)%lag( ia ) = ia ! early lags (first row)
       nfilt%hlx(i)%lag( na+ia) = n1 + ia - na+1 ! late lags (second row)
       nfilt2%hlx(i)%lag( ia ) = ia ! early lags (first row)
       nfilt2%hlx(i)%lag( na+ia) = n1 + ia - na+1 ! late lags (second row)
     end do
   end do
write(0,*) "WHAT THE",nfilt%hlx(1)%lag

  call allocatehelix( lap, 4) ! laplacian filter
  eps=.01
    lap0 = 1. + eps ! zero lag coeff.
    lap%lag = (/1,n1-1,n1,n1+1 /) ! lag(1)= 1; lag(2)=n1 # one side only
  do i=1,90
     dip=-90+(i-1)*2.
     m=cos(dip*atan(1.)/45.)-sin(dip*atan(1.)/45.)
     p=cos(dip*atan(1.)/45.)+sin(dip*atan(1.)/45.)
     m=m/2.
     p=p/2.
     lap%flt = (/2*m*p,-p*p,-2*m*p,-m*m/) ! flt(1)=-1; flt(2)=-1
     call wilson_init( 10 * n1 )
    call wilson_factor( 100, lap0, lap, a0, nfilt%hlx(i))
    call wilson_close()
end do
! probably done already in allocation.


    do is=1,ns 
     do i=1,90
;
s=os+ds*(is-1)
     lap0=1+s*4+eps
     dip=-90+(i-1)*2.
     m=cos(dip*atan(1.)/45.)-sin(dip*atan(1.)/45.)
     p=cos(dip*atan(1.)/45.)+sin(dip*atan(1.)/45.)
     m=m/2.
     p=p/2.
     lap%flt = (/2*m*p-s,-p*p,-2*m*p-s,-m*m/) ! flt(1)=-1; flt(2)=-1
     call wilson_init( 10 * n1 )
    call wilson_factor( 100, lap0, lap, a0, nfilt2%hlx(i))
    call wilson_close()
     
    
!     nfilt%hlx(i)%flt=nfilt%hlx(i)%flt/a0*.95
!call srite("junk.H",nfilt%hlx(i)%flt,60)
  end do
   call nhelicon_init(nfilt)
      call npolydiv_init(n1*n2,nfilt2)
      na=nhelicon_lop(.false.,.false.,sv,out)
     na=nhelicon_lop(.true.,.false.,in,out)
    na=npolydiv_lop(.false.,.false.,in,out)
     na=npolydiv_lop(.true.,.false.,in,out)
       call npolydiv_close()
    call nhelicon_close()
      call srite("out",in,size(out)*4)
  end do
  do i=1,n1*n2
     in(i)=nfilt%pch(i)
   end do
   call srite("out",in,size(out)*4)

end program
