[neutron-mc] Maxwell random deviate generator

[George Apostolopoulos]

Dear all,

I found today a much more elegant solution to this problem. Actually the
following three lines of code:

{
  double u1,u2;
  do u1 = rand01(); while (u1==0.0);
  do u2 = rand01(); while (u2==0.0);
  return sqrt(-log(u1*u2));
}

generate random velocities with the required Maxwell flux distribution
2*exp(-v^2)*v^3. (v is scaled to the thermal velocity v_T = 2kT/m)

What is most interesting for McStas is the following:
Generating velocities with the above routine is slower (about 0.6 times)
than the standard 
McStas method of first generating velocities (or equivalently
wavelengths) uniformly and then calculating the particle weight. But:
the same accuracy
of Monte-Carlo results is reached with fewer particle histories when
sampling the true distribution. I get roughly a performance improvement
of 1.5 to 2 on my computer, depending on compiler.

As an example, the small attached program "maxwell2.c" calculates the
average <v^2> with both methods. 
The output of  "maxwell2.c" on my computer with gcc:

************************************************************************
********
Testing the Maxwell random deviate generator.

The value of <v^2> = 2 will be calculated by 2 methods:
  1. Generating random v distributed uniformly in 0.1 < v <5 and
     multiplying by the distribution weight w = 2*exp(-v^2)*v^3
  2. Generating random v distributed according to 2*exp(-v^2)*v^3

1st method

N          <v^2>      delta<v^2>/<v^2> time(ms)
       100   1.888390       1.463e-001      0.000
      1000   2.157340       4.092e-002      0.000
     10000   2.010312       1.352e-002      0.000
    100000   2.004246       4.283e-003     15.000
   1000000   1.998229       1.355e-003    172.000
  10000000   1.999877       4.282e-004   1690.000

2nd method

N          <v^2>      delta<v^2>/<v^2> time(ms)
       100   2.109436       6.675e-002      0.000
      1000   2.007962       2.170e-002      0.000
     10000   2.007131       7.057e-003      0.000
    100000   2.000076       2.232e-003     16.000
   1000000   1.998449       7.057e-004    235.000
  10000000   1.999998       2.234e-004   2582.000
************************************************************************
*****

The pdf file is a short description of the method. 

Best regards,
George

-----Original Message-----
From: neutron-mc-bounces at risoe.dk [mailto:neutron-mc-bounces at risoe.dk]
On Behalf Of George Apostolopoulos
Sent: Thursday, October 04, 2007 4:56 PM
To: neutron-mc at risoe.dk
Subject: [neutron-mc] Maxwell random deviate generator

Sorry, forgot to paste the c-code.

******************************************************
Dr. Georgios Apostolopoulos
Institute of Nuclear Technology & Radiation Protection
National Center for Scientific Research "Demokritos"
15310 Agia Paraskevi Attikis, Greece

tel +30-210-650-3731 
fax +30-210-653-3431
email: gapost at ipta.demokritos.gr
******************************************************
-------------- next part --------------
A non-text attachment was scrubbed...
Name: maxwell2.pdf
Type: application/octet-stream
Size: 58804 bytes
Desc: maxwell2.pdf
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20071005/de44ced7/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: maxwell2.c
Type: application/octet-stream
Size: 2638 bytes
Desc: maxwell2.c
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20071005/de44ced7/attachment-0001.obj>