RFLU_ModBessel.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Routines to compute Bessel functions, their derivatives, and their
!   zeros.
!
! Description: None.
!
! Notes: 
!   1. From: http://perso.wanadoo.fr/jean-pierre.moreau/f_bessel.html, 
!      labelled: F90 Release 1.0 By J-P Moreau, Paris. The original program
!      header included the following reference:
!      Ref.: www.esrf.fr/computing/expg/libraries/smf/PROGRAMS/MJYZO.FOR
!
! ******************************************************************************
!
! $Id: RFLU_ModBessel.F90,v 1.4 2008/12/06 08:44:20 mtcampbe Exp $
!
! ******************************************************************************

MODULE rflu_modbessel
      
  USE moddatatypes    
      
  IMPLICIT NONE
     
  PRIVATE
  PUBLIC :: rflu_jyndd, & 
            rflu_jyzo, & 
            rflu_jyzom

! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
   
  CHARACTER(CHRLEN) :: RCSIdentString = & 
    '$RCSfile: RFLU_ModBessel.F90,v $ $Revision: 1.4 $'        
             
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS





! ******************************************************************************
!
! Purpose:  Compute Bessel functions Jn(x) and Yn(x) and their first and 
!   second derivatives 
!
! Description: None.
!
! Input:
!   n           Order of Jn(x) and Yn(x)
!   x           Argument of Jn(x) and Yn(x) (x>0)
!
! Output:
!   BJN         Jn(x)
!   DJN         dJn/dx(x)
!   FJN         d2Jn/dx2(x)
!   BYN         Yn(x)
!   DYN         dYn/dx(x)
!   FYN         d2Yn/dx2(x)
!
! Notes:
!   1. Removed labelled loops and CONTINUE statements from original code. 
!
! ******************************************************************************
        
  SUBROUTINE rflu_jyndd(N,X,BJN,DJN,FJN,BYN,DYN,FYN)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================

    INTEGER :: n
    DOUBLE PRECISION, INTENT(IN) :: x
    DOUBLE PRECISION, INTENT(OUT) :: bjn,djn,fjn,byn,dyn,fyn

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: k,m,mt,nt
    DOUBLE PRECISION :: bs,ec,e0,f,f0,f1,su,s1
    DOUBLE PRECISION :: bj(102),by(102)

! ******************************************************************************
!   Start
! ******************************************************************************
                
    DO nt=1,900
       mt=int(0.5*log10(6.28*nt)-nt*log10(1.36*dabs(x)/nt))
       IF (mt.GT.20) THEN 
         EXIT
       END IF ! MT
    END DO ! NT    
    m=nt
    bs=0.0d0
    f0=0.0d0
    f1=1.0d-35
    su=0.0d0    
    DO k=m,0,-1
      f=2.0d0*(k+1.0d0)*f1/x-f0       
      IF (k.LE.n+1) THEN 
        bj(k+1)=f
      END IF ! K       
      IF (k.EQ.2*int(k/2)) THEN
         bs=bs+2.0d0*f
         IF (k.NE.0) THEN
           su=su+(-1)**(k/2)*f/k
         END IF ! K
      END IF       
      f0=f1 
      f1=f
    END DO ! K    
    DO k=0,n+1
      bj(k+1)=bj(k+1)/(bs-f)
    END DO ! K 
    bjn=bj(n+1)
    ec=0.5772156649015329d0
    e0=0.3183098861837907d0
    s1=2.0d0*e0*(dlog(x/2.0d0)+ec)*bj(1)
    f0=s1-8.0d0*e0*su/(bs-f)
    f1=(bj(2)*f0-2.0d0*e0/x)/bj(1)
    by(1)=f0
    by(2)=f1    
    DO k=2,n+1
       f=2.0d0*(k-1.0d0)*f1/x-f0
       by(k+1)=f
       f0=f1
       f1=f
    END DO ! K    
    byn=by(n+1)
    djn=-bj(n+2)+n*bj(n+1)/x
    dyn=-by(n+2)+n*by(n+1)/x
    fjn=(n*n/(x*x)-1.0d0)*bjn-djn/x
    fyn=(n*n/(x*x)-1.0d0)*byn-dyn/x

! ******************************************************************************
!   End
! ******************************************************************************

  END SUBROUTINE rflu_jyndd




! ******************************************************************************
!
! Purpose: Compute the zeros of Bessel functions Jn(x) and Yn(x) and their 
!   first derivatives
!
! Description: None.
!
! Input:
!   N           Order of Bessel functions (0 to 100)
!   NT          Number of zeros
!
! Output:
!   RJ0(L)      L-th zero of Jn(x), L=1,2,...,NT
!   RJ1(L)      L-th zero of dJn/dx(x), L=1,2,...,NT
!   RY0(L)      L-th zero of Yn(x), L=1,2,...,NT
!   RY1(L)      L-th zero of dYn/dx(x), L=1,2,...,NT
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_jyzo(N,NT,RJ0,RJ1,RY0,RY1)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================

    INTEGER :: n,nt
    DOUBLE PRECISION, INTENT(OUT) :: rj0(nt),rj1(nt),ry0(nt),ry1(nt)    

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: l
    DOUBLE PRECISION :: bjn,byn,djn,dyn,fjn,fyn,x,x0

