airyzo Subroutine

public subroutine airyzo(Nt, Kf, Xa, Xb, Xc, Xd)

Compute the first NT zeros of Airy functions Ai(x) and Ai'(x), a and a', and the associated values of Ai(a') and Ai'(a); and the first NT zeros of Airy functions Bi(x) and Bi'(x), b and b', and the associated values of Bi(b') and Bi'(b)

Arguments

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

Total number of zeros
integer, intent(in) :: Kf

Function code:

  • KF=1 for Ai(x) and Ai'(x)
  • KF=2 for Bi(x) and Bi'(x)
real(kind=wp), intent(out) :: Xa(Nt)

a, the m-th zero of Ai(x) or b, the m-th zero of Bi(x)
real(kind=wp), intent(out) :: Xb(Nt)

a', the m-th zero of Ai'(x) or b', the m-th zero of Bi'(x)
real(kind=wp), intent(out) :: Xc(Nt)

Ai(a') or Bi(b')
real(kind=wp), intent(out) :: Xd(Nt)

Ai'(a) or Bi'(b) ( The index is the serial number of zeros )

Calls

proc~~airyzo~~CallsGraph proc~airyzo specfun_module::airyzo proc~airyb specfun_module::airyb proc~airyzo->proc~airyb

Source Code

   subroutine airyzo(Nt,Kf,Xa,Xb,Xc,Xd)

      integer,intent(in) :: Nt !! Total number of zeros
      integer,intent(in) :: Kf !! Function code:
                               !!
                               !!  * KF=1 for Ai(x) and Ai'(x)
                               !!  * KF=2 for Bi(x) and Bi'(x)
      real(wp),intent(out) :: Xa(Nt) !! a, the m-th zero of Ai(x) or b, the m-th zero of Bi(x)
      real(wp),intent(out) :: Xb(Nt) !! a', the m-th zero of Ai'(x) or b', the m-th zero of Bi'(x)
      real(wp),intent(out) :: Xc(Nt) !! Ai(a') or Bi(b')
      real(wp),intent(out) :: Xd(Nt) !! Ai'(a) or Bi'(b) ( The index is the serial number of zeros )

      real(wp) :: ad , ai , bd , bi , err , rt , rt0 , u , u1 , x
      integer :: i

      rt = 0.0_wp
      do i = 1 , Nt
         rt0 = 0.0_wp
         if ( Kf==1 ) then
            u = 3.0_wp*pi*(4.0_wp*i-1)/8.0_wp
            u1 = 1/(u*u)
         elseif ( Kf==2 ) then
            if ( i==1 ) then
               rt0 = -1.17371d0
            else
               u = 3.0_wp*pi*(4.0_wp*i-3.0_wp)/8.0_wp
               u1 = 1.0_wp/(u*u)
            endif
         endif
         ! DLMF 9.9.18
         if ( rt0==0 ) rt0 = -(u*u)**(1.0_wp/3.0_wp) &
                             *(+1.0_wp+u1*(5.0_wp/48.0_wp+u1* &
                             (-5.0_wp/36.0_wp+u1*(77125.0_wp/82944.0_wp+ &
                             u1*(-108056875.0_wp/6967296.0_wp)))))
         do
            x = rt0
            call airyb(x,ai,bi,ad,bd)
            if ( Kf==1 ) rt = rt0 - ai/ad
            if ( Kf==2 ) rt = rt0 - bi/bd
            err = abs((rt-rt0)/rt)
            if ( err>1.0e-12_wp ) then
               rt0 = rt
            else
               Xa(i) = rt
               if ( err>1.0e-14_wp ) call airyb(rt,ai,bi,ad,bd)
               if ( Kf==1 ) Xd(i) = ad
               if ( Kf==2 ) Xd(i) = bd
               exit
            endif
         end do
      enddo
      do i = 1 , Nt
         rt0 = 0.0_wp
         if ( Kf==1 ) then
            if ( i==1 ) then
               rt0 = -1.01879_wp
            else
               u = 3.0_wp*pi*(4.0_wp*i-3.0_wp)/8.0_wp
               u1 = 1/(u*u)
            endif
         elseif ( Kf==2 ) then
            if ( i==1 ) then
               rt0 = -2.29444_wp
            else
               u = 3.0_wp*pi*(4.0_wp*i-1.0_wp)/8.0_wp
               u1 = 1/(u*u)
            endif
         endif
         ! DLMF 9.9.19
         if ( rt0==0 ) rt0 = -(u*u)**(1.0_wp/3.0_wp) &
                             *(+1.0_wp+u1*(-7.0_wp/48.0_wp+u1* &
                             (+35.0_wp/288.0_wp+u1*(-181223.0_wp/207360.0_wp+ &
                             u1*(18683371.0_wp/1244160.0_wp)))))
         do
            x = rt0
            call airyb(x,ai,bi,ad,bd)
            if ( Kf==1 ) rt = rt0 - ad/(ai*x)
            if ( Kf==2 ) rt = rt0 - bd/(bi*x)
            err = abs((rt-rt0)/rt)
            if ( err>1.0e-12_wp ) then
               rt0 = rt
            else
               Xb(i) = rt
               if ( err>1.0e-14_wp ) call airyb(rt,ai,bi,ad,bd)
               if ( Kf==1 ) Xc(i) = ai
               if ( Kf==2 ) Xc(i) = bi
               exit
            endif
         end do
      enddo

   end subroutine airyzo