!>
!> @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