program misso
 use mask1
 use sep
 use steering
 use npolydiv
 use solver_reg_mod
 use solver_prc_mod
 use chain_me_mod
 use cgstep_mod
 implicit none

 integer :: n1,n2,niter
 real :: eps
 real,pointer :: m(:),d(:),p(:)
 integer ,pointer:: mski(:)
 real,pointer :: filts(:,:,:),dips(:)
 logical,pointer :: msk(:)
 integer :: nc,ndip,ns
 integer :: i,j,k
 real :: odip,ddip,os,ds,dp
 integer,pointer :: nsz(:)
 character(len=128) :: label
 type(nfilter),allocatable :: nfilt(:)
 real, allocatable :: r(:,:)

 call sep_init()
 call from_param("niter",niter,10)
 call from_param("eps",eps,.1)
 call from_history("n1",n1)
 call from_history("n2",n2)
 allocate(m(n1*n2),d(n1*n2),mski(n1*n2),msk(n1*n2),dips(n1*n2),p(n1*n2))
 call sreed("in",d,size(d)*4)
 call sreed("mask",mski,size(mski)*4)
 call sreed("dips",dips,size(dips)*4)
 msk=.false.
 where(mski==1) msk=.true.
 call sep_get_data_axis_par("filt",1,nc,os,ds,label)
 call sep_get_data_axis_par("filt",2,ndip,odip,ddip,label)
 call sep_get_data_axis_par("filt",3,ns,os,ds,label)
 do i=1,size(dips)
    ds=atan(-dips(i))*45./atan(1.)
    mski(i)=max(1,min(ndip,nint((ds-odip)/ddip+1.)));
 end do
 allocate(filts(nc,ndip,ns),nsz(ndip))
 call sreed("filt",filts,size(filts)*4)
 allocate(nfilt(ns+1))
 nsz=nc
 allocate(r(n1*n2,niter))
 do k=1,ns
    call nallocate(nfilt(k),nsz,mski)
    do j=1,ndip
       do i=1,nc/2
          nfilt(k)%hlx(j)%lag(i)=i
          nfilt(k)%hlx(j)%lag(i+nc/2)=n1+i-nc/2+1
       end do
       nfilt(k)%hlx(j)%flt=filts(:,j,k)
    end do
 end do
 nsz=5
 call steer_init(.99,n1)
 call nallocate(nfilt(ns+1),nsz,mski)
 do j=1,ndip
    dp=-atan(atan(1.)/45.*(odip+ddip*(j-1)))
    call steering2(dp,nfilt(ns+1)%hlx(j))
 end do

 call to_history("n3",ns+1)
 call nhelicon_init(nfilt(ns))
 do i=1,ns+1
    if(i/=ns+1) then
      call chain_me_init(n1*n2,msk,nfilt(i),nfilt(ns))
    else 
      call chain_me_init(n1*n2,msk,nfilt(i),nfilt(ns+1))
    end if
    niter=3
    call solver_reg(p,d,chain_me_lop,chain_me_reg,cgstep,n1*n2,niter,eps,verb=.true.,rmov=r);
!    call srite("rmov.H",r,size(r)*4)
    call cgstep_close()
    call to_model(p,m)
    call chain_me_close()
    call srite("out",m,size(m)*4)
 end do
end program
