aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2009-11-18 21:23:23 +0100
committerThomas White <taw@bitwiz.org.uk>2009-11-18 21:23:23 +0100
commitdcd717de4415c007cf47faf85ec382549f129797 (patch)
tree3eb4bebe2a70a6e480c169545a007b8b747b9bca
parent7f004b45f39c9207a6a233217f983c9b269bca45 (diff)
Add water calculation
-rw-r--r--src/diffraction.c62
1 files changed, 60 insertions, 2 deletions
diff --git a/src/diffraction.c b/src/diffraction.c
index 25b8c3b4..461063ad 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -24,6 +24,16 @@
#include "sfac.h"
+/* Density of water in kg/m^3 */
+#define WATER_DENSITY (1.0e6)
+
+/* Molar mass of water, in kg/mol */
+#define WATER_MOLAR_MASS (18.01528e3)
+
+/* Avogadro's number */
+#define AVOGADRO (6.022e23)
+
+
static double lattice_factor(struct threevec q, double ax, double ay, double az,
double bx, double by, double bz,
double cx, double cy, double cz)
@@ -85,7 +95,7 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q,
double ph;
- ph= q.u*spec->x[j] + q.v*spec->y[j] + q.w*spec->z[j];
+ ph = q.u*spec->x[j] + q.v*spec->y[j] + q.w*spec->z[j];
/* Conversion from revolutions to radians is required */
contrib += cos(2.0*M_PI*ph) + I*sin(2.0*M_PI*ph);
@@ -101,6 +111,45 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q,
}
+double complex water_factor(struct threevec q, double en)
+{
+ double x;
+ double complex res = 0.0;
+ int n = 0;
+ double s;
+ double molecules_per_m3;
+ double molecules_per_m2;
+ const double rc = 0.5e-6; /* Radius of cylinder */
+ const double rb = 1.5e-6; /* Radius of beam */
+
+ molecules_per_m3 = WATER_DENSITY * (AVOGADRO / WATER_MOLAR_MASS);
+ molecules_per_m2 = molecules_per_m3 / (2*rb*2*rc);
+
+ /* s = sin(theta)/lambda = 1/2d = (1/d)/2.0 */
+ s = modulus(q.u, q.v, q.w) / 2.0;
+
+ for ( x=-rc; x<rc; x+=rc/50.0 ) {
+
+ double ph;
+ double thickness;
+
+ /* The thickness of a cylinder, axis
+ * perpendicular to the beam */
+ thickness = 2.0 * sqrt((rc*rc)-(x*x));
+
+ ph = q.u*x;
+ res += thickness * (cos(2.0*M_PI*ph) + I*sin(2.0*M_PI*ph));
+ n++;
+
+ }
+ res *= (get_sfac("O", s, en) + 2.0*get_sfac("H", s, en))
+ * molecules_per_m2;
+
+ /* Correct for sampling of integration */
+ return (res/n) * (2*rc);
+}
+
+
void get_diffraction(struct image *image, UnitCell *cell)
{
int x, y;
@@ -125,6 +174,7 @@ void get_diffraction(struct image *image, UnitCell *cell)
double f_lattice;
double complex f_molecule;
+ double complex f_water;
struct threevec q;
q = image->qvecs[x + image->width*y];
@@ -133,7 +183,15 @@ void get_diffraction(struct image *image, UnitCell *cell)
f_molecule = molecule_factor(image->molecule, q,
image->xray_energy);
- image->sfacs[x + image->width*y] = f_lattice * f_molecule;
+ /* Nasty approximation follows */
+ if ( y == image->height/2 ) {
+ f_water = water_factor(q, image->xray_energy);
+ } else {
+ f_water = 0.0;
+ }
+
+ image->sfacs[x + image->width*y] = (f_lattice * f_molecule)
+ + f_water;
}
progress_bar(x, image->width-1);