Use LUTs for CPU simulation as well

2 files changed, 57 insertions, 26 deletions

diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 9b96f112..e7449148 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -36,7 +36,6 @@
 
 #define SAMPLING (4)
 #define BWSAMPLING (10)
-
 #define SINC_LUT_ELEMENTS (4096)
 
 
diff --git a/src/diffraction.c b/src/diffraction.c
index a7cf90c7..13a3dd95 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -29,12 +29,53 @@
 
 #define SAMPLING (4)
 #define BWSAMPLING (10)
+#define SINC_LUT_ELEMENTS (4096)
+
+
+static double *get_sinc_lut(int n)
+{
+	int i;
+	double *lut;
+
+	lut = malloc(SINC_LUT_ELEMENTS*sizeof(double));
+	lut[0] = n;
+	if ( n == 1 ) {
+		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
+			lut[i] = 1.0;
+		}
+	} else {
+		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
+			double x, val;
+			x = (double)i/SINC_LUT_ELEMENTS;
+			val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
+			lut[i] = val;
+		}
+	}
+
+	return lut;
+}
+
+
+static double interpolate_lut(double *lut, double val)
+{
+	double i, pos, f;
+	unsigned int low, high;
+
+	pos = SINC_LUT_ELEMENTS * modf(val, &i);
+	low = (int)pos;  /* Discard fractional part */
+	high = low + 1;
+	f = modf(pos, &i);  /* Fraction */
+	if ( high == SINC_LUT_ELEMENTS ) high = 0;
+
+	return (1.0-f)*low + f*high;
+}
 
 
 static double lattice_factor(struct rvec q, double ax, double ay, double az,
                                             double bx, double by, double bz,
                                             double cx, double cy, double cz,
-                                            int na, int nb, int nc)
+                                            double *lut_a, double *lut_b,
+                                            double *lut_c)
 {
 	struct rvec Udotq;
 	double f1, f2, f3;
@@ -43,30 +84,10 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
 	Udotq.v = bx*q.u + by*q.v + bz*q.w;
 	Udotq.w = cx*q.u + cy*q.v + cz*q.w;
 
-	/* At exact Bragg condition, f1 = na */
-	if ( na > 1 ) {
-		f1 = sin(M_PI*(double)na*Udotq.u) / sin(M_PI*Udotq.u);
-	} else {
-		f1 = 1.0;
-	}
+	f1 = interpolate_lut(lut_a, Udotq.u);
+	f2 = interpolate_lut(lut_b, Udotq.u);
+	f3 = interpolate_lut(lut_c, Udotq.u);
 
-	/* At exact Bragg condition, f2 = nb */
-	if ( nb > 1 ) {
-		f2 = sin(M_PI*(double)nb*Udotq.v) / sin(M_PI*Udotq.v);
-	} else {
-		f2 = 1.0;
-	}
-
-	/* At exact Bragg condition, f3 = nc */
-	if ( nc > 1 ) {
-		f3 = sin(M_PI*(double)nc*Udotq.w) / sin(M_PI*Udotq.w);
-	} else {
-		f3 = 1.0;
-	}
-
-	/* At exact Bragg condition, this will multiply the molecular
-	 * part of the structure factor by the number of unit cells,
-	 * as desired (more scattering from bigger crystal!) */
 	return f1 * f2 * f3;
 }
 
@@ -331,6 +352,9 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 	double bx, by, bz;
 	double cx, cy, cz;
 	float klow, khigh, bwstep;
+	double *lut_a;
+	double *lut_b;
+	double *lut_c;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
@@ -344,6 +368,10 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 	khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0));
 	bwstep = (khigh-klow) / BWSAMPLING;
 
+	lut_a = get_sinc_lut(na);
+	lut_b = get_sinc_lut(nb);
+	lut_c = get_sinc_lut(nc);
+
 	for ( fs=0; fs<image->width; fs++ ) {
 	for ( ss=0; ss<image->height; ss++ ) {
 
@@ -371,7 +399,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 			f_lattice = lattice_factor(q, ax, ay, az,
 			                              bx, by, bz,
 			                              cx, cy, cz,
-			                              na, nb, nc);
+			                              lut_a, lut_b, lut_c);
 
 			I_molecule = molecule_factor(intensities,
 			                             phases, flags, q,
@@ -397,4 +425,8 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 	}
 	progress_bar(fs, image->width-1, "Calculating diffraction");
 	}
+
+	free(lut_a);
+	free(lut_b);
+	free(lut_c);
 }