cjyvb Subroutine

   subroutine cjyvb(v,z,Vm,Cbj,Cdj,Cby,Cdy)

!       Input :  z --- Complex argument
!                v --- Order of Jv(z) and Yv(z)
!                      ( v = n+v0, n = 0,1,2,..., 0 ≤ v0 < 1 )
!       Output:  CBJ(n) --- Jn+v0(z)
!                CDJ(n) --- Jn+v0'(z)
!                CBY(n) --- Yn+v0(z)
!                CDY(n) --- Yn+v0'(z)
!                VM --- Highest order computed

      real(wp) a0 , ga , gb , pv0 , rp2 , v , v0 , vg ,    &
                     & Vm , vv , w0
      complex(wp) ca , ca0 , cb , Cbj , Cby , cck , Cdj , Cdy , cec ,    &
               & cf , cf1 , cf2 , cfac0 , ci , cju0 , cjv0 , cjvn ,     &
               & cpz , cqz , cr
      complex(wp) cr0 , crp , crq , cs , cs0 , csk , cyv0 , cyy , z ,    &
               & z1 , z2 , zk
      integer k , k0 , m , n
      dimension Cbj(0:*) , Cdj(0:*) , Cby(0:*) , Cdy(0:*)

      rp2 = .63661977236758d0
      ci = (0.0_wp,1.0_wp)
      a0 = abs(z)
      z1 = z
      z2 = z*z
      n = int(v)
      v0 = v - n
      pv0 = pi*v0
      if ( a0<1.0e-100_wp ) then
         do k = 0 , n
            Cbj(k) = (0.0_wp,0.0_wp)
            Cdj(k) = (0.0_wp,0.0_wp)
            Cby(k) = -(1.0e+300_wp,0.0_wp)
            Cdy(k) = (1.0e+300_wp,0.0_wp)
         enddo
         if ( v0==0.0_wp ) then
            Cbj(0) = (1.0_wp,0.0_wp)
            Cdj(1) = (0.5_wp,0.0_wp)
         else
            Cdj(0) = (1.0e+300_wp,0.0_wp)
         endif
         Vm = v
         return
      endif
      if ( real(z,wp)<0.0_wp ) z1 = -z
      if ( a0<=12.0_wp ) then
         cjv0 = (1.0_wp,0.0_wp)
         cr = (1.0_wp,0.0_wp)
         do k = 1 , 40
            cr = -0.25_wp*cr*z2/(k*(k+v0))
            cjv0 = cjv0 + cr
            if ( abs(cr)<abs(cjv0)*1.0e-15_wp ) exit
         enddo
         vg = 1.0_wp + v0
         call gamma2(vg,ga)
         ca = (0.5_wp*z1)**v0/ga
         cjv0 = cjv0*ca
      else
         k0 = 11
         if ( a0>=35.0_wp ) k0 = 10
         if ( a0>=50.0_wp ) k0 = 8
         vv = 4.0_wp*v0*v0
         cpz = (1.0_wp,0.0_wp)
         crp = (1.0_wp,0.0_wp)
         do k = 1 , k0
            crp = -0.78125d-2*crp*(vv-(4.0_wp*k-3.0_wp)**2.0_wp)                 &
                & *(vv-(4.0_wp*k-1.0_wp)**2.0_wp)/(k*(2.0_wp*k-1.0_wp)*z2)
            cpz = cpz + crp
         enddo
         cqz = (1.0_wp,0.0_wp)
         crq = (1.0_wp,0.0_wp)
         do k = 1 , k0
            crq = -0.78125d-2*crq*(vv-(4.0_wp*k-1.0_wp)**2.0_wp)                 &
                & *(vv-(4.0_wp*k+1.0_wp)**2.0_wp)/(k*(2.0_wp*k+1.0_wp)*z2)
            cqz = cqz + crq
         enddo
         cqz = 0.125_wp*(vv-1.0_wp)*cqz/z1
         zk = z1 - (0.5_wp*v0+0.25_wp)*pi
         ca0 = sqrt(rp2/z1)
         cck = cos(zk)
         csk = sin(zk)
         cjv0 = ca0*(cpz*cck-cqz*csk)
         cyv0 = ca0*(cpz*csk+cqz*cck)
      endif
      if ( a0<=12.0_wp ) then
         if ( v0/=0.0_wp ) then
            cjvn = (1.0_wp,0.0_wp)
            cr = (1.0_wp,0.0_wp)
            do k = 1 , 40
               cr = -0.25_wp*cr*z2/(k*(k-v0))
               cjvn = cjvn + cr
               if ( abs(cr)<abs(cjvn)*1.0e-15_wp ) exit
            enddo
            vg = 1.0_wp - v0
            call gamma2(vg,gb)
            cb = (2.0_wp/z1)**v0/gb
            cju0 = cjvn*cb
            cyv0 = (cjv0*cos(pv0)-cju0)/sin(pv0)
         else
            cec = log(z1/2.0_wp) + gamma
            cs0 = (0.0_wp,0.0_wp)
            w0 = 0.0_wp
            cr0 = (1.0_wp,0.0_wp)
            do k = 1 , 30
               w0 = w0 + 1.0_wp/k
               cr0 = -0.25_wp*cr0/(k*k)*z2
               cs0 = cs0 + cr0*w0
            enddo
            cyv0 = rp2*(cec*cjv0-cs0)
         endif
      endif
      if ( n==0 ) n = 1
      m = msta1(a0,200)
      if ( m<n ) then
         n = m
      else
         m = msta2(a0,n,15)
      endif
      cf2 = (0.0_wp,0.0_wp)
      cf1 = (1.0e-100_wp,0.0_wp)
      do k = m , 0 , -1
         cf = 2.0_wp*(v0+k+1.0_wp)/z1*cf1 - cf2
         if ( k<=n ) Cbj(k) = cf
         cf2 = cf1
         cf1 = cf
      enddo
      cs = cjv0/cf
      do k = 0 , n
         Cbj(k) = cs*Cbj(k)
      enddo
      if ( real(z,wp)<0.0_wp ) then
         cfac0 = exp(pv0*ci)
         if ( aimag(z)<0.0_wp ) then
            cyv0 = cfac0*cyv0 - 2.0_wp*ci*cos(pv0)*cjv0
         elseif ( aimag(z)>0.0_wp ) then
            cyv0 = cyv0/cfac0 + 2.0_wp*ci*cos(pv0)*cjv0
         endif
         do k = 0 , n
            if ( aimag(z)<0.0_wp ) then
               Cbj(k) = exp(-pi*(k+v0)*ci)*Cbj(k)
            elseif ( aimag(z)>0.0_wp ) then
               Cbj(k) = exp(pi*(k+v0)*ci)*Cbj(k)
            endif
         enddo
         z1 = z1
      endif
      Cby(0) = cyv0
      do k = 1 , n
         cyy = (Cbj(k)*Cby(k-1)-2.0_wp/(pi*z))/Cbj(k-1)
         Cby(k) = cyy
      enddo
      Cdj(0) = v0/z*Cbj(0) - Cbj(1)
      do k = 1 , n
         Cdj(k) = -(k+v0)/z*Cbj(k) + Cbj(k-1)
      enddo
      Cdy(0) = v0/z*Cby(0) - Cby(1)
      do k = 1 , n
         Cdy(k) = Cby(k-1) - (k+v0)/z*Cby(k)
      enddo
      Vm = n + v0

   end subroutine cjyvb