From 2c238039c2efda1788ea72c9fb41ff354acc8e97 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 16:04:46 +0100 Subject: Move the "indexed reflection array" thing to where it can't do any harm --- src/diffraction-gpu.c | 530 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/diffraction-gpu.h | 57 ++++++ src/diffraction.c | 464 +++++++++++++++++++++++++++++++++++++++++++ src/diffraction.h | 34 ++++ src/list_tmp.h | 106 ++++++++++ src/pattern_sim.c | 71 +++++++ src/pattern_sim.h | 50 +++++ 7 files changed, 1312 insertions(+) create mode 100644 src/diffraction-gpu.c create mode 100644 src/diffraction-gpu.h create mode 100644 src/diffraction.c create mode 100644 src/diffraction.h create mode 100644 src/list_tmp.h create mode 100644 src/pattern_sim.h (limited to 'src') diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c new file mode 100644 index 00000000..c365cecb --- /dev/null +++ b/src/diffraction-gpu.c @@ -0,0 +1,530 @@ +/* + * diffraction-gpu.c + * + * Calculate diffraction patterns by Fourier methods (GPU version) + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#ifdef HAVE_CL_CL_H +#include +#else +#include +#endif + +#include "image.h" +#include "utils.h" +#include "cell.h" +#include "diffraction.h" +#include "cl-utils.h" +#include "beam-parameters.h" +#include "pattern_sim.h" + + +#define SAMPLING (4) +#define BWSAMPLING (10) +#define DIVSAMPLING (1) +#define SINC_LUT_ELEMENTS (4096) + + +struct gpu_context +{ + cl_context ctx; + cl_command_queue cq; + cl_program prog; + cl_kernel kern; + cl_mem intensities; + cl_mem flags; + + /* Array of sinc LUTs */ + cl_mem *sinc_luts; + cl_float **sinc_lut_ptrs; + int max_sinc_lut; /* Number of LUTs, i.e. one greater than the maximum + * index. This equals the highest allowable "n". */ +}; + + +static void check_sinc_lut(struct gpu_context *gctx, int n) +{ + cl_int err; + cl_image_format fmt; + int i; + + if ( n > gctx->max_sinc_lut ) { + + gctx->sinc_luts = realloc(gctx->sinc_luts, + n*sizeof(*gctx->sinc_luts)); + gctx->sinc_lut_ptrs = realloc(gctx->sinc_lut_ptrs, + n*sizeof(*gctx->sinc_lut_ptrs)); + + for ( i=gctx->max_sinc_lut; isinc_lut_ptrs[i] = NULL; + } + + gctx->max_sinc_lut = n; + } + + if ( gctx->sinc_lut_ptrs[n-1] != NULL ) return; + + /* Create a new sinc LUT */ + gctx->sinc_lut_ptrs[n-1] = malloc(SINC_LUT_ELEMENTS*sizeof(cl_float)); + gctx->sinc_lut_ptrs[n-1][0] = n; + if ( n == 1 ) { + for ( i=1; isinc_lut_ptrs[n-1][i] = 1.0; + } + } else { + for ( i=1; isinc_lut_ptrs[n-1][i] = val; + } + } + + fmt.image_channel_order = CL_INTENSITY; + fmt.image_channel_data_type = CL_FLOAT; + + gctx->sinc_luts[n-1] = clCreateImage2D(gctx->ctx, + CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, + &fmt, SINC_LUT_ELEMENTS, 1, 0, + gctx->sinc_lut_ptrs[n-1], &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't create LUT for %i\n", n); + return; + } +} + + +static int set_arg_float(struct gpu_context *gctx, int idx, float val) +{ + cl_int err; + err = clSetKernelArg(gctx->kern, idx, sizeof(cl_float), &val); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't set kernel argument %i: %s\n", + idx, clError(err)); + return 1; + } + + return 0; +} + + +static int set_arg_int(struct gpu_context *gctx, int idx, int val) +{ + cl_int err; + + err = clSetKernelArg(gctx->kern, idx, sizeof(cl_int), &val); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't set kernel argument %i: %s\n", + idx, clError(err)); + return 1; + } + + return 0; +} + + +static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val) +{ + cl_int err; + + err = clSetKernelArg(gctx->kern, idx, sizeof(cl_mem), &val); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't set kernel argument %i: %s\n", + idx, clError(err)); + return 1; + } + + return 0; +} + + +void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, + int na, int nb, int nc, UnitCell *ucell) +{ + cl_int err; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + float klow, khigh; + int i; + cl_float16 cell; + cl_int4 ncells; + const int sampling = SAMPLING; + cl_float bwstep; + int n_inf = 0; + int n_neg = 0; + cl_float divxlow, divxstep; + cl_float divylow, divystep; + int n_nan = 0; + int sprod; + + if ( gctx == NULL ) { + ERROR("GPU setup failed.\n"); + return; + } + + cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az; + cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz; + cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz; + + /* Calculate wavelength */ + klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0)); + khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0)); + bwstep = (khigh-klow) / BWSAMPLING; + + /* Calculate divergence stuff */ + divxlow = -image->beam->divergence/2.0; + divylow = -image->beam->divergence/2.0; + divxstep = image->beam->divergence / DIVSAMPLING; + divystep = image->beam->divergence / DIVSAMPLING; + + ncells.s[0] = na; + ncells.s[1] = nb; + ncells.s[2] = nc; + ncells.s[3] = 0; /* unused */ + + /* Ensure all required LUTs are available */ + check_sinc_lut(gctx, na); + check_sinc_lut(gctx, nb); + check_sinc_lut(gctx, nc); + + if ( set_arg_float(gctx, 2, klow) ) return; + if ( set_arg_mem(gctx, 9, gctx->intensities) ) return; + if ( set_arg_int(gctx, 12, sampling) ) return; + if ( set_arg_float(gctx, 14, bwstep) ) return; + if ( set_arg_mem(gctx, 15, gctx->sinc_luts[na-1]) ) return; + if ( set_arg_mem(gctx, 16, gctx->sinc_luts[nb-1]) ) return; + if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return; + if ( set_arg_mem(gctx, 18, gctx->flags) ) return; + if ( set_arg_float(gctx, 23, divxlow) ) return; + if ( set_arg_float(gctx, 24, divxstep) ) return; + if ( set_arg_int(gctx, 25, DIVSAMPLING) ) return; + if ( set_arg_float(gctx, 26, divylow) ) return; + if ( set_arg_float(gctx, 27, divystep) ) return; + if ( set_arg_int(gctx, 28, DIVSAMPLING) ) return; + + /* Unit cell */ + err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't set unit cell: %s\n", clError(err)); + return; + } + + /* Local memory for reduction */ + sprod = BWSAMPLING*SAMPLING*SAMPLING*DIVSAMPLING*DIVSAMPLING; + err = clSetKernelArg(gctx->kern, 13, sprod*sizeof(cl_float), NULL); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't set local memory: %s\n", clError(err)); + return; + } + + /* Allocate memory for the result */ + image->data = calloc(image->width * image->height, sizeof(float)); + image->twotheta = calloc(image->width * image->height, sizeof(double)); + + /* Iterate over panels */ + for ( i=0; idet->n_panels; i++ ) { + + size_t dims[3]; + size_t ldims[3] = {SAMPLING, SAMPLING, + BWSAMPLING * DIVSAMPLING * DIVSAMPLING}; + struct panel *p; + cl_mem tt; + size_t tt_size; + cl_mem diff; + size_t diff_size; + float *diff_ptr; + float *tt_ptr; + int pan_width, pan_height; + int fs, ss; + + p = &image->det->panels[i]; + + pan_width = 1 + p->max_fs - p->min_fs; + pan_height = 1 + p->max_ss - p->min_ss; + + /* Buffer for the results of this panel */ + diff_size = pan_width * pan_height * sizeof(cl_float); + diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, + diff_size, NULL, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't allocate diffraction memory\n"); + return; + } + tt_size = pan_width * pan_height * sizeof(cl_float); + tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size, + NULL, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't allocate twotheta memory\n"); + return; + } + + if ( set_arg_mem(gctx, 0, diff) ) return; + if ( set_arg_mem(gctx, 1, tt) ) return; + if ( set_arg_int(gctx, 3, pan_width) ) return; + if ( set_arg_float(gctx, 4, p->cnx) ) return; + if ( set_arg_float(gctx, 5, p->cny) ) return; + if ( set_arg_float(gctx, 6, p->res) ) return; + if ( set_arg_float(gctx, 7, p->clen) ) return; + if ( set_arg_int(gctx, 10, p->min_fs) ) return; + if ( set_arg_int(gctx, 11, p->min_ss) ) return; + if ( set_arg_float(gctx, 19, p->fsx) ) return; + if ( set_arg_float(gctx, 20, p->fsy) ) return; + if ( set_arg_float(gctx, 21, p->ssx) ) return; + if ( set_arg_float(gctx, 22, p->ssy) ) return; + + dims[0] = pan_width * SAMPLING; + dims[1] = pan_height * SAMPLING; + dims[2] = BWSAMPLING * DIVSAMPLING * DIVSAMPLING; + + err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 3, NULL, + dims, ldims, 0, NULL, NULL); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't enqueue diffraction kernel: %s\n", + clError(err)); + return; + } + + clFinish(gctx->cq); + + diff_ptr = clEnqueueMapBuffer(gctx->cq, diff, CL_TRUE, + CL_MAP_READ, 0, diff_size, + 0, NULL, NULL, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't map diffraction buffer: %s\n", + clError(err)); + return; + } + tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ, + 0, tt_size, 0, NULL, NULL, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't map tt buffer\n"); + return; + } + + for ( fs=0; fsmin_fs + fs; + tss = p->min_ss + ss; + image->data[tfs + image->width*tss] = val; + image->twotheta[tfs + image->width*tss] = tt; + + } + } + + clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr, + 0, NULL, NULL); + clEnqueueUnmapMemObject(gctx->cq, tt, tt_ptr, + 0, NULL, NULL); + + clReleaseMemObject(diff); + clReleaseMemObject(tt); + + } + + + if ( n_neg + n_inf + n_nan ) { + ERROR("WARNING: The GPU calculation produced %i negative" + " values, %i infinities and %i NaNs.\n", + n_neg, n_inf, n_nan); + } + +} + + +/* Setup the OpenCL stuff, create buffers, load the structure factor table */ +struct gpu_context *setup_gpu(int no_sfac, + const double *intensities, unsigned char *flags, + const char *sym, int dev_num) +{ + struct gpu_context *gctx; + cl_uint nplat; + cl_platform_id platforms[8]; + cl_context_properties prop[3]; + cl_int err; + cl_device_id dev; + size_t intensities_size; + float *intensities_ptr; + size_t flags_size; + float *flags_ptr; + size_t maxwgsize; + int i; + char cflags[512] = ""; + + STATUS("Setting up GPU...\n"); + + err = clGetPlatformIDs(8, platforms, &nplat); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't get platform IDs: %i\n", err); + return NULL; + } + if ( nplat == 0 ) { + ERROR("Couldn't find at least one platform!\n"); + return NULL; + } + prop[0] = CL_CONTEXT_PLATFORM; + prop[1] = (cl_context_properties)platforms[0]; + prop[2] = 0; + + gctx = malloc(sizeof(*gctx)); + gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU, + NULL, NULL, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't create OpenCL context: %i\n", err); + free(gctx); + return NULL; + } + + dev = get_cl_dev(gctx->ctx, dev_num); + + gctx->cq = clCreateCommandQueue(gctx->ctx, dev, 0, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't create OpenCL command queue\n"); + free(gctx); + return NULL; + } + + /* Create a single-precision version of the scattering factors */ + intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float); + intensities_ptr = malloc(intensities_size); + if ( intensities != NULL ) { + for ( i=0; iintensities = clCreateBuffer(gctx->ctx, + CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, + intensities_size, intensities_ptr, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't allocate intensities memory\n"); + free(gctx); + return NULL; + } + free(intensities_ptr); + + if ( sym != NULL ) { + if ( strcmp(sym, "1") == 0 ) { + strncat(cflags, "-DPG1 ", 511-strlen(cflags)); + } else if ( strcmp(sym, "-1") == 0 ) { + strncat(cflags, "-DPG1BAR ", 511-strlen(cflags)); + } else if ( strcmp(sym, "6/mmm") == 0 ) { + strncat(cflags, "-DPG6MMM ", 511-strlen(cflags)); + } else if ( strcmp(sym, "6") == 0 ) { + strncat(cflags, "-DPG6 ", 511-strlen(cflags)); + } else if ( strcmp(sym, "6/m") == 0 ) { + strncat(cflags, "-DPG6M ", 511-strlen(cflags)); + } else { + ERROR("Sorry! Point group '%s' is not currently" + " supported on the GPU." + " I'm using '1' instead.\n", sym); + strncat(cflags, "-DPG1 ", 511-strlen(cflags)); + } + } else { + if ( intensities != NULL ) { + ERROR("You gave me an intensities file but no point" + " group. I'm assuming '1'.\n"); + strncat(cflags, "-DPG1 ", 511-strlen(cflags)); + } + } + + /* Create a flag array */ + flags_size = IDIM*IDIM*IDIM*sizeof(cl_float); + flags_ptr = malloc(flags_size); + if ( flags != NULL ) { + for ( i=0; iflags = clCreateBuffer(gctx->ctx, + CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, + flags_size, flags_ptr, &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't allocate flag buffer\n"); + free(gctx); + return NULL; + } + free(flags_ptr); + + gctx->prog = load_program(DATADIR"/crystfel/diffraction.cl", gctx->ctx, + dev, &err, cflags); + if ( err != CL_SUCCESS ) { + free(gctx); + return NULL; + } + + gctx->kern = clCreateKernel(gctx->prog, "diffraction", &err); + if ( err != CL_SUCCESS ) { + ERROR("Couldn't create kernel\n"); + free(gctx); + return NULL; + } + + gctx->max_sinc_lut = 0; + gctx->sinc_lut_ptrs = NULL; + gctx->sinc_luts = NULL; + + clGetDeviceInfo(dev, CL_DEVICE_MAX_WORK_GROUP_SIZE, + sizeof(size_t), &maxwgsize, NULL); + STATUS("Maximum work group size = %lli\n", (long long int)maxwgsize); + + return gctx; +} + + +void cleanup_gpu(struct gpu_context *gctx) +{ + int i; + + clReleaseProgram(gctx->prog); + clReleaseMemObject(gctx->intensities); + + /* Release LUTs */ + for ( i=1; i<=gctx->max_sinc_lut; i++ ) { + if ( gctx->sinc_lut_ptrs[i-1] != NULL ) { + clReleaseMemObject(gctx->sinc_luts[i-1]); + free(gctx->sinc_lut_ptrs[i-1]); + } + } + + free(gctx->sinc_luts); + free(gctx->sinc_lut_ptrs); + + clReleaseCommandQueue(gctx->cq); + clReleaseContext(gctx->ctx); + free(gctx); +} diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h new file mode 100644 index 00000000..a3bde4e1 --- /dev/null +++ b/src/diffraction-gpu.h @@ -0,0 +1,57 @@ +/* + * diffraction-gpu.h + * + * Calculate diffraction patterns by Fourier methods (GPU version) + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef DIFFRACTION_GPU_H +#define DIFFRACTION_GPU_H + +#include "image.h" +#include "cell.h" + +struct gpu_context; + +#if HAVE_OPENCL + +extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, + int na, int nb, int nc, UnitCell *ucell); +extern struct gpu_context *setup_gpu(int no_sfac, + const double *intensities, + const unsigned char *flags, + const char *sym, int dev_num); +extern void cleanup_gpu(struct gpu_context *gctx); + +#else + +static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, + int na, int nb, int nc, UnitCell *ucell) +{ + /* Do nothing */ + ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n"); +} + +static struct gpu_context *setup_gpu(int no_sfac, + const double *intensities, + const unsigned char *flags, + const char *sym, int dev_num) +{ + return NULL; +} + +static void cleanup_gpu(struct gpu_context *gctx) +{ +} + +#endif + +#endif /* DIFFRACTION_GPU_H */ diff --git a/src/diffraction.c b/src/diffraction.c new file mode 100644 index 00000000..de994133 --- /dev/null +++ b/src/diffraction.c @@ -0,0 +1,464 @@ +/* + * diffraction.c + * + * Calculate diffraction patterns by Fourier methods + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include +#include +#include +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" +#include "cell.h" +#include "diffraction.h" +#include "beam-parameters.h" +#include "symmetry.h" +#include "pattern_sim.h" + + +#define SAMPLING (4) +#define BWSAMPLING (10) +#define DIVSAMPLING (1) +#define SINC_LUT_ELEMENTS (4096) + + +static double *get_sinc_lut(int n) +{ + int i; + double *lut; + + lut = malloc(SINC_LUT_ELEMENTS*sizeof(double)); + lut[0] = n; + if ( n == 1 ) { + for ( i=1; i= 0.0); + + val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); + val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); + + val1 = val1; + val2 = val2; + + return (1.0-f)*val1 + f*val2; +} + + +static double interpolate_bilinear(const double *ref, + const unsigned char *flags, + const SymOpList *sym, + float hd, float kd, signed int l) +{ + signed int k; + double val1, val2; + float f; + + k = (signed int)kd; + if ( kd < 0.0 ) k -= 1; + f = kd - (float)k; + assert(f >= 0.0); + + val1 = interpolate_linear(ref, flags, sym, hd, k, l); + val2 = interpolate_linear(ref, flags, sym, hd, k+1, l); + + return (1.0-f)*val1 + f*val2; +} + + +static double interpolate_intensity(const double *ref, + const unsigned char *flags, + const SymOpList *sym, + float hd, float kd, float ld) +{ + signed int l; + double val1, val2; + float f; + + l = (signed int)ld; + if ( ld < 0.0 ) l -= 1; + f = ld - (float)l; + assert(f >= 0.0); + + val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l); + val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1); + + return (1.0-f)*val1 + f*val2; +} + + +static double complex interpolate_phased_linear(const double *ref, + const double *phases, + const unsigned char *flags, + const SymOpList *sym, + float hd, + signed int k, signed int l) +{ + signed int h; + double val1, val2; + float f; + double ph1, ph2; + double re1, re2, im1, im2; + double re, im; + + h = (signed int)hd; + if ( hd < 0.0 ) h -= 1; + f = hd - (float)h; + assert(f >= 0.0); + + val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); + val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); + ph1 = sym_lookup_phase(phases, flags, sym, h, k, l); + ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l); + + val1 = val1; + val2 = val2; + + /* Calculate real and imaginary parts */ + re1 = val1 * cos(ph1); + im1 = val1 * sin(ph1); + re2 = val2 * cos(ph2); + im2 = val2 * sin(ph2); + + re = (1.0-f)*re1 + f*re2; + im = (1.0-f)*im1 + f*im2; + + return re + im*I; +} + + +static double complex interpolate_phased_bilinear(const double *ref, + const double *phases, + const unsigned char *flags, + const SymOpList *sym, + float hd, float kd, + signed int l) +{ + signed int k; + double complex val1, val2; + float f; + + k = (signed int)kd; + if ( kd < 0.0 ) k -= 1; + f = kd - (float)k; + assert(f >= 0.0); + + val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l); + val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l); + + return (1.0-f)*val1 + f*val2; +} + + +static double interpolate_phased_intensity(const double *ref, + const double *phases, + const unsigned char *flags, + const SymOpList *sym, + float hd, float kd, float ld) +{ + signed int l; + double complex val1, val2; + float f; + + l = (signed int)ld; + if ( ld < 0.0 ) l -= 1; + f = ld - (float)l; + assert(f >= 0.0); + + val1 = interpolate_phased_bilinear(ref, phases, flags, sym, + hd, kd, l); + val2 = interpolate_phased_bilinear(ref, phases, flags, sym, + hd, kd, l+1); + + return cabs((1.0-f)*val1 + f*val2); +} + + +/* Look up the structure factor for the nearest Bragg condition */ +static double molecule_factor(const double *intensities, const double *phases, + const unsigned char *flags, struct rvec q, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz, + GradientMethod m, const SymOpList *sym) +{ + float hd, kd, ld; + signed int h, k, l; + double r; + + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + + /* No flags -> flat intensity distribution */ + if ( flags == NULL ) return 1.0e5; + + switch ( m ) { + case GRADIENT_MOSAIC : + fesetround(1); /* Round to nearest */ + h = (signed int)rint(hd); + k = (signed int)rint(kd); + l = (signed int)rint(ld); + if ( abs(h) > INDMAX ) r = 0.0; + else if ( abs(k) > INDMAX ) r = 0.0; + else if ( abs(l) > INDMAX ) r = 0.0; + else r = sym_lookup_intensity(intensities, flags, sym, h, k, l); + break; + case GRADIENT_INTERPOLATE : + r = interpolate_intensity(intensities, flags, sym, hd, kd, ld); + break; + case GRADIENT_PHASED : + r = interpolate_phased_intensity(intensities, phases, flags, + sym, hd, kd, ld); + break; + default: + ERROR("This gradient method not implemented yet.\n"); + exit(1); + } + + return r; +} + + +void get_diffraction(struct image *image, int na, int nb, int nc, + const double *intensities, const double *phases, + const unsigned char *flags, UnitCell *cell, + GradientMethod m, const SymOpList *sym) +{ + unsigned int fs, ss; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + float klow, khigh, bwstep; + double *lut_a; + double *lut_b; + double *lut_c; + double divxlow, divylow, divxstep, divystep; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + /* Allocate (and zero) the "diffraction array" */ + image->data = calloc(image->width * image->height, sizeof(float)); + + /* Needed later for Lorentz calculation */ + image->twotheta = malloc(image->width * image->height * sizeof(double)); + + klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0)); + khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0)); + bwstep = (khigh-klow) / BWSAMPLING; + + divxlow = -image->beam->divergence/2.0; + divylow = -image->beam->divergence/2.0; + divxstep = image->beam->divergence / DIVSAMPLING; + divystep = image->beam->divergence / DIVSAMPLING; + + lut_a = get_sinc_lut(na); + lut_b = get_sinc_lut(nb); + lut_c = get_sinc_lut(nc); + + for ( fs=0; fswidth; fs++ ) { + for ( ss=0; ssheight; ss++ ) { + + int fs_step, ss_step, kstep; + int divxval, divyval; + int idx = fs + image->width*ss; + + for ( fs_step=0; fs_stepdata[idx] += intensity; + + if ( fs_step + ss_step + kstep == 0 ) { + image->twotheta[idx] = twotheta; + } + + } + } + } + } + } + + image->data[idx] /= (SAMPLING*SAMPLING*BWSAMPLING + *DIVSAMPLING*DIVSAMPLING); + + + } + progress_bar(fs, image->width-1, "Calculating diffraction"); + } + + free(lut_a); + free(lut_b); + free(lut_c); +} diff --git a/src/diffraction.h b/src/diffraction.h new file mode 100644 index 00000000..f71d3cce --- /dev/null +++ b/src/diffraction.h @@ -0,0 +1,34 @@ +/* + * diffraction.h + * + * Calculate diffraction patterns by Fourier methods + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef DIFFRACTION_H +#define DIFFRACTION_H + +#include "image.h" +#include "cell.h" +#include "symmetry.h" + +typedef enum { + GRADIENT_MOSAIC, + GRADIENT_INTERPOLATE, + GRADIENT_PHASED +} GradientMethod; + +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, + GradientMethod m, const SymOpList *sym); + +#endif /* DIFFRACTION_H */ diff --git a/src/list_tmp.h b/src/list_tmp.h new file mode 100644 index 00000000..a524b2f9 --- /dev/null +++ b/src/list_tmp.h @@ -0,0 +1,106 @@ +/* + * Template for creating indexed 3D lists of a given type, usually indexed + * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX. + * + * These are used, for example, for: + * - a list of 'double complex' values for storing structure factors, + * - a list of 'double' values for storing reflection intensities, + * - a list of 'unsigned int' values for counts of some sort. + * + * When LABEL and TYPE are #defined appropriately, including this header + * defines functions such as: + * - new_list_