cisib Subroutine

   subroutine cisib(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) :: fx , gx , x2

      x2 = x*x
      if ( x==0.0_wp ) then
         Ci = -1.0e+300_wp
         Si = 0.0_wp
      elseif ( x<=1.0_wp ) then
         Ci = ((((-3.0e-8_wp*x2+3.10e-6_wp)*x2-2.3148e-4_wp)*x2+1.041667e-2_wp) &
            & *x2-0.25)*x2 + 0.577215665_wp + log(x)
         Si = ((((3.1e-7_wp*x2-2.834e-5_wp)*x2+1.66667e-003_wp)*x2-5.555556e-002_wp)&
            & *x2+1.0_wp)*x
      else
         fx = ((((x2+38.027264_wp)*x2+265.187033_wp)*x2+335.67732_wp)    &
            & *x2+38.102495_wp)                                          &
            & /((((x2+40.021433_wp)*x2+322.624911_wp)*x2+570.23628_wp)   &
            & *x2+157.105423d0)
         gx = ((((x2+42.242855_wp)*x2+302.757865_wp)*x2+352.018498_wp)   &
            & *x2+21.821899_wp)                                          &
            & /((((x2+48.196927_wp)*x2+482.485984_wp)*x2+1114.978885_wp) &
            & *x2+449.690326_wp)/x
         Ci = fx*sin(x)/x - gx*cos(x)/x
         Si = 1.570796327_wp - fx*cos(x)/x - gx*sin(x)/x
      endif

   end subroutine cisib