program integrate_org

! 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 xgrad_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     :: dat3(:,:,:),mod3(:,:,:),aux3(:,:,:)
  real, pointer     :: dat1(:),mod1(:),aux1(:)
  real, pointer     :: xres3(:,:,:),dip3(:,:,:),coh(:,:,:),res(:,:,:)
  real, pointer     :: xres1(:),dip(:)
  real    :: eps,icycle,epsh
  integer :: n(3),ndim,nt,nwind(3),niter,n3,i3,na
  integer :: maxsize,ierr1,ierr2,stat
  logical :: verb,verb2,v3d


  call sep_init(SOURCE)
  call sep_begin_prog()

  method=" "
  call from_param("verb",verb,.true.) 
  call from_param("v3d",v3d  ,.false.);  
  call from_param('max_size',maxsize,100); 
  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,"xpuck")
  maxsize=maxsize*250000
  write(0,*) "method=",method


  ! === Set smoothing factors
  ndim=sep3d_ndims(data)
  if (ndim>3) then; 
     call erexit("can only calculate the dip on <= 3D data"); 
  end if
  n=1
  n(1:ndim)=data%n(1:ndim)
  nwind(1)=min(n(1),5)
  nwind(2)=min(n(2),5)
  nwind(3)=min(n(3),5)
  call from_param("nwind",nwind,nwind)
  call helderiv_init ( n(1),na,epsh)
  
  ! === Get dimensions 2D,3D , Check mod~aux size
  icycle=.0001
  verb2=.false.
  if (.not. v3d) then; 
     n3=n(3); n(3)=1; 
     if (n3>1 .and. verb) then; 
        verb=.false.; 
        verb2=.true.; 
     end if; 
  else; 
     n3=1; 
  end if
  nt=product(n)
  if (nt/=product(aux%n(1:sep3d_ndims(aux)))) &
       call seperr("Input and auxilary needs to be the same size")








  ! --------- XLINE loop ---------
  do i3=1,n3
     ! Print progress
     if (verb2 .and. real(i3)/real(n3) > icycle ) then; 
        write(0,*) "working on ",nint(real(i3)/real(n3)*100.),"%"; icycle=icycle+.05; 
     end if !3d

     ! === Check Max memory size
     if (nt>maxsize/4) then; 
        write(0,*) "max_size needs to be at least ",nt*16/1000000+1; 
        call seperr("maxsize is too small "); 
     end if

     ! === Allocate in first try
     if (i3.eq.1) &
          allocate(dat3(n(1),n(2),n(3)),aux3(n(1),n(2),n(3)),mod3(n(1),n(2),n(3)),stat=ierr1)
          allocate(dat1(nt),mod1(nt),aux1(nt),stat=ierr2)
          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",dat3,size(dat3)*4)
     aux3=reshape(aux1,(/n(1),n(2),n(3)/))
     dat1= reshape(dat3,(/n(1)*n(2)*n(3)/))
     ! Normalize : Done in xpuck_init     
     !      aux1=aux1/sum(abs(aux1))
     !    ! dat3=dat3/sum(abs(dat3))


     ! === Write out Xpuck on estimated model
     allocate(coh(n(1),n(2),n(3)),dip3(n(1),n(2),n(3)),xres3(n(1),n(2),n(3)),res(n(1),n(2),n(3)))
     if (.not. find_puck_dip_fast(aux3,nwind,dip3,coh,res)) &
          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("coh",data)
     call sep3d_write_description("res",data)
     call sep3d_write_description("dip",data)
     call sep3d_write_description("xres",data)
     call srite("xres",xres3,size(xres3)*4)  
     call srite("coh",coh,size(coh)*4)
     call srite("res",res,size(res)*4)        
     call srite("dip",dip3,size(dip3)*4)  
     !  call srite("dip",dip,size(dip)*4)  
     deallocate(coh,dip3,xres3,res) 
     

     

     if (method(1:5)=='xgrad') then 
        ! ----- XGRAD Method -----
        call xgrad_init(data,aux1)	
        call solver_reg_l2(m=mod1,d=dat1, Fop=I_lop, Aop=xgrad_lop, stepper=cgstep, &
             nAop=nt,niter=niter,eps=eps,m0=dat1,verb=verb)
     
     else
        ! ----- FAST Puck Method -----
        call xpuck_dip_init(data,aux1,nwind)	
        call solver_reg_l2(m=mod1,d=dat1, Fop=helderiv_lop, Aop=xpuck_dip_lop, stepper=cgstep, &
             nAop=nt,niter=niter,eps=eps,m0=dat1,verb=verb)
        mod3= reshape(mod1,(/n(1),n(2),n(3)/))     
        
     end if
        



     call srite("out",mod1,size(mod1)*4)
    
  end do 



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


  call sep_end_prog()
end program 
