cfs Subroutine

   subroutine cfs(z,Zf,Zd)

      complex(wp),intent(in) :: z !! Argument of S(z)
      complex(wp),intent(out) :: Zf !! S(z)
      complex(wp),intent(out) :: Zd !! S'(z)

      complex(wp) :: cf , cf0 , cf1 , cg , cr , d , s , z0 , zp , zp2
      real(wp) :: w0 , wb , wb0
      integer :: k , m

      real(wp),parameter :: eps = 1.0e-14_wp

      w0 = abs(z)
      zp = halfpi*z*z
      zp2 = zp*zp
      z0 = (0.0_wp,0.0_wp)
      if ( z==z0 ) then
         s = z0
      elseif ( w0<=2.5_wp ) then
         s = z*zp/3.0_wp
         cr = s
         wb0 = 0.0_wp
         do k = 1 , 80
            cr = -0.5_wp*cr*(4.0_wp*k-1.0_wp)/k/(2.0_wp*k+1.0_wp) &
                 /(4.0_wp*k+3.0_wp)*zp2
            s = s + cr
            wb = abs(s)
            if ( abs(wb-wb0)<eps .and. k>10 ) exit
            wb0 = wb
         enddo
      elseif ( w0>2.5_wp .and. w0<4.5_wp ) then
         m = 85
         s = z0
         cf1 = z0
         cf0 = (1.0e-100_wp,0.0_wp)
         do k = m , 0 , -1
            cf = (2.0_wp*k+3.0_wp)*cf0/zp - cf1
            if ( k/=int(k/2)*2 ) s = s + cf
            cf1 = cf0
            cf0 = cf
         enddo
         s = 2.0_wp/(pi*z)*sin(zp)/cf*s
      else
         ! Auxiliary functions f(z) and g(z) can be computed using an
         ! asymptotic expansion in the right quadrant |arg(z)| <= pi/4, not pi/2
         ! as sometimes suggested. Use the symmetry S(z) = -iS(-iz).
         ! Interestingly, most of the expansion code is the same across
         ! the quadrants. (The fourth power in Z is the equalizer here.)
         ! Only one constant has to be adapted.
         if ( aimag(z)>-real(z,wp) .and. aimag(z)<=real(z,wp) ) then
            ! right quadrant
            d = cmplx(0.5_wp,0.0_wp,wp)
         elseif ( aimag(z)>real(z,wp) .and. aimag(z)>=-real(z,wp) ) then
            ! upper quadrant
            d = cmplx(0.0_wp,-0.5_wp,wp)
         elseif ( aimag(z)<-real(z,wp) .and. aimag(z)>=real(z,wp) ) then
            ! left quadrant
            d = cmplx(-0.5_wp,0.0_wp,wp)
         else
            ! lower quadrant
            d = cmplx(0.0_wp,0.5_wp,wp)
         endif
         cr = (1.0_wp,0.0_wp)
         cf = (1.0_wp,0.0_wp)
         do k = 1 , 20
            cr = -0.25_wp*cr*(4.0_wp*k-1.0_wp)*(4.0_wp*k-3.0_wp)/zp2
            cf = cf + cr
         enddo
         cr = (1.0_wp,0.0_wp)
         cg = (1.0_wp,0.0_wp)
         do k = 1 , 12
            cr = -0.25_wp*cr*(4.0_wp*k+1.0_wp)*(4.0_wp*k-1.0_wp)/zp2
            cg = cg + cr
         enddo
         cg = cg/(pi*z*z)
         s = d - (cf*cos(zp)+cg*sin(zp))/(pi*z)
      endif
      Zf = s
      Zd = sin(halfpi*z*z)
   end subroutine cfs