cjylv Subroutine

   subroutine cjylv(v,z,Cbjv,Cdjv,Cbyv,Cdyv)

      real(wp),intent(in) :: v !! Order of Jv(z) and Yv(z)
      complex(wp),intent(in) :: z !! Complex argument
      real(wp),intent(out) :: Cbjv !! Jv(z)
      real(wp),intent(out) :: Cdjv !! Jv'(z)
      real(wp),intent(out) :: Cbyv !! Yv(z)
      real(wp),intent(out) :: Cdyv !! Yv'(z)

      real(wp) :: a(91) , v0 , vr
      complex(wp) :: ceta , cf(12) , cfj , cfy , csj , csy , ct , ct2 , cws
      integer :: i , k , km , l , l0 , lf

      km = 12
      call cjk(km,a)
      do l = 1 , 0 , -1
         v0 = v - l
         cws = sqrt(1.0_wp-(z/v0)*(z/v0))
         ceta = cws + log(z/v0/(1.0_wp+cws))
         ct = 1.0_wp/cws
         ct2 = ct*ct
         do k = 1 , km
            l0 = k*(k+1)/2 + 1
            lf = l0 + k
            cf(k) = a(lf)
            do i = lf - 1 , l0 , -1
               cf(k) = cf(k)*ct2 + a(i)
            enddo
            cf(k) = cf(k)*ct**k
         enddo
         vr = 1.0_wp/v0
         csj = (1.0_wp,0.0_wp)
         do k = 1 , km
            csj = csj + cf(k)*vr**k
         enddo
         Cbjv = sqrt(ct/(twopi*v0))*exp(v0*ceta)*csj
         if ( l==1 ) cfj = Cbjv
         csy = (1.0_wp,0.0_wp)
         do k = 1 , km
            csy = csy + (-1)**k*cf(k)*vr**k
         enddo
         Cbyv = -sqrt(2.0_wp*ct/(pi*v0))*exp(-v0*ceta)*csy
         if ( l==1 ) cfy = Cbyv
      enddo
      Cdjv = -v/z*Cbjv + cfj
      Cdyv = -v/z*Cbyv + cfy

   end subroutine cjylv