function gaussian_antideriv ( mu, sigma, x )

!*****************************************************************************80
!
!! gaussian_antideriv() evaluates the antiderivative of a Gaussian function.
!
!  Discussion:
!
!    g_anti(mu,sigmax) = 1/sigma/sqrt(2*pi) * 
!      integral ( 0 <= z <= x ) exp ( - 0.5 * (( z - mu )/sigma)^2 ) dz
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 July 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real MU, the mean.
!
!    real SIGMA: the standard deviation.
!
!    real X(M): evaluation points.
!
!  Output:
!
!    real gaussian_antideriv: values of the antiderivative function.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) gaussian_antideriv
  real ( kind = rk8 ) mu
  real ( kind = rk8 ) r8_erf
  real ( kind = rk8 ) sigma
  real ( kind = rk8 ) x
  real ( kind = rk8 ) z1
  real ( kind = rk8 ) z2

  z1 =   mu       / sqrt ( 2.0D+00 ) / sigma
  z2 = ( mu - x ) / sqrt ( 2.0D+00 ) / sigma

  gaussian_antideriv = r8_erf ( z1 ) - r8_erf ( z2 )
  
  return
end
function gaussian_value ( n, mu, sigma, x )

!*****************************************************************************80
!
!! gaussian_value() evaluates a Gaussian function or a derivative.
!
!  Discussion:
!
!    g(mu,sigmax) = 1/sigma/sqrt(2*pi) * exp ( - 0.5 * (( x - mu )/sigma).^2 )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 July 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the derivative order.
!    0 <= N.
!
!    real MU, the mean.
!
!    real SIGMA: the standard deviation.
!
!    real X(M): evaluation points.
!
!  Output:
!
!    real gaussian_value: the value of the N-th derivative of the 
!    Gaussian function.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) g
  real ( kind = rk8 ) gaussian_value
  real ( kind = rk8 ) hermite_polynomial_value
  real ( kind = rk8 ) hn
  real ( kind = rk8 ) mu
  integer n
  real ( kind = rk8 ) root2
  real ( kind = rk8 ) sigma
  real ( kind = rk8 ) x

  root2 = sqrt ( 2.0D+00 )

  g = exp ( - 0.5D+00 * ( ( x - mu ) / sigma )**2 )

  hn = hermite_polynomial_value ( n, ( x - mu ) / sqrt ( sigma ) / root2 )

  gaussian_value = ( -1.0D+00 )**n / ( sigma * root2 )**n * hn * g
  
  return
end
function hermite_polynomial_value ( n, x )

!*****************************************************************************80
!
!! hermite_polynomial_value() evaluates He(i,x).
!
!  Discussion:
!
!    He(i,x) represents the probabilist's Hermite polynomial.
!
!  Differential equation:
!
!    ( exp ( - 0.5 * x^2 ) * He(n,x)' )' + n * exp ( - 0.5 * x^2 ) * He(n,x) = 0
!
!  First terms:
!
!   1
!   X
!   X^2  -  1
!   X^3  -  3 X
!   X^4  -  6 X^2 +   3
!   X^5  - 10 X^3 +  15 X
!   X^6  - 15 X^4 +  45 X^2 -   15
!   X^7  - 21 X^5 + 105 X^3 -  105 X
!   X^8  - 28 X^6 + 210 X^4 -  420 X^2 +  105
!   X^9  - 36 X^7 + 378 X^5 - 1260 X^3 +  945 X
!   X^10 - 45 X^8 + 630 X^6 - 3150 X^4 + 4725 X^2 - 945
!
!  Recursion:
!
!    He(0,X) = 1,
!    He(1,X) = X,
!    He(N,X) = X * He(N-1,X) - (N-1) * He(N-2,X)
!
!  Orthogonality:
!
!    Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * He(M,X) He(N,X) dX 
!    = sqrt ( 2 * pi ) * N! * delta ( N, M )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    03 July 2025
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    National Bureau of Standards, 1964,
!    ISBN: 0-486-61272-4,
!    LC: QA47.A34.
!
!    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
!    NIST Handbook of Mathematical Functions,
!    Cambridge University Press, 2010,
!    ISBN: 978-0521192255,
!    LC: QA331.N57.
!
!  Input:
!
!    integer N, the order of the polynomial.
!
!    real X, the evaluation point.
!
!  Output:
!
!    real hermite_polynomial_value: the value of the probabilist's 
!    Hermite polynomial of index N at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) hermite_polynomial_value
  integer j
  real ( kind = rk8 ) p(0:n)
  real ( kind = rk8 ) x

  p(0) = 1.0D+00

  if ( 0 < n ) then
    p(1) = x
  end if

  do j = 2, n
    p(j) = x * p(j-1) - real ( j - 1, kind = rk8 ) * p(j-2)
  end do
 
  hermite_polynomial_value = p(n)

  return
end
function r8_erf ( x )

