aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--data/diffraction.cl16
-rw-r--r--src/Makefile.am2
-rw-r--r--src/Makefile.in7
-rw-r--r--src/detector.c72
-rw-r--r--src/detector.h7
-rw-r--r--src/diffraction-gpu.c14
-rw-r--r--src/diffraction.c34
-rw-r--r--src/diffraction.h4
-rw-r--r--src/displaywindow.c4
-rw-r--r--src/image.h2
-rw-r--r--src/index.c31
-rw-r--r--src/index.h4
-rw-r--r--src/indexamajig.c7
-rw-r--r--src/parameters-lcls.tmp11
-rw-r--r--src/pattern_sim.c20
-rw-r--r--src/render.c2
16 files changed, 116 insertions, 121 deletions
diff --git a/data/diffraction.cl b/data/diffraction.cl
index d3fef375..a401cdfa 100644
--- a/data/diffraction.cl
+++ b/data/diffraction.cl
@@ -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
diff --git a/src/Makefile.am b/src/Makefile.am
index 921208f6..3240290a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -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
diff --git a/src/Makefile.in b/src/Makefile.in
index 9f9618b6..0eccb84c 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -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
diff --git a/src/detector.c b/src/detector.c
index d605822b..7a7f5942 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -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;
- }
- }
}
diff --git a/src/detector.h b/src/detector.h
index f71768c6..d6e740a0 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -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);
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 33410701..a2d1c47a 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -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 ) {
diff --git a/src/diffraction.c b/src/diffraction.c
index b24b68b0..f577be45 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -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");
}
}
diff --git a/src/diffraction.h b/src/diffraction.h
index f26d9a4e..503e9d19 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -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 */
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 5831d034..2536c1ed 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -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, "--");
}
diff --git a/src/image.h b/src/image.h
index a8103b92..349e680c 100644
--- a/src/image.h
+++ b/src/image.h
@@ -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;
diff --git a/src/index.c b/src/index.c
index 7439cf66..865715c2 100644
--- a/src/index.c
+++ b/src/index.c
@@ -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)
{
diff --git a/src/index.h b/src/index.h
index 08ffa2fc..429e5ec8 100644
--- a/src/index.h
+++ b/src/index.h
@@ -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 */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 5204b15e..d3cbfcb6 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -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,
diff --git a/src/parameters-lcls.tmp b/src/parameters-lcls.tmp
new file mode 100644
index 00000000..96b2bc33
--- /dev/null
+++ b/src/parameters-lcls.tmp
@@ -0,0 +1,11 @@
+/* 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)
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index a6a2156b..91100426 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -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:
diff --git a/src/render.c b/src/render.c
index b02469f6..c91a2658 100644
--- a/src/render.c
+++ b/src/render.c
@@ -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];
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-10 23:45:16 +0000