Attribute VB_Name = "modStatistics"
Option Compare Database ' declare the default string comparison method based on the sort order of the database
Option Explicit         ' Makes the programmer declare all variables
Option Base 0           ' Explicitly sets the lower bound of arrays to zero

Public Function dGEGMED(indata() As Variant, lngCount As Long) As Double
'
' Purpose:
' To take sorted array and array size find median
' History:
' Version 1.0 January 2004 by G.E. Granato
' indata is a single precision array of input values
' lngcount is a long value indicating the size of the array
'
' NOTE -1 accounts for zero based arrays
' Arguments
Dim lngfi As Long
Dim lngfj As Long
Dim lngfk As Long
Dim lngfin As Long
Dim llbound As Long

On Error GoTo dGEGMED_ERR:

llbound = LBound(indata)

If lngCount = 1 Then
 dGEGMED = indata(0)
 Exit Function
End If

lngfj = lngCount Mod 2 'modulus provides remainder
lngfk = lngCount \ 2 'integer division truncates the -1 is used because VBA uses 0 to N-1 Array
' get median value
If (lngfj = 0) Then  'this is an even number
 dGEGMED = (indata((lngfk - 1)) + indata((lngfk + 1) - 1)) / 2#
Else 'this is an odd number
 dGEGMED = indata((lngfk + 1) - 1)
End If

Exit Function

dGEGMED_ERR:
dGEGMED = -9999#
End Function
Public Function fndCensorUniform01Rescale(dU01 As Double, dFractionCensored As Double, Optional intHiLow As Integer) As Double
' Purpose:
' To take a single uniform random number in the interval between 0 and 1 and
' a proportion censored value and rescale the uncensored values to a U01
' return this U01 value to model the distribution of noncensored values
'
' Note: The user must test if the original U01 is >= the fraction of censoring and
' use this function to rescale values above censoring to generate numbers from a
' different probability.  The original U01 values represent the probability of the whole
' sample.
'
' Algorithm developed by by G.E. Granato
'
' Arguments:
' dU01 As Double input uniform random number in the interval between 0 and 1
' dFractionCensored the dFractionCensored As Double the fraction that are censored
' intHiLow as Integer indicates 0 lower censored 1 upper censored
'
' History:
' Version 1.0 September 19 2008 by G.E. Granato

