Move water calculation to diffraction.c, and work out the consequences

16 files changed, 116 insertions, 121 deletions

diff --git a/data/diffraction.cl b/data/diffraction.cl
index d3fef375..a401cdfa 100644
--- a/data/diffraction.cl
+++ b/data/diffraction.cl
@@ -130,11 +130,11 @@ float2 get_sfac(global float2 *sfacs, float16 cell, float4 q)
 }
 
 
-kernel void diffraction(global float2 *diff, global float *tt, float klow,
+kernel void diffraction(global float *diff, global float *tt, float klow,
                        int w, float cx, float cy,
                        float res, float clen, float16 cell,
                        global float2 *sfacs, float4 z, int4 ncells,
-                       int xmin, int ymin, int sampling, local float2 *tmp,
+                       int xmin, int ymin, int sampling, local float *tmp,
                        float kstep)
 {
 	float ttv;
@@ -149,6 +149,8 @@ kernel void diffraction(global float2 *diff, global float *tt, float klow,
 	float k = klow + kstep * get_local_id(2);
 	const int ax = x / sampling;
 	const int ay = y / sampling;
+	float intensity;
+	float2 val;
 
 	/* Calculate value */
 	q = get_q(x, y, cx, cy, res, clen, k, &ttv, z, sampling);
@@ -156,7 +158,9 @@ kernel void diffraction(global float2 *diff, global float *tt, float klow,
 	f_molecule = get_sfac(sfacs, cell, q);
 
 	/* Write the value to local memory */
-	tmp[lx+sampling*ly+sampling*sampling*lb] = f_molecule * f_lattice;
+	val = f_molecule * f_lattice;
+	intensity = pow(val.x, 2.0f) + pow(val.y, 2.0f);
+	tmp[lx+sampling*ly+sampling*sampling*lb] = intensity;
 
 	/* Memory fence */
 	barrier(CLK_LOCAL_MEM_FENCE);
@@ -165,12 +169,14 @@ kernel void diffraction(global float2 *diff, global float *tt, float klow,
 	if ( lx + ly + lb == 0 ) {
 
 		int i;
-		float2 sum = (float2)(0.0, 0.0);
+		float sum = 0.0;
+		float val;
 
 		for ( i=0; i<sampling*sampling*get_local_size(2); i++ )
 			sum += tmp[i];
 
-		diff[ax+w*ay] = sum / (float)(sampling*sampling*get_local_size(2));
+		val = sum / (float)(sampling*sampling*get_local_size(2));
+		diff[ax+w*ay] = val;
 
 		/* Leader thread also records the 2theta value.
 		 * This should really be averaged across all pixels, but
diff --git a/src/Makefile.am b/src/Makefile.am
index 921208f6..3240290a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,7 +46,7 @@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c
 compare_hkl_LDADD = @LIBS@
 
 powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
-                      detector.c index.c diffraction.c sfac.c
+                      detector.c
 powder_plot_LDADD = @LIBS@
 
 if HAVE_GLIB
diff --git a/src/Makefile.in b/src/Makefile.in
index 9f9618b6..0eccb84c 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -93,11 +93,10 @@ am_pattern_sim_OBJECTS = pattern_sim.$(OBJEXT) diffraction.$(OBJEXT) \
 pattern_sim_OBJECTS = $(am_pattern_sim_OBJECTS)
 pattern_sim_DEPENDENCIES =
 am__powder_plot_SOURCES_DIST = powder_plot.c cell.c utils.c image.c \
-	hdf5-file.c detector.c index.c diffraction.c sfac.c dirax.c
+	hdf5-file.c detector.c dirax.c
 am_powder_plot_OBJECTS = powder_plot.$(OBJEXT) cell.$(OBJEXT) \
 	utils.$(OBJEXT) image.$(OBJEXT) hdf5-file.$(OBJEXT) \
-	detector.$(OBJEXT) index.$(OBJEXT) diffraction.$(OBJEXT) \
-	sfac.$(OBJEXT) $(am__objects_2)
+	detector.$(OBJEXT) $(am__objects_2)
 powder_plot_OBJECTS = $(am_powder_plot_OBJECTS)
 powder_plot_DEPENDENCIES =
 am_process_hkl_OBJECTS = process_hkl.$(OBJEXT) sfac.$(OBJEXT) \
@@ -235,7 +234,7 @@ get_hkl_LDADD = @LIBS@
 compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c
 compare_hkl_LDADD = @LIBS@
 powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
-	detector.c index.c diffraction.c sfac.c $(am__append_5)
+	detector.c $(am__append_5)
 powder_plot_LDADD = @LIBS@
 INCLUDES = "-I$(top_srcdir)/data"
 all: all-am
diff --git a/src/detector.c b/src/detector.c
index d605822b..7a7f5942 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -19,22 +19,42 @@
 #include "utils.h"
 #include "diffraction.h"
 #include "detector.h"
+#include "parameters-lcls.tmp"
 
 
-/* Number of photons in pulse */
-#define FLUENCE (1.0e13)
+/* x,y in pixels relative to central beam */
+int map_position(struct image *image, double dx, double dy,
+                 double *rx, double *ry, double *rz)
+{
+	double d;
+	double twotheta, psi;
+	const double k = 1.0 / image->lambda;
+	struct panel *p;
+	double x = 0.0;
+	double y = 0.0;
+
+	p = find_panel(&image->det, dx, dy);
+
+	x = ((double)dx - p->cx);
+	y = ((double)dy - p->cy);
 
-/* Detector's quantum efficiency (ADU per photon, front detector) */
-#define DQE (167.0)
+	/* Convert pixels to metres */
+	x /= p->res;
+	y /= p->res;	/* Convert pixels to metres */
+	d = sqrt((x*x) + (y*y));
+	twotheta = atan2(d, p->clen);
 
-/* Radius of the water column */
-#define WATER_RADIUS (3.0e-6 / 2.0)
+	psi = atan2(y, x);
 
-/* Radius of X-ray beam */
-#define BEAM_RADIUS (3.0e-6 / 2.0)
+	*rx = k*sin(twotheta)*cos(psi);
+	*ry = k*sin(twotheta)*sin(psi);
+	*rz = k - k*cos(twotheta);
+
+	return 0;
+}
 
 
-void record_image(struct image *image, int do_water, int do_poisson)
+void record_image(struct image *image, int do_poisson)
 {
 	int x, y;
 	double total_energy, energy_density;
@@ -51,39 +71,20 @@ void record_image(struct image *image, int do_water, int do_poisson)
 	       "Total energy = %5.3f microJ\n",
 	       FLUENCE, energy_density/1e7, total_energy*1e6);
 
-	image->hdr = malloc(image->width * image->height * sizeof(double));
-
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 
 		int counts;
 		double cf;
-		double intensity, sa, water;
-		double complex val;
+		double intensity, sa;
 		double pix_area, Lsq;
 		double dsq, proj_area;
 		struct panel *p;
 
-		val = image->sfacs[x + image->width*y];
-		intensity = pow(cabs(val), 2.0);
+		intensity = image->data[x + image->width*y];
 
 		p = find_panel(&image->det, x, y);
 
-		/* FIXME: Move to diffraction.c somehow */
-		if ( do_water ) {
-
-			struct rvec q;
-
-			q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
-
-			/* Add intensity contribution from water */
-			water = water_intensity(q,
-			                        ph_lambda_to_en(image->lambda),
-			                        BEAM_RADIUS, WATER_RADIUS);
-			intensity += water;
-
-		}
-
 		/* Area of one pixel */
 		pix_area = pow(1.0/p->res, 2.0);
 		Lsq = pow(p->clen, 2.0);
@@ -107,20 +108,11 @@ void record_image(struct image *image, int do_water, int do_poisson)
 			counts = (int)rounded;
 		}
 
-		image->hdr[x + image->width*y] = counts;
+		image->data[x + image->width*y] = counts;
 
 	}
 	progress_bar(x, image->width-1, "Post-processing");
 	}
-
-	image->data = malloc(image->width * image->height * sizeof(float));
-	for ( x=0; x<image->width; x++ ) {
-	for ( y=0; y<image->height; y++ ) {
-		int val;
-		val = image->hdr[x + image->width*y];
-		image->data[x + image->width*y] = val;
-	}
-	}
 }
 
 
diff --git a/src/detector.h b/src/detector.h
index f71768c6..d6e740a0 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -38,7 +38,12 @@ struct detector
 	int           n_panels;
 };
 
-extern void record_image(struct image *image, int do_water, int do_poisson);
+
+/* x,y in pixels relative to central beam */
+extern int map_position(struct image *image, double x, double y,
+                        double *rx, double *ry, double *rz);
+
+extern void record_image(struct image *image, int do_poisson);
 
 extern struct panel *find_panel(struct detector *det, int x, int y);
 
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 33410701..a2d1c47a 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -140,7 +140,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	}
 	/* Local memory for reduction */
 	clSetKernelArg(gctx->kern, 15,
-	               BWSAMPLING*SAMPLING*SAMPLING*2*sizeof(cl_float), NULL);
+	               BWSAMPLING*SAMPLING*SAMPLING*sizeof(cl_float), NULL);
 	if ( err != CL_SUCCESS ) {
 		ERROR("Couldn't set arg 15: %s\n", clError(err));
 		return;
@@ -230,20 +230,18 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 
 	free(event);
 
-	image->sfacs = calloc(image->width * image->height,
-	                      sizeof(double complex));
+	image->data = calloc(image->width * image->height, sizeof(float));
 	image->twotheta = calloc(image->width * image->height, sizeof(double));
 
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 
-		float re, im, tt;
+		float val, tt;
 
-		re = diff_ptr[2*(x + image->width*y)+0];
-		im = diff_ptr[2*(x + image->width*y)+1];
+		val = diff_ptr[x + image->width*y];
 		tt = tt_ptr[x + image->width*y];
 
-		image->sfacs[x + image->width*y] = re + I*im;
+		image->data[x + image->width*y] = val;
 		image->twotheta[x + image->width*y] = tt;
 
 	}
