diff options
| -rw-r--r-- | src/partial_sim.c | 29 |
1 files changed, 28 insertions, 1 deletions
diff --git a/src/partial_sim.c b/src/partial_sim.c index 58dc4bef..91aea1aa 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -31,6 +31,25 @@ #include "stream.h" +static int gaussian_noise(double expected, double variance) +{ + double x1, x2, w; + double noise; + + do { + + x1 = 2.0 * ((double)random()/RAND_MAX) - 1.0; + x2 = 2.0 * ((double)random()/RAND_MAX) - 1.0; + w = pow(x1, 2.0) + pow(x2, 2.0); + + } while ( w >= 1.0 ); + + w = sqrt((-2.0*log(w))/w); + noise = w * x1; + + return expected + noise*sqrt(variance); +} + static void mess_up_cell(UnitCell *cell) { double ax, ay, az; @@ -38,7 +57,15 @@ static void mess_up_cell(UnitCell *cell) double cx, cy, cz; cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - ax += 0.05*ax; + ax = gaussian_noise(ax, fabs(ax)/50.0); + ay = gaussian_noise(ay, fabs(ay)/50.0); + az = gaussian_noise(az, fabs(az)/50.0); + bx = gaussian_noise(bx, fabs(bx)/50.0); + by = gaussian_noise(by, fabs(by)/50.0); + bz = gaussian_noise(bz, fabs(bz)/50.0); + cx = gaussian_noise(cx, fabs(cx)/50.0); + cy = gaussian_noise(cy, fabs(cy)/50.0); + cz = gaussian_noise(cz, fabs(cz)/50.0); cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz); } |