!! -------------------------------------------------------------------------------------------
!Subroutine gammaFit (L, gammaParams)
!  
!
!  return
!End Subroutine gammaFit

!! -------------------------------------------------------------------------------------------
!real Function  gammaStandardize (value, params)
!
!  return
!End Function gammaStandardize

!! -------------------------------------------------------------------------------------------
!Subroutine pearsonIIIFit (L, pearsonIIIParams)
!
!  return
!End Subroutine pearsonIIIFit

!! -------------------------------------------------------------------------------------------
!real Function pearsonIIIStandardize (value, params)
!  return
!end Function pearsonIIIStandardize

!! -------------------------------------------------------------------------------------------
!! Estimates the parameters of a Gamma distribution functions
Subroutine logLogisticFit (beta, logLogisticParams)
  real, dimension (0:2) :: beta, logLogisticParams
  real :: g1, g2
  
  !! estimate gamma parameter
  logLogisticParams (2) = (2*beta (1) - beta (0)) / (6*beta(1)-beta(0)-6*beta(2))
  g1 = exp(gammaLn(1+1/logLogisticParams(2)))
  g2 = exp(gammaLn(1-1/logLogisticParams(2)))
  !! estimate alpha parameter
  logLogisticParams(1) = (beta(0)-2*beta(1))*logLogisticParams(2) / (g1*g2)
  !! estimate beta parameter
  logLogisticParams(0) = beta(0) - logLogisticParams(1)*g1*g2

  !!print*,"beta",beta(0), beta(1), beta(2)
  !!print*,"LOG ",logLogisticParams(0), logLogisticParams(1), logLogisticParams(2)
  return
End Subroutine logLogisticFit

!! Gives the cumulative distribution function of 'value', following a LogLogistic distribution
real Function logLogisticCDF (val, params)
  real :: val, params (0:2)
  logLogisticCDF = 1 / (1 + (params(1)/(val-params(0)))**params(2))
  return
end Function logLogisticCDF

!! Gives the inverse cumulative distribution function following a standard Gaussian
!! distribution. I.e., the function gives the standardized value (Z) corresponding to a
!! given probability (0<prob<1) following a standard Gaussian distribution (mean=0 and sd=1).
real Function standardGaussianInvCDF (prob)
  real :: prob
  real :: C(0:2) = (/ 2.515517, 0.802853, 0.010328 /)
  real :: d(0:3) = (/ 0.000000, 1.432788, 0.189269,0.001308/)
  real :: W, WW, WWW, resul

  if (prob .le. 0.5) then
    W = sqrt(-2*log(prob))
  else
    W = sqrt(-2*log(1-prob))
  endif
  WW = W*W
  WWW = WW*W
  resul = W - (C(0) + C(1)*W + C(2)*WW) / (1 + d(1)*W + d(2)*WW + d(3)*WWW)

  if (prob .gt. 0.5) resul = -resul
  !print*, resul

  standardGaussianInvCDF = resul
  return
end Function standardGaussianInvCDF
!!! EOF ###