program integrate

! I_H    : High_F part of image
! I_L    : Low-F  part of image
! P_H(.) : X-Dip operator based on I_H
!   D    = I_L + I_H
! J(m)   = |I_L - m| + e* P_H(I_L) ~0

! Author:  Mohammad Maysami
!          SEP, Stanford University
!          Jul 09



use sep
use sep3d_struct_mod
use helderiv
use identity_mod
use find_xz_dip_mod

use solver_reg_l2_mod
use solver_smp_l2_mod
use solver_prc_l2_mod
use cgstep_mod

  implicit none
  type(sep3d)       :: data,aux
  character(len=20) :: method
  real, pointer, dimension(:,:,:)   :: dat3,aux3,dip3,xres3
  real, pointer, dimension(:)       :: dat1,aux1,mod1,dip1,xres1
  real, allocatable, dimension(:,:) :: gmov,rmov
  real    :: icycle,eps,epsh,eps1,eps2
  integer :: n(3),nwind(3),ndim,nt,n3,i3,niter,na
  integer :: maxsize,stat,ierr1,ierr2
  logical :: verb,verb2,v3d,lo_aux


  call sep_init(SOURCE)
  call sep_begin_prog()

  method=" "
  call from_param("verb",verb,.true.) 
  call from_param("lo_aux",lo_aux,.true.) 
  call from_param("niter" ,niter,50)
  call from_param("eps"   ,eps,0.05)
  call from_param("na"    ,na,8)
  call from_param("epsh"  ,epsh,0.001)
  call init_sep3d("in"    ,data,"INPUT")
  call init_sep3d("aux"   ,aux ,"INPUT")
  call from_param("method",method,"xsergey")
  write(0,*) "method=",method
  write(0,*) "low_aux=",lo_aux


  ndim=sep3d_ndims(data)
  n=1
  n(1:ndim)=data%n(1:ndim)
  nt=product(n)

  ! =================================================================
  ! Set smoothing factors for PUCK
  nwind(1)=min(n(1),5); 
  nwind(2)=min(n(2),5); 
  nwind(3)=min(n(3),5);       
  call from_param("nwind",nwind,nwind); 

  ! Set smoothing factors for SERGEY
  call from_param("niter",niter,30) 
  call from_param("eps1",eps1,0.00002)
  call from_param("eps2",eps2,0.00002)

  ! Init Helderiv
  call helderiv_init ( n(1),na,epsh)
  if (nt/=product(aux%n(1:sep3d_ndims(aux)))) &
       call seperr("Input and auxilary needs to be the same size")


  !  Allocation
  allocate(dat3(n(1),n(2),n(3)),aux3(n(1),n(2),n(3)),stat=ierr1) ! 3D Alloc
  allocate(dat1(nt),mod1(nt),aux1(nt),stat=ierr2)     ! 1D Alloc
  allocate(gmov(nt,niter),rmov(2*nt,niter))         ! 2D Alloc
  if ((ierr1.ne.0) .or. (ierr2.ne.0)) call seperr("trouble allocating data")



  ! ======= Read model and auxiliary data
  call sreed("aux",aux1,size(aux1)*4)
  call sreed("in",dat1,size(dat1)*4)
  ! Normalize : aux Done in xpuck_init     
  !      aux1=aux1/sum(abs(aux1))
  !    ! dat1=dat1/sum(abs(dat1))


     
  ! =================================================================

     ! ----- FAST Puck Method -----
     call xsergey_dip_init(data,aux1,niter,eps1,eps2)
     if (lo_aux) then 
        call solver_reg_l2(m=mod1,d=dat1, Fop=helderiv_lop, Aop=xsergey_dip_lop, stepper=cgstep, &
             nAop=nt,niter=niter,eps=eps,m0=dat1,rmov=rmov,gmov=gmov,verb=verb)
     else
        call solver_reg_l2(m=mod1,d=dat1, Fop=I_lop, Aop=xsergey_dip_lop, stepper=cgstep, &
             nAop=nt,niter=niter,eps=eps,m0=dat1,rmov=rmov,gmov=gmov,verb=verb)

     endif
     
  call srite("out",mod1,size(mod1)*4)

  ! rmov=[rd,rm]   gmov=gradient
  call sep3d_write_description("rmov",data)
  call sep3d_write_description("gmov",data)
  call to_history("n2",2*n(2),"rmov")
  call to_history("n3",niter,"rmov")
  call to_history("n3",niter,"gmov")
  call srite("gmov",gmov,size(gmov)*4)  
  call srite("rmov",rmov,size(rmov)*4)  


  ! =================================================================
  !  Write out XDip on estimated model
  allocate(dip1(nt),xres1(nt))
  call xsergey_dip_init(data,aux1,niter,eps1,eps2)
  stat=xsergey_dip_lop(.false.,.false.,dat1,xres1)
  call get_dip_aux(dip1)

  ! aux3=reshape(aux1,(/n(1),n(2),n(3)/))
  ! dat3= reshape(dat1,(/n(1),n(2),n(3)/))
  ! if (.not. find_puck_dip_fast(aux3,nwind,dip3,coh3,res3)) &
  !     call seperr("trouble finding Puch dip using fast method")
  !   call wavekill_3p(.false.,.false.,1., dip3, dat3, xres3)
  
  !  stat=find_sergey_fast(n,aux1,20,0.05,0.05,dip,.false.,.false.)
  call sep3d_write_description("dip",data)
  call sep3d_write_description("xres",data)
  call srite("dip",dip1,size(dip1)*4)  
  call srite("xres",xres1,size(xres1)*4)  
  deallocate(dip1,xres1) 




  deallocate(aux3,aux1,mod1,dat3,dat1)
  deallocate(rmov,gmov)

  call sep_end_prog()
end program 
