Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/math_util.f90
1a965ce17d3bfeature/refactoring

Last Change
Raw File

math_util.f90
View Options

!>
!> @file math_util.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE math_util
!
! MATH_UTIL: Some math utilities.
!
! T.M. Tran, CRPP-EPFL
! December 2012
!
!
! Notes:
! - Assume the Fortran 2008 intrinsic BESSEL_JN(n,x) exists!
!
IMPLICIT NONE
DOUBLE PRECISION, PARAMETER :: pi=4.0d0*ATAN(1.0d0)
!
CONTAINS
ELEMENTAL FUNCTION bessjp(n,x)
!
! Derivative of J_n
!
DOUBLE PRECISION :: bessjp
INTEGER, INTENT(in) :: n
DOUBLE PRECISION, INTENT(in) :: x
!
IF(n.EQ.0) THEN
bessjp = -bessel_jn(1,x)
ELSE
bessjp = 0.5d0*(bessel_jn(n-1,x)-bessel_jn(n+1,x))
END IF
END FUNCTION bessjp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION root_bessj(n, s, info)
!
! s^th root of j_n
!
DOUBLE PRECISION :: root_bessj
INTEGER, INTENT(in) :: n
INTEGER, INTENT(in) :: s
INTEGER, OPTIONAL, INTENT(out) :: info
!
DOUBLE PRECISION :: b0,b1,b2,b3,b5,b7,t0,t1,t3,t5,t7,fn,fk,f1,f2,f3
DOUBLE PRECISION :: c1=1.8557571d0, c2=1.033150d0, c3=.00397d0, c4=.0908d0,&
& c5=.043d0
DOUBLE PRECISION :: zero
INTEGER :: iter
!
fn = REAL(ABS(n),8)
IF(s.EQ.1) THEN ! first zero
IF(n.EQ.0) THEN
zero = c1+c2-c3-c4+c5
ELSE
f1 = fn**(1.d0/3.d0)
f2 = f1*f1*fn
f3 = f1*fn*fn
zero = fn+c1*f1+(c2/f1)-(c3/fn)-(c4/f2)+(c5/f3)
END IF
ELSE ! Other zeros
t0 = 4.d0*fn*fn
t1 = t0-1.d0
t3 = 4.d0*t1*(7.d0*t0-31.d0)
t5 = 32.d0*t1*((83.d0*t0-982.d0)*t0+3779.d0)
t7 = 64.d0*t1*(((6949.d0*t0-153855.d0)*t0+1585743.d0)*t0 &
-6277237.d0)
fk = REAL(s,8)
!
b0 = (fk+.5d0*fn-.25d0)*pi! mac mahon's series for k>>n
b1 = 8.d0*b0
b2 = b1*b1
b3 = 3.d0*b1*b2
b5 = 5.d0*b3*b2
b7 = 7.d0*b5*b2
zero = b0-(t1/b1)-(t3/b3)-(t5/b5)-(t7/b7)
END IF
CALL newton(iter)
IF(PRESENT(info)) info = iter
root_bessj = zero
CONTAINS
SUBROUTINE newton(iter)
INTEGER, INTENT(out) :: iter
INTEGER :: itermx = 20
DOUBLE PRECISION :: dx, tol
tol = EPSILON(1.0d0)*zero
iter = 0
DO
iter = iter+1
dx = -bessel_jn(n,zero)/bessjp(n,zero)
zero = zero+dx
IF(iter.GE.itermx .OR. ABS(dx).LT.tol) EXIT
END DO
END SUBROUTINE newton
END FUNCTION root_bessj
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION root_bessjp(n, s, info)
!
! s^th root of derivative of j_n
!
DOUBLE PRECISION :: root_bessjp
INTEGER, INTENT(in) :: n
INTEGER, INTENT(in) :: s
INTEGER, OPTIONAL, INTENT(out) :: info
!
DOUBLE PRECISION :: c1=0.8086165D0, c2=0.072490D0, c3=.05097D0, c4=.0094D0
DOUBLE PRECISION :: b0,b1,b2,b3,b5,b7,t0,t1,t3,t5,t7,fn,fk,f1,f2
INTEGER :: iter
DOUBLE PRECISION :: zero
!
IF(n.EQ.0 .AND. s.EQ.1) THEN
root_bessjp = 0.0d0
IF(PRESENT(info)) info = 0
RETURN
END IF
!
fn = REAL(ABS(n),8)
fk = REAL(s,8)
!
IF(s.GT.1) THEN
!
! McMahon's series for s >> n
b0 = (fk+.5d0*fn-.75d0)*pi
b1 = 8.d0*b0
b2 = b1*b1
b3 = 3.d0*b1*b2
b5 = 5.d0*b3*b2
b7 = 7.d0*b5*b2
t0 = 4.d0*fn*fn
t1 = t0+3.d0
t3 = 4.d0*((7.d0*t0+82.d0)*t0-9.d0)
t5 = 32.d0*(((83.d0*t0+2075.d0)*t0-3039.d0)*t0+3537.d0)
t7 = 64.d0*((((6949.d0*t0+296492.d0)*t0-1248002.d0)*t0 &
+7414380.d0)*t0-5853627.d0)
zero = b0-(t1/b1)-(t3/b3)-(t5/b5)-(t7/b7)
ELSE
!
! Tchebychev's series for s <= n
f1 = fn**(1.d0/3.d0)
f2 = f1*f1*fn
zero = fn+c1*f1+(c2/f1)-(c3/fn)+(c4/f2)
END IF
!
CALL newton(iter)
root_bessjp = zero
IF(PRESENT(info)) info = iter
CONTAINS
SUBROUTINE newton(iter)
INTEGER, INTENT(out) :: iter
INTEGER :: itermx = 20
DOUBLE PRECISION :: dx, tol
tol = EPSILON(1.0d0)*zero
iter = 0
DO
iter = iter+1
dx = -bessel_jn(n,zero)/bessjp(n,zero)
dx = -2.0d0 * (bessel_jn(n-1,zero)-bessel_jn(n+1,zero)) / &
& (bessel_jn(n-2,zero)-2.d0*bessel_jn(n,zero)+&
& bessel_jn(n+2,zero))
zero = zero+dx
IF(iter.GE.itermx .OR. ABS(dx).LT.tol) EXIT
END DO
END SUBROUTINE newton
END FUNCTION root_bessjp
END MODULE math_util
c4science · Help