On Error GoTo Calculation_Err:

  intHiLow = Int(Nz(intHiLow, 0)) ' if not specified default is lower censoring
  
  If dU01 <= 0# Or dU01 >= 1# Then fndCensorUniform01Rescale = -999999#
  If dFractionCensored <= 0# Or dFractionCensored >= 1# Then fndCensorUniform01Rescale = -999999#
  If intHiLow = 0 And dFractionCensored <= dU01 Then fndCensorUniform01Rescale = -999999#
  If intHiLow = 1 And dFractionCensored >= dU01 Then fndCensorUniform01Rescale = -999999#
 ' If Uniform01ToExponential = -999999# Then Exit Function
  
  If intHiLow = 0 Then
   fndCensorUniform01Rescale = (dU01 - dFractionCensored) / (1# - dFractionCensored)
  Else
   fndCensorUniform01Rescale = (dU01) / (dFractionCensored)
  End If
  
  Exit Function
  
Calculation_Err:
   'Uniform01ToExponential = -999999#

End Function

Public Sub PrimeAS241(dA0 As Double, dA1 As Double, dA2 As Double, dA3 As Double, dA4 As Double, dA5 As Double, dA6 As Double, dA7 As Double, _
dB1 As Double, dB2 As Double, dB3 As Double, dB4 As Double, dB5 As Double, dB6 As Double, dB7 As Double, _
dC0 As Double, dC1 As Double, dC2 As Double, dC3 As Double, dC4 As Double, dC5 As Double, dC6 As Double, dC7 As Double, _
dD1 As Double, dD2 As Double, dD3 As Double, dD4 As Double, dD5 As Double, dD6 As Double, dD7 As Double, _
dE0 As Double, dE1 As Double, dE2 As Double, dE3 As Double, dE4 As Double, dE5 As Double, dE6 As Double, dE7 As Double, _
dF1 As Double, dF2 As Double, dF3 As Double, dF4 As Double, dF5 As Double, dF6 As Double, dF7 As Double)
'
' Purpose:
' To initialize values for the AS241 Polynomial Values for a uniform to normal distribution function
' Initialize separately to reduce reasignment with each value in an array
'
' History:
' Version 1.0 June 2006 by G.E. Granato ggranato@usgs.gov
' Method developed from algorithm developed and described as Accurate to more than 10^-7 and less than 6x10-16
' Wichura, M.J., 1988, Algorithm AS241--The percentage points of the normal distribution:
' Applied Statistics, Journal of the Royal Statistical Society (Series C)
' v. 37, no. 3, p. 477-484.
'
' dA0-dA7, dC0-dC7, dC0-dC7 polynomial coefficients for numerator
' dB1-dB7, dD1-dD7, dF1-dF7 polynomial coefficients for denominator
'
' Assign values

dA0 = 3.38713287279637
dA1 = 133.141667891784
dA2 = 1971.59095030655
dA3 = 13731.6937655095
dA4 = 45921.9539315499
dA5 = 67265.7709270087
dA6 = 33430.575583588
dA7 = 2509.08092873012
dB1 = 42.3133307016009
dB2 = 687.187007492058
dB3 = 5394.19602142475
dB4 = 21213.7943015866
dB5 = 39307.8958000927
dB6 = 28729.0857357219
dB7 = 5226.49527885285
dC0 = 1.42343711074968
dC1 = 4.63033784615654
dC2 = 5.76949722146069
dC3 = 3.6478483247632
dC4 = 1.27045825245237
dC5 = 0.241780725177451
dC6 = 2.27238449892692E-02
dC7 = 7.74545014278341E-04
dD1 = 2.05319162663776
dD2 = 1.6763848301838
dD3 = 0.6897673349851
dD4 = 0.14810397642748
dD5 = 1.51986665636165E-02
dD6 = 5.47593808499534E-04
dD7 = 1.05075007164442E-09
dE0 = 6.6579046435011
dE1 = 5.46378491116411
dE2 = 1.78482653991729
dE3 = 0.296560571828504
dE4 = 2.65321895265761E-02
dE5 = 1.24266094738807E-03
dE6 = 2.71155556874348E-05
dE7 = 2.01033439929228E-07
dF1 = 0.599832206555887
dF2 = 0.136929880922735
dF3 = 1.48753612908506E-02
dF4 = 7.86869131145613E-04
dF5 = 1.84631831751005E-05
dF6 = 1.42151175831644E-07
dF7 = 2.04426310338993E-15
Exit Sub

End Sub
Public Sub MRG32k3a(dOutUniform01 As Double, dInSeed10 As Double, dInSeed11 As Double, dInSeed12 As Double, _
dInSeed20 As Double, dInSeed21 As Double, dInSeed22 As Double)
'
' Purpose:
' Combined Multiple Recursive Random Number Generator
' Takes input seeds and returns uniform random number and output seed values
' History:
' Version 1.0 February 2004
' Programmed in Visual Basic 6.0 by Gregory E. Granato
' Adapted from data and Design (in C) published in:
' L'Ecuyer, Pierre, 1999, Good parameters and implementations
' for combined multiple recursive random number generators
' Operations Research, v. 47, no. 1, pp. 159-164.
'
' Subroutine arguments
' dInSeed10 -  dInSeed12 input MRRNG seeds for component 1
' dInSeed20 -  dInSeed22 input MRRNG seeds for component 2
' dOutUniform01  Uniform Random Number between 0 and 1
'
' Internal variables
Dim dP1 As Double       ' recursive modulus 1
Dim dP2 As Double       ' recursive modulus 2
Dim dK As Double        ' quotient from integer division
' Multiple Recursive Random Number Generator (MRRNG) constants
Const cdblnorm As Double = 2.32830654929573E-10  'MRNG modulus constant norm =1/(m1+1)
Const cdblm1 As Double = 4294967087#            'MRNG modulus constant 1
Const cdblm2 As Double = 4294944443#            'MRNG modulus constant 2
Const cdbla12 As Double = 1403580#              'MRNG part 1 constant 1
Const cdbla13n As Double = 810728#              'MRNG part 1 constant 2
Const cdbla22 As Double = 527612#               'MRNG part 2 constant 1
Const cdbla23n As Double = 1370589#             'MRNG part 2 constant 2
' Assign values
' set seed
If dInSeed11 = 0 Then dInSeed11 = dInSeed10
If dInSeed12 = 0 Then dInSeed12 = dInSeed11
If dInSeed21 = 0 Then dInSeed21 = dInSeed20
If dInSeed22 = 0 Then dInSeed22 = dInSeed21
' ***********************
' Component 1
' ***********************
 dP1 = (cdbla12 * dInSeed11) - (cdbla13n * dInSeed10)
 ' integer division
 dK = Int(dP1 / cdblm1)
 ' get modulus
 dP1 = dP1 - (dK * cdblm1)
 ' set to ensure nonnegative nonzero value
 If dP1 < 0# Then dP1 = dP1 + cdblm1
 ' increment seeds
 dInSeed10 = dInSeed11
 dInSeed11 = dInSeed12
 dInSeed12 = dP1
' ***********************
' Component 2
' ***********************
 dP2 = (cdbla22 * dInSeed21) - (cdbla23n * dInSeed20)
 ' integer division
 dK = Int(dP2 / cdblm2)
 ' get modulus
 dP2 = dP2 - (dK * cdblm2)
 ' set to ensure nonnegative nonzero value
 If dP2 < 0# Then dP2 = dP2 + cdblm2
 ' increment seeds
  dInSeed20 = dInSeed21
  dInSeed21 = dInSeed22
  dInSeed22 = dP2
' ***********************
' Combime result
' ***********************
 If (dP1 <= dP2) Then
  dOutUniform01 = (((dP1 - dP2) + cdblm1) * cdblnorm)
 Else
  dOutUniform01 = ((dP1 - dP2) * cdblnorm)
 End If
End Sub


Public Sub PrimeWilsonHilfertyKirby(dInskew As Double, dA As Double, dB As Double, dG As Double, dH As Double)
'
' Purpose:
' Computes parameters used in the Wilson-Hilferty transformation for skewed distributions
' Equivalent to WILFR$ by William Kirby, USGS Office of Surface Water 1972
'
' Kirby, W., 1972, Computer-oriented Wilson-Hilferty transformation that preserves the first
' three moments and the lower bound of the Pearson Type 3 Distribution:
' Water Resources Research, v. 8, no.5, p. 1251-1254.
'
' History:
' William Kirby Wilson Hilferty
'      WKIRBY  72-02-25  ORIGINAL WILFRT
'      REVISED 73-02-09   TO ACCEPT ZERO SKEW.
' The method was converted to VBA by Gregory E. Granato 2006 ggranato@usgs.gov
' Updated to consolidate values by Gregory E. Granato September 2008
'
' Arguments
Dim dTSkew(41) As Double ' Skew coefficient from the table, which varies from 0.0 to 9.75
' The array variables are the difference between Kirby Table 1 and:
Dim dDG(41) As Double ' the approximation G= skew-0.063 * Max (0, skew-1)^1.85
Dim dDA(41) As Double ' the approximation A= max (2/skew, 0.4)
Dim dDB(41) As Double ' the approximation B= 1.0 + 0.0144 * Max(0, skew-2.25)^2
'
Dim dSkew As Double    ' internal skew
Dim dISkew As Double   ' Temporary skew-table generator
Dim dP As Double       ' proportional distance between input skew and table skew values
Dim dQ As Double       ' inverse of proportional distance between input skew and table skew values
Dim dAproxA As Double  ' Approximate A
Dim dAproxB As Double  ' Approximate B
Dim dIMVar As Double   ' Intermediate variable for calculations
Dim lngI As Long       ' index variable

' **********************************************
' Populate interpolation-data table with 40 values each
' **********************************************
' Assign 40 skew values from 0 to 9.75
 dISkew = 0#
 dTSkew(1) = dISkew
 For lngI = 2 To 40
  dTSkew(lngI) = dTSkew(lngI - 1) + 1# / 4#
 Next lngI

' Assign 40 dDG() values for skews 0-9.75 with a 0.25 increment
dDG(1) = 0#
dDG(2) = -0.000144
dDG(3) = -0.001137
dDG(4) = -0.003762
dDG(5) = -0.008674
dDG(6) = -0.011555
dDG(7) = -0.010076
dDG(8) = -0.006049
dDG(9) = -0.000921
dDG(10) = 0.004189
dDG(11) = 0.008515
dDG(12) = 0.011584
dDG(13) = 0.013139
dDG(14) = 0.013122
dDG(15) = 0.010945
dDG(16) = 0.007546
dDG(17) = 0.002767
dDG(18) = -0.003181
dDG(19) = -0.010089
dDG(20) = -0.017528
dDG(21) = -0.025476
dDG(22) = -0.033609
dDG(23) = -0.042434
dDG(24) = -0.050525
dDG(25) = -0.058192
dDG(26) = -0.065221
dDG(27) = -0.07141
dDG(28) = -0.076638
dDG(29) = -0.080655
dDG(30) = -0.083349
dDG(31) = -0.084584
dDG(32) = -0.084203
dDG(33) = -0.082089
dDG(34) = -0.078126
dDG(35) = -0.072165
dDG(36) = -0.064188
dDG(37) = -0.054059
dDG(38) = -0.041633
dDG(39) = -0.027005
dDG(40) = -0.010188
'
' Assign 40 dDA() values for skews 0-9.75 with a 0.25 increment
dDA(1) = 0#
dDA(2) = 0.004614
dDA(3) = 0.009159
dDA(4) = 0.013553
dDA(5) = 0.017753
dDA(6) = 0.021764
dDA(7) = 0.025834
dDA(8) = 0.030406
dDA(9) = 0.03571
dDA(10) = 0.04173
dDA(11) = 0.048321
dDA(12) = 0.055309
dDA(13) = 0.062538
dDA(14) = 0.069873
dDA(15) = 0.077334
dDA(16) = 0.084682
dDA(17) = 0.091926
dDA(18) = 0.099028
dDA(19) = 0.105967
dDA(20) = 0.112695
dDA(21) = 0.119245
dDA(22) = 0.106551
dDA(23) = 0.095488
dDA(24) = 0.085671
dDA(25) = 0.07699
dDA(26) = 0.06929
dDA(27) = 0.062443
dDA(28) = 0.056349
dDA(29) = 0.050908
dDA(30) = 0.046047
dDA(31) = 0.041702
dDA(32) = 0.037815
dDA(33) = 0.034339
dDA(34) = 0.031229
dDA(35) = 0.028445
dDA(36) = 0.025964
dDA(37) = 0.023753
dDA(38) = 0.021782
dDA(39) = 0.020043
dDA(40) = 0.018528
'
' Assign 40 dDB() values for skews 0-9.75 with a 0.25 increment
dDB(1) = 0#
dDB(2) = 0#
dDB(3) = -0.000001
dDB(4) = -0.000004
dDB(5) = -0.000021
dDB(6) = -0.000075
dDB(7) = -0.00019
dDB(8) = -0.000326
dDB(9) = -0.000317
dDB(10) = 0.000116
dDB(11) = 0.000434
dDB(12) = 0.000116
dDB(13) = -0.000464
dDB(14) = -0.000981
dDB(15) = -0.001165
dDB(16) = -0.000743
dDB(17) = 0.000435
dDB(18) = 0.002479
dDB(19) = 0.005462
dDB(20) = 0.009353
dDB(21) = 0.014206
dDB(22) = 0.019964
dDB(23) = 0.026829
dDB(24) = 0.034307
dDB(25) = 0.042495
dDB(26) = 0.051293
dDB(27) = 0.060593
dDB(28) = 0.070324
dDB(29) = 0.080332
dDB(30) = 0.090532
dDB(31) = 0.100831
dDB(32) = 0.111114
dDB(33) = 0.121283
dDB(34) = 0.131245
dDB(35) = 0.140853
dDB(36) = 0.15012
dDB(37) = 0.158901
dDB(38) = 0.167085
dDB(39) = 0.174721
dDB(40) = 0.181994
'
' **********************************************
' Compute Approximate values and adjust with the interpolation table
' **********************************************
' Get absolute Value of Skew for table interpolation
  dSkew = Abs(dInskew)

' If the skew is very small then truncate the lower bound
  If dSkew < 0.0005 Then dSkew = 0.0005

' Limit skew coefficient to 9.750 excessive truncate the upper bound
  If dSkew > 9.75 Then dSkew = 9.75

' Find index number for the input skew in the interpolation table
   lngI = 2
  Do While dTSkew(lngI) < dSkew
   lngI = lngI + 1
  Loop
' Calculate interpolation factors
  dP = (dTSkew(lngI) - dSkew) / (1# / 4#)
  dQ = 1# - dP
' Compute sG
  dG = dSkew + (dQ * dDG(lngI) + dP * dDG(lngI - 1))
' Adjust sG
  If (dSkew > 1#) Then dG = dG - 0.063 * (dSkew - 1#) ^ 1.85
' Compute Approximate A where A=Max (2.0/skew, 0.400)
  If ((2# / dSkew) > (4# / 10#)) Then
   dAproxA = (2# / dSkew)
  Else
   dAproxA = (4# / 10#)
  End If
' Compute A
  dA = dAproxA + (dQ * dDA(lngI) + dP * dDA(lngI - 1))
' Compute Approximate B  = 1. + 0.0144 max (0, skew-2.25)
  If (dSkew <= 2.25) Then
   dAproxB = 1#
  Else
   dAproxB = 1# + 0.0144 * (dSkew - 2.25) ^ 2#
  End If
' Compute B
   dB = dAproxB + (dQ * dDB(lngI) + dP * dDB(lngI - 1))
' Compute H = B -((2/skew)/A)^1/3
   dH = ((dB - ((2# / dSkew) / dA)) ^ (1# / 3#))
   
Exit Sub
   
End Sub

Public Function fndAdjustedWilsonHilfertyK(dInputSkew As Double, dInputNormalK As Double, _
Optional dMyA As Double, Optional dMyB As Double, Optional dMyG As Double, Optional dMyH As Double) As Double
'
' Purpose:
' To take the skew and the standard normal variate and return a Modified Wilson Hilferty variate
' Returns normal variate if Skew<0.005 unajusted Wilson Hilferty if Skew < 0.5
' and Kirby's adjusted wilson Hilferty if skew >= 0.5
' If the optional input parameters are input from calling subroutine (for use in large loops) will use those
' Else will get them from interpolation table (PrimeWilsonHilfertyKirby)
'
' Modified from WILFRT -- WILSON-HILFERTY REVISED TRANSFORM QUANTILES
' by William Kirby, USGS Office of Surface Water 1972
' WILFRT = A * MAX(H, (1-(G/6)**2+(G/6)*ZETA)**3 - B)
' ZETA(P) IS THE P-TH STANDARDIZED NORMAL QUANTILE AND PARAMETERS
' A,B,H ARE CHOSEN SUCH THAT WILFRT HAS MEAN 0, STDDEV 1, AND SKEW G.
'
' Differences between fndAdjustedWilsonHilfertyK percentage points and harters tables are of the order of
' a few hundredths of a std. Deviation, except in extreme positive tail (95% or so)
' where error is of order of tenths in magnitude but about 3% in relative mag.
'
' Kirby, W., 1972, Computer-oriented Wilson-Hilferty transformation that preserves the first
' three moments and the lower bound of the Pearson Type 3 Distribution:
' Water Resources Research, v. 8, no.5, p. 1251-1254.
'
' History:
' William Kirby Wilson Hilferty
'      WKIRBY  72-02-25  ORIGINAL WILFRT
'      REVISED 73-02-09   TO ACCEPT ZERO SKEW.
' The method was converted to VBA by Gregory E. Granato 2006 ggranato@usgs.gov
' Updated to consolidate values and handle small skews by Gregory E. Granato September 2008
'
' Uses Subroutine:
' PrimeWilsonHilfertyKirby(dInskew As Double, dA As Double, dB As Double, dG As Double, dH As Double)
' To get interpolation table values
'

Dim dABSkew As Double ' Absolute value of inskew
Dim dSby6 As Double    ' Skew/6 for Unajusted Wilson Hilferty
Dim dMyVal As Double
Dim bCallPrimeSub As Boolean

Dim strMsgStr As String

On Error GoTo fndAdjustedWilsonHilfertyK_Err:

dABSkew = Abs(dInputSkew) ' take absolute value of the skew for processing

' *******************************************************************************
' Ignore very small skew values
' The standard error of the skew coefficient of a normal distribution is sqrt(6/N)
' a skew of 0.005 is within the standard error for sample sizes (N) less than 240,000
If dABSkew < 0.005 Then ' set small skew values to zero set dAdjustedWilsonHilferty K = NormalK
 fndAdjustedWilsonHilfertyK = dInputNormalK
 Exit Function
End If

' *******************************************************************************
' Use unajusted Wilson Hilferty approximation for skews below the second entry in
' Kirby's adjustment table see PrimeWilsonHilfertyKirby(dInputSkew, dMyA, dMyB, dMyG, dMyH)
' Unajusted Wilson Hilferty approximation slightly better for for skews below 0.5
If dABSkew < 0.5 Then
 dSby6 = dInputSkew / 6#
 fndAdjustedWilsonHilfertyK = (2# / dInputSkew) * (((dSby6 * (dInputNormalK - dSby6) + 1#) ^ 3#) - 1)
 Exit Function
End If

' Limit skew coefficient to 9.750 excessive skew truncated
  If dInputSkew > 9.75 Then dInputSkew = 9.75
  If dInputSkew < -9.75 Then dInputSkew = -9.75 ' oops added 11/14/2011

' *******************************************************************************
' Test to see that optional inputs are properly specified
' if not get values from PrimeWilsonHilfertyKirby(dInputSkew, dMyA, dMyB, dMyG, dMyH)
  bCallPrimeSub = False
  If dMyA < 0.4 Then
   bCallPrimeSub = True
  ElseIf dMyB < 0.98 Or dMyB > 2# Then
   bCallPrimeSub = True
  ElseIf dMyG < 0.00045 Or dMyB > 6.5 Then
   bCallPrimeSub = True
  ElseIf dMyH > 1.16 Then
   bCallPrimeSub = True
  End If
  
  If bCallPrimeSub = True Then Call PrimeWilsonHilfertyKirby(dInputSkew, dMyA, dMyB, dMyG, dMyH)
  
' *******************************************************************************
' Calculate Variate

If dInputSkew < 0 Then ' neg skew
 
 dMyVal = (1# - ((dMyG / 6#) * (dMyG / 6#)) - (dMyG / 6#) * dInputNormalK)
 
 If dMyH > dMyVal Then dMyVal = dMyH
 
 fndAdjustedWilsonHilfertyK = -1# * dMyA * ((dMyVal * dMyVal * dMyVal) - dMyB)

ElseIf dInputSkew > 0 Then 'pos skew
 dMyVal = (1# - ((dMyG / 6#) * (dMyG / 6#)) + (dMyG / 6#) * dInputNormalK)
 
 If dMyH > dMyVal Then dMyVal = dMyH
 
 fndAdjustedWilsonHilfertyK = 1# * dMyA * ((dMyVal * dMyVal * dMyVal) - dMyB)

End If
 
Exit Function
' Error
fndAdjustedWilsonHilfertyK_Err:
 fndAdjustedWilsonHilfertyK = dInputNormalK
 Exit Function
End Function

Public Function fndPlottingPosition(lngRank As Long, lngCount As Long, intType As Integer, _
 Optional bSortAscending As Boolean = True, Optional bPercentage As Boolean = False) As Double
'
' Purpose: to generate plotting positions from a rank
' History: Version 1.0 January 25 2009 by Gregory E. Granato
'
' Function arguments
' lngRank As Double the rank of the value in the sample
' lngCount As Long the total number of values in the sample
' intType As Integer the type of plotting poisition formula
' bSortAscending As Boolean True (default) = Sort ascending False = Sort descending
' bPercentage As Boolean True (default) = percentage (0<PP<100); False = P (0<PP<1)
'
' Arguments
Dim dA As Double  ' Plotting position formula a value

On Error GoTo fndPlottingPosition_Err:

 bSortAscending = CBool(Nz(bSortAscending, "True"))
 bPercentage = CBool(Nz(bPercentage, "True"))

Select Case intType
 Case 0 'Blom
  dA = 0.375
 Case 1 'Cunnane
  dA = 0.4
 Case 2 'Gringorten
  dA = 0.44
 Case 3 'Hazen
  dA = 0.5
 Case 4 'Median
  dA = 0.3175
 Case 5 'Weibull
  dA = 0#
End Select

 fndPlottingPosition = (CDbl(lngRank) - dA) / (CDbl(lngCount) + 1# - 2# * dA)
 ' False = Sort descending
 If bSortAscending = False Then fndPlottingPosition = 1# - fndPlottingPosition
 '
 If bPercentage = True Then fndPlottingPosition = 100# * fndPlottingPosition
 
 Exit Function

fndPlottingPosition_Err:
 fndPlottingPosition = -9999#
  Exit Function
End Function
Public Function fnlngRankFromPP(dInputPP As Double, lngCount As Long, intType As Integer, _
 Optional bSortAscending As Boolean, Optional bPercentage As Boolean) As Long
'
' Purpose: to generate a rank from a plotting position
' History: Version 1.0 April 09 2009 by Gregory E. Granato
'
' Function arguments
' dInputPP As Double the input plotting position (or percentile) value
' lngCount As Long the total number of values in the sample
' intType As Integer the type of plotting poisition formula
' bSortAscending As Boolean True (default) = Sort ascending False = Sort descending
' bPercentage As Boolean False (default) False (0<dInputPP<1); True percentage (0< dInputPP <100);
'
' Arguments
Dim dA As Double  ' Plotting position formula a value
Dim lngTempRank As Long

On Error GoTo fnlngRankFromPP_Err:

 bSortAscending = CBool(Nz(bSortAscending, "True"))
 bPercentage = CBool(Nz(bPercentage, "False"))

Select Case intType
 Case 0 'Blom
  dA = 0.375
 Case 1 'Cunnane
  dA = 0.4
 Case 2 'Gringorten
  dA = 0.44
 Case 3 'Hazen
  dA = 0.5
 Case 4 'Median
  dA = 0.3175
 Case 5 'Weibull
  dA = 0#
End Select

 '
 If bPercentage = True Then dInputPP = dInputPP / 100#

 ' False = Sort descending
 If bSortAscending = False Then dInputPP = 1# - dInputPP
 
 lngTempRank = CLng(Round((dInputPP * (lngCount + 1# - (2# * dA)) + dA), 0))
 
 If lngTempRank < 1 Then lngTempRank = 1
 If lngTempRank > lngCount Then lngTempRank = lngCount
 
 fnlngRankFromPP = lngTempRank
 
 Exit Function

fnlngRankFromPP_Err:
 fnlngRankFromPP = -9999#
  Exit Function
End Function

Public Function fndUniform01ToExponential(dU01 As Double, dMean As Double, dMin As Double) As Double
'
' Purpose:
' To take a single uniform random number in the interval between 0 and 1 and
' return the equivalent exponential number using the distribution function (CDF inversion)
'
' Algorithm from:
' Saucier, Richard, 2000, Computer generation of statistical distributions:
' U.S. Army Research Laboratory Report ARL-TR-2168, 105 p.
'
' Arguments:
' dU01 As Double input uniform random number in the interval between 0 and 1
' dMean As Double the mean of the exponential distribution
' dMin As Double the minimum of the exponential distribution
'
' History:
' Version 1.0 September 18 2008 by G.E. Granato

On Error GoTo Calculation_Err:
' check input values
  fndUniform01ToExponential = dU01
  If dU01 <= 0# Then dU01 = 0.0000000002328
  If dU01 >= 1# Then dU01 = 1# - 0.0000000002328
  If dMin < 0# Then fndUniform01ToExponential = -999999#
  If dMean < dMin Then fndUniform01ToExponential = -999999#
  If fndUniform01ToExponential = -999999# Then Exit Function
  ' Note use (1# - dU01) so exponential values go up with increasing U
  fndUniform01ToExponential = dMin - ((dMean - dMin) * Log((1# - dU01)))
  
  Exit Function
  
Calculation_Err:
   fndUniform01ToExponential = -999999#
   Exit Function
End Function

Public Function fndEncapsulatedUniform01ToNormalAS241(dMyInputUniform As Double) As Double
'
' Purpose:
' To convert a uniform random number (0<U<1) to the equivalent standard normal score value with Algorithm AS241
' Encapsulates the function and the primeing function to be used when only one or a few normals are needed
'
' History:
' Version 1.0 Feb 07 2009 by G.E. Granato ggranato@usgs.gov
'
' Get the equation coefficients
Dim dA0 As Double, dA1 As Double, dA2 As Double, dA3 As Double, dA4 As Double, dA5 As Double, dA6 As Double, dA7 As Double
Dim dB1 As Double, dB2 As Double, dB3 As Double, dB4 As Double, dB5 As Double, dB6 As Double, dB7 As Double
Dim dC0 As Double, dC1 As Double, dC2 As Double, dC3 As Double, dC4 As Double, dC5 As Double, dC6 As Double, dC7 As Double
Dim dD1 As Double, dD2 As Double, dD3 As Double, dD4 As Double, dD5 As Double, dD6 As Double, dD7 As Double
Dim dE0 As Double, dE1 As Double, dE2 As Double, dE3 As Double, dE4 As Double, dE5 As Double, dE6 As Double, dE7 As Double
Dim dF1 As Double, dF2 As Double, dF3 As Double, dF4 As Double, dF5 As Double, dF6 As Double, dF7 As Double

Call PrimeAS241(dA0, dA1, dA2, dA3, dA4, dA5, dA6, dA7, _
dB1, dB2, dB3, dB4, dB5, dB6, dB7, _
dC0, dC1, dC2, dC3, dC4, dC5, dC6, dC7, _
dD1, dD2, dD3, dD4, dD5, dD6, dD7, _
dE0, dE1, dE2, dE3, dE4, dE5, dE6, dE7, _
dF1, dF2, dF3, dF4, dF5, dF6, dF7)

' Get the Normal Random Variate value
fndEncapsulatedUniform01ToNormalAS241 = fndUniform01ToNormalAS241(dMyInputUniform, _
dA0, dA1, dA2, dA3, dA4, dA5, dA6, dA7, _
dB1, dB2, dB3, dB4, dB5, dB6, dB7, _
dC0, dC1, dC2, dC3, dC4, dC5, dC6, dC7, _
dD1, dD2, dD3, dD4, dD5, dD6, dD7, _
dE0, dE1, dE2, dE3, dE4, dE5, dE6, dE7, _
dF1, dF2, dF3, dF4, dF5, dF6, dF7)

End Function


Public Function fndUniform01ToNormalAS241(dMyUniform As Double, _
dA0 As Double, dA1 As Double, dA2 As Double, dA3 As Double, dA4 As Double, dA5 As Double, dA6 As Double, dA7 As Double, _
dB1 As Double, dB2 As Double, dB3 As Double, dB4 As Double, dB5 As Double, dB6 As Double, dB7 As Double, _
dC0 As Double, dC1 As Double, dC2 As Double, dC3 As Double, dC4 As Double, dC5 As Double, dC6 As Double, dC7 As Double, _
dD1 As Double, dD2 As Double, dD3 As Double, dD4 As Double, dD5 As Double, dD6 As Double, dD7 As Double, _
dE0 As Double, dE1 As Double, dE2 As Double, dE3 As Double, dE4 As Double, dE5 As Double, dE6 As Double, dE7 As Double, _
dF1 As Double, dF2 As Double, dF3 As Double, dF4 As Double, dF5 As Double, dF6 As Double, dF7 As Double) As Double
'
' Purpose:
' To convert a uniform random number (0<U<1) to the equivalent standard normal score value with Algorithm AS241
'
' History:
' Version 1.0 June 2006 by G.E. Granato ggranato@usgs.gov
' Method developed from algorithm developed and described as Accurate to more than 10^-7 and less than 6x10-16
' Wichura, M.J., 1988, Algorithm AS241--The percentage points of the normal distribution:
' Applied Statistics, Journal of the Royal Statistical Society (Series C)
' v. 37, no. 3, p. 477-484.
'
' dMyUniform input uniform random numbers in interval 1E-19 < dMyUniform <1
' dMyNormal output standard normal score values
' dA0-dA7, dC0-dC7, dC0-dC7 polynomial coefficients for numerator
' dB1-dB7, dD1-dD7, dF1-dF7 polynomial coefficients for denominator
'
' Assign P values
' Arguments
Dim dMyP As Double
Dim dY As Double
Dim dR As Double
Dim dQ As Double
Dim dMyNormal As Double

On Error GoTo fndUniform01ToNormalAS241_Err:

' Short cut if near zero
If (dMyUniform > 0.499999999999999) And (dMyUniform < 0.500000000000001) Then  ' if U = 0.5 then Z=0
 fndUniform01ToNormalAS241 = 0#
 Exit Function
End If

' Condition Input
If dMyUniform <= 1E-21 Then dMyUniform = 1E-21           ' Keep within bounds
If dMyUniform >= 1# Then dMyUniform = 0.999999999999999  ' Keep within bounds

dMyP = dMyUniform

dQ = dMyP - 0.5

If Abs(dQ) <= 0.425 Then ' Is the value close to one-half use first polygon

 dR = 0.180625 - (dQ * dQ)
 dMyNormal = dQ * (((((((dA7 * dR + dA6) * dR + dA5) * dR + dA4) * dR + dA3) * dR + dA2) * dR + dA1) * dR + dA0) / _
                (((((((dB7 * dR + dB6) * dR + dB5) * dR + dB4) * dR + dB3) * dR + dB2) * dR + dB1) * dR + 1#)


Else  ' use other polygons
 If dQ < 0# Then
  dR = dMyP
 Else
  dR = 1# - dMyP
 End If
 
 dR = Sqr(-1# * Log(dR))
 
 If dR <= 5# Then
  dR = dR - 1.6
  dMyNormal = (((((((dC7 * dR + dC6) * dR + dC5) * dR + dC4) * dR + dC3) * dR + dC2) * dR + dC1) * dR + dC0) / _
                (((((((dD7 * dR + dD6) * dR + dD5) * dR + dD4) * dR + dD3) * dR + dD2) * dR + dD1) * dR + 1#)
 Else
  dR = dR - 5#
  dMyNormal = (((((((dE7 * dR + dE6) * dR + dE5) * dR + dE4) * dR + dE3) * dR + dE2) * dR + dE1) * dR + dE0) / _
                (((((((dF7 * dR + dF6) * dR + dF5) * dR + dF4) * dR + dF3) * dR + dF2) * dR + dF1) * dR + 1#)
 End If
 
 If (dQ < 0) Then dMyNormal = dMyNormal * -1#
 
End If

fndUniform01ToNormalAS241 = dMyNormal  ' Assign final value

Exit Function
' Error
fndUniform01ToNormalAS241_Err:
 fndUniform01ToNormalAS241 = 0#
 Exit Function
End Function

Public Function fndUniform01ToTrapezoid(dU01 As Double, dMin As Double, dLower As Double, _
dUpper As Double, dMax As Double) As Double
'
' Purpose:
' Trapezoidal Random Number Generator
' Takes input uniform random number and returns values in a trapezoidal distribution
' that has a lower triangle (dMin to dLower) , a middle plateau (dLower to dUpper), and
' an upper triangle (dUpper to dMax)
' Special cases:
' if (dMin = dLower) and (dUpper = dMax) it is rectangular
' if (dLower = dUpper) it is a triangle
'
' History:
' Version 1.0 October 18 2008
' Programmed in Visual Basic 6.0 by Gregory E. Granato
' Adapted from data inverse cumulative distribution function published in:
' Kacker, R.N., and Lawrence, J.F., 2007, Trapezoid and triangular distributions for Type B evaluation of
' standard uncertainty: Metrologia, v. 44, p. 117-127.
'
' Function arguments
' dU01 As Double Uniform Random Number between 0 and 1
' dMin As Double Lower limit of the trapezoid
' dLower As Double Lower limit of the plateau
' dUpper As Double Upper limit of the plateau
' dMax As Double Upper limit of the trapezoid
' Internal variables
Dim dH As Double ' Height of the trapezoid

On Error GoTo fndUniform01ToTrapezoid_Err:

' check for errors
If (dMin > dLower) Or (dMin > dUpper) Or (dMin >= dMax) Then ' error return U01
 fndUniform01ToTrapezoid = dU01
 Exit Function
ElseIf (dLower > dUpper) Or (dLower > dMax) Then ' error return U01
 fndUniform01ToTrapezoid = dU01
 Exit Function
ElseIf (dUpper > dMax) Then ' error return U01
 fndUniform01ToTrapezoid = dU01
 Exit Function
End If
 
If (dMin = dLower) And (dUpper = dMax) Then ' Uniform/rectangle
 fndUniform01ToTrapezoid = dMin * (dMax - dMin) * dU01
 Exit Function
End If

dH = 2# / ((dMax - dMin) + (dUpper - dLower))

If dU01 >= 0# And dU01 <= ((dH / 2#) * (dLower - dMin)) Then
 fndUniform01ToTrapezoid = dMin + (Sqr(2# * ((dLower - dMin) / dH))) * Sqr(dU01)
 
ElseIf dU01 > ((dH / 2#) * (dLower - dMin)) And dU01 <= (1 - ((dH / 2#) * (dMax - dUpper))) Then
 fndUniform01ToTrapezoid = ((dMin + dLower) / 2#) + (dU01 / dH)
 
ElseIf dU01 > (1 - ((dH / 2#) * (dMax - dUpper))) And dU01 <= 1 Then
 fndUniform01ToTrapezoid = dMax - (Sqr((2# * (dMax - dUpper)) / dH) * Sqr(1# - dU01))
End If

Exit Function
' Error
fndUniform01ToTrapezoid_Err:
 fndUniform01ToTrapezoid = dU01
 Exit Function
End Function

Public Function fndUniform01ToTriangular(dU01 As Double, dMode As Double, dMin As Double, dMax As Double) As Double
'
' Purpose:
' To take a single uniform random number in the interval between 0 and 1 and
' return the equivalent triangular number using the distribution function (CDF inversion)
'
' Algorithm from:
' Saucier, Richard, 2000, Computer generation of statistical distributions:
' U.S. Army Research Laboratory Report ARL-TR-2168, 105 p.
'
' Arguments:
' dU01 As Double input uniform random number in the interval between 0 and 1
' dMode As Double the mode of the triangular distribution
' dMin As Double the minimum of the triangular distribution
' dMax As Double the maximum of the triangular distribution
'
' History:
' Version 1.0 September 18 2008 by G.E. Granato

On Error GoTo Calculation_Err:

' check input values
  fndUniform01ToTriangular = dU01
  If dU01 <= 0# Or dU01 >= 1# Then fndUniform01ToTriangular = -999999#
  If dMin >= dMax Then fndUniform01ToTriangular = -999999#
  If dMode < dMin Or dMode > dMax Then fndUniform01ToTriangular = -999999#
  If fndUniform01ToTriangular = -999999# Then Exit Function
  ' transform
  If dU01 <= ((dMode - dMin) / (dMax - dMin)) Then
   fndUniform01ToTriangular = dMin + Sqr(((dMax - dMin) * (dMode - dMin)) * dU01)
  Else
   fndUniform01ToTriangular = dMax - Sqr(((dMax - dMin) * (dMax - dMode)) * (1# - dU01))
  End If
  
  Exit Function
  
Calculation_Err:
   fndUniform01ToTriangular = -999999#
End Function
Public Sub SeedTableStatistics(lngMaxSeed As Long, lngMinSeed As Long, lngCountSeed As Long, _
Optional intIOErr As Integer)
'
' Purpose: To get the minimum, maximum, and seed count from tblURNSeeds to scale seed selections
'
' History: Version 1.0 Feburary 7 2009 by G.E. Granato
'
Dim bHasdata As Boolean                 ' test data
Dim rstURNSeedRecordset As ADODB.Recordset   ' recordset
Dim strSQL As String                    ' Query text
Dim intOutErr As Integer

On Error GoTo SeedTableStatistics_Err:
intOutErr = 0
intIOErr = 0

 strSQL = "SELECT Count(tblURNSeeds.URNSeed_ID) AS CountOfURNSeed_ID, " & _
 "Max(tblURNSeeds.URNSeed_ID) AS MaxOfURNSeed_ID, " & _
 "Min(tblURNSeeds.URNSeed_ID) AS MinOfURNSeed_ID " & _
 "FROM tblURNSeeds;"

'Reference an ADO Recordset
 Set rstURNSeedRecordset = New ADODB.Recordset
  
'Populate Recordset
 rstURNSeedRecordset.Open strSQL, CurrentProject.Connection, adOpenStatic, adLockReadOnly

' Check values
 bHasdata = False
 If rstURNSeedRecordset.BOF And rstURNSeedRecordset.EOF Then  ' There is no record
  bHasdata = False
 ElseIf IsNull(rstURNSeedRecordset.Fields(0).Value) = False Then
  If IsEmpty(rstURNSeedRecordset.Fields(0).Value) = False Then bHasdata = True ' There are data
 End If
 
 If rstURNSeedRecordset.Fields(0).Value = 0 Then bHasdata = False
 
 If bHasdata = False Then
  lngCountSeed = 0
  lngMaxSeed = 0
  lngMinSeed = 0
  intIOErr = -1
 Else
  lngCountSeed = rstURNSeedRecordset.Fields("CountOfURNSeed_ID").Value
  lngMaxSeed = rstURNSeedRecordset.Fields("MaxOfURNSeed_ID").Value
  lngMinSeed = rstURNSeedRecordset.Fields("MinOfURNSeed_ID").Value
  intIOErr = 0
 End If
 
CleanUp:
 If rstURNSeedRecordset.State = adStateOpen Then rstURNSeedRecordset.Close  'IF Recordset is open, close it
 Set rstURNSeedRecordset = Nothing   'Free memory
 
Exit Sub
SeedTableStatistics_Err:
  intIOErr = -1
  intOutErr = intOutErr + 1
  If intOutErr < 2 Then Resume CleanUp:
  Exit Sub
End Sub

Public Sub GetRankCorrelation(dInRho As Double, dInU01 As Double, dOutU01 As Double, _
dU2Seed10 As Double, dU2Seed20 As Double, _
Optional dU2Seed11 As Double, Optional dU2Seed12 As Double, _
Optional dU2Seed21 As Double, Optional dU2Seed22 As Double)
' Purpose:
' To take a correlation value, the primary uniform random number (0-1),
' seed values for the correlated random number and return a uniform random number (0-1)
' with the correct correlation.  This is done by starting with a 1:1 relation between ranks
' (direction of relation set by the sign of the correlation), using the correlation coefficient
' to estimate the noise around the 1:1 line, and rejecting values outside the 0-1 range
'
' Algorithm from:
' Mykytka E.F., and Cheng, C.Y., 1994,
' Generating correlated random variates based on an analogy between correlation and force
' in Tew, J.D., Manivannan, S., Sadowski, D.A., and Seila, A.F. eds.,
' Proceedings of the 1994 Winter Simulation Conference of the Association for Computing Machinery,
' December 11-14, 1994, Lake Buena Vista, FL, p. 1413-1416.
'
'
' History:
' Version 1.0 October 22, 2012 by G.E. Granato
'
Dim dA As Double ' The first factor rho' (the absolute value of the adjusted rho)
Dim dB As Double ' The second factor (1-(rho')^2)^0.5
Dim dC As Double ' The third factor 0.5 * (1 - A - B)
Dim dMyRho As Double
Dim dRho As Double

On Error GoTo GetRankCorrelation_Err:

dMyRho = Abs(dInRho)

' Adjust the absolute value of Rho to acount for the Bias in average Rho in Mykytka's algorithm
If dMyRho <= 0.2 Then
 dRho = dMyRho + (0.0578 * dMyRho) - 0.0012
ElseIf dMyRho <= 0.7 Then
 dRho = dMyRho - (0.3245 * (dMyRho ^ 2)) + (0.3155 * dMyRho) - 0.0527
ElseIf dMyRho <= 0.77 Then
 dRho = dMyRho - (0.126 * dMyRho) + 0.0974
ElseIf dMyRho <= 0.97 Then
 dRho = dMyRho - (0.6814 * (dMyRho ^ 3)) + (2.2569 * (dMyRho ^ 2)) - (2.3823 * dMyRho) + 0.8078
Else
 dRho = dMyRho
End If

' Calculate the factors with the adjusted absolute value of Rho
dA = dRho
dB = (1# - (dRho ^ 2#)) ^ 0.5
dC = (1# - dA - dB) / 2#

' Get the uniform random number
Call MykytkaRho(dRho, dInU01, dOutU01, dA, dB, dC, dU2Seed10, dU2Seed20, _
                dU2Seed11, dU2Seed12, dU2Seed21, dU2Seed22)

' Adjust results for negative correlation
If dInRho < 0# Then dOutU01 = 1# - dOutU01

Exit Sub

GetRankCorrelation_Err:

 dOutU01 = dInU01 ' out = in user can tell if paying attention

Exit Sub
End Sub

Public Sub MykytkaRho(dInR As Double, dInU01 As Double, dOutU01 As Double, _
dA As Double, dB As Double, dC As Double, _
dU2Seed10 As Double, dU2Seed20 As Double, _
Optional dU2Seed11 As Double, Optional dU2Seed12 As Double, _
Optional dU2Seed21 As Double, Optional dU2Seed22 As Double)
' Purpose:
' To take a correlation value, the primary uniform random number (0-1),
' seed values for the correlated random number and return a uniform random number (0-1)
' with the correct correlation.  This is done by starting with a 1:1 relation between ranks
' (direction of relation set by the sign of the correlation), using the correlation coefficient
' to estimate the noise around the 1:1 line, and rejecting values outside the 0-1 range
' ****Note this sub does not handle negative correlation it is handled by the calling function
' Algorithm from:
' Mykytka E.F., and Cheng, C.Y., 1994,
' Generating correlated random variates based on an analogy between correlation and force
' in Tew, J.D., Manivannan, S., Sadowski, D.A., and Seila, A.F. eds.,
' Proceedings of the 1994 Winter Simulation Conference of the Association for Computing Machinery,
' December 11-14, 1994, Lake Buena Vista, FL, p. 1413-1416.
'
'
' History:
' Version 1.0 October 22, 2012 by G.E. Granato
'
' dInR As Double     The input absolute value of the adjusted correlation coefficient (rho')
' dInU01 As Double   The input uniform random number
' dOutU01 As Double  The output uniform random number
' dA As Double       The first factor rho' (the absolute value of the adjusted rho)
' dB As Double       The second factor (1-(rho')^2)^0.5
' dC As Double       The third factor 0.5 * (1 - A - B)
' dU2Seed10 As Double - dU2Seed22 As Double the random number seeds for the second value
'
' Arguments
 Dim dY2 As Double
 Dim dU3 As Double
 Dim dTempR As Double       ' Temporary value of R (absolute value)

On Error GoTo Mykytka_Err:

dTempR = Abs(dInR)
' First test for Very High Rho
' Assign values and Exit sub if essentially the same or different
 If dTempR > 0.9975 Then  ' for all intents and purposes they are on the 1:1 line
  dOutU01 = dInU01
  Exit Sub
 End If
 
' If Input seeds not set return input value
  If Nz(dU2Seed10, 0#) = 0# And Nz(dU2Seed20, 0#) = 0# Then
   dU2Seed10 = CDbl(Round(999999999# * Rnd(), 0))
   dU2Seed20 = CDbl(Round(7777777# * Rnd(), 0))
  End If
 ' Prime seed if necessary
 If Nz(dU2Seed11, 0#) = 0# Or Nz(dU2Seed12, 0#) = 0# Or Nz(dU2Seed21, 0#) = 0# Or Nz(dU2Seed22, 0#) = 0# Then
  Call MRG32k3a(dU3, dU2Seed10, dU2Seed11, dU2Seed12, dU2Seed20, dU2Seed21, dU2Seed22)
  Call MRG32k3a(dU3, dU2Seed10, dU2Seed11, dU2Seed12, dU2Seed20, dU2Seed21, dU2Seed22)
  Call MRG32k3a(dU3, dU2Seed10, dU2Seed11, dU2Seed12, dU2Seed20, dU2Seed21, dU2Seed22)
 End If

 ' Assign values and if essentially random
 If dTempR < 0.08 Then  ' for all intents and purposes they are independent
  Call MRG32k3a(dOutU01, dU2Seed10, dU2Seed11, dU2Seed12, dU2Seed20, dU2Seed21, dU2Seed22)
  Exit Sub
 End If
  
 Call MRG32k3a(dU3, dU2Seed10, dU2Seed11, dU2Seed12, dU2Seed20, dU2Seed21, dU2Seed22)
 
 ' Initial Calculation:
 dY2 = dA * dInU01 + dB * dU3 + dC

 If dTempR <= (1# / (2# ^ 0.5)) Then
  If dY2 < dC Then
   'This should not happen but...
   dOutU01 = dInU01
  ElseIf dY2 <= (dA + dC) Then
   dOutU01 = ((dY2 - dC) ^ 2#) / (2# * dA * dB)
  ElseIf dY2 <= (dB + dC) Then
   dOutU01 = (dY2 - (dA / 2#) - dC) / dB
  ElseIf dY2 <= (dA + dB + dC) Then
   dOutU01 = 1# - (((dA + dB + dC - dY2) ^ 2#) / (2# * dA * dB))
  Else
   'this should not happen but...
   dOutU01 = dInU01
  End If
 Else
  If dY2 < dC Then
   'this should not happen but...
    dOutU01 = dInU01
  ElseIf dY2 <= (dB + dC) Then
   dOutU01 = ((dY2 - dC) ^ 2#) / (2# * dA * dB)
  ElseIf dY2 <= (dA + dC) Then
   dOutU01 = (dY2 - (dB / 2#) - dC) / dA
  ElseIf dY2 <= (dA + dB + dC) Then
   dOutU01 = 1# - (((dA + dB + dC - dY2) ^ 2#) / (2# * dA * dB))
  Else
   'this should not happen but...
   dOutU01 = dInU01
  End If
 End If
 
Exit Sub

Mykytka_Err:
 dOutU01 = dInU01 ' out = in user can tell if paying attention
Exit Sub
End Sub