Proper water intensity calculation

Intensity is calculated at the detector stage now

3 files changed, 44 insertions, 8 deletions

diff --git a/src/detector.c b/src/detector.c
index f4de05ff..dc378480 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -156,14 +156,16 @@ void record_image(struct image *image)
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 
-		double counts;
-		double intensity;
-		double sa;
+		double counts,  intensity, sa;
 		double complex val;
 
 		val = image->sfacs[x + image->width*y];
 		intensity = (double)(val * conj(val));
 
+		/* Add intensity contribution from water */
+		intensity += water_intensity(image->qvecs[x + image->width*y],
+		                            image->xray_energy);
+
 		/* What solid angle is subtended by this pixel? */
 		sa = sa_per_pixel * cos(image->twotheta[x + image->width*y]);
 
diff --git a/src/diffraction.c b/src/diffraction.c
index 9a9e31b6..07749ac7 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -99,9 +99,44 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q,
 }
 
 
-static double complex water_factor(struct threevec q, double en)
+double water_intensity(struct threevec q, double en)
 {
-	return 0.0;
+	complex double fH, fO;
+	double s, modq;
+	double intensity;
+
+	/* Interatomic distances in water molecule */
+	const double rOH = 0.09584e-9;
+	const double rHH = 0.1515e-9;
+
+	/* Dimensions of water column */
+	const double water_r = 0.5e-6;
+	const double beam_r = 1.5e-6;
+
+	/* Volume of water column */
+	const double water_v = M_PI*pow(water_r, 2.0) * 2.0 * beam_r;
+
+	/* Number of water molecules */
+	const double n_water = water_v * WATER_DENSITY
+	                                        * (AVOGADRO / WATER_MOLAR_MASS);
+
+	/* s = sin(theta)/lambda = 1/2d = |q|/2 */
+	modq = modulus(q.u, q.v, q.w);
+	s = modq / 2.0;
+
+	fH = get_sfac("H", s, en);
+	fO = get_sfac("O", s, en);
+
+	/* Four O-H cross terms */
+	intensity = 4.0*fH*fO * sin(modq*rOH)/(modq*rOH);
+
+	/* Three H-H cross terms */
+	intensity += 3.0*fH*fH * sin(modq*rHH)/(modq*rHH);
+
+	/* Three diagonal terms */
+	intensity += 2.0*fH*fH + fO*fO;
+
+	return intensity * n_water;
 }
 
 
@@ -137,7 +172,6 @@ void get_diffraction(struct image *image)
 
 		double f_lattice;
 		double complex f_molecule;
-		double complex f_water;
 		struct threevec q;
 		double complex val;
 
@@ -147,8 +181,7 @@ void get_diffraction(struct image *image)
 		f_molecule = molecule_factor(image->molecule, q,
 		                             ax,ay,az,bx,by,bz,cx,cy,cz);
 
-		f_water = water_factor(q, image->xray_energy);
-		val = (f_molecule * f_lattice) + f_water;
+		val = f_molecule * f_lattice;
 		image->sfacs[x + image->width*y] = val;
 
 	}
diff --git a/src/diffraction.h b/src/diffraction.h
index 112818bc..e2994a0b 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -20,5 +20,6 @@
 #include "cell.h"
 
 extern void get_diffraction(struct image *image);
+extern double water_intensity(struct threevec q, double en);
 
 #endif	/* DIFFRACTION_H */