! user_gamma 2 REAL REAL .TRUE.
! user_normal 2 REAL REAL .TRUE.

MODULE user_dist
  IMPLICIT NONE

CONTAINS

  SUBROUTINE user_normal_gen(mu, s, buffer)
    ! A wrapper subroutine for the G05FDF NAG normal generation.
    IMPLICIT NONE
    INTEGER :: param
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), DIMENSION(:), INTENT(OUT) :: buffer
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(IN) :: mu, s
    EXTERNAL G05FDF

    CALL G05FDF (mu, s, SIZE(buffer), buffer)

  END SUBROUTINE user_normal_gen

  SUBROUTINE user_gamma_gen(a, b, buffer)
    ! UNDER CONSTRUCTION
    IMPLICIT NONE
    REAL, INTENT(IN) :: a, b
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), DIMENSION(:), INTENT(OUT) :: buffer
    INTEGER :: ifail
    EXTERNAL G05FFF

    ifail = 0
    CALL g05FFF(a, b, SIZE(buffer), buffer, ifail)

    IF (ifail /= 0) THEN
       WRITE (*, *) 'Error when trying to generate proposals of the gamma distribution, parameters: a = ', a, ' b = ', b 
       STOP
    END IF

  END SUBROUTINE user_gamma_gen

  SUBROUTINE user_normal_freq(my, sigma2, point, freq)
    ! palauttaa x:n tiheyden logaritmin Norm(my,sigma2)-jakaumassa
    IMPLICIT NONE
    INTRINSIC SQRT
    INTRINSIC LOG
    INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(8,0)
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN) :: my
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN) :: sigma2
    REAL (KIND=SELECTED_REAL_KIND(8,0)) :: sigma
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN) :: point
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(OUT) :: freq

    sigma = 1.0 / sigma2
    sigma = SQRT(sigma)
!    IF (sigma == 0) THEN
!       sigma = sigma + 1.0D-99
!    END IF
    freq = -LOG(sigma) + user_z_freq((point-my)/sigma)

  END SUBROUTINE user_normal_freq

  FUNCTION user_z_freq(x) RESULT (result)
    ! palauttaa x:n tiheyden logaritmin Norm(0,1)-jakaumassa 
    IMPLICIT NONE
    INTRINSIC SQRT
    INTRINSIC LOG
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN):: x
    REAL (KIND=SELECTED_REAL_KIND(8,0)) :: result


    result =  -0.91893853 - 0.5*x*x

  END FUNCTION user_z_freq

  SUBROUTINE user_gamma_freq(alpha, lambda, point, freq)
    ! palauttaa x:n tiheyden logaritmin Gamma(alpha,lambda)-jakaumassa 
    IMPLICIT NONE
    INTRINSIC LOG
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN):: alpha, lambda, point
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(OUT):: freq

    freq =  (alpha * LOG(lambda) + (alpha-1) * LOG(point) - lambda * point - user_gammaln (alpha))

  END SUBROUTINE user_gamma_freq

  FUNCTION user_gammaln(xx) RESULT (result)
    ! palauttaa xx:n gammafunktion arvon, tarvitaan laskettaessa 
    ! gamma-jakauman tiheyksiä
    IMPLICIT NONE
    INTRINSIC LOG
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN):: xx
    REAL (KIND=SELECTED_REAL_KIND(8,0)) :: x, tmp, ser, result
    REAL (KIND=SELECTED_REAL_KIND(8,0)), DIMENSION(6) :: cof
    INTEGER :: j

    cof(1) =  76.18009173
    cof(2) = -86.50532033
    cof(3) =  24.01409822
    cof(4) = -1.231739516
    cof(5) =  0.120858003e-2
    cof(6) = -0.536382e-5

    x = xx - 1.0
    tmp = x + 5.5
    tmp = tmp - (x + 0.5) * LOG (tmp)
    ser = 1.0

    DO j = 1, 6
       x = x + 1.0;
       ser = ser + cof(j) / x;
    END DO

    result = -tmp + LOG (2.50662827465 * ser)

  END FUNCTION user_gammaln

END MODULE user_dist

