| 55 | } |
| 56 | |
| 57 | double gaussian_random(double mean, double stddev) { |
| 58 | init_rand(); |
| 59 | // Box-Muller algorithm to generate Gaussian from uniform |
| 60 | // see Knuth vol II algorithm P, sec. 3.4.1 |
| 61 | double v1, v2, s; |
| 62 | do { |
| 63 | v1 = 2 * meep_mt_genrand_res53() - 1; |
| 64 | v2 = 2 * meep_mt_genrand_res53() - 1; |
| 65 | s = v1 * v1 + v2 * v2; |
| 66 | } while (s >= 1.0); |
| 67 | if (s == 0) { return mean; } |
| 68 | else { return mean + v1 * sqrt(-2 * log(s) / s) * stddev; } |
| 69 | } |
| 70 | |
| 71 | } // namespace meep |
no test coverage detected