functions.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: functions.f90                                                     #
!#                                                                           #
!# Copyright (C) 2006-2013                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Björn Sperling <sperling@astrophysik.uni-kiel.de>                         #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program 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, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE functions
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  REAL, PARAMETER    :: pi = 3.14159265358979323846
  REAL, PARAMETER    :: sqrt_two = 1.41421356237309504880
  REAL, PARAMETER    :: eps_agm = epsilon(eps_agm)   ! precision of the AGM
  INTEGER, PARAMETER :: max_agm_iterations = 20      ! limit iteration steps
  ! store the first 20 coefficients used in computing the Legendre polynomials
  INTEGER            :: iii
  INTEGER, PARAMETER :: max_pl_coeff = 20
  REAL, DIMENSION(MAX_PL_COEFF), PARAMETER &
                     :: pl_coeff = (/ (iii/(iii+1.0), iii=1,max_pl_coeff) /)
  !--------------------------------------------------------------------------!
  ! exclude interface block from doxygen processing
  INTERFACE legendrepolynomial
     MODULE PROCEDURE legendrepolynomial_one, legendrepolynomial_all
  END INTERFACE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       asinh, acosh, atanh, &
       heaviside, &
       barrier, &
       normalrand, &
       ellipticintegrals, &
       ellipticintegral_k, &
       ellipticintegral_e, &
       incompleteellipticintegral_f, &
       incompleteellipticintegral_e, &
       errorfunc,&
       comperrorfunc,&
       inverf,&
       legendrepolynomials, &
       legendrepolynomial, &
       legendrefunction_qminhalf, &
       bessel_i0, &
       bessel_i1, &
       bessel_k0, &
       bessel_k0e, &
       bessel_k1, &
       lngamma, &
       ei, &
       gamma
  !--------------------------------------------------------------------------!

CONTAINS

  ! not included in FORTRAN before 2008 standard
  ELEMENTAL FUNCTION asinh(x) RESULT(fx)
    IMPLICIT NONE
    REAL, INTENT(IN) :: x
    REAL :: fx
    fx = log(x+sqrt(x*x+1.0))
  END FUNCTION asinh


  ! not included in FORTRAN before 2008 standard
  ELEMENTAL FUNCTION acosh(x) RESULT(fx)
    IMPLICIT NONE
    REAL, INTENT(IN) :: x
    REAL :: fx
    fx = log(x+sqrt(x*x-1.0))
  END FUNCTION acosh


  ! not included in FORTRAN before 2008 standard
  ELEMENTAL FUNCTION atanh(x) RESULT(fx)
    IMPLICIT NONE
    REAL, INTENT(IN) :: x
    REAL :: fx
    fx = 0.5*log((1.0+x)/(1.0-x))
  END FUNCTION atanh


  !      returns:   0   for x<a
  !                 1/2 for x=a
  !                 1   for x>a
  ELEMENTAL FUNCTION heaviside(x,a) RESULT(fx)
    IMPLICIT NONE
    REAL, INTENT(IN) :: x,a
    REAL :: fx
    fx = 0.5*(sign(1.0,x-a)+1.0)
  END FUNCTION heaviside

  !         returns    0   for x<a and x>b
  !                    1/2 for x=a and x=b
  !                    1   for a<x<b
  ELEMENTAL FUNCTION barrier(x,a,b) RESULT(fx)
    IMPLICIT NONE
    REAL, INTENT(IN) :: x,a,b
    REAL :: fx
    fx = heaviside(x,a)-heaviside(x,b)
  END FUNCTION barrier

  ELEMENTAL SUBROUTINE normalrand(u1,u2,x,y)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: u1,u2
    REAL, INTENT(OUT) :: x,y
    !------------------------------------------------------------------------!
    REAL    :: r,t
    !------------------------------------------------------------------------!
    r = sqrt(-2.0*log(u1))
    t = 2.0*pi*u2
    x = r*cos(t)
    y = r*sin(t)
  END SUBROUTINE normalrand


  ELEMENTAL SUBROUTINE ellipticintegrals(k,Kell,Eell)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: k ! elliptic modulus
    REAL, INTENT(OUT) :: eell,kell
    !------------------------------------------------------------------------!
    INTEGER           :: i,n
    REAL              :: an,bn,cn,tmp
    !------------------------------------------------------------------------!
    tmp = 0.0
    IF (abs(k).LT.1.0) THEN
       i   = 1
       an  = 1.0
       bn  = sqrt(1.0-k*k)
       cn  = k
!NEC$ UNROLL(20)
       DO n=1,max_agm_iterations ! do not loop forever
          tmp = tmp + cn*cn*i
          IF (abs(cn).GT.0.5*eps_agm*abs(an)) THEN
             cn = 0.5*(an-bn)
             bn = sqrt(an*bn)
             an = an-cn     ! = 0.5*(an+bn_old)
             i = ishft(i,1) ! = 2*i
#if !defined(NECSXAURORA)
          ELSE
             EXIT  ! exit prohibits vectorization
#endif
          END IF
       END DO
       ! Kell = 0.5*PI / an better: Kell = 0.5*PI/0.5*(an+bn)
       kell = pi / (an + bn + tiny(kell)) ! avoid division by 0
    ELSE
       kell = sqrt(-1.0*abs(k)) ! return NaN
    END IF
    eell = (1.0-0.5*tmp)*kell
  END SUBROUTINE ellipticintegrals


  ELEMENTAL FUNCTION ellipticintegral_k(k) RESULT(Kell)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: k ! elliptic modulus
    REAL :: kell
    !------------------------------------------------------------------------!
    REAL :: an,bn,tmp
    INTEGER :: i
    !------------------------------------------------------------------------!
    IF (abs(k).LT.1.0) THEN
       an = 1.0
       bn = sqrt(1.0-k*k) ! = NaN if |k| > 1
!NEC$ UNROLL(20)
       DO i=1,max_agm_iterations ! do not loop forever
          IF (abs(an-bn).GT.eps_agm*abs(an)) THEN
             tmp = 0.5*(an + bn)
             bn  = sqrt(an*bn)
             an  = tmp
#if !defined(NECSXAURORA)
          ELSE
             EXIT  ! exit prohibits vectorization
#endif
          END IF
       END DO
       ! Kell = 0.5*PI / an -> next iteration would compute an = 0.5*(an+bn)
       ! hence return an even better approximation:
       kell = pi / ( an + bn + tiny(kell) ) ! avoid division by 0
    ELSE
       kell = sqrt(-1.0*abs(k)) ! return NaN
    END IF
  END FUNCTION ellipticintegral_k


  ELEMENTAL FUNCTION ellipticintegral_e(k) RESULT(Eell)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: k ! elliptic modulus
    REAL             :: eell
    !------------------------------------------------------------------------!
    REAL :: dummy
    !------------------------------------------------------------------------!
    CALL ellipticintegrals(k,dummy,eell)
  END FUNCTION ellipticintegral_e


  ELEMENTAL FUNCTION rf(x,y,z) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL            :: x,y,z,res
    !------------------------------------------------------------------------!
    REAL, PARAMETER :: errtol = .0025
    REAL, PARAMETER :: third = 1./3.
    REAL, PARAMETER :: tin = 1.5e-38
    REAL, PARAMETER :: big = 3.e37
    REAL, PARAMETER :: c1 = 1./24.
    REAL, PARAMETER :: c2 = 0.1
    REAL, PARAMETER :: c3 = 3./44.
    REAL, PARAMETER :: c4 = 1./14.
    REAL            :: alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,&
                       yt,zt
    !------------------------------------------------------------------------!
    INTENT(IN)      :: x,y,z
    !------------------------------------------------------------------------!
    IF((min(x,y,z).LT.0.).OR.&
       (min(x+y,x+z,y+z).LT.tin).OR.&
       (max(x,y,z).GT.big)) THEN
      res = sqrt(-1.*abs(x)) ! return NAN
    ELSE
      xt = x
      yt = y
      zt = z
      DO
        sqrtx = sqrt(xt)
        sqrty = sqrt(yt)
        sqrtz = sqrt(zt)
        alamb = sqrtx * (sqrty+sqrtz) + sqrty*sqrtz
        xt = .25*(xt+alamb)
        yt = .25*(yt+alamb)
        zt = .25*(zt+alamb)
        ave = third*(xt+yt+zt)
        delx = (ave-xt)/ave
        dely = (ave-yt)/ave
        delz = (ave-zt)/ave
        IF(max(abs(delx),abs(dely),abs(delz)).LE.errtol) &
          EXIT
      END DO

      e2 = delx*dely-delz**2
      e3 = delx*dely*delz
      res = (1. + (c1*e2-c2-c3*e3)*e2 + c4*e3)/sqrt(ave)
    END IF
  END FUNCTION rf


  ELEMENTAL FUNCTION rd(x,y,z) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL            :: x,y,z,res
    !------------------------------------------------------------------------!
    REAL, PARAMETER :: errtol = .0015
    REAL, PARAMETER :: third = 1./3.
    REAL, PARAMETER :: tin = 1.5e-25
    REAL, PARAMETER :: big = 4.5e+21
    REAL, PARAMETER :: c1 = 3./14.
    REAL, PARAMETER :: c2 = 1./6.
    REAL, PARAMETER :: c3 = 9./22.
    REAL, PARAMETER :: c4 = 3./26.
    REAL, PARAMETER :: c5 = .25*c3
    REAL, PARAMETER :: c6 = 1.5*c4
    REAL            :: alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,&
                       sqrty,sqrtz,summ,xt,yt,zt
    !------------------------------------------------------------------------!
    INTENT(IN)      :: x,y,z
    !------------------------------------------------------------------------!
    IF((min(x,y).LT.0.).OR.&
       (min(x+y,z).LT.tin).OR.&
       (max(x,y,z).GT.big)) THEN
      res = sqrt(-1.*abs(x)) ! return NAN
    ELSE
      xt = x
      yt = y
      zt = z
      summ = 0.
      fac = 1.
      DO
        sqrtx = sqrt(xt)
        sqrty = sqrt(yt)
        sqrtz = sqrt(zt)
        alamb = sqrtx * (sqrty+sqrtz) + sqrty*sqrtz
        summ = summ + fac/(sqrtz*(zt+alamb))
        fac = .25*fac
        xt = .25*(xt+alamb)
        yt = .25*(yt+alamb)
        zt = .25*(zt+alamb)
        ave = .2*(xt+yt+3.*zt)
        delx = (ave-xt)/ave
        dely = (ave-yt)/ave
        delz = (ave-zt)/ave
        IF(max(abs(delx),abs(dely),abs(delz)).LE.errtol) &
          EXIT
      END DO
      ea = delx*dely
      eb = delz*delz
      ec = ea - eb
      ed = ea - 6. * eb
      ee = ed + ec + ec
      res = 3. * summ + fac*(1.+ed*(-c1+c5*ed-c6*delz*ee) &
        + delz*(c2*ee+delz*(-c3*ec+delz*c4*ea)))/(ave*sqrt(ave))
    END IF
  END FUNCTION rd

  ! Legendre elliptic integral of the first kind F(phi,k), evaluated using
  ! Carlson's function R_F. The argument ranges are 0 <= phi <= pi/2,
  ! 0 <= k*sin(phi) <= 1
  ! similar to Numerical recipes 2nd edition
  ELEMENTAL FUNCTION incompleteellipticintegral_f(phi,ak) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL            :: phi,ak,res
    !------------------------------------------------------------------------!
    REAL            :: s
    !------------------------------------------------------------------------!
    INTENT(IN)      :: phi,ak
    !------------------------------------------------------------------------!
    s = sin(phi)
    res = s*rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.)
  END FUNCTION incompleteellipticintegral_f


  ELEMENTAL FUNCTION incompleteellipticintegral_e(phi,ak) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL            :: phi,ak,res
    !------------------------------------------------------------------------!
    REAL            :: cc,q,s
    !------------------------------------------------------------------------!
    INTENT(IN)      :: phi,ak
    !------------------------------------------------------------------------!
    s = sin(phi)
    cc = cos(phi)**2
    q = (1.-s*ak)*(1.+s*ak)
    res = s*(rf(cc,q,1.)-((s*ak)**2)*rd(cc,q,1.)/3.)
  END FUNCTION incompleteellipticintegral_e


  PURE SUBROUTINE legendrepolynomials(l,x,P)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    INTEGER              :: l
    REAL                 :: x
    REAL, DIMENSION(0:l) :: p
    !------------------------------------------------------------------------!
    INTEGER              :: i
    !------------------------------------------------------------------------!
    INTENT(IN)           :: l,x
    INTENT(OUT)          :: p
    !------------------------------------------------------------------------!
    IF(l.LT.0) THEN
       p(:) = sqrt(1.0*l) ! return NaN
    ELSE
       p(0) = 1.0
       IF (l.GE.1) THEN
          p(1) = x
          ! use the constant coefficients for the first
          ! MAX_PL_COEFF Legendre Polynomials
          DO i=2,min(l,max_pl_coeff)
             p(i) = (1.0+pl_coeff(i))*x*p(i-1) - pl_coeff(i)*p(i-2)
          END DO
          ! from MAX_PL_COEFF+1 to l compute the coefficients i/(i+1);
          ! this is probably slower, because of the division
          DO i=max_pl_coeff+1,l
             ! recurrence formula for the Legendre polynomials:
             ! P(i) = ( (2*i+1)*x*P(i-1) - i*P(i-2) ) / (i+1)
             !      = (1+i/(i+1))*x*P(i-1) - i/(i+1)*P(i-2)
             p(i) = i/(i+1.0) ! temporary (is a number close to 1)
             p(i) = (1.0+p(i))*x*p(i-1) - p(i)*p(i-2)
          END DO
       END IF
    END IF
  END SUBROUTINE legendrepolynomials


  ELEMENTAL FUNCTION legendrepolynomial_one(l,x,Plminus1,Plminus2) RESULT (Pl)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    INTEGER    :: l
    REAL       :: x,plminus1,plminus2,pl
    !------------------------------------------------------------------------!
    REAL       :: c
    !------------------------------------------------------------------------!
    INTENT(IN) :: l,x,plminus1,plminus2
    !------------------------------------------------------------------------!
    IF(l.LT.0) THEN
       pl = sqrt(1.0*l) ! return NaN
    ELSE
       SELECT CASE(l)
       CASE(0)
          pl = 1.0
       CASE(1)
          pl = x
       CASE DEFAULT
          ! speed up things a little with predefined coefficients
          IF (l.LE.max_pl_coeff) THEN
             c = pl_coeff(l)
          ELSE
             c = l/(l+1.0) ! temporary storage
          END IF
          ! do recursion step
          pl = (1.0+c)*x*plminus1 - c*plminus2
       END SELECT
    END IF
  END FUNCTION legendrepolynomial_one


  ELEMENTAL FUNCTION legendrepolynomial_all(l,x) RESULT (Pl)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    INTEGER    :: l
    REAL       :: x,pl
    !------------------------------------------------------------------------!
    INTEGER    :: i
    REAL       :: plminus1,plminus2
    !------------------------------------------------------------------------!
    INTENT(IN) :: l,x
    !------------------------------------------------------------------------!
    IF(l.LT.0) THEN
        pl = sqrt(1.0*l) ! return NaN
    ELSE
