; program that demonstrates Monte Carlo error propagation
; written by C. Federrath, 2017
pro mc

    ; define two arrays with random numbers (Gaussian distributions)
    nrand = long(1d6)
    x1 = randomn(33333, nrand, /double)
    x2 = randomn(55555, nrand, /double)

    ; define variables a and b with uncertainties
    a_val = 10d0
    a_err = 1d0
    b_val = 1d0
    b_err = 0.1d0
    a = a_val + a_err * x1
    b = b_val + b_err * x2

    ; get the distribution of a
    get_pdf_astr4004, a, $
      PDF=pdf_a, LOC=loc_a, BINSIZE=binsize

    get_pdf_astr4004, b, $
      PDF=pdf_b, LOC=loc_b, BINSIZE=binsize

    setcolors
    thick = 2
    xrange = [min([loc_a, loc_b]), max([loc_a, loc_b])]
    yrange = [min([pdf_a, pdf_b]), max([pdf_a, pdf_b])]
    plot, loc_a, pdf_a, xr=xrange, yr=yrange, charsize=2.0, $
        xtitle='a or b', ytitle='PDF(a) or PDF(b)', thick=thick
    oplot, loc_b, pdf_b, color=3, linestyle=2, thick=thick

    ; define derived variable
    c = a^2/b^2

    ; Gaussian error propagation
    c_val = a_val^2/b_val^2
    c_err = c_val * sqrt( (2d0*a_err/a_val)^2 + (2d0*b_err/b_val)^2 )
    print, 'Gaussian error propagation: c = a^2/b^2 = '
    print, c_val, ' +/- ', c_err

    ; get PDF of c
    get_pdf_astr4004, c, $
      PDF=pdf_c, LOC=loc_c, BINSIZE=binsize, NBINS=1000, CDF=cdf_c
    plot, loc_c, pdf_c, xrange=[0,400], yrange=[1d-5, 1d-1], charsize=2.0, $
        xtitle='c', ytitle='PDF(c)', thick=thick, psym=10, /ylog

    ; get statistical moments of c
    c_mean    = total(pdf_c * (loc_c+binsize/2d0)   * binsize)
    c_mean_sq = total(pdf_c * (loc_c+binsize/2d0)^2 * binsize)
    c_stddev  = sqrt(c_mean_sq - c_mean^2)
    print, 'Monte Carlo error propagation: c = a^2/b^2 = '
    print, c_mean, ' +/- ', c_stddev
    
    ; compute the mode of c
    ind = where(pdf_c eq max(pdf_c))
    c_mode = loc_c[ind]
    print, 'mode = ', c_mode
    
    ; compute 16th and 84th percentile values of c
    ind = where(cdf_c ge 0.16)
    sixteenth = loc_c[ind[0]]
    ind = where(cdf_c ge 0.84)
    eightyfourth = loc_c[ind[0]]
    print, '16th percentile value = ', sixteenth
    print, '84th percentile value = ', eightyfourth

    ; compute min, max asymmetric error bars
    c_min_diff = c_mean - sixteenth
    c_max_diff = eightyfourth - c_mean
    print, 'c = ' + string(c_mean,     format='(F5.1)') +' -' $
                  + string(c_min_diff, format='(F4.1)') +' +' $
                  + string(c_max_diff, format='(F4.1)') +' (' $
                  + string(c_stddev,   format='(F4.1)') +')'

    stop
  
end