module wei_cigs
  !! compute CIGs
  !!        offset domain (weihcig)
  !! ray parameter domain (weipcig)

#ifndef USE_FFTW
use fftw_sep_fake_mod
#endif


  use wei_image_mod
  use wei_wavespace_mod
!  use wei_util
!  use wei_external
  use wei_flags
  use wei_tds
  use wei_attr
  !---------------
  implicit none

  type, public        :: axis
     integer          :: n
     real             :: o,d
     character(len=2) :: l    !! label (i.e. axis name)
  end type axis


  real,               private :: pcigeps

contains
  !----------------------------------------------------------------  
!!$=head1 NAME
!!$
!!$weipcig_init -initialize offset ray parameter imaging
!!$
!!$=head1 SYNOPSIS
!!$
!!$ierr= weihcig_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initializes migration operator whose image (R) domain is p_h
!!$
!!$ Grabs pcig_eps (default=1.0)
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weihcig_init> L<weipcig>
!!$
!!$=cut
!!$

  !----------------------------------------------------------------  
  !! p-domain CIG
  integer function weipcigs_init() result(st)
    integer, dimension(3) :: ni
    logical               :: fftmeasure
    integer               :: fftsetup
    integer               :: er,i
    integer               :: ipx,ipy, jpx,jpy, ihx,ihy
    real                  :: arg
#ifdef DEBUG
    character(len=128), parameter :: unit='weipcig_init'
!    if(run%debug) call in(unit)
#endif



#ifdef DEBUG
!    if(run%debug) call out(unit)
#endif
    st=0
  end function weipcigs_init
!!$=head1 NAME
!!$
!!$weipcig - Perform offset ray parameter imaging 
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weipcig(wfld,iws,ith)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item iws - integer
!!$
!!$  current frequency slice
!!$
!!$=item ith - integer
!!$
!!$  current thread
!!$
!!$=item wfld - complex(:,:,:,:,:)
!!$
!!$  wavefield
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Imaging operation when in the offset ray parameter (p_h) domain
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weipcig_init> L<weihcig>
!!$
!!$=cut
  !----------------------------------------------------------------
  integer function weipcigs(wfld,iws,ith) result(st)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    complex, dimension(:,:,:,:),allocatable :: midin
    complex, dimension(:,:),allocatable :: exp_multx,exp_multy
    integer, intent(in) ::iws,ith
    logical :: adj
    type(sep3d) :: wsep,rsep
		real :: w,px,py,hx,hy
		complex :: shift
		integer :: ipx,ipy,imx,imy,ihy,ihx
#ifdef DEBUG
    character(len=128), parameter :: unit='weipcig'
!    if(run%depar) call in(unit)
#endif

   call grab_current_image_space(rsep)
   call grab_current_wave_space_pad(wsep)
   adj= grab_current_adj()
		allocate(midin(wsep%n(1),wsep%n(2),wsep%n(3),wsep%n(4)),&
 exp_multx(wsep%n(3),rsep%n(3)),exp_multy(wsep%n(4),rsep%n(4)))
		midin=0.0


    call grab_current_omega(w)

exp_multx=cmplx(0.,0.)
		do ipx=1,rsep%n(3)
			px=(ipx-1)*rsep%d(3)+rsep%o(3)
			do ihx=1,wsep%n(3)
				hx=(ihx-1)*wsep%d(3)+wsep%o(3)
        exp_multx(ihx,ipx)=cexp(cmplx(0,-w*px*hx))
			end do
		end do


    if(adj) then !! MIGRATION
		
			midin=wfld(:,:,:,:,ith)
			wfld(:,:,:,:,ith)=cmplx(0.0,0.)
			!! get the p_hx axis done
					do ipx=1,rsep%n(3)
						do ihy=1,wsep%n(4)
							do ihx=1,wsep%n(3)
			do imy=1,wsep%n(2)
				do imx=1,wsep%n(1)
wfld(imx,imy,ipx,ihy,ith)=wfld(imx,imy,ipx,ihy,ith)+midin(imx,imy,ihx,ihy)*exp_multx(ihx,ipx)
							end do
						end do
					end do
				end do
			end do

!	if(run%m_ith==0) write(0,*)"sum wfld2",sum(wfld)

!if(1==3) then
!	if(wsep%n(4) > 1) then
!			do ipy=1,rsep%n(4)
!				py=(ipy-1)*rsep%d(4)+rsep%o(4)
!				do ihy=1,wsep%n(4)
!					hy=(ihy-1)*wsep%d(4)+wsep%o(4)
!					exp_multy(ihy,ipy)=exp(cmplx(0,-w*py*hy))
!				end do
!			end do
!
!			!! get the p_hy axis done -- NEVER TESTED
!			do imy=1,wsep%n(2)
!				do imx=1,wsep%n(1)
!					midin(:,:)=wfld(imx,imy,:,:,ith)
!					wfld(imx,imy,:,:,ith)=0.0
!					do ipy=1,rsep%n(4)
!						do ihx=1,wsep%n(3)
!							do ihy=1,wsep%n(4)
!								wfld(imx,imy,ihx,ipy,ith)=wfld(imx,imy,ihx,ipy,ith)+midin(ihy,ipy)*exp_multy(ihy,ipy)
!							end do
!						end do
!					end do
!				end do
!			end do
!	end if
!end if
!	if(run%m_ith==0) write(0,*)"sum wfld3",sum(wfld)

    else             !! MODELING
					midin=wfld(:,:,:,:,ith)
					wfld(:,:,:,:,ith)=0.0

			!! get the p_hx axis done
					do ipx=1,rsep%n(3)
						do ihy=1,wsep%n(4)
							do ihx=1,wsep%n(3)
			do imy=1,wsep%n(2)
				do imx=1,wsep%n(1)
wfld(imx,imy,ihx,ihy,ith)=wfld(imx,imy,ihx,ihy,ith)+midin(imx,imy,ipx,ihy)*conjg(exp_multx(ihx,ipx))
							end do
						end do
					end do
				end do
			end do

!if(1==3) then
!	if(wsep%n(4) > 1) then
!			do ipy=1,rsep%n(4)
!				py=(ipy-1)*rsep%d(4)+rsep%o(4)
!				do ihy=1,wsep%n(4)
!					hy=(ihy-1)*wsep%d(4)+wsep%o(4)
!					exp_multy(ihy,ipy)=exp(cmplx(0,w*py*hy))
!				end do
!			end do
!
!			!! get the p_hy axis done -- NEVER TESTED
!			do imy=1,wsep%n(2)
!				do imx=1,wsep%n(1)
!					midin(:,:)=wfld(imx,imy,:,:,ith)
!					wfld(imx,imy,:,:,ith)=0.0
!					do ipy=1,rsep%n(4)
!						do ihx=1,wsep%n(3)
!							do ihy=1,wsep%n(4)
!								wfld(imx,imy,ihx,ipy,ith)=wfld(imx,imy,ihx,py,ith)+midin(ihx,ipy)*exp_multy(ihy,ipy)
!							end do
!						end do
!					end do
!				end do
!			end do
!		end if
!end if

    end if
 
!  if(run%m_ith==0) write(0,*) "CHECK OUT",sum(cabs(wfld)),sum(exp_multx)
		deallocate(midin,exp_multx,exp_multy)

#ifdef DEBUG
!    if(run%depar) call out(unit)
#endif
    st=0
  end function weipcigs
  !----------------------------------------------------------------  

  !----------------------------------------------------------------
end module wei_cigs

