cjyva Subroutine

public subroutine cjyva(v, z, Vm, Cbj, Cdj, Cby, Cdy)

Compute Bessel functions Jv(z), Yv(z) and their derivatives for a complex argument

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: v

Order of Jv(z) and Yv(z) ( v = n+v0, n = 0,1,2,..., 0 ≤ v0 < 1)
complex(kind=wp), intent(in) :: z

Complex argument
integer, intent(out) :: Vm

Highest order computed
complex(kind=wp), intent(out), dimension(0:*) :: Cbj

CBJ(n) --- Jn+v0(z)
complex(kind=wp), intent(out), dimension(0:*) :: Cdj

CDJ(n) --- Jn+v0'(z)
complex(kind=wp), intent(out), dimension(0:*) :: Cby

CBY(n) --- Yn+v0(z)
complex(kind=wp), intent(out), dimension(0:*) :: Cdy

CDY(n) --- Yn+v0'(z)

Calls

proc~~cjyva~~CallsGraph proc~cjyva specfun_module::cjyva proc~gamma2 specfun_module::gamma2 proc~cjyva->proc~gamma2 proc~msta1 specfun_module::msta1 proc~cjyva->proc~msta1 proc~msta2 specfun_module::msta2 proc~cjyva->proc~msta2 proc~envj specfun_module::envj proc~msta1->proc~envj proc~msta2->proc~envj

