From 7d8662ffe897dc2438141ecc8848863bad9b9d92 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 26 Feb 2010 16:47:27 +0100 Subject: Move water calculation to diffraction.c, and work out the consequences --- src/Makefile.am | 2 +- src/Makefile.in | 7 +++-- src/detector.c | 72 ++++++++++++++++++++++--------------------------- src/detector.h | 7 ++++- src/diffraction-gpu.c | 14 +++++----- src/diffraction.c | 34 +++++++++++++++++------ src/diffraction.h | 4 +-- src/displaywindow.c | 4 +-- src/image.h | 2 -- src/index.c | 31 --------------------- src/index.h | 4 --- src/indexamajig.c | 7 +++-- src/parameters-lcls.tmp | 11 ++++++++ src/pattern_sim.c | 20 +++++++++----- src/render.c | 2 +- 15 files changed, 105 insertions(+), 116 deletions(-) create mode 100644 src/parameters-lcls.tmp (limited to 'src') 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; xwidth; x++ ) { for ( y=0; yheight; 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; xwidth; x++ ) { - for ( y=0; yheight; 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; xwidth; x++ ) { for ( y=0; yheight; 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]; } -- cgit v1.2.3