MODULE GlobalMath
  IMPLICIT NONE

  REAL(8) :: pi , machine_precision , machine_zero , machine_infinity

  REAL(8), PRIVATE :: minlog,maxlog,minexp,maxexp
  REAL(8), PRIVATE :: minlogarg,maxlogarg,minexparg,maxexparg

CONTAINS

  !****************************************************
  SUBROUTINE Init_GlobalConstants

    INTEGER :: i
    REAL(8)    :: tmp , a1,a2,a3

    ! Calculate machine accuracy
    machine_precision = 0
    a1 = 4.d0/3.d0
    DO WHILE (machine_precision == 0.d0)
       a2 = a1 - 1.d0
       a3 = a2 + a2 + a2
       machine_precision = ABS(a3 - 1.d0)
    ENDDO

    machine_zero= machine_precision**4
    machine_infinity = 1.d0/machine_zero

    pi = ACOS(-1.d0)

    minlogarg=machine_precision; minlog=LOG(minlogarg)
    maxlogarg=1.d0/machine_precision; maxlog=LOG(maxlogarg)
    minexparg=LOG(machine_precision);  minexp=0.d0
    maxexparg=-LOG(machine_precision);  maxexp=EXP(maxexparg)


    RETURN
  END SUBROUTINE Init_GlobalConstants

  !****************************************************
  FUNCTION ddlog(arg)
    REAL(8) :: arg,ddlog

    IF (arg>maxlogarg) THEN
       ddlog=maxlog
    ELSE IF (arg<minlogarg) THEN
       ddlog=minlog
    ELSE
       ddlog=LOG(arg)
    ENDIF

    RETURN
  END FUNCTION ddlog

  !****************************************************
  FUNCTION ddexp(arg)
    REAL(8) :: arg,ddexp

    IF (arg>maxexparg) THEN
       ddexp=maxexp
    ELSE IF (arg<minexparg) THEN
       ddexp=minexp
    ELSE
       ddexp=EXP(arg)
    ENDIF

    RETURN
  END FUNCTION ddexp

  !*****************************************************

  SUBROUTINE nderiv(h,y,z,ndim,ierr)
    INTEGER, INTENT(IN) :: ndim
    INTEGER, INTENT(INOUT) :: ierr
    REAL(8) , INTENT(IN) :: h,y(ndim)
    REAL(8) , INTENT(INOUT) :: z(ndim)
    !     subroutine ddet5(h,y,z,ndim)
    !      ssp routine modified by nawh 6/8/76

    REAL(8) :: hh,yy,a,b,c
    INTEGER :: i
    ierr=-1
    IF (ndim.LT.5) RETURN
    !        prepare differentiation loop
    hh=.08333333333333333d0/h
    yy=y(ndim-4)
    b=hh*(-25.d0*y(1)+48.d0*y(2)-36.d0*y(3)+16.d0*y(4)-3.d0*y(5))
    c=hh*(-3.d0*y(1)-10.d0*y(2)+18.d0*y(3)-6.d0*y(4)+y(5))
    !
    !        start differentiation loop
    DO  i=5,ndim
       a=b
       b=c
       c=hh*(y(i-4)-y(i)+8.d0*(y(i-1)-y(i-3)))
       z(i-4)=a
    ENDDO
    !        end of differentiation loop
    !
    !        normal exit
    a=hh*(-yy+6.d0*y(ndim-3)-18.d0*y(ndim-2)+10.d0*y(ndim-1)          &
         +3.d0*y(ndim))                                               
    z(ndim)=hh*(3.d0*yy-16.d0*y(ndim-3)+36.d0*y(ndim-2)               &
         -48.d0*y(ndim-1)+25.d0*y(ndim))
    z(ndim-1)=a
    z(ndim-2)=c
    z(ndim-3)=b
    !
    ierr=0
    RETURN
  END SUBROUTINE nderiv

  !**********************************************
  SUBROUTINE sphbes(lval,ndim,x)
    INTEGER, INTENT(IN)    :: lval
    INTEGER, INTENT(IN)    :: ndim
    REAL,    INTENT(INOUT) :: x(ndim)

    INTEGER :: it,lll,l,n,il
    REAL :: z,bs,bs0,bs1,xx,z2,arg,sum
    REAL, PARAMETER :: tol=1.e-11

    DO it=1,ndim
       z=x(it)

       IF(z.LE.0.5) THEN ! series expansion for small arguments
          IF (z.LE.tol) THEN ! - special treatment if z = 0
             bs=0.d0
             IF(lval.EQ.0) bs=1.d0

          ELSE
             lll=2*lval+1
             xx=1.d0
             DO l=1,lll,2
                xx=xx*l
             END DO

             z2=0.5d0*z*z
             lll=lll+2
             arg=-z2/lll
             sum=1.d0+arg

             DO  n=2,500
                lll=lll+2
                arg=-arg*z2/(n*lll)
                sum=sum+arg
                IF (ABS(arg).LE.tol) EXIT
             END DO

             bs=(z**lval)*sum/xx
          END IF
       ELSE        ! -- trigometric form
          bs=SIN(z)/z

          IF (lval.GT.0) THEN
             bs0=bs
             bs=(bs0-COS(z))/z
             IF (lval.GT.1) THEN
                bs1=bs

                DO  il=2,lval
                   bs=(2*il-1)*bs1/z-bs0
                   bs0=bs1
                   bs1=bs
                END DO
             END IF
          END IF
       END IF

       IF (ABS(bs)<1D-50) bs = 0
       x(it)=bs

    END DO

    RETURN
  END  SUBROUTINE sphbes

  !**********************************************************************
  FUNCTION overint(n,h,f1)
    !   function to calculate the integral of one vectors f1
    !      using simpsons rule, assuming f1(0)=0
    !      assumes f1 tabulated at points r=i*h; f1(1)=f1(i*h)
    REAL(8) :: overint
    REAL(8), INTENT(IN) :: h,f1(n)
    INTEGER, INTENT(IN) :: n

    INTEGER :: i,j
    REAL(8) :: fac

    IF (n==1) THEN
       overint=h*f1(1)/2
       RETURN
    ENDIF
    fac=h/3
    overint=4*f1(1)+f1(2)
    j=(n/2)*2
    IF (j>=4) THEN
       DO i=4,j,2
          overint=overint+f1(i-2)+4*f1(i-1)+f1(i)
       ENDDO
       overint=fac*overint
    ENDIF
    IF (n>j) THEN
       overint=overint+(f1(j)+f1(n))*h/2
    ENDIF
    RETURN
  END FUNCTION overint
  !*****************************************************************

  FUNCTION overlap(n,h,f1,f2)
    !   function to calculate the overlap between two vectors f1 and f2
    !      using simpsons rule, assuming f1(0)*f2(0)=0
    !      assumes f1, f2 tabulated at points r=i*h
    REAL(8) :: overlap
    INTEGER, INTENT(IN) :: n
    REAL(8), INTENT(IN) :: h, f1(n),f2(n)

    INTEGER :: i,j
    REAL(8) :: fac

    IF (n==1) THEN
       overlap=h*f1(1)*f2(1)/2
       RETURN
    ENDIF
    fac=h/3
    overlap=4*f1(1)*f2(1)+f1(2)*f2(2)
    j=(n/2)*2
    IF(j>=4) THEN
       DO i=4,j,2
          overlap=overlap+f1(i-2)*f2(i-2)+4*f1(i-1)*f2(i-1)+f1(i)*f2(i)
       ENDDO
       overlap=fac*overlap
    ENDIF
    IF (n>j) THEN
       overlap=overlap+h*(f1(j)*f2(j)+f1(n)*f2(n))/2
    ENDIF
    RETURN
  END FUNCTION overlap


  SUBROUTINE thomas(n,a,b,c,d)
    !  subroutine to use Thomas's algorithm to solve tridiagonal matrix
    !   Dale U. von Rosenberg, "Methods for the Numerical Solution of
    !     Partial Differential Equations, Am. Elsevier Pub., 1969, pg. 113
    !      a(i)*u(i-1)+b(i)*u(i)+c(i)*u(i+1)=d(i)
    !      a(1)=c(n)==0
    !      upon return, d(i)=u(i), a,b, & c modified
    INTEGER, INTENT(IN) :: n
    REAL(8), INTENT(INOUT) :: a(n),b(n),c(n),d(n)

    INTEGER :: i
    IF (n.LT.1) THEN
       WRITE(6,*) '***error in thomas -- n = ',n
       STOP
    ELSE IF (n.EQ.1) THEN
       d(1)=d(1)/b(1)
       RETURN
    ELSE
       DO i=2,n
          b(i)=b(i)-a(i)*c(i-1)/b(i-1)
       ENDDO
       d(1)=d(1)/b(1)
       DO i=2,n
          d(i)=(d(i)-a(i)*d(i-1))/b(i)
       ENDDO
       DO i=n-1,1,-1
          d(i)=d(i)-c(i)*d(i+1)/b(i)
       ENDDO
    ENDIF
    RETURN
  END SUBROUTINE thomas
  !**********************************************************************

  FUNCTION ranx()
    REAL(8) :: ranx
    INTEGER, PARAMETER :: konst=125
    INTEGER  :: m=100001
    m=m*konst
    m=m-2796203*(m/2796203)
    ranx=m/2796203.d0
    RETURN
  END FUNCTION ranx

  SUBROUTINE shift4(v1,v2,v3,v4,NEW)
    REAL(8), INTENT(IN) :: NEW
    REAL(8), INTENT(INOUT) :: v1,v2,v3,v4

    v1=v2
    v2=v3
    v3=v4
    v4=NEW

    RETURN
  END SUBROUTINE shift4

  SUBROUTINE mkname(i,stuff)
    !  subroutine to take an integer (.le. 4 digits) and return it
    !   in the form of a character string
    CHARACTER(4) stuff
    INTEGER i,i1,i10,i100,i1000
    stuff='?'
    IF (i.GT.10000) i=MOD(i,10000)
    i1000=i/1000
    i100=(i-1000*i1000)/100
    i10=(i-1000*i1000-100*i100)/10
    i1=(i-1000*i1000-100*i100-10*i10)
    IF (i.GE.1000) THEN
       stuff=CHAR(i1000+48)//CHAR(i100+48)//CHAR(i10+48)//CHAR(i1+48)
       RETURN
    ENDIF
    IF (i.GE.100) THEN
       stuff=CHAR(i100+48)//CHAR(i10+48)//CHAR(i1+48)
       RETURN
    ENDIF
    IF (i.GE.10) THEN
       stuff=CHAR(i10+48)//CHAR(i1+48)
       RETURN
    ENDIF
    IF (i.GE.0) stuff=CHAR(i1+48)
    RETURN
  END SUBROUTINE mkname
END MODULE GlobalMath