fgauss.f90

Go to the documentation of this file.

function fgauss(x, nv, equidistant, s_acc, weights, sig) result(f_fgauss)

!----------------------------------------------------------------------------!
! computes a sum of weighted Gauss functions                                 !
!  at position x                                                             !
!  s_acc(1..nv)     = centers of Gauss functions                             !
!  weigths(1..nv)   = weighting factors                                      !
!  sig(1..nv)       = widths of Gauss functions                              !
!  equidistant = 0..1: 1 produces equidistant grid, 0 produces pure Gaussian !
!----------------------------------------------------------------------------!

  use itm_types

  implicit none

  integer(itm_i4) :: nv
  real(r8), dimension(nv) :: s_acc, weights, sig
  real(r8) :: x, equidistant
  real(r8) :: f_fgauss

  real(r8) :: gauss_sum, weight
  integer(itm_i4) :: i

  weight = 0._r8
  gauss_sum = 0._r8
  do i = 1, nv
    weight = weight + weights(i)
! to avoid floating underflow: exp(-large_number) = 0
    if (((x - s_acc(i))**2 / sig(i)**2) < 100) then
      gauss_sum = gauss_sum + weights(i) * (0.39894_r8 / sig(i)) &
       * exp(-0.5_r8 * (x - s_acc(i))**2 / sig(i)**2)
    end if
  end do
  f_fgauss = equidistant + (1._r8 - equidistant) * gauss_sum / weight

end function fgauss