scale_specfunc.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
module scale_specfunc
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: sf_gamma

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  ! Lanczos approximation
  ! g = 7.0, N = 8
  function sf_gamma(x)
    implicit none

    real(RP), intent(in) :: x
    real(RP)             :: sf_gamma

    real(DP), parameter :: lrtp = 0.9189385332046727_dp ! log(sqrt(2*PI))
    real(DP)            :: ag, dx, work
    !---------------------------------------------------------------------------

    dx = real(x-1.0_rp, kind=dp)

    ag = 0.999999999999809932276847004735_dp                     &
       + 676.520368121885098567009190444_dp    / ( dx + 1.0_dp ) &
       - 1259.13921672240287047156078755_dp    / ( dx + 2.0_dp ) &
       + 771.323428777653078848652825889_dp    / ( dx + 3.0_dp ) &
       - 176.615029162140599065845513540_dp    / ( dx + 4.0_dp ) &
       + 12.5073432786869048144589368533_dp    / ( dx + 5.0_dp ) &
       - 0.138571095265720116895547069851_dp   / ( dx + 6.0_dp ) &
       + 9.98436957801957085956266899569e-6_dp / ( dx + 7.0_dp ) &
       + 1.50563273514931155833835577539e-7_dp / ( dx + 8.0_dp )

    work = exp( lrtp + (dx+0.5_dp)*log(dx+7.5_dp) - (dx+7.5_dp) + log(ag) )
    sf_gamma = real(work, kind=rp)

    return
  end function sf_gamma

end module scale_specfunc