Commit 7d8662ff authored by Thomas White's avatar Thomas White
Browse files

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

parent 86dd71e8
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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;
}
}
}
......
......@@ -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);
......
......@@ -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 ) {
......
......@@ -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");
}
}
......@@ -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 */
......@@ -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, "--");
}
......
......@@ -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;
......
......@@ -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)
{
......
......@@ -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 */
......@@ -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,
......
/* 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)
......@@ -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:
......
......@@ -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];
}
......
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