Commit 43337144 authored by Thomas White's avatar Thomas White
Browse files

Remove water simulation stuff

It's not really useful - our dominant background is from elsewhere
parent 243e7220
......@@ -41,7 +41,6 @@ struct beam_params *get_beam_parameters(const char *filename)
b->divergence = -1.0;
b->dqe = -1.0;
b->adu_per_photon = -1.0;
b->water_radius = -1.0;
do {
......@@ -81,8 +80,6 @@ struct beam_params *get_beam_parameters(const char *filename)
b->dqe = atof(bits[2]);
} else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) {
b->adu_per_photon = atof(bits[2]);
} else if ( strcmp(bits[0], "jet/radius") == 0 ) {
b->water_radius = atof(bits[2]);
} else {
ERROR("Unrecognised field '%s'\n", bits[0]);
}
......@@ -125,10 +122,6 @@ struct beam_params *get_beam_parameters(const char *filename)
" 'detector/adu_per_photon'.\n");
reject = 1;
}
if ( b->water_radius < 0.0 ) {
ERROR("Invalid or unspecified value for 'jet/radius'.\n");
reject = 1;
}
if ( reject ) {
ERROR("Please fix the above problems with the beam"
......
......@@ -31,8 +31,6 @@ struct beam_params
double dqe; /* Detector DQE (fraction) */
double adu_per_photon; /* Detector "gain" */
double water_radius; /* metres */
};
......
......@@ -319,55 +319,10 @@ static double molecule_factor(const double *intensities, const double *phases,
}
double water_diffraction(struct rvec q, double en,
double beam_r, double water_r)
{
double complex fH, fO;
double s, modq;
double width;
double complex ifac;
/* Interatomic distances in water molecule */
const double rOH = 0.09584e-9;
const double rHH = 0.1515e-9;
/* Volume of water column, approximated as:
* (2water_r) * (2beam_r) * smallest(2beam_r, 2water_r)
* neglecting the curvature of the faces of the volume */
if ( beam_r > water_r ) {
width = 2.0 * water_r;
} else {
width = 2.0 * beam_r;
}
const double water_v = 2.0*beam_r * 2.0*water_r * width;
/* 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 */
ifac = 4.0*fH*fO * sin(2.0*M_PI*modq*rOH)/(2.0*M_PI*modq*rOH);
/* Three H-H cross terms */
ifac += 3.0*fH*fH * sin(2.0*M_PI*modq*rHH)/(2.0*M_PI*modq*rHH);
/* Three diagonal terms */
ifac += 2.0*fH*fH + fO*fO;
return cabs(ifac) * n_water;
}
void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities, const double *phases,
const unsigned char *flags, UnitCell *cell, int do_water,
const unsigned char *flags, UnitCell *cell,
GradientMethod m, const char *sym)
{
unsigned int xs, ys;
......@@ -436,21 +391,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
}
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),
image->beam->beam_radius,
image->beam->water_radius) * sw;
}
}
progress_bar(xs, SAMPLING*image->width-1, "Calculating diffraction");
......
......@@ -28,6 +28,6 @@ typedef enum {
extern void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities, const double *phases,
const unsigned char *flags, UnitCell *cell,
int do_water, GradientMethod m, const char *sym);
GradientMethod m, const char *sym);
#endif /* DIFFRACTION_H */
......@@ -97,7 +97,6 @@ static void show_help(const char *s)
"speed, or for testing, you can choose to disable certain things using the\n"
"following options.\n"
"\n"
" --no-water Do not simulate water background.\n"
" --no-noise Do not calculate Poisson noise.\n"
);
}
......@@ -124,10 +123,6 @@ static void show_details()
"This can be approximated to varying levels of accuracy by the methods given by\n"
"the '--gradients' option.\n"
"\n"
"Intensity from water is added according to the first term of equation 5\n"
"from Phys Chem Chem Phys 2003 (5) 1981--1991. This simulates the\n"
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
......@@ -208,7 +203,6 @@ int main(int argc, char *argv[])
int config_nearbragg = 0;
int config_randomquat = 0;
int config_noimages = 0;
int config_nowater = 0;
int config_nonoise = 0;
int config_nosfac = 0;
int config_gpu = 0;
......@@ -241,7 +235,6 @@ int main(int argc, char *argv[])
{"random-orientation", 0, NULL, 'r'},
{"number", 1, NULL, 'n'},
{"no-images", 0, &config_noimages, 1},
{"no-water", 0, &config_nowater, 1},
{"no-noise", 0, &config_nonoise, 1},
{"intensities", 1, NULL, 'i'},
{"symmetry", 1, NULL, 'y'},
......@@ -375,12 +368,6 @@ 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;
}
if ( grad_str == NULL ) {
STATUS("You didn't specify a gradient calculation method, so"
" I'm using the 'mosaic' method, which is fastest.\n");
......@@ -553,8 +540,7 @@ int main(int argc, char *argv[])
get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
} else {
get_diffraction(&image, na, nb, nc, intensities, phases,
flags, cell, !config_nowater, grad,
sym);
flags, cell, grad, sym);
}
if ( image.data == NULL ) {
ERROR("Diffraction calculation failed.\n");
......
......@@ -240,7 +240,7 @@ static double get_waas_kirf(const char *n, double s)
/* Get complex scattering factors for element 'n' at energy 'en' (J/photon),
* at resolution 's' = sin(theta)/lambda (in m^-1) */
double complex get_sfac(const char *n, double s, double en)
static double complex get_sfac(const char *n, double s, double en)
{
double complex f1f2;
double fq;
......
......@@ -55,10 +55,6 @@ struct molecule
double zc;
};
/* This is so that the water background calculation can use it */
extern double complex get_sfac(const char *n, double s, double en);
extern struct molecule *load_molecule(const char *filename);
extern void free_molecule(struct molecule *mol);
extern double *get_reflections(struct molecule *mol, double en, double res,
......
......@@ -42,15 +42,6 @@
/* Thomson scattering length (m) */
#define THOMSON_LENGTH (2.81794e-15)
/* 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)
/* ------------------------------ Quaternions ------------------------------- */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment