ittjyb Subroutine

   subroutine ittjyb(x,Ttj,Tty)

      real(wp),intent(in) :: x !! Variable in the limits  ( x ≥ 0 )
      real(wp),intent(out) :: Ttj !! Integration of [1-J0(t)]/t from 0 to x
      real(wp),intent(out) :: Tty !! Integration of Y0(t)/t from x to ∞

      real(wp) :: e0 , f0 , g0 , t , t1 , x1 , xt

      if ( x==0.0_wp ) then
         Ttj = 0.0_wp
         Tty = -1.0e+300_wp
      elseif ( x<=4.0_wp ) then
         x1 = x/4.0_wp
         t = x1*x1
         Ttj = ((((((.35817e-4_wp*t-.639765e-3_wp)*t+.7092535e-2_wp)*t- &
               .055544803_wp)*t+.296292677_wp)*t-.999999326_wp) &
               *t+1.999999936_wp)*t
         Tty = (((((((-.3546e-5_wp*t+.76217e-4_wp)*t-.1059499e-2_wp)*t+ &
               .010787555_wp)*t-.07810271_wp)*t+.377255736_wp) &
               *t-1.114084491_wp)*t+1.909859297_wp)*t
         e0 = gamma + log(x/2.0_wp)
         Tty = pi/6.0_wp + e0/pi*(2.0_wp*Ttj-e0) - Tty
      elseif ( x<=8.0_wp ) then
         xt = x + .25_wp*pi
         t1 = 4.0_wp/x
         t = t1*t1
         f0 = (((((.0145369_wp*t-.0666297_wp)*t+.1341551_wp)*t-.1647797_wp) &
              *t+.1608874_wp)*t-.2021547_wp)*t + .7977506_wp
         g0 = ((((((.0160672_wp*t-.0759339_wp)*t+.1576116_wp)*t-.1960154_wp) &
              *t+.1797457_wp)*t-.1702778_wp)*t+.3235819_wp)*t1
         Ttj = (f0*cos(xt)+g0*sin(xt))/(sqrt(x)*x)
         Ttj = Ttj + gamma + log(x/2.0_wp)
         Tty = (f0*sin(xt)-g0*cos(xt))/(sqrt(x)*x)
      else
         t = 8.0_wp/x
         xt = x + .25_wp*pi
         f0 = (((((.18118e-2_wp*t-.91909e-2_wp)*t+.017033_wp)*t-.9394e-3_wp) &
              *t-.051445_wp)*t-.11e-5_wp)*t + .7978846_wp
         g0 = (((((-.23731e-2_wp*t+.59842e-2_wp)*t+.24437e-2_wp)*t-.0233178_wp) &
              *t+.595e-4_wp)*t+.1620695_wp)*t
         Ttj = (f0*cos(xt)+g0*sin(xt))/(sqrt(x)*x) &
               + gamma + log(x/2.0_wp)
         Tty = (f0*sin(xt)-g0*cos(xt))/(sqrt(x)*x)
      endif

   end subroutine ittjyb