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

!****s* WaveOpticsPropagator/ropp_pp_wopt_fresnel *
!
! NAME
!   ropp_pp_wopt_fresnel
!
! SYNOPSIS
!   CALL ropp_pp_wopt_fresnel(x, s, c)
!
! DESCRIPTION
!   This routine computes the Fresnel integrals (S and C)
!   at x by means of a rational approximation.
!
! INPUTS
!   REAL   :: x          ! arguments
!
! OUTPUT
!   REAL   :: S, C       ! The Fresnel integrals S and C at x.
!
! REFERENCES
!   Rational polynomial approximation given by
!   M.A Heald, 1985, Rational Approximations for the Fresnel Integrals,
!   Mathematics of Computation, vol 44, no. 170, pp 459-461.
!   See Table 1, Row 5.
!
! 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) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE ropp_pp_wopt_fresnel(x, s, c)

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

  IMPLICIT NONE

! Input/output

  REAL(wp), INTENT(IN)  :: x(:)
  REAL(wp), INTENT(OUT) :: s(SIZE(x)), c(SIZE(x))

!Local

! See Heald 1985, Table 1, row 3 (ie A_24 and R_34)

!  INTEGER, PARAMETER  :: na = 3
!  REAL(wp), PARAMETER :: aval(na) = &
!    (/1.0_wp,0.08218_wp,0.15108_wp/)

!  INTEGER, PARAMETER  :: nb = 5
!  REAL(wp), PARAMETER :: bval(nb) = &
!    (/2.0_wp,2.7097_wp,2.3185_wp,1.2389_wp,0.6561_wp/)

!  INTEGER, PARAMETER  :: nc = 4
!  REAL(wp), PARAMETER :: cval(nc) = &
!    (/1.0_wp,0.60427_wp,0.41159_wp,0.1917_wp/)

!  INTEGER, PARAMETER  :: nd = 5
!  REAL(wp), PARAMETER :: dval(nd) = &
!    (/1.414213562_wp,2.26794_wp,2.15594_wp,1.26057_wp,0.60353_wp/)

! See Heald 1985, Table 1, row 5 (ie A_46 and R_56)

  INTEGER, PARAMETER  :: na = 5
  REAL(wp), PARAMETER :: aval(na) = &
    (/1.0_wp,0.1945161_wp,0.2363641_wp,0.068324_wp,0.0241212_wp/)

  INTEGER, PARAMETER  :: nb = 7
  REAL(wp), PARAMETER :: bval(nb) = &
    (/2.0_wp,2.9355041_wp,2.7570246_wp,1.875721_wp,0.978113_wp,0.356681_wp,0.118247_wp/)

  INTEGER, PARAMETER  :: nc = 6
  REAL(wp), PARAMETER :: cval(nc) = &
    (/1.0_wp,0.7769507_wp,0.6460117_wp,0.3460509_wp,0.1339259_wp,0.0433995_wp/)

  INTEGER, PARAMETER  :: nd = 7
  REAL(wp), PARAMETER :: dval(nd) = &
    (/1.414213562_wp,2.5129806_wp,2.7196741_wp,1.9840524_wp,1.0917325_wp,0.4205217_wp,0.13634704_wp/)

  REAL(wp), PARAMETER :: Piby2 = pi1*0.5_wp

  REAL(wp)            :: asum(SIZE(x)),bsum(SIZE(x)),csum(SIZE(x)),dsum(SIZE(x))

  REAL(wp)            :: Rterm(SIZE(x)),Aterm(SIZE(x)),x1(SIZE(x))

  INTEGER             :: i


  x1 = ABS(x)

! compute the polynomials by nested multiplication

  asum = aval(na)
  DO i = na-1, 1, -1
    asum = asum * x1 + aval(i)
  END DO

  bsum = bval(nb)
  DO i = nb-1, 1, -1
    bsum = bsum * x1 + bval(i)
  END DO

  Aterm = asum / bsum

  csum = cval(nc)
  DO i = nc-1, 1, -1
    csum = csum * x1 + cval(i)
  END DO

  dsum = dval(nd)
  DO i = nd-1, 1, -1
    dsum = dsum * x1 + dval(i)
  END DO

  Rterm = csum / dsum

! compute the rational approximations

  c = 0.5_wp - Rterm*SIN(Piby2*(Aterm-x1**2))

  s = 0.5_wp - Rterm*COS(Piby2*(Aterm-x1**2))

! correct the sign where x < 0 (S and C are both odd)

  WHERE ( x < 0.0_wp )

    c = -1.0_wp * c
    s = -1.0_wp * s

  END WHERE

END SUBROUTINE ropp_pp_wopt_fresnel