!NEC$ UNROLL(20)
        DO i=0,l
           pl = legendrepolynomial_one(i,x,plminus1,plminus2)
           plminus2 = plminus1
           plminus1 = pl
        END DO
    END IF
  END FUNCTION legendrepolynomial_all

  ELEMENTAL FUNCTION legendrefunction_qminhalf(x) RESULT (q)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL       :: x,q
    !------------------------------------------------------------------------!
    INTEGER    :: i
    REAL       :: an,bn,tmp
    !------------------------------------------------------------------------!
    INTENT(IN) :: x
    !------------------------------------------------------------------------!
    IF (x.GE.1.0) THEN
       an = sqrt(x-1.0)
       bn = sqrt(x+1.0) ! if x < 1 bn = NaN, hence q = NaN
!NEC$ UNROLL(20)
       DO i=1,max_agm_iterations ! do not loop forever
          IF (abs(an-bn).GT.eps_agm*abs(an)) THEN
             tmp = 0.5*(an + bn)
             bn = sqrt(an*bn)
             an = tmp
#if !defined(NECSXAURORA)
          ELSE
             EXIT  ! exit prohibits vectorization
#endif
          END IF
       END DO
!       q = 0.5*PI*SQRT_TWO / an
       q = pi*sqrt_two / (an+bn+tiny(q))
    ELSE
       q = sqrt(-1.0*abs(x)) ! return NaN
    END IF
  END FUNCTION legendrefunction_qminhalf



  ELEMENTAL FUNCTION bessel_i0(x) RESULT(I0)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL                :: x
    REAL                :: i0
    !------------------------------------------------------------------------!
    INTENT(IN)          :: x
    !------------------------------------------------------------------------!
    REAL, DIMENSION(7), PARAMETER &
                        :: p = (/ 1.0d0, 3.5156229d0, 3.0899424d0, &
                                  1.2067492d0, 0.2659732d0, 0.360768d-1, &
                                  0.45813d-2 /)
    REAL, DIMENSION(9), PARAMETER &
                        :: q = (/ 0.39894228d0, 0.1328592d-1, 0.225319d-2, &
                                  -0.157565d-2, 0.916281d-2, -0.2057706d-1, &
                                  0.2635537d-1, -0.1647633d-1, 0.392377d-2 /)
    REAL                :: t, absx
    !------------------------------------------------------------------------!

    IF(abs(x).LT.3.75) THEN
        t = (x/3.75)**2
        i0 = p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7))))))
    ELSE
        absx = abs(x)
        t = 3.75/absx
        i0 = (exp(absx)/sqrt(absx)) &
             * (q(1)+t*(q(2)+t*(q(3)+t*(q(4) &
                + t*(q(5)+t*(q(6)+t*(q(7)+t*(q(8)+t*q(9)))))))))
    ENDIF

    END FUNCTION bessel_i0



  ELEMENTAL FUNCTION bessel_i1(x) RESULT(I1)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL                :: x
    REAL                :: i1
    !------------------------------------------------------------------------!
    INTENT(IN)          :: x
    !------------------------------------------------------------------------!
    REAL, DIMENSION(7), PARAMETER &
                        :: p = (/ 0.5d0, 0.87890594d0, 0.51498869d0, &
                                  0.15084934d0, 0.2658733d-1, 0.301532d-2, &
                                  0.32411d-3 /)
    REAL, DIMENSION(9), PARAMETER &
                        :: q = (/ 0.39894228d0, -0.3988024d-1, -0.362018d-2, &
                                  0.163801d-2, -0.1031555d-1, 0.2282967d-1, &
                                  -0.2895312d-1, 0.1787654d-1, -0.420059d-2 /)
    REAL                :: t, absx
    !------------------------------------------------------------------------!

    IF(abs(x).LT.3.75) THEN
        t = (x/3.75)**2
        i1 = x*(p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
    ELSE
        absx = abs(x)
        t = 3.75/absx
        i1 = (exp(absx)/sqrt(absx)) &
             * (q(1)+t*(q(2)+t*(q(3)+t*(q(4) &
                + t*(q(5)+t*(q(6)+t*(q(7)+t*(q(8)+t*q(9)))))))))
        IF(x.LT.0.) THEN
            i1 = -i1
        END IF

    ENDIF

    END FUNCTION bessel_i1



  ELEMENTAL FUNCTION bessel_k0(x) RESULT(K0)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL                :: x
    REAL                :: k0
    !------------------------------------------------------------------------!
    INTENT(IN)          :: x
    !------------------------------------------------------------------------!
    REAL, DIMENSION(7), PARAMETER &
                        :: p = (/ -0.57721566d0, 0.42278420d0, 0.23069756d0, &
                                  0.3488590d-1, 0.262698d-2, 0.10750d-3, &
                                  0.74d-5 /)
    REAL, DIMENSION(7), PARAMETER &
                        :: q = (/ 1.25331414d0, -0.7832358d-1, 0.2189568d-1, &
                                  -0.1062446d-1, 0.587872d-2, -0.251540d-2, &
                                  0.53208d-3 /)
    REAL                :: t
    !------------------------------------------------------------------------!

    IF(x.LE.2.0) THEN
            t = x*x / 4.0
            k0 = (-log(x/2.0)*bessel_i0(x)) &
                 + (p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
    ELSE
            t = (2.0/x)
            k0 = (exp(-x)/sqrt(x)) &
                 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4)+t*(q(5)+t*(q(6)+t*q(7)))))))
    ENDIF

    END FUNCTION bessel_k0

  ELEMENTAL FUNCTION bessel_k0e(x) RESULT(K0e)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL                :: x
    REAL                :: k0e
    !------------------------------------------------------------------------!
    INTENT(IN)          :: x
    !------------------------------------------------------------------------!
    REAL, DIMENSION(7), PARAMETER &
                        :: q = (/ 1.25331414, -0.07832358, 0.02189568,&
                                  -0.01062446, 0.00587872, -0.00251540,&
                                  0.00053208 /)
    REAL                :: t
    !------------------------------------------------------------------------!

    IF(x.LT.2.0) THEN
            k0e = exp(x) * bessel_k0(x)
    ELSE
            t = (2.0/x)
            k0e = (1.0/sqrt(x)) &
                 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4)+t*(q(5)+t*(q(6)+t*q(7)))))))
    ENDIF

    END FUNCTION bessel_k0e

  ELEMENTAL FUNCTION bessel_k1(x) RESULT(K1)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL                :: x
    REAL                :: k1
    !------------------------------------------------------------------------!
    INTENT(IN)          :: x
    !------------------------------------------------------------------------!
    REAL, DIMENSION(7), PARAMETER &
                        :: p = (/ 1.0d0, 0.15443144d0, -0.67278579d0, &
                                  -0.18156897d0, -0.1919402d-1, -0.110404d-2, &
                                  -0.4686d-4 /)
    REAL, DIMENSION(7), PARAMETER &
                        :: q = (/ 1.25331414d0, 0.23498619d0, -0.3655620d-1, &
                                  0.1504268d-1, -0.780353d-2, 0.325614d-2, &
                                  -0.68245d-3 /)
    REAL                :: t
    !------------------------------------------------------------------------!

    IF(x.LE.2.0) THEN
            t = x*x / 4.0
            k1 = (log(x/2.0)*bessel_i1(x)) &
                 + (1.0/x) &
                    * (p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
    ELSE
            t = (2.0/x)
            k1 = (exp(-x)/sqrt(x)) &
                 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4)+t*(q(5)+t*(q(6)+t*q(7)))))))
    ENDIF

    END FUNCTION bessel_k1


  ELEMENTAL FUNCTION ei(x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN):: x
    REAL            :: res
    !------------------------------------------------------------------------!
    REAL, PARAMETER :: eps=epsilon(x), euler=.577215664902, fpmin=tiny(x)/eps
    !------------------------------------------------------------------------!
    IF(x.LT.-5.) THEN
      res = continued_fraction_ei(x)
    ELSE IF(x.EQ.0.) THEN
      res = -huge(res)
    ELSE IF(x.LT.6.8) THEN
      res = power_series_ei(x)
    ELSE IF(x.LT.50.) THEN
      res = argument_addition_series_ei(x)
    ELSE
      res = continued_fraction_ei(x)
    END IF
    END FUNCTION ei


  ELEMENTAL FUNCTION continued_fraction_ei(x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN):: x
    REAL            :: res
    !------------------------------------------------------------------------!
    REAL            :: am1,a0,bm1,b0,a,b,ap1,bp1,eps
    INTEGER         :: j
    !------------------------------------------------------------------------!
    am1 = 1.
    a0 = 0.
    bm1 = 0.
    b0 = 1.
    a = exp(x)
    b = -x + 1.
    ap1 = b* a0 + a * am1
    bp1 = b* b0 + a * bm1
    eps = 10. * epsilon(x)
    j = 1
    a = 1.
    DO WHILE(abs(ap1 * b0 - a0 * bp1).GT.eps*abs(a0*bp1))
      IF(abs(bp1).GT.1.) THEN
        am1 = a0 / bp1
        a0 = ap1 / bp1
        bm1 = b0 / bp1
        b0 = 1.
      ELSE
        am1 = a0
        a0 = ap1
        bm1 = b0
        b0 = bp1
      END IF
      a = -j*j
      b = b + 2.
      ap1 = b * a0 + a * am1
      bp1 = b * b0 + a * bm1
      j = j + 1
    END DO
    res = -ap1 / bp1
  END FUNCTION continued_fraction_ei


  ELEMENTAL FUNCTION power_series_ei(x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN):: x
    REAL            :: res
    !------------------------------------------------------------------------!
    REAL            :: xn,sn,sm1,hsum,g,y,factorial,eps
    !------------------------------------------------------------------------!
    xn = -x
    sn = -x
    sm1 = 0.
    hsum = 1.
    g = 0.5772156649015328606065121
    y = 1.
    factorial = 1.
    eps = 10.*epsilon(x)
    IF(x.EQ.0.) &
      res = -huge(res)
    DO WHILE(abs(sn-sm1).GT.eps*abs(sm1))
      sm1 = sn
      y = y + 1.
      xn = xn * (-x)
      factorial = factorial * y
      hsum = hsum + 1. / y
      sn = sn + hsum * xn / factorial
    END DO
    res = g + log(abs(x)) - exp(x) * sn
  END FUNCTION power_series_ei

  ELEMENTAL FUNCTION argument_addition_series_ei(x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN):: x
    REAL            :: res
    !------------------------------------------------------------------------!
    REAL, DIMENSION(44), PARAMETER :: eic = &
    (/1.915047433355013959531e2,  4.403798995348382689974e2, &
      1.037878290717089587658e3,  2.492228976241877759138e3, &
      6.071406374098611507965e3,  1.495953266639752885229e4, &
      3.719768849068903560439e4,  9.319251363396537129882e4, &
      2.349558524907683035782e5,  5.955609986708370018502e5, &
      1.516637894042516884433e6,  3.877904330597443502996e6, &
      9.950907251046844760026e6,  2.561565266405658882048e7, &
      6.612718635548492136250e7,  1.711446713003636684975e8, &
      4.439663698302712208698e8,  1.154115391849182948287e9, &
      3.005950906525548689841e9,  7.842940991898186370453e9, &
      2.049649711988081236484e10, 5.364511859231469415605e10,&
      1.405991957584069047340e11, 3.689732094072741970640e11,&
      9.694555759683939661662e11, 2.550043566357786926147e12,&
      6.714640184076497558707e12, 1.769803724411626854310e13,&
      4.669055014466159544500e13, 1.232852079912097685431e14,&
      3.257988998672263996790e14, 8.616388199965786544948e14,&
      2.280446200301902595341e15, 6.039718263611241578359e15,&
      1.600664914324504111070e16, 4.244796092136850759368e16,&
      1.126348290166966760275e17, 2.990444718632336675058e17,&
      7.943916035704453771510e17, 2.111342388647824195000e18,&
      5.614329680810343111535e18, 1.493630213112993142255e19,&
      3.975442747903744836007e19, 1.058563689713169096306e20 /)
    INTEGER :: k, j
    REAL    :: xx,dx,xxj,edx,sm,sn,term,factorial,dxj,eps
    !------------------------------------------------------------------------!
    k = nint(x + 0.5)
    j = 0
    xx = k
    dx = x - xx
    xxj = xx
    edx = exp(dx)
    sm = 1.
    sn = (edx - 1.) / xxj
    term = huge(x)
    factorial = 1.
    dxj = 1.
    eps = 10.*epsilon(x)

    DO WHILE(abs(term).GT.eps*abs(sn))
      j = j + 1
      factorial = factorial * j
      xxj = xxj * xx
      dxj = dxj * (-dx)
      sm = sm + dxj / factorial
      term = (factorial * (edx * sm - 1.)) / xxj
      sn = sn + term
    END DO
    res = eic(k-6) + sn * exp(xx)
  END FUNCTION argument_addition_series_ei


   ELEMENTAL FUNCTION errorfunc(x) RESULT(erf)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: x
    REAL             :: erf
    !------------------------------------------------------------------------!
    REAL, PARAMETER  :: xerf  = 0.46875
    REAL, PARAMETER  :: xerfc = 4.0
    !------------------------------------------------------------------------!
    !Cody error function approximation
    IF (abs(x) .LT. tiny(1.) ) THEN !!0 besser abfangen
            erf = 0.0
    ELSE IF(x.GT.0.0)THEN
           IF (x .LT. xerf) THEN
               erf = codyerfapprox(x)
           ELSE IF (x .LT. xerfc) THEN
               erf = 1. - codycerf1approx(x)
           ELSE
               erf = 1. - codycerf2approx(x)
           END IF
    !use erf,erfc symmetry for negative arguments
    ELSE
           IF (-x .LT. xerf) THEN
               erf = - codyerfapprox(-x)
           ELSE IF (-x .LT. xerfc) THEN
               erf = codycerf1approx(-x) - 1.
           ELSE
               erf = codycerf2approx(-x) - 1.
           END IF
    END IF

  END FUNCTION errorfunc

  ! returns the complementary error function
  ELEMENTAL FUNCTION comperrorfunc(x) RESULT(cerf)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: x
    REAL             :: cerf
    !------------------------------------------------------------------------!
    REAL, PARAMETER  :: xerf  = 0.46875
    REAL, PARAMETER  :: xerfc = 4.0
    !------------------------------------------------------------------------!
    !Cody complementary error function approximation
    IF(x.GE.0.0) THEN
        IF (x .LT.xerf) THEN
            cerf = 1. - codyerfapprox(x)
        ELSE IF (x .LT. xerfc) THEN
            cerf = codycerf1approx(x)
        ELSE
            cerf = codycerf2approx(x)
        END IF
    !use erf,erfc symmetry for negative arguments
    ELSE
        IF (-x .LT. xerf) THEN
            cerf = 1. + codyerfapprox(-x)
        ELSE IF (-x .LT. xerfc) THEN
            cerf = 2. - codycerf1approx(-x)
        ELSE
            cerf = 2. - codycerf2approx(-x)
        END IF
    END IF

  END FUNCTION comperrorfunc


  ELEMENTAL FUNCTION codyerfapprox(x) RESULT(erf)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: x
    REAL              :: erf
    !------------------------------------------------------------------------!
    REAL :: denom, nom, x2, x4, x6
    REAL,DIMENSION(4) :: p_denom
    REAL,DIMENSION(5) :: p_nom
    !------------------------------------------------------------------------!
    p_nom   = (/ 3.20937758913846947e03 , 3.77485237685302021e02 ,&
                 1.13864154151050156e02 , 3.16112374387056560    ,&
                 1.85777706184603153e-1   /)
    p_denom = (/ 2.84423683343917062e03 , 1.28261652607737228e03 ,&
                 2.44024637934444173e02 , 2.36012909523441209e01 /)

    x2 = x  * x
    x4 = x2 * x2
    x6 = x4 * x2

    nom   = p_nom(1)  + x2 * (p_nom(2)   + p_nom(3)* x2   + p_nom(4) *   x4 &
            + p_nom(5)*x6)
    denom = p_denom(1)+ x2 * (p_denom(2) + p_denom(3)* x2 + p_denom(4) * x4 &
            + x6)

    erf   = x * nom / denom

  END FUNCTION codyerfapprox

  ELEMENTAL FUNCTION codycerf1approx(x) RESULT(cerf)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: x
    REAL              :: cerf
    !------------------------------------------------------------------------!
    INTEGER           :: i
    REAL              :: denom, nom
    REAL,DIMENSION(9) :: p_nom , p_denom
    !------------------------------------------------------------------------!
     p_nom   = (/ 1.23033935479799725e03 , 2.05107837782607147e03 ,&
                  1.71204761263407058e03 , 8.81952221241769090e02 ,&
                  2.98635138197400131e02 , 6.61191906371416295e01 ,&
                  8.88314979438837594e0  , 5.64188496988670089e-1 ,&
                  2.15311535474403846e-8 /)
     p_denom = (/ 1.23033935480374942e03 , 3.43936767414372164e03 ,&
                  4.36261909014324716e03 , 3.29079923573345963e03 ,&
                  1.62138957456669019e03 , 5.37181101862009858e02 ,&
                  1.17693950891312499e02 , 1.57449261107098347e01 ,&
                  1.0 /)

     nom   = 0.0
     denom = 0.0

     DO i= 0,7
        nom   = (nom   + p_nom(9-i))   * x
        denom = (denom + p_denom(9-i)) * x
     END DO

     nom   = nom   + p_nom(1)
     denom = denom + p_denom(1)

     cerf = exp(-x*x)* nom / denom

  END FUNCTION codycerf1approx

  ELEMENTAL FUNCTION codycerf2approx(x) RESULT(cerf)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: x
    REAL              :: cerf
    !------------------------------------------------------------------------!
    REAL              :: denom, nom, x2, x4, x6, x8
    REAL,DIMENSION(5) :: p_denom
    REAL,DIMENSION(6) :: p_nom
    !------------------------------------------------------------------------!
    p_nom   = (/ 6.58749161529837803e-4 , 1.60837851487422766e-2 ,&
                 1.25781726111229246e-1 , 3.60344899949804439e-1 ,&
                 3.05326634961232344e-1 , 1.63153871373020978e-2 /)
    p_denom = (/ 2.33520497626869185e-3 , 6.05183413124413191e-2 ,&
                 5.27905102951428412e-1 , 1.87295284992346047e00 ,&
                 2.56852019228982242e00 /)

    x2 = x  * x
    x4 = x2 * x2
    x6 = x4 * x2
    x8 = x4 * x4

    nom   = -x2 * p_nom(1) - p_nom(2) - p_nom(3)/x2 - p_nom(4)/(x4) - &
            p_nom(5)/(x6) - p_nom(6)/(x8)
    denom = x2 * p_denom(1) + p_denom(2) + p_denom(3)/x2 + p_denom(4)/(x4) + &
            p_denom(5)/(x6) + 1./(x8)

    cerf  = exp(-x*x)/x * (1./sqrt(pi) +1./x2 * nom/denom)

  END FUNCTION codycerf2approx

  ELEMENTAL FUNCTION inverf(x) RESULT(ierf)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: x
    REAL             :: ierf
    !------------------------------------------------------------------------!
