ittikb Subroutine

   subroutine ittikb(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) :: e0 , t , t1 , x1

      if ( x==0.0_wp ) then
         Tti = 0.0_wp
      elseif ( x<=5.0_wp ) then
         x1 = x/5.0_wp
         t = x1*x1
         Tti = (((((((.1263e-3_wp*t+.96442e-3_wp)*t+.968217e-2_wp)*t+.06615507_wp)&
               *t+.33116853_wp)*t+1.13027241_wp)*t+2.44140746_wp) &
               *t+3.12499991_wp)*t
      else
         t = 5.0_wp/x
         Tti = (((((((((2.1945464_wp*t-3.5195009_wp)*t-11.9094395_wp)*t+   &
               40.394734_wp)*t-48.0524115_wp)*t+28.1221478_wp)             &
               *t-8.6556013_wp)*t+1.4780044_wp)*t-.0493843_wp)             &
               *t+.1332055_wp)*t + .3989314_wp
         Tti = Tti*exp(x)/(sqrt(x)*x)
      endif
      if ( x==0.0_wp ) then
         Ttk = 1.0e+300_wp
      elseif ( x<=2.0_wp ) then
         t1 = x/2.0_wp
         t = t1*t1
         Ttk = (((((.77e-6_wp*t+.1544e-4_wp)*t+.48077e-3_wp)*t+.925821e-2_wp) &
             & *t+.10937537_wp)*t+.74999993_wp)*t
         e0 = gamma + log(x/2.0_wp)
         Ttk = pi*pi/24.0_wp + e0*(0.5_wp*e0+Tti) - Ttk
      elseif ( x<=4.0_wp ) then
         t = 2.0_wp/x
         Ttk = (((.06084_wp*t-.280367_wp)*t+.590944_wp)*t-.850013_wp) &
             & *t + 1.234684_wp
         Ttk = Ttk*exp(-x)/(sqrt(x)*x)
      else
         t = 4.0_wp/x
         Ttk = (((((.02724_wp*t-.1110396_wp)*t+.2060126_wp)*t-.2621446_wp)  &
             & *t+.3219184_wp)*t-.5091339_wp)*t + 1.2533141_wp
         Ttk = Ttk*exp(-x)/(sqrt(x)*x)
      endif

   end subroutine ittikb