@@ -312,7 +310,7 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
 	}
 
 	/* Create buffer for the picture */
-	gctx->diff_size = image->width*image->height*sizeof(cl_float)*2;
+	gctx->diff_size = image->width*image->height*sizeof(cl_float);
 	gctx->diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY,
 	                            gctx->diff_size, NULL, &err);
 	if ( err != CL_SUCCESS ) {
diff --git a/src/diffraction.c b/src/diffraction.c
index b24b68b0..f577be45 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -21,6 +21,7 @@
 #include "cell.h"
 #include "diffraction.h"
 #include "sfac.h"
+#include "parameters-lcls.tmp"
 
 
 #define SAMPLING (4)
@@ -91,8 +92,8 @@ static double complex molecule_factor(struct molecule *mol, struct rvec q,
 }
 
 
-double water_intensity(struct rvec q, double en,
-                       double beam_r, double water_r)
+double water_diffraction(struct rvec q, double en,
+                                 double beam_r, double water_r)
 {
 	double complex fH, fO;
 	double s, modq;
@@ -169,7 +170,8 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 }
 
 
-void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
+void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac,
+                     int do_water)
 {
 	unsigned int xs, ys;
 	double ax, ay, az;
@@ -184,8 +186,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 		                                  &cx, &cy, &cz);
 
 	/* Allocate (and zero) the "diffraction array" */
