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

#ifndef USE_FFTW
use fftw_sep_fake_mod
#endif
use wei_image_mod
  implicit none
  real,               private :: pcigeps
contains
  !----------------------------------------------------------------  
!!$=head1 NAME
!!$
!!$weihcig_init -initialize offset imaging operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$ierr= weihcig_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initializes migration operator whose image (R) domain is offset
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weipcig_init> L<weihcig>
!!$
!!$=cut
!!$

  !! h-domain CIG
  integer function weihcig_init() result(st)
#ifdef DEBUG
    character(len=128), parameter :: unit='weihcig_init'
    if(wei_debug()) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(wei_debug()) call out(unit)
#endif
    st=OK
  end function weihcig_init


!!$=head1 NAME
!!$
!!$weihcig - Perform offset domain imaging 
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weihcig(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 domain
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weihcig_init> L<weipcig>
!!$
!!$=cut




  !----------------------------------------------------------------
  integer function weihcig(wfld,iws,ith) result(st)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) ::iws,ith
#ifdef DEBUG
    character(len=128), parameter :: unit='weihcig'
    if(wei_debug()) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(wei_debug()) call out(unit)
#endif
    st=OK
  end function weihcig
  !----------------------------------------------------------------  
!!$=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
    type(sep3d) :: win
#ifdef DEBUG
    character(len=128), parameter :: unit='weipcig_init'
    if(my_debug()) call in(unit)
#endif


#ifdef DEBUG
    if(my_debug()) call out(unit)
#endif
    st=OK
  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(:,:) :: midin
    complex, dimension(:,:) :: exp_multx,exp_multy
    integer, intent(in) ::iws,ith
		real :: w,px,py,hx,hy
		complex :: shift
		integer :: ipx,ipy,imx,imy,ihy,ihx
    type(sep3d) :: wsep,rsep
    logical :: adj
#ifdef DEBUG
    character(len=128), parameter :: unit='weipcig'
    if(wei_debug()) call in(unit)
#endif

    call grab_current_wave_space(wsep)
    call grab_current_image_space(rsep)
    call grab_current_omega(w)
    call grab_current_adj(adj)
   
		allocate(midin(wsep%n(3),wsep%n(4)),exp_multx(wsep%n(3),rsep%n(3)))
    call add_mem(size(midin))
    call add_mem(size(exp_multx))
    allocate(exp_multy(wsep%n(4),rsep%n(4)))
    call add_mem(size(exp_multy))
		midin=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)=exp(cmplx(0,-w*px*hx))
			end do
		end do


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

	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
else             !! MODELING
	if(wsep%n(3) > 1) then
			!! get the p_hx axis done
		do imy=1,wsep%n(2)
			do imx=1,wsep%n(1)
				midin(:,:)=wfld(imx,imy,:,:,ith)
				wfld(imx,imy,:,:,ith)=0.0
				do ipx=1,rsep%n(3)
					do ihy=1,wsep%n(4)
						do ihx=1,wsep%n(3)
							wfld(imx,imy,ihx,ihy,ith)=wfld(imx,imy,ihx,ihy,ith)+&
               midin(ipx,ihy)*conjg(exp_multx(ihx,ipx))
						end do
					end do
				end do
			end do
		end do
  end if

	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(ihx,ipy)*exp_multy(ihy,ipy)
						end do
					end do
				end do
			end do
		end do
	end if
end if

 
call del_mem(size(midin)+size(exp_multx)+size(exp_multy))
deallocate(midin,exp_multx,exp_multy)

#ifdef DEBUG
    if(wei_debug()) call out(unit)
#endif
    st=OK
  end function weipcigs
  !----------------------------------------------------------------  

  !----------------------------------------------------------------
end module wei_cig

