ittika Subroutine

   subroutine ittika(x,Tti,Ttk)

      real(wp),intent(in) :: x !! Variable in the limits  ( x ≥ 0 )
      real(wp),intent(out) :: Tti !! Integration of [I0(t)-1]/t from 0 to x
      real(wp),intent(out) :: Ttk !! Integration of K0(t)/t from x to ∞

      real(wp) :: b1 , e0 , r , r2 , rc , rs
      integer :: k

      real(wp),dimension(8),parameter :: c = [1.625_wp , &
                                              4.1328125_wp , &
                                              1.45380859375e+1_wp , &
                                              6.553353881835e+1_wp , &
                                              3.6066157150269e+2_wp , &
                                              2.3448727161884e+3_wp , &
                                              1.7588273098916e+4_wp , &
                                              1.4950639538279e+5_wp]

      if ( x==0.0_wp ) then
         Tti = 0.0_wp
         Ttk = 1.0e+300_wp
         return
      endif
      if ( x<40.0_wp ) then
         Tti = 1.0_wp
         r = 1.0_wp
         do k = 2 , 50
            r = .25_wp*r*(k-1.0_wp)/(k*k*k)*x*x
            Tti = Tti + r
            if ( abs(r/Tti)<1.0e-12_wp ) exit
         enddo
         Tti = Tti*.125_wp*x*x
      else
         Tti = 1.0_wp
         r = 1.0_wp
         do k = 1 , 8
            r = r/x
            Tti = Tti + c(k)*r
         enddo
         rc = x*sqrt(twopi*x)
         Tti = Tti*exp(x)/rc
      endif
      if ( x<=12.0_wp ) then
         e0 = (0.5_wp*log(x/2.0_wp)+gamma)*log(x/2.0_wp) + pi*pi/24.0_wp + &
              0.5_wp*gamma*gamma
         b1 = 1.5_wp - (gamma+log(x/2.0_wp))
         rs = 1.0_wp
         r = 1.0_wp
         do k = 2 , 50
            r = .25_wp*r*(k-1.0_wp)/(k*k*k)*x*x
            rs = rs + 1.0_wp/k
            r2 = r*(rs+1.0_wp/(2.0_wp*k)-(gamma+log(x/2.0_wp)))
            b1 = b1 + r2
            if ( abs(r2/b1)<1.0e-12_wp ) exit
         enddo
         Ttk = e0 - .125_wp*x*x*b1
      else
         Ttk = 1.0_wp
         r = 1.0_wp
         do k = 1 , 8
            r = -r/x
            Ttk = Ttk + c(k)*r
         enddo
         rc = x*sqrt(2.0_wp/pi*x)
         Ttk = Ttk*exp(-x)/rc
      endif

   end subroutine ittika