-	image->sfacs = calloc(image->width * image->height,
-	                      sizeof(double complex));
+	image->data = calloc(image->width * image->height, sizeof(float));
 
 	if ( !no_sfac ) {
 		if ( image->molecule->reflections == NULL ) {
@@ -220,6 +221,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 			float k;
 			double kw = 1.0/BWSAMPLING;
 			double complex val;
+			double intensity;
 
 			/* Calculate k this time round */
 			k = klow + kstep * bwstep;
@@ -238,12 +240,28 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 			                            ax,ay,az,bx,by,bz,cx,cy,cz);
 			}
 
-			val = sw * kw * f_molecule * f_lattice;
-			image->sfacs[x + image->width*y] += val;
+			val = f_molecule * f_lattice;
+			intensity = sw * kw * pow(cabs(val), 2.0);
+			image->data[x + image->width*y] += intensity;
 
 		}
 
+		if ( do_water ) {
+
+			/* Bandwidth not simulated for water */
+			struct rvec q;
+
+			q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
+
+			/* Add intensity contribution from water */
+			image->data[x + image->width*y] += water_diffraction(q,
+			                        ph_lambda_to_en(image->lambda),
+			                        BEAM_RADIUS, WATER_RADIUS);
+
+		}
+
+
 	}
-	progress_bar(xs, SAMPLING*image->width-1, "Calculating lattice factors");
+	progress_bar(xs, SAMPLING*image->width-1, "Calculating diffraction");
 	}
 }
