Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/math_util.f90
93a1e9712d04feature/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/>.
!>
!> @authors
!> (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
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!!$ PURE FUNCTION BESSEL_JN(n,x) RESULT(bessj)
!!$ DOUBLE PRECISION, INTENT(in) :: x
!!$ DOUBLE PRECISION :: bessj
!!$ DOUBLE PRECISION BIGNO,BIGNI
!!$ INTEGER n,IACC
!!$ PARAMETER (IACC=40,BIGNO=1.d10,BIGNI=1.d-10)
!!$ INTEGER j,jsum,m
!!$ DOUBLE PRECISION ax,bj,bjm,bjp,sum,tox,bessj0,bessj1
!!$ IF( n.EQ.0 ) THEN
!!$ bessj = bessj0(x)
!!$ RETURN
!!$ ELSE IF( n.EQ.1 ) THEN
!!$ bessj = bessj1(x)
!!$ RETURN
!!$ ENDIF
!!$ ax=ABS(x)
!!$ IF(ax.EQ.0.d0)THEN
!!$ bessj=0.d0
!!$ ELSE IF(ax.GT.float(n))THEN
!!$ tox=2./ax
!!$ bjm=bessj0(ax)
!!$ bj=bessj1(ax)
!!$ DO j=1,n-1
!!$ bjp=j*tox*bj-bjm
!!$ bjm=bj
!!$ bj=bjp
!!$ END DO
!!$ bessj=bj
!!$ ELSE
!!$ tox=2./ax
!!$ m=2*((n+INT(SQRT(float(IACC*n))))/2)
!!$ bessj=0.d0
!!$ jsum=0
!!$ sum=0.d0
!!$ bjp=0.d0
!!$ bj=1.
!!$ DO j=m,1,-1
!!$ bjm=j*tox*bj-bjp
!!$ bjp=bj
!!$ bj=bjm
!!$ IF(ABS(bj).GT.BIGNO)THEN
!!$ bj=bj*BIGNI
!!$ bjp=bjp*BIGNI
!!$ bessj=bessj*BIGNI
!!$ sum=sum*BIGNI
!!$ ENDIF
!!$ IF(jsum.NE.0)sum=sum+bj
!!$ jsum=1-jsum
!!$ IF(j.EQ.n)bessj=bjp
!!$ END DO
!!$ sum=2.*sum-bj
!!$ bessj=bessj/sum
!!$ ENDIF
!!$ IF(x.LT.0.d0.AND.MOD(n,2).EQ.1)bessj=-bessj
!!$ RETURN
!!$ END FUNCTION bessel_jn
!!$
!!$ PURE FUNCTION bessj0(x)
!!$ DOUBLE PRECISION, INTENT(in) :: x
!!$ DOUBLE PRECISION bessj0
!!$ DOUBLE PRECISION ax,xx,z
!!$ DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y
!!$ SAVE 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.)THEN
!!$ y=x**2
!!$ bessj0=(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
!!$ bessj0=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
!!$
!!$ PURE FUNCTION bessj1(x)
!!$ DOUBLE PRECISION, INTENT(in) :: x
!!$ DOUBLE PRECISION bessj1
!!$ DOUBLE PRECISION ax,xx,z
!!$ DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y
!!$ SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6
!!$ 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
!!$ bessj1=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.d0/ax
!!$ y=z**2
!!$ xx=ax-2.356194491d0
!!$ bessj1=SQRT(.636619772d0/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(1.d0,x)
!!$ ENDIF
!!$ RETURN
!!$ END FUNCTION bessj1
END MODULE math_util
c4science · Help