jdzo Subroutine

public subroutine jdzo(Nt, n, m, p, Zo)

Compute the zeros of Bessel functions Jn(x) and Jn'(x), and arrange them in the order of their magnitudes

Note

SciPy: Changed P from a character array to an integer array.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Nt

Number of total zeros ( NT ≤ 1200 )
integer, intent(out) :: n(1400)

n, order of Jn(x) or Jn'(x) associated with the L-th zero
integer, intent(out) :: m(1400)

m, serial number of the zeros of Jn(x) or Jn'(x) associated with the L-th zero ( L is the serial number of all the zeros of Jn(x) and Jn'(x) )
integer, intent(out) :: p(1400)

0 (TM) or 1 (TE), a code for designating the zeros of Jn(x) or Jn'(x). In the waveguide applications, the zeros of Jn(x) correspond to TM modes and those of Jn'(x) correspond to TE modes
real(kind=wp), intent(out) :: Zo(0:1400)

Value of the L-th zero of Jn(x) and Jn'(x)

Calls

proc~~jdzo~~CallsGraph proc~jdzo specfun_module::jdzo proc~bjndd specfun_module::bjndd proc~jdzo->proc~bjndd

Source Code

   subroutine jdzo(Nt,n,m,p,Zo)


      integer,intent(in) :: Nt !! Number of total zeros ( NT ≤ 1200 )
      integer,intent(out) :: n(1400) !! n, order of Jn(x) or Jn'(x) associated
                                     !! with the L-th zero
      integer,intent(out) :: m(1400) !! m, serial number of the zeros of Jn(x)
                                     !! or Jn'(x) associated with the L-th zero
                                     !! ( L is the serial number of all the
                                     !! zeros of Jn(x) and Jn'(x) )
      integer,intent(out) :: p(1400) !! 0 (TM) or 1 (TE), a code for designating the
                                     !! zeros of Jn(x) or Jn'(x).
                                     !! In the waveguide applications, the zeros
                                     !! of Jn(x) correspond to TM modes and
                                     !! those of Jn'(x) correspond to TE modes
      real(wp),intent(out) :: Zo(0:1400) !! Value of the L-th zero of Jn(x)
                                         !! and Jn'(x)

      real(wp) :: x , x0 , x1 , x2 , xm
      integer :: i , j , k , l , l0 , l1 , l2 , mm , nm
      integer,dimension(70) :: n1, m1, p1
      real(wp),dimension(0:70) :: zoc
      real(wp),dimension(101) :: bj , dj , fj

      real(wp),parameter :: tol = 1.0e-10_wp

      x = 0
      zoc(0) = 0
      if ( Nt<600 ) then
         xm = -1.0_wp + 2.248485_wp*Nt**0.5_wp - .0159382_wp*Nt + 3.208775e-4_wp*Nt**1.5_wp
         nm = int(14.5_wp+.05875_wp*Nt)
         mm = int(.02_wp*Nt) + 6
      else
         xm = 5.0_wp + 1.445389_wp*Nt**.5_wp + .01889876_wp*Nt - 2.147763e-4_wp*Nt**1.5_wp
         nm = int(27.8_wp+.0327_wp*Nt)
         mm = int(.01088_wp*Nt) + 10
      endif
      l0 = 0
      do i = 1 , nm
         x1 = .407658_wp + .4795504_wp*(i-1)**.5_wp + .983618_wp*(i-1)
         x2 = 1.99535_wp + .8333883_wp*(i-1)**.5_wp + .984584_wp*(i-1)
         l1 = 0
         do j = 1 , mm
            main : block
               if ( i/=1 .or. j/=1 ) then
                  x = x1
                  do
                     call bjndd(i,x,bj,dj,fj)
                     x0 = x
                     x = x - dj(i)/fj(i)
                     if ( x1>xm ) exit main
                     if ( abs(x-x0)<=tol) exit
                  end do
               endif
               l1 = l1 + 1
               n1(l1) = i - 1
               m1(l1) = j
               if ( i==1 ) m1(l1) = j - 1
               p1(l1) = 1
               zoc(l1) = x
               if ( i<=15 ) then
                  x1 = x + 3.057_wp + .0122_wp*(i-1) + (1.555_wp+.41575_wp*(i-1))/(j+1)**2
               else
                  x1 = x + 2.918_wp + .01924_wp*(i-1) + (6.26_wp+.13205_wp*(i-1))/(j+1)**2
               endif
            end block main
            x = x2
            do
               call bjndd(i,x,bj,dj,fj)
               x0 = x
               x = x - bj(i)/dj(i)
               if ( x<=xm ) exit
               if ( abs(x-x0)<=tol ) then
                  l1 = l1 + 1
                  n1(l1) = i - 1
                  m1(l1) = j
                  p1(l1) = 0
                  zoc(l1) = x
                  if ( i<=15 ) then
                     x2 = x + 3.11_wp + .0138_wp*(i-1) + (.04832_wp+.2804_wp*(i-1))/(j+1)**2
                  else
                     x2 = x + 3.001_wp + .0105_wp*(i-1) + (11.52_wp+.48525_wp*(i-1))/(j+3)**2
                  endif
                  exit
               end if
            end do
         enddo
         l = l0 + l1
         l2 = l
         do
            if ( l0==0 ) then
               do k = 1 , l
                  Zo(k) = zoc(k)
                  n(k) = n1(k)
                  m(k) = m1(k)
                  p(k) = p1(k)
               enddo
               l1 = 0
            elseif ( l0/=0 ) then
               if ( Zo(l0)>=zoc(l1) ) then
                  Zo(l0+l1) = Zo(l0)
                  n(l0+l1) = n(l0)
                  m(l0+l1) = m(l0)
                  p(l0+l1) = p(l0)
                  l0 = l0 - 1
               else
                  Zo(l0+l1) = zoc(l1)
                  n(l0+l1) = n1(l1)
                  m(l0+l1) = m1(l1)
                  p(l0+l1) = p1(l1)
                  l1 = l1 - 1
               endif
            endif
            if ( l1==0 ) exit
         end do
         l0 = l2
      enddo
   end subroutine jdzo