diff --git a/src/diffraction.h b/src/diffraction.h
index f26d9a4e..503e9d19 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -20,9 +20,7 @@
 #include "cell.h"
 
 extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            int nosfac);
-extern double water_intensity(struct rvec q, double en,
-                              double beam_r, double water_r);
+                            int nosfac, int do_water);
 extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
                          unsigned int sampling, float *ttp, float k);
 #endif	/* DIFFRACTION_H */
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 5831d034..2536c1ed 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -618,7 +618,7 @@ static void numbers_update(DisplayWindow *dw)
 	for ( py=0; py<17; py++ ) {
 
 		char s[32];
-		int16_t val;
+		float val;
 		GtkWidget *l;
 		int x, y;
 		int valid;
@@ -635,7 +635,7 @@ static void numbers_update(DisplayWindow *dw)
 		}
 
 		if ( (x>0) && (y>0) && valid ) {
-			snprintf(s, 31, "%i", val);
+			snprintf(s, 31, "%.0f", val);
 		} else {
 			strcpy(s, "--");
 		}
diff --git a/src/image.h b/src/image.h
index a8103b92..349e680c 100644
--- a/src/image.h
+++ b/src/image.h
@@ -69,9 +69,7 @@ struct rvec
 /* Structure describing an image */
 struct image {
 
-	int			*hdr;      /* Actual counts */
 	float			*data;     /* Integer counts after bloom */
-	double complex		*sfacs;
 	double			*twotheta;
 	struct molecule		*molecule;
 	UnitCell		*indexed_cell;
diff --git a/src/index.c b/src/index.c
index 7439cf66..865715c2 100644
--- a/src/index.c
+++ b/src/index.c
@@ -29,37 +29,6 @@
 #include "index.h"
 
 
-/* x,y in pixels relative to central beam */
-int map_position(struct image *image, double dx, double dy,
-                 double *rx, double *ry, double *rz)
-{
-	double d;
-	double twotheta, psi;
-	const double k = 1.0 / image->lambda;
-	struct panel *p;
-	double x = 0.0;
-	double y = 0.0;
-
-	p = find_panel(&image->det, dx, dy);
-
-	x = ((double)dx - p->cx);
-	y = ((double)dy - p->cy);
-
-	/* Convert pixels to metres */
-	x /= p->res;
-	y /= p->res;	/* Convert pixels to metres */
-	d = sqrt((x*x) + (y*y));
-	twotheta = atan2(d, p->clen);
-
-	psi = atan2(y, x);
-
-	*rx = k*sin(twotheta)*cos(psi);
-	*ry = k*sin(twotheta)*sin(psi);
-	*rz = k - k*cos(twotheta);
-
-	return 0;
-}
-
 
 static void write_drx(struct image *image)
 {
diff --git a/src/index.h b/src/index.h
index 08ffa2fc..429e5ec8 100644
--- a/src/index.h
+++ b/src/index.h
@@ -27,9 +27,5 @@ typedef enum {
 
 extern void index_pattern(struct image *image, IndexingMethod indm,
                           int no_match);
-/* x,y in pixels relative to central beam */
-extern int map_position(struct image *image, double x, double y,
-                        double *rx, double *ry, double *rz);
-
 
 #endif	/* INDEX_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 5204b15e..d3cbfcb6 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -223,9 +223,8 @@ int main(int argc, char *argv[])
 			if ( config_nearbragg || config_simulate ) {
 
 				/* Simulate a diffraction pattern */
-				image.sfacs = NULL;
 				image.twotheta = NULL;
-				image.hdr = NULL;
+				image.data = NULL;
 
 				/* View head-on (unit cell is tilted) */
 				image.orientation.w = 1.0;
@@ -255,13 +254,13 @@ int main(int argc, char *argv[])
 					get_diffraction_gpu(gctx, &image,
 					                    8, 8, 8);
 				} else {
-					get_diffraction(&image, 8, 8, 8, 0);
+					get_diffraction(&image, 8, 8, 8, 0, 0);
 				}
 				if ( image.molecule == NULL ) {
 					ERROR("Couldn't open molecule.pdb\n");
 					return 1;
 				}
-				record_image(&image, 0, 0);
+				record_image(&image, 0);
 
 				hdf5_write("simulated.h5", image.data,
 				           image.width, image.height,
diff --git a/src/parameters-lcls.tmp b/src/parameters-lcls.tmp
new file mode 100644
index 00000000..96b2bc33
--- /dev/null
+++ b/src/parameters-lcls.tmp
@@ -0,0 +1,11 @@
+/* Number of photons in pulse */
+#define FLUENCE (1.0e13)
+
+/* Detector's quantum efficiency (ADU per photon, front detector) */
+#define DQE (167.0)
+
+/* Radius of the water column */
+#define WATER_RADIUS (3.0e-6 / 2.0)
+
+/* Radius of X-ray beam */
+#define BEAM_RADIUS (3.0e-6 / 2.0)
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index a6a2156b..91100426 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -106,6 +106,9 @@ static void show_details()
 "algorithm.  When the intensity is sufficiently high that Knuth's algorithm\n"
 "would result in machine precision problems, a normal distribution with\n"
 "standard deviation sqrt(I) is used instead.\n"
+"\n"
+"The coherent, elastic part of the diffuse scattering from the water jet can\n"
+"be simulated.\n"
 );
 }
 
@@ -215,6 +218,12 @@ int main(int argc, char *argv[])
 		return 0;
 	}
 
+	if ( (!config_nowater) && config_gpu ) {
+		ERROR("Cannot simulate water scattering on the GPU.\n");
+		ERROR("Please try again with the --no-water option.\n");
+		return 1;
+	}
+
 	/* Define image parameters */
 	image.width = 1024;
 	image.height = 1024;
@@ -257,10 +266,8 @@ int main(int argc, char *argv[])
 		}
 
 		/* Ensure no residual information */
