NormDist.f90

Go to the documentation of this file.

!===================================================================== 
!===============================================================
 
SUBROUTINE normdist(val, mu, s, x)
 
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
 
    IMPLICIT NONE
 
    ! <Passed variables>--------------------------------------------------------------------------------
    REAL, INTENT(in)  :: mu, &       ! mean                                                         [=mu]
                         s, &        ! standard deviation                                            [=s]
                         x           ! value of rand.var of which distrib. it to be eval.            [=x]
    REAL, INTENT(out) :: val         ! distribution value (output)                              [=height]
 
    ! <Local variables>---------------------------------------------------------------------------------
    REAL             :: xtransf, &   ! local variable                                          [=xtransf]
                        rval, &      ! local variable                                             [=rval]
                        v1, &        ! local variable                                               [=v1]
                        v2, &        ! local variable                                               [=v2]
                        nd(0:201)    ! local variable                           [=NormDistFuncArr(index)]
    INTEGER          :: index        ! local variable                                            [=index]
 
    ! <Here the SUBROUTINE starts>**********************************************************************
    !----- First, initialize the Standard Normal Distribution array nd(..) 202 values, start at -5, end at +5
 
    nd = (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0000106885, 0.0000133457, &
           0.0000166238, 0.0000206575, 0.0000256088, 0.0000316712, 0.0000390756, 0.0000480963, 0.0000590589, &
           0.000072348, 0.0000884173, 0.0001078, 0.00013112, 0.000159109, 0.000192616, 0.000232629, &
           0.000280293, 0.000336929, 0.000404058, 0.000483424, 0.000577025, 0.000687138, 0.000816352, &
           0.000967603, 0.00114421, 0.0013499, 0.00158887, 0.00186581, 0.00218596, 0.00255513, 0.00297976, &
           0.00346697, 0.00402459, 0.00466119, 0.00538615, 0.00620967, 0.00714281, 0.00819754, 0.00938671, &
           0.0107241, 0.0122245, 0.0139034, 0.0157776, 0.0178644, 0.0201822, 0.0227501, 0.0255881, 0.0287166, &
           0.0321568, 0.0359303, 0.0400592, 0.0445655, 0.0494715, 0.0547993, 0.0605708, 0.0668072, 0.0735293, &
           0.0807567, 0.088508, 0.0968005, 0.10565, 0.11507, 0.125072, 0.135666, 0.146859, 0.158655, 0.171056, &
           0.18406, 0.197663, 0.211855, 0.226627, 0.241964, 0.257846, 0.274253, 0.29116, 0.308538, 0.326355, &
           0.344578, 0.363169, 0.382089, 0.401294, 0.42074, 0.440382, 0.460172, 0.480061, 0.5, 0.519939, &
           0.539828, 0.559618, 0.57926, 0.598706, 0.617911, 0.636831, 0.655422, 0.673645, 0.691462, 0.70884, &
           0.725747, 0.742154, 0.758036, 0.773373, 0.788145, 0.802337, 0.81594, 0.828944, 0.841345, 0.853141, &
           0.864334, 0.874928, 0.88493, 0.89435, 0.9032, 0.911492, 0.919243, 0.926471, 0.933193, 0.939429, &
           0.945201, 0.950529, 0.955435, 0.959941, 0.96407, 0.967843, 0.971283, 0.974412, 0.97725, 0.979818, &
           0.982136, 0.984222, 0.986097, 0.987776, 0.989276, 0.990613, 0.991802, 0.992857, 0.99379, 0.994614, &
           0.995339, 0.995975, 0.996533, 0.99702, 0.997445, 0.997814, 0.998134, 0.998411, 0.99865, 0.998856, &
           0.999032, 0.999184, 0.999313, 0.999423, 0.999517, 0.999596, 0.999663, 0.99972, 0.999767, 0.999807, &
           0.999841, 0.999869, 0.999892, 0.999912, 0.999928, 0.999941, 0.999952, 0.999961, 0.999968, 0.999974, &
           0.999979, 0.999983, 0.999987, 0.999989, 0.999991, 0.999993, 0.999995, 0.999996, 0.999997, 0.999997, &
           0.999998, 0.999998, 0.999999, 0.999999, 0.999999, 0.999999, 1.0, 1.0, 1.0, 1.0/)
 
    !----- Transform the value where the Normal Distr. distribution function is to be evaluated
    xtransf = (x - mu)/s
    !----- If transformed evaluation point is outside critical values, set the result...
    if (xtransf <= -5.0) then
        val = 0.0
    else if (xtransf >= 5.0) then
        val = 1.0
        !----- ... otherwise use the values from the array nd by transformation and interpolation
    else
        rval = 20.0*(xtransf + 5.0)
        index = max(0, min(201, int(rval)))
        !----- Interpolate
        v1 = nd(index)
        v2 = nd(index + 1)
        val = v1 + (rval - float(index))*(v2 - v1)
    end if
 
end SUBROUTINE normdist
!done