! FIXME : find global error variable, depending on precision(single,double,..)
    REAL, PARAMETER  :: error = 1e-14
    REAL             :: q,y ,sqrtpi,sqrtpi_0_5,fy
    !------------------------------------------------------------------------!

     IF(x.LT.1.0)THEN
       sqrtpi = sqrt(pi)
       sqrtpi_0_5 = sqrtpi*0.5

       ! initial guess:
       y  = sqrtpi*(0.5*x + pi/24.*x**3 + 7.*pi**2/960.*x**5)
       q  = sqrtpi_0_5 * exp(y*y)*(errorfunc(y)-x)   !f(y)/f'(y)
       fy = errorfunc(y) - x

       DO WHILE (abs(fy) .GT. error)
          y = y - q/(1.+q*y)
          fy = errorfunc(y) - x
          q = sqrtpi_0_5 * exp(y*y)*(errorfunc(y)-x)
       END DO
     ierf = y
     ELSE
     ! infinity
        ierf = huge(1.)
     END IF

   END FUNCTION inverf

  ELEMENTAL FUNCTION gamma(x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: x
    REAL             :: res
    !------------------------------------------------------------------------!
    REAL, DIMENSION(26), PARAMETER :: g = &
    (/ 1.0e0,                 0.5772156649015329e0,&
      -0.6558780715202538e0, -0.420026350340952e-1, &
       0.1665386113822915e0, -0.421977345555443e-1, &
      -0.96219715278770e-2,   0.72189432466630e-2, &
      -0.11651675918591e-2,  -0.2152416741149e-3, &
       0.1280502823882e-3,   -0.201348547807e-4, &
      -0.12504934821e-5,      0.11330272320e-5, &
      -0.2056338417e-6,       0.61160950e-8, &
       0.50020075e-8,        -0.11812746e-8, &
       0.1043427e-9,          0.77823e-11, &
      -0.36968e-11,           0.51e-12, &
      -0.206e-13,            -0.54e-14, &
       0.14e-14,              0.1e-15/)
    INTEGER         :: m1,k,m
    REAL            :: z,r,gr
    !------------------------------------------------------------------------!
    r=1 ! this surpresses the 'maybe-uninitialized' compiler warning
    IF (x.EQ.int(x)) THEN
      IF (x.GT.0.) THEN
        res=1.0
        m1=x-1
        DO k=2,m1
          res=res*k
        END DO
      ELSE
        res=1.0e+300
      END IF
    ELSE
      IF (abs(x).GT.1.) THEN
        z=abs(x)
        m=int(z)
        !R=1.
        DO k=1,m
          r=r*(z-k)
        END DO
        z=z-m
      ELSE
        z=x
      ENDIF
      gr=g(26)
      DO k=25,1,-1
        gr=gr*z+g(k)
      END DO
      res=1.0/(gr*z)
      IF (abs(x).GT.1.) THEN
        res=res*r
        IF (x.LT.0.) &
          res=-pi/(x*res*sin(pi*x))
      ENDIF
    ENDIF
  END FUNCTION gamma

 ELEMENTAL  FUNCTION lngamma(x) RESULT(lg)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    DOUBLE PRECISION :: x
    REAL             :: lg
    DOUBLE PRECISION :: alpha,beta,ag,con
    INTEGER          :: i
    !------------------------------------------------------------------------!
    INTENT(IN)       :: x
    !------------------------------------------------------------------------!
    DOUBLE PRECISION, DIMENSION(15), PARAMETER &
                     :: p = &
                 (/ .99999999999999709182d0,      57.156235665862923517d0,&
                   -59.597960355475491248d0,      14.136097974741747174d0,&
                   -.49191381609762019978d0,     .33994649984811888699d-4,&
                   .46523628927048575665d-4,    -.98374475304879564677d-4,&
                   .15808870322491248884d-3,    -.21026444172410488319d-3,&
                   .21743961811521264320d-3,    -.16431810653676389022d-3,&
                   .84418223983852743293d-4,    -.26190838401581408670d-4,&
                   .36899182659531622704d-5                                /)

    ! prefactor
    alpha = 2.506628274631000502415d0
    beta  = x + 5.24218750d0
    con = log(alpha)+log(beta)*(x+0.5d0) - beta
    lg = 0

    ! sum
    IF (x.GE.0.0) THEN
      ag = p(1)
      DO i=2,15
        ag = ag + p(i)/(x+i-1.d0)
      END DO
      lg  = con + log(ag)
    ELSE
!      CALL Error(this,"LnGamma","bad value for for gammafunction")
    END IF
  END FUNCTION lngamma
END MODULE functions