-		image.sfacs = NULL;
 		image.data = NULL;
 		image.twotheta = NULL;
-		image.hdr = NULL;
 
 		cell_get_parameters(image.molecule->cell, &a, &b, &c, &d, &d, &d);
 		STATUS("Particle size = %i x %i x %i (=%5.2f x %5.2f x %5.2f nm)\n",
@@ -273,18 +280,19 @@ int main(int argc, char *argv[])
 			}
 			get_diffraction_gpu(gctx, &image, na, nb, nc);
 		} else {
-			get_diffraction(&image, na, nb, nc, config_nosfac);
+			get_diffraction(&image, na, nb, nc, config_nosfac,
+			                !config_nowater);
 		}
 		if ( image.molecule == NULL ) {
 			ERROR("Couldn't open molecule.pdb\n");
 			return 1;
 		}
-		if ( image.sfacs == NULL ) {
+		if ( image.data == NULL ) {
 			ERROR("Diffraction calculation failed.\n");
 			goto skip;
 		}
 
-		record_image(&image, !config_nowater, !config_nonoise);
+		record_image(&image, !config_nonoise);
 
 		if ( config_nearbragg ) {
 			output_intensities(&image, image.molecule->cell);
@@ -324,8 +332,6 @@ int main(int argc, char *argv[])
 
 		/* Clean up */
 		free(image.data);
-		free(image.hdr);
-		free(image.sfacs);
 		free(image.twotheta);
 
 skip:
diff --git a/src/render.c b/src/render.c
index b02469f6..c91a2658 100644
--- a/src/render.c
+++ b/src/render.c
@@ -52,7 +52,7 @@ static void *render_bin(float *in, int inw, int inh, int binning, float *maxp)
 		}
 		}
 
-		data[x+w*y] = total / (binning * binning);
+		data[x+w*y] = total / ((double)binning * (double)binning);
 		if ( data[x+w*y] > max ) max = data[x+w*y];
 
 	}