bessel.f90

Go to the documentation of this file.

!
!Authors: Jean-Luc Guermond, Raphael Laguerre, Caroline Nore, Copyright 2005
!
MODULE bessel

  PUBLIC :: bessk, bessk0, bessk1
  PUBLIC :: bessi, bessi0, bessi1
  PUBLIC :: bessj0, bessj1
  PRIVATE

CONTAINS
  ! ----------------------------------------------------------------------

  FUNCTION bessk(N,X) RESULT(vv)
    REAL *8 x,vv,tox,bk,bkm,bkp
    INTEGER :: n, j
    ! ------------------------------------------------------------------------
    !     CE SOUS-PROGRAMME CALCULE LA FONCTION BESSEL MODIFIFIEE 3E ESPECE
    !     D'ORDRE N ENTIER POUR TOUT X REEL POSITIF > 0.  ON UTILISE ICI LA
    !     FORMULE DE RECURRENCE CLASSIQUE EN PARTANT DE BESSK0 ET BESSK1.
    !
    !     THIS ROUTINE CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD
    !     KIND OF INTEGER ORDER, N FOR ANY POSITIVE REAL ARGUMENT, X. THE
    !     CLASSICAL RECURSION FORMULA IS USED, STARTING FROM BESSK0 AND BESSK1.
    ! ------------------------------------------------------------------------ 
    !     REFERENCE:
    !     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
    !     MATHEMATICAL TABLES, VOL.5, 1962.
    ! ------------------------------------------------------------------------
    IF (n.EQ.0) THEN
       vv = bessk0(x)
       RETURN
    ENDIF
    IF (n.EQ.1) THEN
       vv = bessk1(x)
       RETURN
    ENDIF
    IF (x.EQ.0.d0) THEN
       vv = 1.d30
       RETURN
    ENDIF
    tox = 2.d0/x
    bk  = bessk1(x)
    bkm = bessk0(x)
    DO 11 j=1,n-1
       bkp = bkm+dfloat(j)*tox*bk
       bkm = bk
       bk  = bkp
11     CONTINUE
       vv = bk
       RETURN
     END FUNCTION bessk

     ! ----------------------------------------------------------------------
     FUNCTION bessk0(X) RESULT(vv)
       !     CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0
       !     POUR TOUT X REEL NON NUL.
       !
       !     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
       !     ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X.
       ! ----------------------------------------------------------------------
       REAL*8 x,vv,y,ax,p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
       DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0, &
            0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/
       DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1, & 
            -0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/
       IF(x.EQ.0.d0) THEN
          vv=1.d30
          RETURN
       ENDIF
       IF(x.LE.2.d0) THEN
          y=x*x/4.d0
          ax=-log(x/2.d0)*bessi0(x)
          vv=ax+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
       ELSE
          y=(2.d0/x)
          ax=exp(-x)/sqrt(x)
          vv=ax*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))))
       ENDIF
       RETURN
     END FUNCTION bessk0
     ! ----------------------------------------------------------------------
     FUNCTION bessk1(X) RESULT(vv)
       !     CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1
       !     POUR TOUT X REEL POSITF NON NUL.
       !
       !     CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
       !     ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X.
       ! ----------------------------------------------------------------------
       REAL*8 x,vv,y,ax,p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
       DATA p1,p2,p3,p4,p5,p6,p7/1.d0,0.15443144d0,-0.67278579d0,  &
            -0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
       DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1, &
            0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
       IF(x.EQ.0.d0) THEN
          vv=1.d32
          RETURN
       ENDIF
       IF(x.LE.2.d0) THEN
          y=x*x/4.d0
          ax=log(x/2.d0)*bessi1(x)
          vv=ax+(1.d0/x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
       ELSE
          y=(2.d0/x)
          ax=exp(-x)/sqrt(x)
          vv=ax*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))))
       ENDIF
       RETURN
     END FUNCTION bessk1


     FUNCTION bessi(N,X) RESULT(vv)
       !
       !     This subroutine calculates the first kind modified Bessel function
       !     of integer order N, for any REAL X. We use here the classical
       !     recursion formula, when X > N. For X < N, the Miller's algorithm
       !     is used to avoid overflows. 
       !     REFERENCE:
       !     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
       !     MATHEMATICAL TABLES, VOL.5, 1962.
       !
       IMPLICIT NONE
       INTEGER, PARAMETER      :: iacc = 40
       REAL(KIND=8), PARAMETER :: bigno = 1.d10, bigni = 1.d-10
       REAL(KIND=8) ::  x,vv,tox,bim,bi,bip
       INTEGER :: n, m, j
       IF (n.EQ.0) THEN
          vv = bessi0(x)
          RETURN
       ENDIF
       IF (n.EQ.1) THEN
          vv = bessi1(x)
          RETURN
       ENDIF
       IF(x.EQ.0.d0) THEN
          vv=0.d0
          RETURN
       ENDIF
       tox = 2.d0/x
       bip = 0.d0
       bi  = 1.d0
       vv = 0.d0
       m = 2*((n+int(sqrt(float(iacc*n)))))
       DO 12 j = m,1,-1
          bim = bip+dfloat(j)*tox*bi
          bip = bi
          bi  = bim
          IF (abs(bi).GT.bigno) THEN
             bi  = bi*bigni
             bip = bip*bigni
             vv = vv*bigni
          ENDIF
          IF (j.EQ.n) vv = bip
