ittjya Subroutine

   subroutine ittjya(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) :: a0 , b1 , bj0 , bj1 , by0 , by1 , e0 , g0 , &
                  g1 , px , qx , r , r0 , r1 , r2 , rs , t , vt , xk
      integer :: k , l

      if ( x==0.0_wp ) then
         Ttj = 0.0_wp
         Tty = -1.0e+300_wp
      elseif ( x<=20.0_wp ) then
         Ttj = 1.0_wp
         r = 1.0_wp
         do k = 2 , 100
            r = -0.25_wp*r*(k-1.0_wp)/(k*k*k)*x*x
            Ttj = Ttj + r
            if ( abs(r)<abs(Ttj)*1.0e-12_wp ) exit
         enddo
         Ttj = Ttj*.125_wp*x*x
         e0 = 0.5_wp*(pi*pi/6.0_wp-gamma*gamma) - (0.5_wp*log(x/2.0_wp)+gamma) &
              *log(x/2.0_wp)
         b1 = gamma + log(x/2.0_wp) - 1.5_wp
         rs = 1.0_wp
         r = -1.0_wp
         do k = 2 , 100
            r = -0.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)<abs(b1)*1.0e-12_wp ) exit
         enddo
         Tty = 2.0_wp/pi*(e0+.125_wp*x*x*b1)
      else
         a0 = sqrt(2.0_wp/(pi*x))
         bj0 = 0.0_wp
         by0 = 0.0_wp
         bj1 = 0.0_wp
         do l = 0 , 1
            vt = 4.0_wp*l*l
            px = 1.0_wp
            r = 1.0_wp
            do k = 1 , 14
               r = -.0078125_wp*r*(vt-(4.0_wp*k-3.0_wp)**2)/(x*k) &
                   *(vt-(4.0_wp*k-1.0_wp)**2)/((2.0_wp*k-1.0_wp)*x)
               px = px + r
               if ( abs(r)<abs(px)*1.0e-12_wp ) exit
            enddo
            qx = 1.0_wp
            r = 1.0_wp
            do k = 1 , 14
               r = -.0078125_wp*r*(vt-(4.0_wp*k-1.0_wp)**2)/(x*k) &
                   *(vt-(4.0_wp*k+1.0_wp)**2)/(2.0_wp*k+1.0_wp)/x
               qx = qx + r
               if ( abs(r)<abs(qx)*1.0e-12_wp ) exit
            enddo
            qx = .125_wp*(vt-1.0_wp)/x*qx
            xk = x - (.25_wp+.5_wp*l)*pi
            bj1 = a0*(px*cos(xk)-qx*sin(xk))
            by1 = a0*(px*sin(xk)+qx*cos(xk))
            if ( l==0 ) then
               bj0 = bj1
               by0 = by1
            endif
         enddo
         t = 2.0_wp/x
         g0 = 1.0_wp
         r0 = 1.0_wp
         do k = 1 , 10
            r0 = -k*k*t*t*r0
            g0 = g0 + r0
         enddo
         g1 = 1.0_wp
         r1 = 1.0_wp
         do k = 1 , 10
            r1 = -k*(k+1.0_wp)*t*t*r1
            g1 = g1 + r1
         enddo
         Ttj = 2.0_wp*g1*bj0/(x*x) - g0*bj1/x + gamma + log(x/2.0_wp)
         Tty = 2.0_wp*g1*by0/(x*x) - g0*by1/x
      endif

   end subroutine ittjya