Source Code

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

      real(wp),intent(in) :: v !! Order of `Jv(z)` and `Yv(z)`
                               !! ( `v = n+v0`, `n = 0,1,2,...`, `0 ≤ v0 < 1 `)
      complex(wp),intent(in) :: z !! Complex argument
      integer,intent(out) :: Vm !! Highest order computed
      complex(wp),dimension(0:*),intent(out) :: Cbj !! `CBJ(n)` --- `Jn+v0(z)`
      complex(wp),dimension(0:*),intent(out) :: Cdj !! `CDJ(n)` --- `Jn+v0'(z)`
      complex(wp),dimension(0:*),intent(out) :: Cby !! `CBY(n)` --- `Yn+v0(z)`
      complex(wp),dimension(0:*),intent(out) :: Cdy !! `CDY(n)` --- `Yn+v0'(z)`

      real(wp) :: a0 , ga , gb , pv0 , pv1 , v0 , vg , vl , &
                  vv , w0 , w1 , wa , ya0 , ya1 , yak
      complex(wp) :: ca , ca0 , cb , cck , cec , &
                     cf , cf0 , cf1 , cf2 , cfac0 , cfac1 , cg0 , cg1 , &
                     ch0 , ch1 , ch2 , ci , cju0 , cju1 , cjv0 , cjv1 , &
                     cjvl , cp11 , cp12 , cp21 , cp22 , cpz , cqz , cr , &
                     cr0 , cr1 , crp , crq , cs , cs0 , cs1 , csk , cyk , &
                     cyl1 , cyl2 , cylk , cyv0 , cyv1 , z1 , z2 , zk
      integer :: j , k , k0 , l , lb , lb0 , m , n

      real(wp),parameter :: rp2 = 2.0_wp / pi ! 0.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
      pv1 = pi*(1.0_wp+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
      lb0 = 0.0_wp
      if ( real(z,wp)<0.0_wp ) z1 = -z
      if ( a0<=12.0_wp ) then
         do l = 0 , 1
            vl = v0 + l
            cjvl = (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+vl))
               cjvl = cjvl + cr
               if ( abs(cr)<abs(cjvl)*1.0e-15_wp ) exit
            enddo
            vg = 1.0_wp + vl
            call gamma2(vg,ga)
            ca = (0.5_wp*z1)**vl/ga
            if ( l==0 ) cjv0 = cjvl*ca
            if ( l==1 ) cjv1 = cjvl*ca
         enddo
      else
         k0 = 11
         if ( a0>=35.0_wp ) k0 = 10
         if ( a0>=50.0_wp ) k0 = 8
         do j = 0 , 1
            vv = 4.0_wp*(j+v0)*(j+v0)
            cpz = (1.0_wp,0.0_wp)
            crp = (1.0_wp,0.0_wp)
            do k = 1 , k0
               crp = -0.78125e-2_wp*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.78125e-2_wp*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*(j+v0)+0.25_wp)*pi
            ca0 = sqrt(rp2/z1)
            cck = cos(zk)
            csk = sin(zk)
            if ( j==0 ) then
               cjv0 = ca0*(cpz*cck-cqz*csk)
               cyv0 = ca0*(cpz*csk+cqz*cck)
            elseif ( j==1 ) then
               cjv1 = ca0*(cpz*cck-cqz*csk)
               cyv1 = ca0*(cpz*csk+cqz*cck)
            endif
         enddo
      endif
      if ( a0<=12.0_wp ) then
         if ( v0/=0.0_wp ) then
            do l = 0 , 1
               vl = v0 + l
               cjvl = (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-vl))
                  cjvl = cjvl + cr
                  if ( abs(cr)<abs(cjvl)*1.0e-15_wp ) exit
               enddo
               vg = 1.0_wp - vl
               call gamma2(vg,gb)
               cb = (2.0_wp/z1)**vl/gb
               if ( l==0 ) cju0 = cjvl*cb
               if ( l==1 ) cju1 = cjvl*cb
            enddo
            cyv0 = (cjv0*cos(pv0)-cju0)/sin(pv0)
            cyv1 = (cjv1*cos(pv1)-cju1)/sin(pv1)
         else
            cec = log(z1/2.0_wp) + gamma ! 0.5772156649015329
            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)
            cs1 = (1.0_wp,0.0_wp)
            w1 = 0.0_wp
            cr1 = (1.0_wp,0.0_wp)
            do k = 1 , 30
               w1 = w1 + 1.0_wp/k
               cr1 = -0.25_wp*cr1/(k*(k+1))*z2
               cs1 = cs1 + cr1*(2.0_wp*w1+1.0_wp/(k+1.0_wp))
            enddo
            cyv1 = rp2*(cec*cjv1-1.0_wp/z1-0.25_wp*z1*cs1)
         endif
      endif
      if ( real(z,wp)<0.0_wp ) then
         cfac0 = exp(pv0*ci)
         cfac1 = exp(pv1*ci)
         if ( aimag(z)<0.0_wp ) then
            cyv0 = cfac0*cyv0 - 2.0_wp*ci*cos(pv0)*cjv0
            cyv1 = cfac1*cyv1 - 2.0_wp*ci*cos(pv1)*cjv1
            cjv0 = cjv0/cfac0
            cjv1 = cjv1/cfac1
         elseif ( aimag(z)>0.0_wp ) then
            cyv0 = cyv0/cfac0 + 2.0_wp*ci*cos(pv0)*cjv0
            cyv1 = cyv1/cfac1 + 2.0_wp*ci*cos(pv1)*cjv1
            cjv0 = cfac0*cjv0
            cjv1 = cfac1*cjv1
         endif
      endif
      Cbj(0) = cjv0
      Cbj(1) = cjv1
      if ( n>=2 .and. n<=int(0.25_wp*a0) ) then
         cf0 = cjv0
         cf1 = cjv1
         do k = 2 , n
            cf = 2.0_wp*(k+v0-1.0_wp)/z*cf1 - cf0
            Cbj(k) = cf
            cf0 = cf1
            cf1 = cf
         enddo
      elseif ( n>=2 ) then
         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)/z*cf1 - cf2
            if ( k<=n ) Cbj(k) = cf
            cf2 = cf1
            cf1 = cf
         enddo
         if ( abs(cjv0)>abs(cjv1) ) cs = cjv0/cf
         if ( abs(cjv0)<=abs(cjv1) ) cs = cjv1/cf2
         do k = 0 , n
            Cbj(k) = cs*Cbj(k)
         enddo
      endif
      Cdj(0) = v0/z*Cbj(0) - Cbj(1)
      do k = 1 , n
         Cdj(k) = -(k+v0)/z*Cbj(k) + Cbj(k-1)
      enddo
      Cby(0) = cyv0
      Cby(1) = cyv1
      ya0 = abs(cyv0)
      lb = 0
      cg0 = cyv0
      cg1 = cyv1
      do k = 2 , n
         cyk = 2.0_wp*(v0+k-1.0_wp)/z*cg1 - cg0
         if ( abs(cyk)<=1.0e+290_wp ) then
            yak = abs(cyk)
            ya1 = abs(cg0)
            if ( yak<ya0 .and. yak<ya1 ) lb = k
            Cby(k) = cyk
            cg0 = cg1
            cg1 = cyk
         endif
      enddo
      if ( lb>4 .and. aimag(z)/=0.0_wp ) then
         do
           if ( lb/=lb0 ) then
               ch2 = (1.0_wp,0.0_wp)
               ch1 = (0.0_wp,0.0_wp)
               lb0 = lb
               do k = lb , 1 , -1
                  ch0 = 2.0_wp*(k+v0)/z*ch1 - ch2
                  ch2 = ch1
                  ch1 = ch0
               enddo
               cp12 = ch0
               cp22 = ch2
               ch2 = (0.0_wp,0.0_wp)
               ch1 = (1.0_wp,0.0_wp)
               do k = lb , 1 , -1
                  ch0 = 2.0_wp*(k+v0)/z*ch1 - ch2
                  ch2 = ch1
                  ch1 = ch0
               enddo
               cp11 = ch0
               cp21 = ch2
               if ( lb==n ) Cbj(lb+1) = 2.0_wp*(lb+v0)/z*Cbj(lb) - Cbj(lb-1)
               if ( abs(Cbj(0))>abs(Cbj(1)) ) then
                  Cby(lb+1) = (Cbj(lb+1)*cyv0-2.0_wp*cp11/(pi*z))/Cbj(0)
                  Cby(lb) = (Cbj(lb)*cyv0+2.0_wp*cp12/(pi*z))/Cbj(0)
               else
                  Cby(lb+1) = (Cbj(lb+1)*cyv1-2.0_wp*cp21/(pi*z))/Cbj(1)
                  Cby(lb) = (Cbj(lb)*cyv1+2.0_wp*cp22/(pi*z))/Cbj(1)
               endif
               cyl2 = Cby(lb+1)
               cyl1 = Cby(lb)
               do k = lb - 1 , 0 , -1
                  cylk = 2.0_wp*(k+v0+1.0_wp)/z*cyl1 - cyl2
                  Cby(k) = cylk
                  cyl2 = cyl1
                  cyl1 = cylk
               enddo
               cyl1 = Cby(lb)
               cyl2 = Cby(lb+1)
               do k = lb + 1 , n - 1
                  cylk = 2.0_wp*(k+v0)/z*cyl2 - cyl1
                  Cby(k+1) = cylk
                  cyl1 = cyl2
                  cyl2 = cylk
               enddo
               do k = 2 , n
                  wa = abs(Cby(k))
                  if ( wa<abs(Cby(k-1)) ) lb = k
               enddo
            else
               exit
            endif
         end do
      endif
      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 cjyva