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(8),    INTENT(INOUT) :: x(ndim)

    INTEGER :: it,lll,l,n,il
    REAL(8) :: z,bs,bs0,bs1,xx,z2,arg,sum
    REAL(8), 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
    INTEGER, INTENT(IN) :: n
    REAL(8), INTENT(IN) :: h,f1(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

  SUBROUTINE linsol(a,b,kk)
    REAL(8), INTENT(INOUT) :: a(:,:),b(:)
    INTEGER, INTENT(IN) :: kk

    REAL(8) :: d,s,r
    INTEGER :: kkm,i,j,k,l,ipo,n,kmo

    d = 1.00000d0 
    kkm=kk-1
    IF (kkm == 0) THEN
       b(1)=b(1)/a(1,1)
    ELSE IF (kkm > 0) THEN
       DO i=1, kkm
          s = 0.0d0
          DO j=i,kk
             r=ABS(a(j,i))
             IF(r >  s) THEN
                s=r
                l=j
             ENDIF
          ENDDO
          IF(l /= i) THEN
             DO j=i,kk
                s=a(i,j)
                a(i,j)=a(l,j)
                a(l,j)=s
             ENDDO
             s=b(i)
             b(i)=b(l)
             b(l)=s
             d = -d
          ENDIF
          IF (a(i,i) /= 0.0d0) THEN
             ipo=i+1
             DO j=ipo,kk
                IF (a(j,i) /= 0.0d0) THEN
                   s=a(j,i)/a(i,i)
                   a(j,i) = 0.0d0
                   DO k=ipo,kk
                      a(j,k)=a(j,k)-a(i,k)*s
                   ENDDO
                   b(j)=b(j)-b(i)*s
                ENDIF
             ENDDO
          ENDIF
       ENDDO
       DO i=1,kk
          d=d*a(i,i)
       ENDDO
       kmo=kk-1
       b(kk)=b(kk)/a(kk,kk)
       DO i=1,kmo
          n=kk-i
          DO j=n,kmo
             b(n)=b(n)-a(n,j+1)*b(j+1)
          ENDDO
          b(n)=b(n)/a(n,n)
       ENDDO
    ENDIF
    IF(ABS(d).LT.1.d-10) WRITE(6,*) '**warning from linsol --',&
         'determinant too small --',d
  END SUBROUTINE linsol

  SUBROUTINE minverse(a,kk)
    REAL(8), INTENT(INOUT) :: a(:,:)
    INTEGER, INTENT(IN) :: kk

    REAL(8) :: d,s,r
    REAL(8), allocatable :: ai(:,:)
    INTEGER :: kkm,i,j,k,l,ipo,n,kmo

    allocate(ai(kk,kk),stat=i)
         if (i/=0) then
             write(6,*) 'Allocation error in minverse',i,kk
             stop
         endif

    ai=0
    do i=1,kk
       ai(i,i)=1
    enddo

    d = 1.00000d0 
    kkm=kk-1
    IF (kkm == 0) THEN
       ai(1,1)=ai(1,1)/a(1,1)
    ELSE IF (kkm > 0) THEN
       DO i=1, kkm
          s = 0.0d0
          DO j=i,kk
             r=ABS(a(j,i))
             IF(r >  s) THEN
                s=r
                l=j
             ENDIF
          ENDDO
          IF(l /= i) THEN
             DO j=i,kk
                s=a(i,j)
                a(i,j)=a(l,j)
                a(l,j)=s
             ENDDO
             Do k=1,kk
                s=ai(i,k)
                ai(i,k)=ai(l,k)
                ai(l,k)=s
             Enddo
             d = -d
          ENDIF
          IF (a(i,i) /= 0.0d0) THEN
             ipo=i+1
             DO j=ipo,kk
                IF (a(j,i) /= 0.0d0) THEN
                   s=a(j,i)/a(i,i)
                   a(j,i) = 0.0d0
                   DO k=ipo,kk
                      a(j,k)=a(j,k)-a(i,k)*s
                   ENDDO
                   Do k=1,kk
                      ai(j,k)=ai(j,k)-ai(i,k)*s
                   enddo
                ENDIF
             ENDDO
          ENDIF
       ENDDO
       DO i=1,kk
          d=d*a(i,i)
       ENDDO
       kmo=kk-1
       do k=1,kk
          ai(kk,k)=ai(kk,k)/a(kk,kk)
          DO i=1,kmo
             n=kk-i
             DO j=n,kmo
                ai(n,k)=ai(n,k)-a(n,j+1)*ai(j+1,k)
             ENDDO
             ai(n,k)=ai(n,k)/a(n,n)
          ENDDO
       enddo
    ENDIF
    a=ai
    deallocate(ai)

    IF(ABS(d).LT.1.d-10) WRITE(6,*) '**warning from linsol --',&
         'determinant too small --',d
  END SUBROUTINE minverse

  FUNCTION factorial(n)
    REAL(8) :: factorial
    INTEGER, INTENT(IN) :: n
    INTEGER :: i
    factorial=1
    IF (n.LT.2) RETURN
    DO i=2,n
       factorial=factorial*i
    ENDDO
  END FUNCTION factorial

  FUNCTION hwfn(nz,np,l,r)
    ! function to calculate the radial H wfn for nuclear charge nz
    !                                            principal qn   np
    !                                            orbital qn     l
    !   r*(radial H wfn) is returned
    REAL(8) :: hwfn
    REAL(8), INTENT(IN) :: r
    INTEGER, INTENT(IN) :: nz,np,l
    INTEGER :: node,k
    REAL(8) :: scale,rho,pref,term,sum
    node=np-l-1
    scale=2.d0*nz/np
    rho=scale*r
    pref=scale*SQRT(scale*factorial(np+l)/(2*np*factorial(node)))
    term=(rho**l)/factorial(2*l+1)
    sum=term
    IF (node.GT.0) THEN
       DO k=1,node
          term=-term*(node-k+1)*rho/(k*(2*l+1+k))
          sum=sum+term
       ENDDO
    ENDIF
    hwfn=r*pref*EXP(-0.5d0*rho)*sum
    !     write(6,*) 'r,hwfn=',r,hwfn
  END FUNCTION hwfn

    !***************************************************************************
    ! subroutine filter(n,func,small)
    !***************************************************************************
    SUBROUTINE filter(n,func,small)
      INTEGER, INTENT(IN) :: n
      REAL(8), INTENT(INOUT) :: func(:)
      REAL(8), INTENT(IN) :: small

      INTEGER :: i
      DO i=1,n
         IF (ABS(func(i)).LT.small) func(i)=0.d0
      ENDDO
    END SUBROUTINE filter

!******************************************************************************
!  From:
!  File : misc.f90
!    by : Alan Tackett
!    on : 10/17/95
!   for : Misc general purpose functions
!
!  This module contains misc. general purpose routines that don't fit
!  in another library. Below is a list of routines contained :
!
!  PrintDate(Unit, Text)
!      Prints the date and time to the specified unit along with the TEXT 
!
!*****************************************************************************

!******************************************************************************
!
!  PrintDate - Prints the date to the specified unit with TEXT prepended.
!
!    Unit - Output unit
!    Text - Text for prepending. All trailing blanks are removed. 
!
!******************************************************************************

Subroutine PrintDate(Unit, Text)
  Integer,       Intent(IN) :: Unit
  Character*(*), Intent(IN) :: Text

  Character*10 :: DateStr, TimeStr
  Character*50 :: FmtStr

  Call Date_And_Time(DateStr, TimeStr)

  FmtStr=' ' // DateStr(5:6) // '/' // DateStr(7:8) // '/' // DateStr(1:4) // &
        ', ' // TimeStr(1:2) // ':' // TimeStr(3:4) // ':' // &
          TimeStr(5:10)

  Write(Unit, '(a,a)')  Trim(Text), Trim(FmtStr)

  Return
End Subroutine


Function stripchar(inputchar)
    CHARACTER(132) :: stripchar
    CHARACTER*(*), INTENT(IN) :: inputchar

    INTEGER :: i,j,n
   
    n=LEN(inputchar)
   
    do i=1,132
      stripchar(i:i)=''
    enddo

    j=0
    do i=1,n
      if (inputchar(i:i) /= '') then
         j=j+1
         stripchar(j:j)=inputchar(i:i)
      endif
    enddo

End function stripchar

 Subroutine ConvertChar(inchar,outn)
    CHARACTER*(*),INTENT(IN) :: inchar
    INTEGER, INTENT(OUT) :: outn

    INTEGER :: i, j, k, n, fac

    n=LEN(inchar)

    fac=1; outn=0
    do i=n,1,-1
       If (inchar(i:i)==''.or.inchar(i:i)=="-") exit
       j=ichar(inchar(i:i))-48
       outn=outn+fac*j
       fac=fac*10
    enddo
 
    outn=outn
    write(6,*) 'exiting ConvertChar', inchar,outn
 end subroutine ConvertChar



END MODULE GlobalMath