ffk Subroutine

subroutine ffk ( ks, x, fr, fi, fm, fa, gr, gi, gm, ga )

  !*****************************************************************************80
  !
  !! FFK computes modified Fresnel integrals F+/-(x) and K+/-(x).
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    23 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) KS, the sign code.
  !    0, to calculate F+(x) and K+(x);
  !    1, to calculate F_(x) and K_(x).
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) FR, FI, FM, FA, the values of
  !    Re[F+/-(x)], Im[F+/-(x)], |F+/-(x)|, Arg[F+/-(x)]  (Degs.).
  !
  !    Output, real ( kind = 8 ) GR, GI, GM, GA, the values of
  !    Re[K+/-(x)], Im[K+/-(x)], |K+/-(x)|, Arg[K+/-(x)]  (Degs.).
  !       
  implicit none

  real ( kind = 8 ) c1
  real ( kind = 8 ) cs
  real ( kind = 8 ) eps
  real ( kind = 8 ) fa
  real ( kind = 8 ) fi
  real ( kind = 8 ) fi0
  real ( kind = 8 ) fm
  real ( kind = 8 ) fr
  real ( kind = 8 ) ga
  real ( kind = 8 ) gi
  real ( kind = 8 ) gm
  real ( kind = 8 ) gr
  integer ( kind = 4 ) k
  integer ( kind = 4 ) ks
  integer ( kind = 4 ) m
  real ( kind = 8 ) p2p
  real ( kind = 8 ) pi
  real ( kind = 8 ) pp2
  real ( kind = 8 ) s1
  real ( kind = 8 ) srd
  real ( kind = 8 ) ss
  real ( kind = 8 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) x4
  real ( kind = 8 ) xa
  real ( kind = 8 ) xc
  real ( kind = 8 ) xf
  real ( kind = 8 ) xf0
  real ( kind = 8 ) xf1
  real ( kind = 8 ) xg
  real ( kind = 8 ) xp
  real ( kind = 8 ) xq
  real ( kind = 8 ) xq2
  real ( kind = 8 ) xr
  real ( kind = 8 ) xs
  real ( kind = 8 ) xsu
  real ( kind = 8 ) xw

  srd = 57.29577951308233D+00
  eps = 1.0D-15
  pi = 3.141592653589793D+00
  pp2 = 1.2533141373155D+00
  p2p = 0.7978845608028654D+00
  xa = abs ( x )
  x2 = x * x
  x4 = x2 * x2

  if ( x == 0.0D+00 ) then

     fr = 0.5D+00 * sqrt ( 0.5D+00 * pi )
     fi = ( -1.0D+00 ) ** ks * fr
     fm = sqrt ( 0.25D+00 * pi )
     fa = ( -1.0D+00 ) ** ks * 45.0D+00
     gr = 0.5D+00
     gi = 0.0D+00
     gm = 0.5D+00
     ga = 0.0D+00

  else

     if ( xa <= 2.5D+00 ) then

        xr = p2p * xa
        c1 = xr
        do k = 1, 50
           xr = -0.5D+00 * xr * ( 4.0D+00 * k - 3.0D+00 ) / k &
                / ( 2.0D+00 * k - 1.0D+00 ) &
                / ( 4.0D+00 * k + 1.0D+00 ) * x4
           c1 = c1 + xr
           if ( abs ( xr / c1 ) < eps ) then
              exit
           end if
        end do

        s1 = p2p * xa * xa * xa / 3.0D+00
        xr = s1
        do k = 1, 50
           xr = -0.5D+00 * xr * ( 4.0D+00 * k - 1.0D+00 ) &
                / k / ( 2.0D+00 * k + 1.0D+00 ) &
                / ( 4.0D+00 * k + 3.0D+00 ) * x4
           s1 = s1 + xr
           if ( abs ( xr / s1 ) < eps ) then
              exit
           end if
        end do

     else if ( xa < 5.5D+00 ) then

        m = int ( 42.0D+00 + 1.75D+00 * x2 )
        xsu = 0.0D+00
        xc = 0.0D+00
        xs = 0.0D+00
        xf1 = 0.0D+00
        xf0 = 1.0D-100
        do k = m, 0, -1
           xf = ( 2.0D+00 * k + 3.0D+00 ) * xf0 / x2 - xf1
           if ( k == 2 * int ( k / 2 ) )  then
              xc = xc + xf
           else
              xs = xs + xf
           end if
           xsu = xsu + ( 2.0D+00 * k + 1.0D+00 ) * xf * xf
           xf1 = xf0
           xf0 = xf
        end do
        xq = sqrt ( xsu )
        xw = p2p * xa / xq
        c1 = xc * xw
        s1 = xs * xw

     else

        xr = 1.0D+00
        xf = 1.0D+00
        do k = 1, 12
           xr = -0.25D+00 * xr * ( 4.0D+00 * k - 1.0D+00 ) &
                * ( 4.0D+00 * k - 3.0D+00 ) / x4
           xf = xf + xr
        end do
        xr = 1.0D+00 / ( 2.0D+00 * xa * xa )
        xg = xr
        do k = 1, 12
           xr = -0.25D+00 * xr * ( 4.0D+00 * k + 1.0D+00 ) &
                * ( 4.0D+00 * k - 1.0D+00 ) / x4
           xg = xg + xr
        end do
        c1 = 0.5D+00 + ( xf * sin ( x2 ) - xg * cos ( x2 ) ) &
             / sqrt ( 2.0D+00 * pi ) / xa
        s1 = 0.5D+00 - ( xf * cos ( x2 ) + xg * sin ( x2 ) ) &
             / sqrt ( 2.0D+00 * pi ) / xa

     end if

     fr = pp2 * ( 0.5D+00 - c1 )
     fi0 = pp2 * ( 0.5D+00 - s1 )
     fi = ( -1.0D+00 ) ** ks * fi0
     fm = sqrt ( fr * fr + fi * fi )

     if ( 0.0D+00 <= fr ) then
        fa = srd * atan ( fi / fr )
     else if ( 0.0D+00 < fi ) then
        fa = srd * ( atan ( fi / fr ) + pi )
     else if ( fi < 0.0D+00 ) then
        fa = srd * ( atan ( fi / fr ) - pi )
     end if

     xp = x * x + pi / 4.0D+00
     cs = cos ( xp )
     ss = sin ( xp )
     xq2 = 1.0D+00 / sqrt ( pi )
     gr = xq2 * ( fr * cs + fi0 * ss )
     gi = ( -1.0D+00 ) ** ks * xq2 * ( fi0 * cs - fr * ss )
     gm = sqrt ( gr * gr + gi * gi )

     if ( 0.0D+00 <= gr ) then
        ga = srd * atan ( gi / gr )
     else if ( 0.0D+00 < gi ) then
        ga = srd * ( atan ( gi / gr ) + pi )
     else if ( gi < 0.0D+00 ) then
        ga = srd * ( atan ( gi / gr ) - pi )
     end if

     if ( x < 0.0D+00 ) then
        fr = pp2 - fr
        fi = ( -1.0D+00 ) ** ks * pp2 - fi
        fm = sqrt ( fr * fr + fi * fi )
        fa = srd * atan ( fi / fr )
        gr = cos ( x * x ) - gr
        gi = - ( -1.0D+00 ) ** ks * sin ( x * x ) - gi
        gm = sqrt ( gr * gr + gi * gi )
        ga = srd * atan ( gi / gr )
     end if

  end if

  return
end subroutine ffk