program dipit
  use sep
  use wilson
  use helix
use steering
use polydiv
  use nhelicon
  use autocorr
  use helicon
  use npolydiv
  implicit none
  integer :: n1,n2,na,i,ia,ip,ns,idip,loc,is
  real :: a0, lap0,eps,dip,s0,os,scale,ds,dipa
 real :: u(2),m,p,s
logical :: notch
  type(filter) :: filt,lap,ss,bb,st
  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 from_param("notch",notch,.false.)
  call to_history("n1",n1)
  call to_history("n2",n2)
  call to_history("n3",ns+1)
  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

  call steer_init(.99,n1)

dips=1.
write(0,*) "DIE 1",dips(1)
call allocatehelix(st,5)
 call steering2(-dips(1),st) 
write(0,*) "DIE 1",st%flt,st%lag
   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
pch=11
   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=.0016
    lap%lag = (/1,n1-1,n1,n1+1 /) ! lag(1)= 1; lag(2)=n1 # one side only
! probably done already in allocation.


    do is=1,ns 
     s=os+ds*(is-1)
     if(notch) then
       lap0=1.+eps+s
     else 
       lap0=1+eps+s*4
     end if 
     do i=1,90
     dip=-90+(i-1)*2.
     dipa=atan(1.)*dip/45.
     u(1)=cos(dipa)
     u(2)=sin(dipa)

                        m=(u(1)-u(2))/2.;
                        p=(u(1)+u(2))/2.;
                        lap%flt = 0.;
      if(.not. notch ) then
                   lap%flt = (/2*m*p-s,-p*p,-2*m*p+s,-m*m/)
      else
                   lap%flt = (/2*m*p,-p*p,-2*m*p,-m*m/)
     end if


if(i==11) then
     call wilson_init( 10 * n1 )
call wilson_factor( 100, lap0, lap, a0, nfilt2%hlx(i))
    nfilt2%hlx(i)%flt= nfilt2%hlx(i)%flt!/a0
    call wilson_close()
end if
     
    
!     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
call helicon_init(st)
    na=helicon_lop(.false.,.false.,sv,out);
    na=helicon_lop(.true.,.false.,in,out);
 
   call srite("out",in,size(out)*4)

end program
