! This is Sergey Fomels dip finder and Jon Claerbout Puck
module find_xz_dip_mod 

  use OW_type
  use OW_parms_areal
  use tri_smooth_mod
  use grad3_dip_mod !good
  use allpass_3d    !good
  use allder_3d
  use divide

  implicit none
  integer, private       :: pn1,pn2,pn3, pniter
  real, private          :: pd1,pd2,pd3, peps1,peps2
  real, private, pointer :: paux1(:),pdip1(:)
  integer, private       :: pnwind(3)


contains



  subroutine xsergey_dip_init(aux1_in,niter,eps1,eps2)
   ! Initialize Xpuck function
    ! aux_s types 3D argument passes the dimensions 
    ! aux_in passes the soft data

    real, dimension(:), target :: aux1_in
    integer                    :: smooth(3),niter
    real                       :: eps1,eps2

    ! Set Sergeys iter and eps (30,0.00002,0.00002)
    pniter=niter
    peps1=eps1
    peps2=eps2

    ! Set size of aux
    pn1=img%xm%n
    pn2=img%z%n
    pd1=img%xm%d
    pd2=img%z%d
! change for 3D !!!!!!!!!!!!!!!!!!!!!
    pn3=1
    pd3=1.
    
    allocate(pdip1(pn1*pn2*pn3))

    ! Assumes 3D array as input
    nullify(paux1)
    paux1 => aux1_in
    paux1 = paux1/sum(abs(paux1))
  end subroutine

  ! --------------------------------------------------------------
  function xsergey_dip_lop(adj,add,mm,dd) result(stat)
    ! Find 3d dip(p) and coh field by sergey method
    logical, intent(in)         :: adj,add
    real                        :: mm(:),dd(:)
    real, allocatable, dimension(:,:,:) :: aux3,field3,xres3,dip3
    integer                     :: nwind(3), ierr,stat

    stat=1
    allocate(dip3(pn1,pn2,pn3))
    allocate(field3(pn1,pn2,pn3),xres3(pn1,pn2,pn3),aux3(pn1,pn2,pn3))


    ! Get SERGEYs slopes from auxiliary 
    aux3 = reshape(paux1,(/pn1,pn2,pn3/))
    ierr=find_sergey_fast((/pn1,pn2,pn3/),paux1,pniter,peps1,peps2,pdip1,.false.,.false.)

    ! Reshape 3D <=> 1D
    dip3   = reshape(pdip1,(/pn1,pn2,pn3/))
    field3 = reshape(mm,(/pn1,pn2,pn3/))
    xres3  = reshape(dd,(/pn1,pn2,pn3/))

    ! Apply dip to model
    call wavekill_3p(adj,add,1., dip3, field3, xres3)
    mm= reshape(field3,(/pn1*pn2*pn3/))
    dd= reshape(xres3,(/pn1*pn2*pn3/))

    if (exist_file("dip3")) then
      call to_history("n3",nz,"dip3") ; call to_history("d3",img%z%d,"dip3") ; call to_history("o3",img%z%o,"dip3")
      call to_history("n2",1,"dip3") ; call to_history("d2",1.,"dip3") ; call to_history("o2",0.,"dip3")
      call to_history("n1",img%xm%n,"dip3"); call to_history("d1",img%xm%d,"dip3"); call to_history("o1",img%xm%o,"dip3")  
      call to_history("n4",1,"dip3"); call to_history("d4",1,"dip3"); call to_history("o4",0,"dip3")
      call to_history("n5",1,"dip3"); call to_history("d5",1,"dip3"); call to_history("o5",0,"dip3")  
      call to_history("esize",4,"dip3")
      call srite("dip3",dip3,4*img%z%n*img%xm%n)
    end if
    deallocate(aux3,field3,xres3,dip3)
    stat=0
  end function
  ! ==============================================================
  ! --------------------------------------------------------------
  ! ==============================================================


  subroutine wavekill( adj,add, aa, pp, uu, vv )
    ! Apply (aa.d/dx + pp.d/dt) uu = vv
    logical, intent(in)         :: adj,add
    real    :: aa, pp, s11, s12, s21, s22, uu(:,:,:), vv(:,:,:)
    integer :: i1,i2,n1,n2,n3,i3

    ! Adjoint Null
    if (.not. add) then
       if (adj) then
          uu=0.
       else
          vv=0.
       endif
    endif

    ! Compute and print Sij
    s11 = -aa-pp; s12 =  aa-pp; 
    s21 = -aa+pp; s22 =  aa+pp
    
    ! Get dimensions
    n1=size(uu,1); 
    n2=size(uu,2); 
    n3=size(uu,3)
    
    ! Loop over axis and apply filter
    if (adj) then
       do i3=1, n3 !xline
          do i1= 1, n1-1 !time 
             do i2= 1, n2-1  !inline 
                uu(i1  ,i2  ,i3) = uu(i1  ,i2  ,i3)+ s11*vv(i1,i2,i3)
                uu(i1  ,i2+1,i3) = uu(i1  ,i2+1,i3)+ s12*vv(i1,i2,i3)
                uu(i1+1,i2  ,i3) = uu(i1+1,i2  ,i3)+ s21*vv(i1,i2,i3)
                uu(i1+1,i2+1,i3) = uu(i1+1,i2+1,i3)+ s22*vv(i1,i2,i3)
             end do
          end do
       end do
    else
       do i3=1, n3 !xline
          do i1= 1, n1-1 !time 
             do i2= 1, n2-1  !inline 
                vv(i1,i2,i3) = vv(i1,i2,i3)+&
                     uu(i1  ,i2,i3)*s11 + uu(i1  ,i2+1,i3)*s12 +&
                     uu(i1+1,i2,i3)*s21 + uu(i1+1,i2+1,i3)*s22
             end do
          end do
       end do
       ! Set boundaries
       vv(n1,:,:)=0. !vv(n1-1,:,:)
       vv(:,n2,:)=0. !vv(:,n2-1,:)
    end if

  end subroutine 

  ! --------------------------------------------------------------
  subroutine wavekill_3p(adj,add,aa, pp, uu, res)
    ! Find res = (aa.d/dx + pp.d/dt) with 3D pp field
    logical, intent(in)         :: adj,add
    integer :: i1,i2,n1,n2,n3,i3
    real    :: aa,uu(:,:,:), res(:,:,:), pp(:,:,:)
    real, dimension(:,:,:), allocatable :: s11, s12, s21, s22


    ! Adjoint Null
    if (.not. add) then
       if (adj) then
          uu=0.
       else
          res=0.
       endif
    endif
   
    ! Get dimensions
    n1=size(uu,1); 
    n2=size(uu,2); 
    n3=size(uu,3)


    
    ! Compute and print Sij
    allocate(s11(n1,n2,n3),s12(n1,n2,n3),s21(n1,n2,n3),s22(n1,n2,n3))
    s11 = -aa-pp; s12 =  aa-pp; 
    s21 = -aa+pp; s22 =  aa+pp

    

    ! Loop over axis and apply filter
    if (adj) then
       do i3=1, n3 !xline
          do i1= 1, n1-1 !time 
             do i2= 1, n2-1  !inline 
                uu(i1  ,i2  ,i3) = uu(i1  ,i2  ,i3)+ s11(i1,i2,i3)*res(i1,i2,i3)
                uu(i1  ,i2+1,i3) = uu(i1  ,i2+1,i3)+ s12(i1,i2,i3)*res(i1,i2,i3)
                uu(i1+1,i2  ,i3) = uu(i1+1,i2  ,i3)+ s21(i1,i2,i3)*res(i1,i2,i3)
                uu(i1+1,i2+1,i3) = uu(i1+1,i2+1,i3)+ s22(i1,i2,i3)*res(i1,i2,i3)
             end do
          end do
       end do
    else
       do i3=1, n3 !xline
          do i1= 1, n1-1 !time 
             do i2= 1, n2-1  !inline 
                res(i1,i2,i3) = res(i1,i2,i3)+&
                     uu(i1  ,i2,i3)*s11(i1,i2,i3) + uu(i1  ,i2+1,i3)*s12(i1,i2,i3) +&
                     uu(i1+1,i2,i3)*s21(i1,i2,i3) + uu(i1+1,i2+1,i3)*s22(i1,i2,i3)
             end do
          end do
       end do
       ! Set boundaries
       res(n1,:,:)=0. !res(n1-1,:,:)
       res(:,n2,:)=0. !res(:,n2-1,:)
    end if

    deallocate(s11,s12,s21,s22)

  end subroutine 


  ! --------------------------------------------------------------
  logical function find_sergey_fast(n,field,niter,eps1,eps2,dip,verb,dave)
        logical, intent(in) :: verb
        real,pointer :: field(:)
        real :: eps1,eps2,mean
        real,pointer ::  dip(:) 
        integer, intent(in) :: niter
        integer :: n(3),ierr,n123,ndim,iter
        real, allocatable  :: u1(:),u2(:),dp(:)
        logical            :: dave
                find_sergey_fast=.false.
                n123=product(n)
                ! write(0,*) "dave",dave
                allocate(u1(n123),u2(n123),dp(n123),stat=ierr)
                if (ierr.ne.0)then; write(0,*) "trouble allocating temp arrays"; return; end if
                call  allpass_3d_init(n(1),n(2),n(3),dip,dave)
                ! write(0,*) "dave",dave
                call   allder_3d_init(n(1),n(2),n(3),dip,dave)
                ! write(0,*) "dave",dave
                dip=0.

                do iter=1,niter 
                        if (verb) write(0,*) "working on non-linear iteration111",iter
                        ! write(0,*) "dave",dave
                        ierr=allpass_3d_lop(.false.,.false.,field,u2)
                        ! write(0,*) "dave",dave
                        ierr=allder_3d_lop(.false.,.false.,field,u1)
                        mean=sqrt(sum(u1*u1))
                        u1=u1/mean
                        u2=u2/mean
                        call divide_apply_tri(n,u2,u1,dp,eps1,eps2)
                        dip=dip+dp
                end do 
                deallocate(dp,u1,u2)
                call allpass_3d_close()
                call allder_3d_close()
                find_sergey_fast=.true.
  end function 
end module 
