diff options
| -rw-r--r-- | src/diffraction.c | 62 |
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); |