! ******************************************************************************
!   Start
! ******************************************************************************
        
    IF (n.LE.20) THEN
       x=2.82141+1.15859*n
    ELSE
       x=n+1.85576*n**0.33333+1.03315/n**0.33333
    ENDIF
    l=0
10  x0=x
    CALL rflu_jyndd(n,x,bjn,djn,fjn,byn,dyn,fyn)
    x=x-bjn/djn
    IF (dabs(x-x0).GT.1.0d-9) go to 10
    l=l+1
    rj0(l)=x
    x=x+3.1416+(0.0972+0.0679*n-0.000354*n**2)/l
    IF (l.LT.nt) go to 10
    IF (n.LE.20) THEN
       x=0.961587+1.07703*n
    ELSE
       x=n+0.80861*n**0.33333+0.07249/n**0.33333
    ENDIF
    IF (n.EQ.0) x=3.8317
    l=0
15  x0=x
    CALL rflu_jyndd(n,x,bjn,djn,fjn,byn,dyn,fyn)
    x=x-djn/fjn
    IF (dabs(x-x0).GT.1.0d-9) go to 15
    l=l+1
    rj1(l)=x
    x=x+3.1416+(0.4955+0.0915*n-0.000435*n**2)/l
    IF (l.LT.nt) go to 15
    IF (n.LE.20) THEN
       x=1.19477+1.08933*n
    ELSE
       x=n+0.93158*n**0.33333+0.26035/n**0.33333
    ENDIF           
    l=0
20  x0=x
    CALL rflu_jyndd(n,x,bjn,djn,fjn,byn,dyn,fyn)
    x=x-byn/dyn
    IF (dabs(x-x0).GT.1.0d-9) go to 20
    l=l+1
    ry0(l)=x
    x=x+3.1416+(0.312+0.0852*n-0.000403*n**2)/l
    IF (l.LT.nt) go to 20
    IF (n.LE.20) THEN
       x=2.67257+1.16099*n
    ELSE
       x=n+1.8211*n**0.33333+0.94001/n**0.33333
    ENDIF  
    l=0
25  x0=x
    CALL rflu_jyndd(n,x,bjn,djn,fjn,byn,dyn,fyn)
    x=x-dyn/fyn
    IF (dabs(x-x0).GT.1.0d-9) go to 25
    l=l+1
    ry1(l)=x
    x=x+3.1416+(0.197+0.0643*n-0.000286*n**2)/l 
    IF (l.LT.nt) go to 25

! ******************************************************************************
!   End
! ******************************************************************************

  END SUBROUTINE rflu_jyzo
        
        
        
        
        
! ******************************************************************************
!
! Purpose: Compute Mth zero of Bessel functions Jn(x) and Yn(x) and their 
!   first derivatives
!
! Description: None.
!
! Input:
!   N           Order of Bessel functions (0 to 100)
!   M           Index of zero
!
! Output:
!   RJ0M        M-th zero of Jn(x)
!   RJ1M        M-th zero of dJn/dx(x)
!   RY0M        M-th zero of Yn(x)
!   RY1M        M-th zero of dYn/dx(x)
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_jyzom(N,M,RJ0M,RJ1M,RY0M,RY1M)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================

    INTEGER :: n,m
    DOUBLE PRECISION, INTENT(OUT) :: rj0m,rj1m,ry0m,ry1m    

! ==============================================================================
!   Locals
! ==============================================================================

    DOUBLE PRECISION :: bjn,byn,djn,dyn,fjn,fyn,x,x0
    DOUBLE PRECISION :: rj0(m),rj1(m),ry0(m),ry1(m) 

! ******************************************************************************
!   Start
! ******************************************************************************
        
    CALL rflu_jyzo(n,m,rj0,rj1,ry0,ry1)
        
    rj0m = rj0(m)
    rj1m = rj1(m)
    ry0m = ry0(m)
    ry1m = ry1(m)

! ******************************************************************************
!   End
! ******************************************************************************

  END SUBROUTINE rflu_jyzom       
        
     


! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modbessel


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModBessel.F90,v $
! Revision 1.4  2008/12/06 08:44:20  mtcampbe
! Updated license.
!
! Revision 1.3  2008/11/19 22:17:31  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.2  2006/04/07 15:19:18  haselbac
! Removed tabs
!
! Revision 1.1  2005/03/15 20:42:53  haselbac
! Initial revision
!
! ******************************************************************************