12        CONTINUE
          vv = vv*bessi0(x)/bi
          RETURN
        END FUNCTION bessi
        ! ----------------------------------------------------------------------
        ! ----------------------------------------------------------------------
        ! Auxiliary Bessel functions for N=0, N=1
        FUNCTION bessi0(X) RESULT(vv)
          IMPLICIT NONE
          REAL(KIND=8) ::  x,vv,y,p1,p2,p3,p4,p5,p6,p7,  &
               q1,q2,q3,q4,q5,q6,q7,q8,q9,ax,bx
          DATA p1,p2,p3,p4,p5,p6,p7/1.d0,3.5156229d0,3.0899424d0,1.2067429d0,  &
               0.2659732d0,0.360768d-1,0.45813d-2/
          DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, &
               0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,  &
               0.2635537d-1,-0.1647633d-1,0.392377d-2/
          IF(abs(x).LT.3.75d0) THEN
             y=(x/3.75d0)**2
             vv=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
          ELSE
             ax=abs(x)
             y=3.75d0/ax
             bx=exp(ax)/sqrt(ax)
             ax=q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))
             vv=ax*bx
          ENDIF
          RETURN
        END FUNCTION bessi0
        ! ----------------------------------------------------------------------
        FUNCTION bessi1(X) RESULT(vv)
          IMPLICIT NONE
          REAL(KIND=8) :: x,vv,y,p1,p2,p3,p4,p5,p6,p7,  &
               q1,q2,q3,q4,q5,q6,q7,q8,q9,ax,bx
          DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,  &
               0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
          DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1, &
               -0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1, &
               -0.2895312d-1,0.1787654d-1,-0.420059d-2/
          IF(abs(x).LT.3.75d0) THEN
             y=(x/3.75d0)**2
             vv=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
          ELSE
             ax=abs(x)
             y=3.75d0/ax
             bx=exp(ax)/sqrt(ax)
             ax=q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))
             vv=ax*bx
          ENDIF
          RETURN
        END FUNCTION bessi1
        !-----------------------------------------------------------------------
        FUNCTION bessj0(x) RESULT(vv)
          IMPLICIT NONE
          REAL(KIND=8) :: x, y, ax, xx,vv , z
          REAL(KIND=8) ::  p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6
          DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, &
               -.2073370639d-5,.2093887211d-6/
          DATA q1,q2,q3,q4,q5/-.1562499995d-1, &
               .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/

          DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0, &
               651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/
          DATA s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,  &
               59272.64853d0,267.8532712d0,1.d0/
          !
          IF(abs(x).LT.8.d0)THEN
             y=x**2
             vv=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*  &
                  (s4+y*(s5+y*s6)))))
          ELSE
             ax=abs(x)
             z=8./ax
             y=z**2
             xx=ax-.785398164
             vv=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*  &
                  p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
          ENDIF
          RETURN
        END FUNCTION bessj0
        !-----------------------------------------------------------------------

        FUNCTION bessj1(x) RESULT(vv)
          IMPLICIT NONE
          REAL(KIND=8) :: x, y, ax, xx,vv , z
          REAL(KIND=8) ::  p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6
          !April 18 2008, JLG
          REAL(KIND=8) :: one
          DATA one/1.d0/
          !April 18 2008
          DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0, &
               242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/
          DATA s1,s2,s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0, &
               99447.43394d0,376.9991397d0,1.d0/
          DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, &
               .2457520174d-5,-.240337019d-6/ 
          DATA  q1,q2,q3,q4,q5/.04687499995d0, &
               -.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
          !
          IF(abs(x).LT.8.)THEN
             y=x**2
             vv=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ &
                  y*(s4+y*(s5+y*s6)))))
          ELSE
             ax=abs(x)
             z=8./ax
             y=z**2
             xx=ax-2.356194491
             !April 18 2008: 1.d0 -> one
             vv=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* &
                  p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(one,x)
             !April 18 2008
          ENDIF
          RETURN
        END FUNCTION bessj1


      END MODULE bessel