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
  type(filter) :: filt,lap,ss,bb
  type(nfilter) :: nfilt,nfilt2
  real,allocatable :: in(:),out(:),sv(:),dips(:)
  real :: s
   
   integer ,allocatable :: pch(:),n_s(:)
  
  call sep_init();
  call from_param("os",os)
  call from_param("ds",ds)
  call from_param("ns",ns)
   na=7
 
 n1=100; n2=100
   allocate(pch(n1*n2),n_s(180))
    
   n_s=na*2
   call nallocate(nfilt,n_s,pch)
   call nallocate(nfilt2,n_s,pch)
   do i=1,180
     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,180
     dip=-90+(i-1)
     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.



call sep_put_data_axis_par("out",1,14,0.,1.,"coef")
call sep_put_data_axis_par("out",2,180,-90.,1.,"dip")
call sep_put_data_axis_par("out",3,ns,os,ds,"sm")

    do is=1,ns
    s=os+ds*(is-1)
     lap0=1+eps+s
     do i=1,180
     dip=-90+(i-1)
     
write(0,*) "working on",dip,is,ns,s
     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
    ! lap%flt=lap%flt*(os+ds*(is-1))
     call wilson_init( 10 * n1 )
    call wilson_factor( 100, lap0, lap, a0, nfilt2%hlx(i))
    call wilson_close()
     call srite("out",nfilt2%hlx(i)%flt,14*4)
  end do
  end do

end program
