[continued from previous message]   
      
   : set-atomic-mass ( F: m_amu -- )  kg/amu f* m_atom_kg f! ;   
      
   m_He set-atomic-mass  \ default value of mass of atom   
      
   fvariable a   
      
   : set-T ( F: Tk -- )   
        kB f* 2.0e0 f* m_atom_kg f@ fswap f/ a f! ;   
      
   \ Return probability density, p(v), for the M-B distribution   
   \ with parameter "a"   
   : pdf ( F: v -- p )   
        fsquare fdup a f@ f* fnegate fexp f*   
        4pi f* a f@ pi f/ 1.5e0 f** f* ;   
      
   \ Return cumulative probability, C(v1) i.e. P(0 <= v < v1),   
   \ the integral of p(v) from 0 to v1 for the M-B distribution.   
   : cdf ( F: v1 -- P )   
        fdup fsquare a f@ f* fdup fsqrt erf   
        fswap fnegate fexp frot f* 2e f* a f@ pi f/ fsqrt f* f- ;   
      
   \ Invert cdf to obtain velocity from value of cdf   
   fvariable pv1   
   : cdf^-1 ( F: P -- v )   
       pv1 f!   
       [: cdf pv1 f@ f- ;] 0.0e0  DEF_VMAX  DEF_TOL )root ;   
      
   \ N random trials of velocities following the M-B distribution   
   \ are sampled and the moments v, v^2, and v^3 are computed   
   \ from the samples and compared to theory.   
   fvariable vsum   
   fvariable v2sum   
   fvariable v3sum   
   : mb-prng ( Ntrials -- ) ( F: --    )   
        DEF_TEMPERATURE set-T   
        0.0e0 fdup fdup vsum f! v2sum f! v3sum f!   
        dup   
        0 ?DO   
          randf   \ random draw of cdf value   
          cdf^-1  \ invert cdf to get random v in M-B distribution   
          fdup         vsum f+!   
          fdup fsquare v2sum f+!   
               fcube   v3sum f+!   
        LOOP   
        s>f   
        vsum  f@ fover f/ fswap   
        v2sum f@ fover f/ fswap   
        v3sum f@ fswap f/ ;   
      
   \ Theoretical moments of v for the M-B distribution at given   
   \ temperature and mass, specified by parameter "a".   
   : _th   ( F: --  )   4.0e0 pi    a f@ f* f/ fsqrt ;   
   : _th ( F: --  ) 3.0e0 2.0e0 a f@ f* f/ ;   
   : _th ( F: --  ) 16.0e0 pi   a f@ fcube f* f/ fsqrt ;   
      
   \ Compute relative errors of randomly sampled averages with   
   \ respect to the theoretical moments.   
   : moment-errors ( F:    -- e1 e2 e3 )   
        frot _th   rel-error   
        frot _th rel-error   
        frot _th rel-error ;   
      
      
   9 constant MAX-POW   
      
   MAX-POW 1+  3 float matrix v_mom{{   
   MAX-POW 1+  3 float matrix m_rel{{   
      
   : print-tables ( -- )   
        cr ." Moments of speed"   
        cr ."  N        (m/s)     (m/s)^2     (m/s)^3" cr   
        MAX-POW 1- 2 DO   
          v_mom{{ I 2 }} f@ v_mom{{ I 1 }} f@ v_mom{{ I 0 }} f@   
          ." 10^" I 1 .r 2 spaces   
          12 4 f.rd 2 spaces   
          12 1 f.rd 6 spaces   
          12 0 f.rd cr   
        LOOP   
        cr ." Relative Errors"   
        cr ."  N       |e1|       |e2|       |e3|"  cr   
        precision  \ save precision   
        3 set-precision   
        MAX-POW 1- 2 DO   
          m_rel{{ I 2 }} f@ m_rel{{ I 1 }} f@ m_rel{{ I 0 }} f@   
          ." 10^" I 1 .r 2 spaces   
          fabs fs. 2 spaces   
          fabs fs. 2 spaces   
          fabs fs. cr   
        LOOP   
        set-precision \ restore precision   
        cr   
      
      
   : test-prng ( xt -- )   
        is prngU   
        prngU-init   
        MAX-POW 1- 2 DO   
          I s>f falog fround>s   
          mb-prng   
          f3dup   
          v_mom{{ I 2 }} f!  v_mom{{ I 1 }} f! v_mom{{ I 0 }} f!   
          moment-errors   
          m_rel{{ I 2 }} f!  m_rel{{ I 1 }} f!  m_rel{{ I 0 }} f!   
        LOOP   
        print-tables   
      
      
   cr cr .( Usage: '  test-prng )   
   cr cr .(   e.g. ' random test-prng ) cr cr   
      
   === end test-prng-mb.4th ===   
      
   --- SoupGate-Win32 v1.05   
    * Origin: you cannot sedate... all the things you hate (1:229/2)