aboutsummaryrefslogtreecommitdiff
path: root/src/partial_sim.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/partial_sim.c')
-rw-r--r--src/partial_sim.c29
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);
}