module elementary_functions_m !! Elementary mathematical functions use precision_m, only : FP, FP_EPS, FP_NAN, FP_SMALLEST_EXP implicit none public :: factorial public :: gaussian public :: step_tanh public :: step_hermite public :: heaviside public :: binomial_coefficient_r contains pure real(FP) function factorial(n) result(res) !! Factorial of integer n; wrapper around intrinsic gamma function. !! Output is real to increase overflow threshold. For overflow case, !! result = Infinity; for n < 0, result = NaN. integer, intent(in) :: n !! Integer to take the factorial of ! For n < 0: forcing n to be an integer ensures that gamma is evaluated ! exactly at its poles, returning the expected NaN ! Special case: gamma(0) = Infinity, whereas factorial(-1) should be ! undefined if (n == -1) then res = FP_NAN return end if res = gamma(n + 1.0_FP) end function factorial pure real(FP) function gaussian(x0, wx, x) !! Gaussian exp[-(x-x0)^2/wx^2] real(FP), intent(in) :: x0 !! Center of Gaussion real(FP), intent(in) :: wx !! Width of Gaussion real(FP), intent(in) :: x !! Position where Gaussian is evaluated real(FP) :: expon if (wx < FP_EPS) then gaussian = 0.0_FP else expon = -1.0_FP * (x - x0) * (x - x0) / (wx * wx) ! avoid exponential underflow if(expon < FP_SMALLEST_EXP) then gaussian = 0.0_FP else gaussian = exp(expon) endif endif end function gaussian pure real(FP) function step_tanh(x0, x, wx) !! Smooth step function based on hyperbolic tangent real(FP), intent(in) :: x0 !! Position of center of step real(FP), intent(in) :: x !! Position where function is evaluated real(FP), intent(in), optional :: wx !! Width of step, if set to zero yields heaviside function step_tanh = heaviside(x0, x) if (present(wx)) then if (wx > FP_EPS) then step_tanh = (tanh((x - x0) / wx) + 1.0_FP) / 2.0_FP endif endif end function step_tanh pure real(FP) function heaviside(x0, x) !! Heaviside function, i.e. discontinuous step function real(FP), intent(in) :: x0 !! Position of step real(FP), intent(in) :: x !! Position where function is evaluated if (x >= x0) then heaviside = 1.0_FP else heaviside = 0.0_FP endif end function heaviside pure real(FP) function step_hermite(x0, x, wx, order) !! Smooth step function based on Hermite interpolation !! exactly zero for x < wx/2 and one for x > wx/2 !! see https://en.wikipedia.org/wiki/Smoothstep real(FP), intent(in) :: x0 !! Position of step real(FP), intent(in) :: x !! Position where function is evaluated real(FP), intent(in), optional :: wx !! Width of step integer, intent(in), optional :: order !! Order of smoothstep function !! Default value 2 integer :: n, k real(FP) :: xn, sh if (present(order)) then n = order else n = 2 endif step_hermite = heaviside(x0, x) if (present(wx)) then if (wx > FP_EPS) then xn = (x - x0) / wx + 0.5_FP; if (xn <= 0.0_FP) then step_hermite = 0.0_FP elseif (xn >= 1) then step_hermite = 1.0_FP else sh = 0.0_FP do k= 0, n sh = sh + binomial_coefficient_r(n+k, k) * binomial_coefficient_r(2*n+1, n-k) * (-xn)**k enddo step_hermite = xn**(n+1) * sh endif endif endif end function pure real(FP) function binomial_coefficient_r(n, k) !! Computes binomial coefficient n over k !! Returns -1 for invalid input, i.e. n or k < 0 or k > n ! Algorithm is safe for large numbers (therefore it returns a real number) integer, intent(in) :: n !! Top number of binomal coefficient integer, intent(in) :: k !! Bottom number of binomal coefficient integer :: i if ((n < 0) .or. (k < 0) .or. (k > n)) then binomial_coefficient_r = -1 return endif binomial_coefficient_r = 1.0_FP do i = 1, k binomial_coefficient_r = binomial_coefficient_r * (n-k+i)/(1.0_FP*i) enddo end function end module