cisia Subroutine

   subroutine cisia(x,Ci,Si)

      real(wp),intent(in) :: x !! Argument of `Ci(x)` and `Si(x)`
      real(wp),intent(out) :: Ci !! `Ci(x)`
      real(wp),intent(out) :: Si !! `Si(x)`

      real(wp) :: bj(101) , eps, x2 , xa , &
                  xa0 , xa1 , xcs , xf , xg , xg1 , xg2 , xr , xs ,&
                  xss
      integer :: k , m

      real(wp),parameter :: p2 = pi / 2.0_wp !1.570796326794897d0

      eps = 1.0e-15_wp
      x2 = x*x
      if ( x==0.0_wp ) then
         Ci = -1.0e+300_wp
         Si = 0.0_wp
      elseif ( x<=16.0_wp ) then
         xr = -0.25_wp*x2
         Ci = gamma + log(x) + xr
         do k = 2 , 40
            xr = -0.5_wp*xr*(k-1)/(k*k*(2*k-1))*x2
            Ci = Ci + xr
            if ( abs(xr)<abs(Ci)*eps ) exit
         enddo
         xr = x
         Si = x
         do k = 1 , 40
            xr = -0.5_wp*xr*(2*k-1)/k/(4*k*k+4*k+1)*x2
            Si = Si + xr
            if ( abs(xr)<abs(Si)*eps ) return
         enddo
      elseif ( x<=32.0_wp ) then
         m = int(47.2_wp+.82_wp*x)
         xa1 = 0.0_wp
         xa0 = 1.0e-100_wp
         do k = m , 1 , -1
            xa = 4.0_wp*k*xa0/x - xa1
            bj(k) = xa
            xa1 = xa0
            xa0 = xa
         enddo
         xs = bj(1)
         do k = 3 , m , 2
            xs = xs + 2.0_wp*bj(k)
         enddo
         bj(1) = bj(1)/xs
         do k = 2 , m
            bj(k) = bj(k)/xs
         enddo
         xr = 1.0_wp
         xg1 = bj(1)
         do k = 2 , m
            xr = .25_wp*xr*(2.0_wp*k-3.0_wp)**2/((k-1.0_wp)*(2.0_wp*k-1.0_wp)**2)*x
            xg1 = xg1 + bj(k)*xr
         enddo
         xr = 1.0_wp
         xg2 = bj(1)
         do k = 2 , m
            xr = .25_wp*xr*(2.0_wp*k-5.0_wp)**2/((k-1.0_wp)*(2.0_wp*k-3.0_wp)**2)*x
            xg2 = xg2 + bj(k)*xr
         enddo
         xcs = cos(x/2.0_wp)
         xss = sin(x/2.0_wp)
         Ci = gamma + log(x) - x*xss*xg1 + 2*xcs*xg2 - 2*xcs*xcs
         Si = x*xcs*xg1 + 2*xss*xg2 - sin(x)
      else
         xr = 1.0_wp
         xf = 1.0_wp
         do k = 1 , 9
            xr = -2.0_wp*xr*k*(2*k-1)/x2
            xf = xf + xr
         enddo
         xr = 1.0_wp/x
         xg = xr
         do k = 1 , 8
            xr = -2.0_wp*xr*(2*k+1)*k/x2
            xg = xg + xr
         enddo
         Ci = xf*sin(x)/x - xg*cos(x)/x
         Si = p2 - xf*cos(x)/x - xg*sin(x)/x
      endif

   end subroutine cisia