! $Id: ropp_pp_wopt.f90 r4887 2016-05-25 15:48:28Z sti $

!****s* WaveOpticsPropagator/ropp_pp_wopt_propagate *
!
! NAME
!   ropp_pp_wopt_propagate
!
! SYNOPSIS
!
!   CALL ropp_pp_wopt_propagate(U0, k, Dy, Dx, U)
!
! DESCRIPTION
!
!   Propagate complex signal (phase and amplitude) from one screen to another,
!   by means of an FFT decomposition, phase adjustment, and 
!   FFT recombination.
!
! INPUTS
!   COMPLEX(wp), INTENT(IN)  :: U0
!   REAL(wp), INTENT(IN)     :: k
!   REAL(wp), INTENT(IN)     :: Dy
!   REAL(wp), INTENT(IN)     :: Dx
!
! OUTPUT
!   COMPLEX(wp), INTENT(OUT) :: U
!
! REFERENCES
!   See Zverev, Radiooptics, Moscow, Nauka, 1976) for details
!   Plane wave in form exp(i*kx*x + i*ky*y) is accepted.
!   Inhomogeneous waves are excluded.
!
! NOTES
!   See RSR 28 for details
!
! AUTHOR
!   ECMWF, Reading, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) Copyright 1998, M. E. Gorbunov.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE ropp_pp_wopt_propagate(U0, k, Dy, Dx, U)

!-------------------------------------------------------------------------------
! 1. Declarations
!-------------------------------------------------------------------------------

  USE typesizes, ONLY: wp => EightByteReal
  USE ropp_pp,   ONLY: ropp_pp_fft
  USE ropp_pp_constants, ONLY: pi1 => pi

  IMPLICIT NONE

! Input
  COMPLEX(wp), INTENT(IN)  :: U0(:)   ! Complex field in source plane
                                      ! (Size(U0) must be a power of 2)
  REAL(wp), INTENT(IN)     :: k       ! Wave vector (2*Pi/Wavelength)
  REAL(wp), INTENT(IN)     :: Dy      ! Discretization step in source plane
  REAL(wp), INTENT(IN)     :: Dx      ! Distance from source to observation plane

! Output
  COMPLEX(wp), INTENT(OUT) :: U(SIZE(U0))  ! Propagated field in observation plane

! Local Parameters:
  COMPLEX(wp), PARAMETER ::  Ci = (0.0_wp, 1.0_wp)  ! Sqrt(-1)

! Local Scalars:
  INTEGER      :: N     ! Number of data points
  INTEGER      :: NP    ! Number of homegeneous plane waves
  INTEGER      :: i     ! Array index
  INTEGER      :: Sgn   ! Fourier transform direction
  REAL(wp)     :: L     ! Plane length
  REAL(wp)     :: twopi ! 2*pi

  COMPLEX(wp), ALLOCATABLE  :: U_tmp(:)     ! for storing the arrays
  COMPLEX(wp)               :: phase_factor ! to save time

!----------------------------------------------------------
! 2. Plane wave expansion
!----------------------------------------------------------

  twopi = 2.0_wp * pi1

  N = SIZE(U0)

  U = U0

  Sgn = -1

  CALL ropp_pp_fft(U, sgn)

  U(:) = U(:) / N      ! division needed because of FFT definition

  ALLOCATE (U_tmp(N))

  U_tmp(:) = U(:)

!----------------------------------------------------------
! 3. Propagation of plane waves
!----------------------------------------------------------

  L  = N*ABS(Dy)

  NP = MIN(N/2, INT(L*k/twopi))

  U(1) = EXP( Ci * MODULO(k*Dx, twopi) ) * U_tmp(1)

  DO i=1,NP
    phase_factor = EXP( Ci * MODULO(SQRT(k**2 - (twopi*i/L)**2)*Dx, twopi) )
    U(i+1)     = phase_factor * U_tmp(i+1)
    U(N+1-i)   = phase_factor * U_tmp(N+1-i)
  END DO

  DEALLOCATE (U_tmp)

!----------------------------------------------------------
! 4. Summation of propagated plane waves
!----------------------------------------------------------

  Sgn = 1

  CALL ropp_pp_fft(U, Sgn)

END SUBROUTINE ropp_pp_wopt_propagate