aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2020-08-04 17:56:07 +0200
committerThomas White <taw@physics.org>2020-08-04 17:56:07 +0200
commitee5e7f49c4c582651d1c3760f204ac861864c9d5 (patch)
tree4afd939ebe2a739fd02b183a68803086d655a958 /src
parentf9d454ae2a5c2c3ac2d7fc7161a1f02fc87da612 (diff)
Convert pattern_sim to DataTemplate
Diffstat (limited to 'src')
-rw-r--r--src/diffraction-gpu.c17
-rw-r--r--src/diffraction.c13
-rw-r--r--src/pattern_sim.c235
3 files changed, 91 insertions, 174 deletions
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 419b5957..74f03ebf 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -38,6 +38,7 @@
#include <string.h>
#include <complex.h>
+#define CL_TARGET_OPENCL_VERSION 220
#ifdef HAVE_CL_CL_H
#include <CL/cl.h>
#else
@@ -187,18 +188,18 @@ static int do_panels(struct gpu_context *gctx, struct image *image,
if ( set_arg_float(gctx, 2, weight) ) return 1;
/* Iterate over panels */
- for ( i=0; i<image->det->n_panels; i++ ) {
+ for ( i=0; i<image->detgeom->n_panels; i++ ) {
size_t dims[2];
size_t ldims[2];
- struct panel *p;
+ struct detgeom_panel *p;
cl_mem diff;
size_t diff_size;
float *diff_ptr;
int fs, ss;
cl_int err;
- p = &image->det->panels[i];
+ p = &image->detgeom->panels[i];
/* Buffer for the results of this panel */
diff_size = p->w * p->h * sizeof(cl_float);
@@ -220,8 +221,8 @@ static int do_panels(struct gpu_context *gctx, struct image *image,
if ( set_arg_float(gctx, 9, p->ssx) ) return 1;
if ( set_arg_float(gctx, 10, p->ssy) ) return 1;
if ( set_arg_float(gctx, 11, p->ssz) ) return 1;
- if ( set_arg_float(gctx, 12, p->res) ) return 1;
- if ( set_arg_float(gctx, 13, p->clen) ) return 1;
+ if ( set_arg_float(gctx, 12, 1.0/p->pixel_pitch) ) return 1;
+ if ( set_arg_float(gctx, 13, p->cnz*p->pixel_pitch) ) return 1;
dims[0] = p->w * sampling;
dims[1] = p->h * sampling;
@@ -325,13 +326,13 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
if ( set_arg_mem(gctx, 19, gctx->sinc_luts[nc-1]) ) return 1;
/* Allocate memory for the result */
- image->dp = malloc(image->det->n_panels * sizeof(float *));
+ image->dp = malloc(image->detgeom->n_panels * sizeof(float *));
if ( image->dp == NULL ) {
ERROR("Couldn't allocate memory for result.\n");
return 1;
}
- for ( i=0; i<image->det->n_panels; i++ ) {
- struct panel *p = &image->det->panels[i];
+ for ( i=0; i<image->detgeom->n_panels; i++ ) {
+ struct detgeom_panel *p = &image->detgeom->panels[i];
image->dp[i] = calloc(p->w * p->h, sizeof(float));
if ( image->dp[i] == NULL ) {
ERROR("Couldn't allocate memory for panel %i\n", i);
diff --git a/src/diffraction.c b/src/diffraction.c
index 8762b0af..4f17a27d 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -376,7 +376,7 @@ static void diffraction_panel(struct image *image, const double *intensities,
int fs, ss;
const int nxs = 4;
const int nys = 4;
- struct panel *p = &image->det->panels[pn];
+ struct detgeom_panel *p = &image->detgeom->panels[pn];
weight /= nxs*nys;
@@ -386,17 +386,22 @@ static void diffraction_panel(struct image *image, const double *intensities,
int idx;
double f_lattice, I_lattice;
double I_molecule;
- struct rvec q;
int xs, ys;
float xo, yo;
for ( xs=0; xs<nxs; xs++ ) {
for ( ys=0; ys<nys; ys++ ) {
+ double qv[3];
+ struct rvec q;
+
xo = (1.0/nxs) * xs;
yo = (1.0/nys) * ys;
- q = get_q_for_panel(p, fs+xo, ss+yo, NULL, k);
+ detgeom_transform_coords(p, fs+xo, ss+yo,
+ image->lambda, qv);
+
+ q.u = qv[0]; q.v = qv[1]; q.w = qv[2];
f_lattice = lattice_factor(q, ax, ay, az,
bx, by, bz,
@@ -435,7 +440,7 @@ static void diffraction_at_k(struct image *image, const double *intensities,
{
int i;
- for ( i=0; i<image->det->n_panels; i++ ) {
+ for ( i=0; i<image->detgeom->n_panels; i++ ) {
diffraction_panel(image, intensities, phases, flags, cell, m,
sym, k, ax, ay, az, bx, by, bz, cx, cy, cz,
lut_a, lut_b, lut_c, i, weight);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 916257e3..4d74589b 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -46,12 +46,12 @@
#include <cell.h>
#include <cell-utils.h>
#include <utils.h>
-#include <detector.h>
#include <peaks.h>
#include <symmetry.h>
#include <reflist.h>
#include <reflist-utils.h>
#include <stream.h>
+#include <datatemplate.h>
#include "diffraction.h"
#include "diffraction-gpu.h"
@@ -117,7 +117,8 @@ static void show_help(const char *s)
}
-static void record_panel(struct panel *p, float *dp, int do_poisson,
+static void record_panel(struct detgeom_panel *p, float *dp,
+ int do_poisson,
gsl_rng *rng, double ph_per_e, double background,
double lambda,
int *n_neg1, int *n_inf1, int *n_nan1,
@@ -144,16 +145,16 @@ static void record_panel(struct panel *p, float *dp, int do_poisson,
if ( isnan(intensity) ) (*n_nan1)++;
/* Area of one pixel */
- pix_area = pow(1.0/p->res, 2.0);
- Lsq = pow(p->clen, 2.0);
+ pix_area = pow(p->pixel_pitch, 2.0);
+ Lsq = pow(p->cnz*p->pixel_pitch, 2.0);
/* Calculate distance from crystal to pixel */
xs = fs*p->fsx + ss*p->ssx;
ys = ss*p->fsy + ss*p->ssy;
- rx = (xs + p->cnx) / p->res;
- ry = (ys + p->cny) / p->res;
+ rx = (xs + p->cnx) * p->pixel_pitch;
+ ry = (ys + p->cny) * p->pixel_pitch;
dsq = pow(rx, 2.0) + pow(ry, 2.0);
- twotheta = atan2(sqrt(dsq), p->clen);
+ twotheta = atan2(sqrt(dsq), p->cnz*p->pixel_pitch);
/* Area of pixel as seen from crystal (approximate) */
proj_area = pix_area * cos(twotheta);
@@ -216,11 +217,9 @@ void record_image(struct image *image, int do_poisson, double background,
"Total energy = %5.3f microJ\n",
nphotons, energy_density/1e7, total_energy*1e6);
- fill_in_adu(image);
+ for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
- for ( pn=0; pn<image->det->n_panels; pn++ ) {
-
- record_panel(&image->det->panels[pn], image->dp[pn],
+ record_panel(&image->detgeom->panels[pn], image->dp[pn],
do_poisson, rng, ph_per_e, background,
image->lambda,
&n_neg1, &n_inf1, &n_nan1,
@@ -364,28 +363,8 @@ static int random_ncells(double len, double min_m, double max_m)
}
-static void fixup_geom(struct detector *det)
-{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
- det->panels[i].clen += det->panels[i].coffset;
- }
-}
-
-
-static int geom_contains_references(struct detector *det)
-{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
- if ( det->panels[i].clen_from != NULL ) return 1;
- }
- return 0;
-}
-
-
-static void add_metadata(const char *filename, struct quaternion q,
+static void add_metadata(const char *filename,
+ struct quaternion q,
UnitCell *cell)
{
hid_t fh, gh, sh, dh;
@@ -488,9 +467,10 @@ static void add_metadata(const char *filename, struct quaternion q,
int main(int argc, char *argv[])
{
int c;
- struct image image;
+ struct image *image;
+ DataTemplate *dtempl;
struct gpu_context *gctx = NULL;
- struct image powder;
+ struct image *powder;
char *intfile = NULL;
double *intensities;
char *rval;
@@ -536,8 +516,6 @@ int main(int argc, char *argv[])
double photon_energy = 9000.0;
double sase_spike_width = 2e6; /* m^-1 */
double twocol_sep = 8e6; /* m^-1 */
- struct beam_params beam;
- int i;
/* Long options */
const struct option longopts[] = {
@@ -802,7 +780,7 @@ int main(int argc, char *argv[])
ERROR("You cannot use -r and --template together.\n");
return 1;
}
- st = open_stream_for_read(template_file);
+ st = stream_open_for_read(template_file);
if ( st == NULL ) {
ERROR("Failed to open stream.\n");
return 1;
@@ -836,26 +814,13 @@ int main(int argc, char *argv[])
ERROR("You need to specify a geometry file with --geometry\n");
return 1;
}
- image.beam = &beam;
- image.det = get_detector_geometry(geometry, image.beam);
- if ( image.det == NULL ) {
- ERROR("Failed to read detector geometry from '%s'\n", geometry);
+ dtempl = data_template_new_from_file(geometry);
+ if ( dtempl == NULL ) {
+ ERROR("Failed to read detector geometry from '%s'\n",
+ geometry);
return 1;
}
free(geometry);
- if ( (beam.photon_energy > 0.0) && (beam.photon_energy_from == NULL) ) {
- ERROR("WARNING: An explicit photon energy was found in the "
- "geometry file. It will be ignored!\n");
- ERROR("The value given on the command line "
- "(with --photon-energy) will be used instead.\n");
- }
- if ( geom_contains_references(image.det) ) {
- ERROR("Geometry file contains a reference to an HDF5 location"
- " for the camera length. Change it to a numerical value "
- " and try again.\n");
- return 1;
- }
- fixup_geom(image.det);
if ( spectrum_str == NULL ) {
STATUS("You didn't specify a spectrum type, so"
@@ -952,28 +917,6 @@ int main(int argc, char *argv[])
}
- image.lambda = ph_en_to_lambda(eV_to_J(photon_energy));
- image.bw = bandwidth;
-
- /* Initialise stuff */
- image.filename = NULL;
- image.features = NULL;
- image.bad = NULL;
- image.dp = malloc(image.det->n_panels*sizeof(float *));
- if ( image.dp == NULL ) {
- ERROR("Failed to allocate data\n");
- return 1;
- }
- for ( i=0; i<image.det->n_panels; i++ ) {
- image.dp[i] = calloc(image.det->panels[i].w *
- image.det->panels[i].h,
- sizeof(float));
- if ( image.dp[i] == NULL ) {
- ERROR("Failed to allocate data\n");
- return 1;
- }
- }
-
rng = gsl_rng_alloc(gsl_rng_mt19937);
if ( config_random ) {
FILE *fh;
@@ -987,30 +930,16 @@ int main(int argc, char *argv[])
fclose(fh);
}
- powder.det = image.det;
- powder.beam = NULL;
- powder.lambda = 0.0;
- powder.spectrum = NULL;
- powder.dp = malloc(image.det->n_panels*sizeof(float *));
- if ( powder.dp == NULL ) {
- ERROR("Failed to allocate powder data\n");
- return 1;
- }
- for ( i=0; i<image.det->n_panels; i++ ) {
- powder.dp[i] = calloc(image.det->panels[i].w *
- image.det->panels[i].h,
- sizeof(float));
- if ( powder.dp[i] == NULL ) {
- ERROR("Failed to allocate powder data\n");
- return 1;
- }
- }
+ image = image_create_for_simulation(dtempl);
+ powder = image_create_for_simulation(dtempl);
/* Splurge a few useful numbers */
STATUS("Simulation parameters:\n");
- STATUS(" Photon energy: %.2f eV (wavelength %.5f A)\n",
- photon_energy, image.lambda*1e10);
- STATUS(" Number of photons: %.0f (%.2f mJ)\n", nphotons,
+ STATUS(" Wavelength: %.5f A (photon energy %.2f eV)\n",
+ image->lambda*1e10,
+ ph_lambda_to_eV(image->lambda));
+ STATUS(" Number of photons: %.0f (%.2f mJ)\n",
+ nphotons,
eV_to_J(photon_energy)*nphotons*1e3);
STATUS(" Beam divergence: not simulated\n");
STATUS(" Beam radius: %.2f microns\n",
@@ -1021,22 +950,22 @@ int main(int argc, char *argv[])
case SPECTRUM_TOPHAT:
STATUS(" X-ray spectrum: top hat, "
- "width %.5f %%\n", image.bw*100.0);
+ "width %.5f %%\n", image->bw*100.0);
break;
case SPECTRUM_GAUSSIAN:
STATUS(" X-ray spectrum: Gaussian, "
- "bandwidth %.5f %%\n", image.bw*100.0);
+ "bandwidth %.5f %%\n", image->bw*100.0);
break;
case SPECTRUM_SASE:
STATUS(" X-ray spectrum: SASE, "
- "bandwidth %.5f %%\n", image.bw*100.0);
+ "bandwidth %.5f %%\n", image->bw*100.0);
break;
case SPECTRUM_TWOCOLOUR:
STATUS(" X-ray spectrum: two colour, "
- "separation %.5f %%\n", image.bw*100.0);
+ "separation %.5f %%\n", image->bw*100.0);
break;
case SPECTRUM_FROMFILE:
@@ -1066,12 +995,13 @@ int main(int argc, char *argv[])
int err = 0;
int pi;
- for ( pi=0; pi<image.det->n_panels; pi++ ) {
+ /* Reset image data to zero for new pattern */
+ for ( pi=0; pi<image->detgeom->n_panels; pi++ ) {
long j;
- long np = image.det->panels[pi].w
- * image.det->panels[pi].h;
+ long np = image->detgeom->panels[pi].w
+ * image->detgeom->panels[pi].h;
for ( j=0; j<np; j++ ) {
- image.dp[pi][j] = 0.0;
+ image->dp[pi][j] = 0.0;
}
}
@@ -1115,70 +1045,54 @@ int main(int argc, char *argv[])
} else {
- struct image image;
- int i;
+ struct image *templ_image;
Crystal *cr;
- image.det = NULL;
-
/* Get data from next chunk */
- if ( read_chunk(st, &image) ) break;
+ templ_image = stream_read_chunk(st, dtempl,
+ STREAM_CRYSTALS);
+ if ( templ_image == NULL ) break;
+ if ( templ_image->n_crystals == 0 ) continue;
- free(image.filename);
- image_feature_list_free(image.features);
+ cr = templ_image->crystals[0];
+ cell = cell_new_from_cell(crystal_get_cell(cr));
- if ( image.n_crystals == 0 ) continue;
-
- cr = image.crystals[0];
- cell = crystal_get_cell(cr);
-
- if ( image.n_crystals > 1 ) {
+ if ( templ_image->n_crystals > 1 ) {
ERROR("Using the first crystal only.\n");
}
- for ( i=1; i<image.n_crystals; i++ ) {
-
- Crystal *cr = image.crystals[i];
- cell = crystal_get_cell(cr);
-
- reflist_free(crystal_get_reflections(cr));
- cell_free(crystal_get_cell(cr));
- crystal_free(cr);
-
- }
-
- free(image.crystals);
+ image_free(templ_image);
}
switch ( spectrum_type ) {
case SPECTRUM_TOPHAT :
- image.spectrum = spectrum_generate_tophat(image.lambda,
- image.bw);
+ image->spectrum = spectrum_generate_tophat(image->lambda,
+ image->bw);
break;
case SPECTRUM_GAUSSIAN :
- image.spectrum = spectrum_generate_gaussian(image.lambda,
- image.bw);
+ image->spectrum = spectrum_generate_gaussian(image->lambda,
+ image->bw);
break;
case SPECTRUM_SASE :
- image.spectrum = spectrum_generate_sase(image.lambda,
- image.bw,
- sase_spike_width,
- rng);
+ image->spectrum = spectrum_generate_sase(image->lambda,
+ image->bw,
+ sase_spike_width,
+ rng);
break;
- case SPECTRUM_TWOCOLOUR :
- image.spectrum = spectrum_generate_twocolour(image.lambda,
- image.bw,
- twocol_sep);
+ case SPECTRUM_TWOCOLOUR :
+ image->spectrum = spectrum_generate_twocolour(image->lambda,
+ image->bw,
+ twocol_sep);
break;
- case SPECTRUM_FROMFILE :
- image.spectrum = spectrum_load(spectrum_fn);
+ case SPECTRUM_FROMFILE :
+ image->spectrum = spectrum_load(spectrum_fn);
break;
}
@@ -1195,12 +1109,12 @@ int main(int argc, char *argv[])
intensities, flags, sym_str,
gpu_dev);
}
- err = get_diffraction_gpu(gctx, &image, na, nb, nc,
+ err = get_diffraction_gpu(gctx, image, na, nb, nc,
cell, no_fringes, flat,
nsamples);
} else {
- get_diffraction(&image, na, nb, nc, intensities, phases,
+ get_diffraction(image, na, nb, nc, intensities, phases,
flags, cell, grad, sym, no_fringes, flat,
nsamples);
}
@@ -1210,28 +1124,28 @@ int main(int argc, char *argv[])
goto skip;
}
- record_image(&image, !config_nonoise, background, rng,
+ record_image(image, !config_nonoise, background, rng,
beam_radius, nphotons);
if ( powder_fn != NULL ) {
int pn;
- for ( pn=0; pn<image.det->n_panels; pn++ ) {
+ for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
size_t w, i;
- w = image.det->panels[pn].w
- * image.det->panels[pn].h;
+ w = image->detgeom->panels[pn].w
+ * image->detgeom->panels[pn].h;
for ( i=0; i<w; i++ ) {
- powder.dp[pn][i] += (double)image.dp[pn][i];
+ powder->dp[pn][i] += (double)image->dp[pn][i];
}
}
if ( !(ndone % 10) ) {
- hdf5_write_image(powder_fn, &powder, NULL);
+ image_write(powder, dtempl, powder_fn);
}
}
@@ -1249,7 +1163,9 @@ int main(int argc, char *argv[])
number++;
/* Write the output file */
- hdf5_write_image(filename, &image, NULL);
+ image_write(image, dtempl, filename);
+
+ /* Add some pattern_sim-specific metadata */
add_metadata(filename, orientation, cell);
}
@@ -1264,21 +1180,16 @@ skip:
} while ( !done );
if ( powder_fn != NULL ) {
- hdf5_write_image(powder_fn, &powder, NULL);
+ image_write(powder, dtempl, powder_fn);
}
if ( gctx != NULL ) {
cleanup_gpu(gctx);
}
- for ( i=0; i<image.det->n_panels; i++ ) {
- free(image.dp[i]);
- free(powder.dp[i]);
- }
- free(image.dp);
- free(powder.dp);
+ image_free(image);
+ image_free(powder);
free(intfile);
- free(image.det);
cell_free(input_cell);
free(intensities);
free(outfile);