!*****************************************************************************80
!
!! r8_erf() evaluates the error function.
!
!  Discussion:
!
!    Since some compilers already supply a routine named ERF which evaluates
!    the error function, this routine has been given a distinct, if
!    somewhat unnatural, name.
!    
!    The function is defined by:
!
!      ERF(X) = ( 2 / sqrt ( PI ) ) * Integral ( 0 <= t <= X ) exp ( - t^2 ) dt
!
!    Properties of the function include:
!
!      Limit ( X -> -oo ) ERF(X) =          -1.0;
!                         ERF(0) =           0.0;
!                         ERF(0.476936...) = 0.5;
!      Limit ( X -> +oo ) ERF(X) =          +1.0.
!
!      0.5 * ( ERF(X/sqrt(2)) + 1 ) = Normal_01_CDF(X)
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 November 2006
!
!  Author:
!
!    Original Fortran77 version by William Cody.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Cody,
!    Rational Chebyshev Approximations for the Error Function,
!    Mathematics of Computation, 
!    1969, pages 631-638.
!
!  Input:
!
!    real ( kind = rk8 ) X, the argument of the error function.
!
!  Output:
!
!    real ( kind = rk8 ) R8_ERF, the value of the error function.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), parameter, dimension ( 5 ) :: a = (/ &
    3.16112374387056560D+00, &
    1.13864154151050156D+02, &
    3.77485237685302021D+02, &
    3.20937758913846947D+03, &
    1.85777706184603153D-01 /)
  real ( kind = rk8 ), parameter, dimension ( 4 ) :: b = (/ &
    2.36012909523441209D+01, &
    2.44024637934444173D+02, &
    1.28261652607737228D+03, &
    2.84423683343917062D+03 /)
  real ( kind = rk8 ), parameter, dimension ( 9 ) :: c = (/ &
    5.64188496988670089D-01, &
    8.88314979438837594D+00, &
    6.61191906371416295D+01, &
    2.98635138197400131D+02, &
    8.81952221241769090D+02, &
    1.71204761263407058D+03, &
    2.05107837782607147D+03, &
    1.23033935479799725D+03, &
    2.15311535474403846D-08 /)
  real ( kind = rk8 ), parameter, dimension ( 8 ) :: d = (/ &
    1.57449261107098347D+01, &
    1.17693950891312499D+02, &
    5.37181101862009858D+02, &
    1.62138957456669019D+03, &
    3.29079923573345963D+03, &
    4.36261909014324716D+03, &
    3.43936767414372164D+03, &
    1.23033935480374942D+03 /)
  real ( kind = rk8 ) del
  integer i
  real ( kind = rk8 ), parameter, dimension ( 6 ) :: p = (/ &
    3.05326634961232344D-01, &
    3.60344899949804439D-01, &
    1.25781726111229246D-01, &
    1.60837851487422766D-02, &
    6.58749161529837803D-04, &
    1.63153871373020978D-02 /)
  real ( kind = rk8 ), parameter, dimension ( 5 ) :: q = (/ &
    2.56852019228982242D+00, &
    1.87295284992346047D+00, &
    5.27905102951428412D-01, &
    6.05183413124413191D-02, &
    2.33520497626869185D-03 /)
  real ( kind = rk8 ) r8_erf
  real ( kind = rk8 ), parameter :: sqrpi = 0.56418958354775628695D+00
  real ( kind = rk8 ), parameter :: thresh = 0.46875D+00
  real ( kind = rk8 ) x
  real ( kind = rk8 ) xabs
  real ( kind = rk8 ), parameter :: xbig = 26.543D+00
  real ( kind = rk8 ) xden
  real ( kind = rk8 ) xnum
  real ( kind = rk8 ) xsq

  xabs = abs ( ( x ) )
!
!  Evaluate ERF(X) for |X| <= 0.46875.
!
  if ( xabs <= thresh ) then

    if ( epsilon ( xabs ) < xabs ) then
      xsq = xabs * xabs
    else
      xsq = 0.0D+00
    end if

    xnum = a(5) * xsq
    xden = xsq
    do i = 1, 3
      xnum = ( xnum + a(i) ) * xsq
      xden = ( xden + b(i) ) * xsq
    end do

    r8_erf = x * ( xnum + a(4) ) / ( xden + b(4) )
!
!  Evaluate ERFC(X) for 0.46875 <= |X| <= 4.0.
!
  else if ( xabs <= 4.0D+00 ) then

    xnum = c(9) * xabs
    xden = xabs
    do i = 1, 7
      xnum = ( xnum + c(i) ) * xabs
      xden = ( xden + d(i) ) * xabs
    end do

    r8_erf = ( xnum + c(8) ) / ( xden + d(8) )
    xsq = real ( int ( xabs * 16.0D+00 ), kind = rk8 ) / 16.0D+00
    del = ( xabs - xsq ) * ( xabs + xsq )
    r8_erf = exp ( - xsq * xsq ) * exp ( - del ) * r8_erf

    r8_erf = ( 0.5D+00 - r8_erf ) + 0.5D+00

    if ( x < 0.0D+00 ) then
      r8_erf = - r8_erf
    end if
!
!  Evaluate ERFC(X) for 4.0D+00 < |X|.
!
  else

    if ( xbig <= xabs ) then

      if ( 0.0D+00 < x ) then
        r8_erf = 1.0D+00
      else
        r8_erf = - 1.0D+00
      end if

    else

      xsq = 1.0D+00 / ( xabs * xabs )

      xnum = p(6) * xsq
      xden = xsq
      do i = 1, 4
        xnum = ( xnum + p(i) ) * xsq
        xden = ( xden + q(i) ) * xsq
      end do

      r8_erf = xsq * ( xnum + p(5) ) / ( xden + q(5) )
      r8_erf = ( sqrpi - r8_erf ) / xabs
      xsq = real ( int ( xabs * 16.0D+00 ), kind = rk8 ) / 16.0D+00
      del = ( xabs - xsq ) * ( xabs + xsq )
      r8_erf = exp ( - xsq * xsq ) * exp ( - del ) * r8_erf

      r8_erf = ( 0.5D+00 - r8_erf ) + 0.5D+00

      if ( x < 0.0D+00 ) then
        r8_erf = - r8_erf
      end if

    end if

  end if

  return
end