aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--INSTALL.md13
-rw-r--r--config.h.in2
-rw-r--r--data/diffraction.cl213
-rwxr-xr-xdata/gen-resources20
-rw-r--r--doc/man/crystfel.714
-rw-r--r--doc/man/crystfel_geometry.58
-rw-r--r--doc/man/get_hkl.13
-rw-r--r--doc/man/partial_sim.1158
-rw-r--r--doc/man/partialator.13
-rw-r--r--doc/man/pattern_sim.1250
-rw-r--r--meson.build32
-rw-r--r--src/cl-utils.c354
-rw-r--r--src/cl-utils.h55
-rw-r--r--src/diffraction-gpu.c596
-rw-r--r--src/diffraction-gpu.h77
-rw-r--r--src/diffraction.c502
-rw-r--r--src/diffraction.cl.h513
-rw-r--r--src/diffraction.h66
-rw-r--r--src/list_tmp.h104
-rw-r--r--src/partial_sim.c1034
-rw-r--r--src/pattern_sim.c1183
-rw-r--r--src/pattern_sim.h67
-rw-r--r--tests/gpu_sim_check.c198
-rw-r--r--tests/gpu_sim_check.geom28
-rw-r--r--tests/meson.build10
25 files changed, 9 insertions, 5494 deletions
diff --git a/INSTALL.md b/INSTALL.md
index 87346ab3..d167328d 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -97,7 +97,6 @@ roughly in order of importance:
* [msgpack-c](https://github.com/msgpack/msgpack-c) (for streaming data in MsgPack format)
* [Seedee](https://gitlab.desy.de/fs-sc/seedee) (for streaming data serialised with Seedee)
* [cJSON](https://github.com/DaveGamble/cJSON/) \[\*\] (extra dependency if Seedee is found)
-* [OpenCL](https://www.khronos.org/opencl/) plus appropriate GPU drivers (for GPU-accelerated pattern simulation)
Most of the dependencies mentioned above should be available from your Linux
distribution's package manager, or from [Homebrew](https://brew.sh/) on Mac OS.
@@ -167,18 +166,6 @@ as follows, replacing the `meson build` step:
meson build -Dhdf5=disabled
```
-Using OpenCL for GPU acceleration
----------------------------------
-
-OpenCL is disabled by default when building with Meson, to keep things simple
-for the majority of people who don't need it. To enable it, add
-`-Dopencl=enabled` to the `meson` command line:
-```
-meson build -Dopencl=enabled
-```
-When building with CMake, OpenCL will always be used if the libraries are
-found.
-
Fedora 22 or later
------------------
diff --git a/config.h.in b/config.h.in
index 2cc2a916..c3f6d2c3 100644
--- a/config.h.in
+++ b/config.h.in
@@ -2,8 +2,6 @@
#mesondefine HAVE_GTK
#mesondefine HAVE_CAIRO
-#mesondefine HAVE_OPENCL
-#mesondefine HAVE_CL_CL_H
#mesondefine HAVE_CLOCK_GETTIME
#mesondefine HAVE_ZMQ
#mesondefine HAVE_HDF5
diff --git a/data/diffraction.cl b/data/diffraction.cl
deleted file mode 100644
index 6e3b0f69..00000000
--- a/data/diffraction.cl
+++ /dev/null
@@ -1,213 +0,0 @@
-/*
- * diffraction.cl
- *
- * GPU calculation kernel for truncated lattice diffraction
- *
- * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2009-2014 Thomas White <taw@physics.org>
- * 2013 Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-/* Maxmimum index to hold values up to (can be increased if necessary)
- * WARNING: Altering this value constitutes an ABI change, and means you must
- * update src/pattern_sim.h then recompile and reinstall everything. */
-#define INDMAX 130
-#define IDIM (INDMAX*2 +1)
-
-#ifndef M_PI
-#define M_PI ((float)(3.14159265))
-#endif
-
-const sampler_t sampler_a = CLK_NORMALIZED_COORDS_TRUE
- | CLK_ADDRESS_REPEAT
- | CLK_FILTER_LINEAR;
-const sampler_t sampler_b = CLK_NORMALIZED_COORDS_TRUE
- | CLK_ADDRESS_REPEAT
- | CLK_FILTER_LINEAR;
-const sampler_t sampler_c = CLK_NORMALIZED_COORDS_TRUE
- | CLK_ADDRESS_REPEAT
- | CLK_FILTER_LINEAR;
-
-
-float4 get_q(float fs, float ss, float res, float clen, float k,
- float corner_x, float corner_y,
- float fsx, float fsy, float fsz, float ssx, float ssy, float ssz)
-{
- float rx, ry, rz;
- float az, tt;
- float4 q;
- float xs, ys;
- float kx, ky, kz;
- float ctt;
-
- /* Calculate 3D position of given position, in m */
- rx = (corner_x + fs*fsx + ss*ssx) / res;
- ry = (corner_y + fs*fsy + ss*ssy) / res;
- rz = clen + (fs*fsz + ss*ssz)/res;
-
- ctt = rz / sqrt(rx*rx + ry*ry + rz*rz); /* cos(2theta) */
- tt = acos(ctt);
- az = atan2(ry, rx);
-
- kx = k*native_sin(tt)*native_cos(az);
- ky = k*native_sin(tt)*native_sin(az);
- kz = k*(ctt - 1.0);
-
- q = (float4)(kx, ky, kz, 0.0);
-
- return q;
-}
-
-
-float lattice_factor(float16 cell, float4 q,
- read_only image2d_t func_a,
- read_only image2d_t func_b,
- read_only image2d_t func_c)
-{
- float f1, f2, f3, v;
- float4 Udotq;
-
- Udotq.x = cell.s0*q.x + cell.s1*q.y + cell.s2*q.z;
- Udotq.y = cell.s3*q.x + cell.s4*q.y + cell.s5*q.z;
- Udotq.z = cell.s6*q.x + cell.s7*q.y + cell.s8*q.z;
-
- /* Look up values from precalculated sinc() table */
- f1 = read_imagef(func_a, sampler_a, (float2)(Udotq.x, 0.0)).s0;
- f2 = read_imagef(func_b, sampler_b, (float2)(Udotq.y, 0.0)).s0;
- f3 = read_imagef(func_c, sampler_c, (float2)(Udotq.z, 0.0)).s0;
-
- return f1 * f2 * f3;
-}
-
-
-float lookup_intensity(global float *intensities,
- signed int h, signed int k, signed int l)
-{
- int idx;
-
- /* Out of range? */
- if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
- return 0.0;
- }
-
- h = (h>=0) ? h : h+IDIM;
- k = (k>=0) ? k : k+IDIM;
- l = (l>=0) ? l : l+IDIM;
-
- idx = h + (IDIM*k) + (IDIM*IDIM*l);
-
- return intensities[idx];
-}
-
-
-float lookup_flagged_intensity(global float *intensities, global float *flags,
- signed int h, signed int k, signed int l)
-{
- return lookup_intensity(intensities, h, k, l)
- * lookup_intensity(flags, h, k, l);
-}
-
-
-float molecule_factor(global float *intensities, global float *flags,
- float16 cell, float4 q)
-{
- float hf, kf, lf;
- int h, k, l;
- float val = 0.0;
-
- #ifdef FLAT_INTENSITIES
- return 100.0;
- #else
-
- hf = cell.s0*q.x + cell.s1*q.y + cell.s2*q.z; /* h */
- kf = cell.s3*q.x + cell.s4*q.y + cell.s5*q.z; /* k */
- lf = cell.s6*q.x + cell.s7*q.y + cell.s8*q.z; /* l */
-
- h = round(hf);
- k = round(kf);
- l = round(lf);
-
- /* Symmetry stuff goes here */
- INSERT_HERE
-
- return val;
- #endif /* FLAT_INTENSITIIES */
-}
-
-
-kernel void diffraction(global float *diff, float k, float weight,
- int w, float corner_x, float corner_y,
- float fsx, float fsy, float fsz,
- float ssx, float ssy, float ssz,
- float res, float clen, float16 cell,
- global float *intensities, global float *flags,
- read_only image2d_t func_a,
- read_only image2d_t func_b,
- read_only image2d_t func_c,
- local float *tmp)
-{
- float fs, ss;
- float f_lattice, I_lattice;
- float I_molecule;
- float4 q;
- const int ls0 = get_local_size(0);
- const int ls1 = get_local_size(1);
- const int li0 = get_local_id(0);
- const int li1 = get_local_id(1);
- const int ls = ls0 * ls1;
-
- /* Calculate fractional coordinates in fs/ss */
- fs = convert_float(get_global_id(0)) / convert_float(ls0);
- ss = convert_float(get_global_id(1)) / convert_float(ls1);
-
- /* Get the scattering vector */
- q = get_q(fs, ss, res, clen, k,
- corner_x, corner_y, fsx, fsy, fsz, ssx, ssy, ssz);
-
- /* Calculate the diffraction */
- f_lattice = lattice_factor(cell, q, func_a, func_b, func_c);
- I_molecule = molecule_factor(intensities, flags, cell, q);
- I_lattice = pow(f_lattice, 2.0f);
-
- tmp[li0 + ls0*li1] = I_molecule * I_lattice;
-
- barrier(CLK_LOCAL_MEM_FENCE);
-
- /* First thread in group sums the samples */
- if ( li0 + li1 == 0 ) {
-
- int i;
- float sum = 0.0;
- float val;
- int idx;
-
- idx = convert_int_rtz(fs) + w*convert_int_rtz(ss);
-
- for ( i=0; i<ls; i++ ) sum += tmp[i];
-
- val = weight * sum / convert_float(ls);
- diff[idx] = val;
-
- }
-
-}
diff --git a/data/gen-resources b/data/gen-resources
deleted file mode 100755
index 273f7db7..00000000
--- a/data/gen-resources
+++ /dev/null
@@ -1,20 +0,0 @@
-#!/bin/sh
-
-function gen_resource()
-{
- INFILE=$1
- OUTFILE=$2
-
- echo "/*" > $OUTFILE
- echo " * This file was generated from $INFILE" >> $OUTFILE
- echo " * using the following command:" >> $OUTFILE
- echo " * xxd -i $INFILE $OUTFILE" >> $OUTFILE
- echo " *" >> $OUTFILE
- echo " * If you have 'xxd' installed, you can run the script" >> $OUTFILE
- echo " * data/gen-resources to re-create this file." >> $OUTFILE
- echo " */" >> $OUTFILE
- echo >> $OUTFILE
- xxd -i $INFILE >> $OUTFILE
-}
-
-gen_resource data/diffraction.cl src/diffraction.cl.h
diff --git a/doc/man/crystfel.7 b/doc/man/crystfel.7
index e04b7644..3ff25c16 100644
--- a/doc/man/crystfel.7
+++ b/doc/man/crystfel.7
@@ -1,7 +1,7 @@
.\"
.\" CrystFEL main man page
.\"
-.\" Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
+.\" Copyright © 2012-2023 Deutsches Elektronen-Synchrotron DESY,
.\" a research centre of the Helmholtz Association.
.\"
.\" Part of CrystFEL - crystallography with a FEL
@@ -24,15 +24,12 @@ Many patterns, numbering tens of thousands or more, are required, so high throug
The crystal orientations in each pattern are random and uncorrelated, which leads to special considerations during scaling and merging.
.RE
-CrystFEL includes programs for simulating and processing patterns subject to the
-above characteristics. The programs are:
+CrystFEL includes programs processing patterns subject to the above
+characteristics. The programs are:
.IP \fBindexamajig\fR
Batch indexing, integration and data reduction program, which produces a "stream" containing the indexing and integration results for each diffraction pattern.
-.IP \fBpattern_sim\fR
-A diffraction pattern simulation tool.
-
.IP \fBprocess_hkl\fR
A tool merging intensities from many patterns into a single reflection list, via the Monte Carlo method.
@@ -54,9 +51,6 @@ A tool for examining the distributions of unit cell parameters.
.IP "\fBcompare_hkl\fR and \fBcheck_hkl\fR"
Tools for calculating figures of merit, such as completeness and R-factors.
-.IP \fBpartial_sim\fB
-A tool for calculating partial reflection intensities, perhaps for testing the convergence of Monte Carlo merging.
-
.IP \fBrender_hkl\fR
A tool for rendering slices of reciprocal space in two dimensions.
@@ -160,8 +154,6 @@ You should have received a copy of the GNU General Public License along with Cry
.BR indexamajig (1),
.BR process_hkl (1),
.BR partialator (1),
-.BR pattern_sim (1),
-.BR partial_sim (1),
.BR compare_hkl (1),
.BR check_hkl (1),
.BR render_hkl (1),
diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5
index 6056c463..998aa213 100644
--- a/doc/man/crystfel_geometry.5
+++ b/doc/man/crystfel_geometry.5
@@ -14,8 +14,9 @@ CRYSTFEL DETECTOR GEOMETRY DESCRIPTION FILES
.SH CRYSTFEL DETECTOR GEOMETRY FILES
The detector geometry is taken from a text file rather than hardcoded into the
program. Programs which care about the geometry (particularly
-\fBindexamajig\fR and \fBpattern_sim\fR) take an argument
-\fB--geometry=\fR\fIfilename\fR, where \fIfilename\fR contains the geometry.
+\fBindexamajig\fR) take an argument \fB--geometry=\fR\fIfilename\fR, where
+\fIfilename\fR contains the geometry.
+
.PP
A flexible (and pedantic) representation of the detector has been developed to
avoid all possible sources of ambiguity. CrystFEL's representation of a
@@ -175,7 +176,7 @@ The range of pixels in the data block specified by the 'data' property that corr
.IP \fBadu_per_eV\fR
.IP \fBadu_per_photon\fR
.PD
-The number of detector intensity units (ADU) which will arise from either one electron-Volt of photon energy, or one photon. This is used to estimate Poisson errors. Note that setting different values for this parameter for different panels does \fBnot\fR result in the intensities being scaled accordingly when integrating data, but it does affect the intensities calculated by \fBpattern_sim\fR. You should only specify one out of \fBadu_per_eV\fR and \fBadu_per_photon\fR.
+The number of detector intensity units (ADU) which will arise from either one electron-Volt of photon energy, or one photon. This is used to estimate Poisson errors. Note that setting different values for this parameter for different panels does \fBnot\fR result in the intensities being scaled accordingly when integrating data. You should only specify one out of \fBadu_per_eV\fR and \fBadu_per_photon\fR.
.PD 0
.IP \fBres\fR
@@ -339,5 +340,4 @@ You should have received a copy of the GNU General Public License along with Cry
.SH SEE ALSO
.BR crystfel (7),
-.BR pattern_sim (1),
.BR indexamajig (1)
diff --git a/doc/man/get_hkl.1 b/doc/man/get_hkl.1
index b15a33df..bdc7df35 100644
--- a/doc/man/get_hkl.1
+++ b/doc/man/get_hkl.1
@@ -126,5 +126,4 @@ You should have received a copy of the GNU General Public License along with Cry
.SH SEE ALSO
.BR crystfel (7),
.BR check_hkl (1),
-.BR compare_hkl (1),
-.BR pattern_sim (1)
+.BR compare_hkl (1)
diff --git a/doc/man/partial_sim.1 b/doc/man/partial_sim.1
deleted file mode 100644
index 102e3cdf..00000000
--- a/doc/man/partial_sim.1
+++ /dev/null
@@ -1,158 +0,0 @@
-.\"
-.\" partial_sim man page
-.\"
-.\" Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
-.\" a research centre of the Helmholtz Association.
-.\"
-.\" Part of CrystFEL - crystallography with a FEL
-.\"
-
-.TH PARTIAL_SIM 1
-.SH NAME
-partial_sim \- calculate partial reflections
-.SH SYNOPSIS
-.PP
-.BR partial_sim
-\fB-o\fR \fIsimulated.stream\fR
-\fB-g\fR \fIgeometry.geom\fR
-\fB-p\fR \fIunitcell.pdb\fR
-[\fIoptions\fR] \fB...\fR
-
-.BR partial_sim
-\fB--help\fR
-
-.SH DESCRIPTION
-partial_sim calculates the intensities of idealised partial reflections from crystals in random orientations, which is useful for testing the convergence of Monte Carlo integration or scaling/post-refinement techniques.
-
-.P
-You need to provide a CrystFEL geometry file (with \fB--geometry=\fR\fImy.geom\fR or \fB-g\fR \fImy.geom\fR), a file containing the unit cell to use for the simulation (with \fB--pdb=\fR\fImy.pdb\fR or \fB-p\fR \fImy.pdb\fR), and an output filename with \fB--output=\fR\fImy.stream\fR or \fB-o\fR \fImy.stream\fR.
-
-For each randomly generated orientation, partial_sim calculates which reflections would appear on the detector with the specified beam parameters. It calculates the partiality for each reflection and multiplies it by the fully integrated intensity to produce a partial intensity. The fully integrated intensities can be taken from a file you provide (see below), otherwise they will be randomly generated (by taking the absolute value of a Gaussian random number, mean zero and standard deviation 1000). All the partial intensities for the orientation are multiplied by an overall scaling factor, which is randomly generated with a Gaussian distribution with mean 1 and standard deviation 0.3. The partial intensities are written to the output stream, and the process repeated for as many different orientations as you ask for (see below, default: 2).
-
-Set the incident radiation wavelength and bandwidth in the geometry file. Set the reciprocal space profile radius on the command line with \fI--profile-radius\fR. Beam divergence is not taken into account by the \fIxsphere\fR partiality model used by partial_sim. See \fBman crystfel_geometry\fR for information about CrystFEL geometry description files.
-
-.SH OPTIONS
-.PD 0
-.B
-.IP "-i \fIfile.hkl\fR"
-.B
-.IP --input=\fIfile.hkl\fR
-.PD
-Take the fully integrated reflection intensities from \fIfile.hkl\fR, instead of generating them randomly.
-
-.B
-.IP "\fB-n\fR \fIn\fR"
-Specify the number of different orientations to simulate. Default: 2.
-
-.PD 0
-.IP "\fB-p\fR \fIunitcell.cell\fR"
-.IP "\fB-p\fR \fIunitcell.pdb\fR"
-.IP \fB--pdb=\fR\fIunitcell.pdb\fR
-.PD
-Specify the name of the file containing unit cell information, in PDB or CrystFEL format.
-
-.PD 0
-.B
-.IP "-r \fIrandom.hkl\fR"
-.B
-.IP --save-random=\fIrandom.hkl\fR
-.PD
-If you did not provide your own fully integrated reflection intensities, they will be generated randomly for you. Use this option to save the random intensities for future comparisons.
-
-.PD 0
-.B
-.IP "\fB-y\fR \fIpointgroup\fR"
-.B
-.IP "\fB--symmetry=\fR\fIpointgroup\fR"
-.PD
-When combined with with \fB-i\fR, specifies the symmetry of the input reflection list. Otherwise, specifies the symmetry of the randomly generated intensities.
-
-.PD 0
-.B
-.IP "\fB-c\fR \fIval\fR"
-.B
-.IP "\fB--cnoise=\fR\fIval\fR"
-.PD
-Add random values with a flat distribution to the components of the reciprocal lattice vectors written to the stream, simulating indexing errors. The maximum value that will be added is +/- \fIval\fR percent.
-
-.PD 0
-.B
-.IP "\fB--pgraph=\fR\fIfilename\fR"
-.PD
-Save a table of values to \fIfilename\fR containing, in resolution shells, the following columns: resolution shell centre in nm^-1, number of reflections in shell, mean partiality, maximum partiality.
-
-.PD 0
-.B
-.IP "\fB--osf-stddev=\fR\fIval\fR"
-.PD
-Set the standard deviation of the distribution of overall scaling factors to \fIval\fR. The distribution will be cut at zero, i.e. negative or zero scaling factors are not allowed. The distribution will be Gaussian centered on 1. The default is \fB--osf-stddev=2.0\fR.
-
-.PD 0
-.B
-.IP "\fB--full-stddev=\fR\fIval\fR"
-.PD
-Set the standard deviation of the distribution of randomly generated full intensities to \fIval\fR. The distribution will be Gaussian, centered on zero, and the absolute value will be taken (i.e. there will be no negative values). The default is \fB--full-stddev=1000.0\fR. This option has no effect if you also use \fB-i\fR or \fB--input\fR.
-
-.PD 0
-.B
-.IP "\fB--noise-stddev=\fR\fIval\fR"
-.PD
-Set the standard deviation of the noise added to the partial intensities to \fIval\fR. The noise will be Gaussian, and the same for all reflections. The default is \fB--noise-stddev=20.0\fR.
-
-.PD 0
-.B
-.IP "\fB-j\fR \fIn\fR"
-.PD
-Use \fIn\fR threads for simulation. Default: 1.
-
-.PD 0
-.B
-.IP "\fB--images=\fR\fIprefix\fR"
-.PD
-For each chunk in the output stream, write a 'sketch' image in HDF5 format to \fIprefix\fR\fB/sim-\fR\fINNN\fR\fB.h5\fR, where \fINNN\fR is the sequence number of the chunk in the output stream. This option is incompatible with \fB-j\fR. The intensities in the peaks in the sketches will be equal to the partial intensities in the stream, including noise and overall scaling factors. The images will also contain a random Poisson-distributed background according to \fB--background\fR.
-
-.PD 0
-.B
-.IP "\fB--background=\fIval\fR"
-.PD
-Add a Poisson-distributed background with \fIval\fR photons to the sketches (see \fB--images\fR). The default is \fB--background=3000\fR.\fR.
-
-.PD 0
-.B
-.IP "\fB--profile-radius=\fIval\fR"
-.PD
-Set the radius of the scattering density surrounding each reciprocal lattice point, in m^-1. The default is \fB--profile-radius=0.001e9\fR m^-1.
-
-.PD 0
-.IP \fB--really-random\fR
-.PD
-Seed the random number generator using the kernel random number generator (/dev/urandom). This means that truly random (although not "cryptographically random") numbers will be used for the orientation and crystal size, instead of the same sequence being used for each new run.
-
-.IP "\fB--template-stream=\fImy.stream\fR"
-.PD
-Get the crystal cell parameters, orientations and the reflections to calculate from \fImy.stream\fR. This allows you to re-calculate partial intensities using new beam parameters. There must only be one crystal per chunk in the template. If more than one thread is used (see \fB-j\fR), note that there is no guarantee that the order of crystals in the output stream will match that of \fImy.stream\fR.
-
-
-.SH AUTHOR
-This page was written by Thomas White.
-
-.SH REPORTING BUGS
-Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
-
-.SH COPYRIGHT AND DISCLAIMER
-Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
-.P
-partial_sim, and this manaul, are part of CrystFEL.
-.P
-CrystFEL is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-.P
-CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-.P
-You should have received a copy of the GNU General Public License along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
-
-.SH SEE ALSO
-.BR crystfel (7),
-.BR process_hkl (1),
-.BR partialator (1),
-.BR pattern_sim (1),
-.BR crystfel_geometry (5).
diff --git a/doc/man/partialator.1 b/doc/man/partialator.1
index 6fa67593..d900476d 100644
--- a/doc/man/partialator.1
+++ b/doc/man/partialator.1
@@ -281,5 +281,4 @@ You should have received a copy of the GNU General Public License along with Cry
.SH SEE ALSO
.BR crystfel (7),
.BR indexamajig (1),
-.BR process_hkl (1),
-.BR partial_sim (1)
+.BR process_hkl (1)
diff --git a/doc/man/pattern_sim.1 b/doc/man/pattern_sim.1
deleted file mode 100644
index 36287162..00000000
--- a/doc/man/pattern_sim.1
+++ /dev/null
@@ -1,250 +0,0 @@
-.\"
-.\" pattern_sim man page
-.\"
-.\" Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
-.\" a research centre of the Helmholtz Association.
-.\"
-.\" Part of CrystFEL - crystallography with a FEL
-.\"
-
-.TH PATTERN_SIM 1
-.SH NAME
-pattern_sim \- Simulation of nanocrystal diffraction patterns
-.SH SYNOPSIS
-.PP
-.B pattern_sim
-\fB-g\fR \fIdetector.geom\fR \fB-p\fR \fImy.pdb\fR
-[\fBoptions\fR] \fB...\fR
-.PP
-.B pattern_sim
-\fB--help\fR
-
-.SH DESCRIPTION
-
-pattern_sim simulates diffraction patterns from small crystals probed with femtosecond pulses of X-rays from a free electron laser. Typical use might be of the form:
-
-pattern_sim -g mydetector.geom -p my.pdb -r -i myintensities.hkl
-
-The unit cell geometry will be taken from the unit cell file you provide, and the intensities of the reflections will be interpolated from the reflection list file you provide. The reflection list format is the same as that output by process_hkl and handled by get_hkl. You also need a geometry description file (-g). See `man crystfel_geometry' for details of how to create a geometry file. Examples of both files can be found in the installation directory, which is normally /usr/local/share/doc/crystfel.
-
-Set the incident radiation wavelength and bandwidth in the geometry file. Beam divergence is not taken into account. See \fBman crystfel_geometry\fR for information about CrystFEL geometry description files.
-
-The result will be written to an HDF5 file in the current directory with the name `sim.h5'.
-
-.SH OPTIONS
-
-.PD 0
-.IP "\fB-p\fR \fIunitcell.cell\fR"
-.IP "\fB-p\fR \fIunitcell.pdb\fR"
-.IP \fB--pdb=\fR\fIunitcell.pdb\fR
-.PD
-Specify the name of the file containing unit cell information, in PDB or CrystFEL format.
-
-.PD 0
-.IP \fB--gpu\fR
-.PD
-Use the GPU to speed up the calculation. Requires that OpenCL libraries and drivers are available, and that CrystFEL was compiled with OpenCL enabled.
-
-.PD 0
-.IP \fB--gpu-dev=\fRIn\fR
-.PD
-Use GPU device number \fIn\fR. If you omit this option, the list of GPU devices will be shown when you run pattern_sim.
-
-.PD 0
-.IP "\fB-g\fR \fIfilename\fR"
-.IP \fB--geometry=\fR\fIfilename\fR
-.PD
-Read the detector geometry description from \fIfilename\fR. See \fBman crystfel_geometry\fR for more information.
-
-.PD 0
-.IP "\fB-n\fR \fn\fR"
-.IP \fB--number=\fR\fIn\fR
-.PD
-Simulate \fIn\fR patterns. Default: \fB-n 1\fR.
-
-.PD 0
-.IP \fB--no-images\fR
-.PD
-Do not save any HDF5 files apart from the powder pattern (if requested).
-
-.PD 0
-.IP "\fB-o\fR \fIfilename\fR"
-.IP \fB--output=\fR\fIfilename\fR
-.PD
-Write the pattern to \fIfilename\fR. The default is \fB--output=sim.5\fR. If more than one pattern is to be simulated (see \fB--number\fR), the filename will be postfixed with a hyphen, the image number and then '.h5'.
-
-.PD 0
-.IP \fB-r\fR
-.IP \fB--random-orientation\fR
-.PD
-Make up a random orientation for each pattern simulated.
-
-.PD 0
-.IP \fB--powder=\fR\fIfilename\fR
-.PD
-Write the sum of all patterns to \fIfilename\fR.
-
-.PD 0
-.IP "\fB-i\fR \fIfilename\fR"
-.IP \fB--intensities=\fR\fIfilename\fR
-.PD
-Get the intensities and phases at the reciprocal lattice points from \fIfilename\fR.
-
-.PD 0
-.IP "\fB-y\fR \fIpointgroup\fR"
-.IP \fB--symmetry=\fR\fIpointgroup\fR
-.PD
-Use \fIpointgroup\fR as the symmetry of the intensity list (see \fB--intensities\fR).
-
-.PD 0
-.IP "\fB-t\fR \fImethod\fR"
-.IP \fB--gradients=\fR\fImethod\fR
-.PD
-Use \fImethod\fR as way of calculating the molecular transform between reciprocal lattice points. See the section \fBGRADIENT METHODS\fR below.
-
-.PD 0
-.IP \fB--really-random\fR
-.PD
-Seed the random number generator using the kernel random number generator (/dev/urandom). This means that truly random numbers for the orientation and crystal size, instead of the same sequence being used for each new run.
-
-.PD 0
-.IP \fB--min-size=\fR\fImin\fR
-.IP \fB--min-size=\fR\fImax\fR
-.PD
-Generate random crystal sizes between \fImin\fR and \fImax\fR nanometres. These options must be used together.
-
-.PD 0
-.IP \fB--no-noise\fR
-.PD
-Do not calculate Poisson noise.
-
-.PD 0
-.IP "\fB-s\fR \fIn\fR"
-.IP \fB--sample-spectrum=\fR\fIn\fR
-.PD
-Include \fIn\fR samples from the spectrum in the calculation.
-
-.PD 0
-.IP "\fB-x\fR \fItype\fR"
-.IP \fB--spectrum=\fR\fItype\fR
-.PD
-Use \fItype\fR of spectrum. \fItype\fR can be one of \fBtophat\fR, \fBsase\fR or \fBtwocolour\fR. See the section \fBSPECTRUM TYPES\fR below.
-
-.PD 0
-.IP \fB--background=\fR\fIn\fR
-.PD
-Add \fIn\fR photons of Poisson-distributed background uniformly over the detector surface.
-
-.PD 0
-.IP \fB--no-fringes\fR
-.PD
-Suppress the subsidiary maxima of the shape transforms by setting I_latt(q) to zero beyond the first minimum of the function.
-
-.PD 0
-.B
-.IP "\fB--nphotons=\fIval\fR"
-.PD
-Set the number of photons per X-ray pulse. The default is \fB--nphotons=1e12\fR. A physically reasonable value is such that the pulse energy (number of photons multiplied by photon energy) is about 1 mJ.
-
-.IP "\fB--beam-radius=\fIval\fR"
-.PD
-Set the radius of the X-ray beam, in metres. The default is \fB--beam-radius=1e-6\fR, i.e. a beam of 2 microns' diameter.
-
-.SH REFLECTION LISTS
-
-You'll need to create a file containing the intensities of the reflections. The normal way to do this is to use CCP4 via the "gen-sfs" script in CrystFEL's script folder. Run it like this:
-
-$ \fBgen-sfs\fR \fImymodel.pdb\fR \fB"\fR\fIP6\fR\fB"\fR \fI3\fR
-
-You need to give the PDB model, the symmetry of the output reflections (use the lowest symmetry space group with the right point group), and optionally the maximum resolution in Angstroms. If you don't specify the resolution, it'll use 3 Angstroms.
-
-The reflections will be output as \fImymodel.pdb\fR.hkl ready for input to pattern_sim. You'll need to give the Laue class of the symmetry you gave to gen-sfs, "6/m" in this case, to pattern_sim with \fB-y\fR. By default, \fBgen-sfs\fR calculates the values for CuKa radiation (8.3 keV, 1.5 A). It will not calculate the anomalous contribution to scattering, i.e. the differences in intensities between Bijoet pairs. Both of these are the default behaviour for "sfall" in CCP4, so read the manual for that for further details. If you need something different, get the "ano_sfall.com" script from James Holton and use the
-\fBgen-sfs-ano\fR script instead of \fBgen-sfs\fR.
-
-.SH CALCULATION DETAILS
-
-The lattice transform from the specified number of unit cells is calculated
-using the closed-form solution for a truncated lattice faceted on the
-(001), (010) and (100) planes:
-.IP
-I_latt(q) = sin^2(pi*na*g.a)/sin^2(pi*g.a)
- * sin^2(pi*nb*g.b)/sin^2(pi*g.b)
- * sin^2(pi*nc*g.c)/sin^2(pi*g.c)
-.IP
-na = number of unit cells in 'a' direction (likewise nb, nc)
-.br
- g = reciprocal vector (1/d convention, not 2pi/d)
-.PP
-This is multiplied by a model of the underlying molecular transform, I_mol(g).
-This can be approximated to varying levels of accuracy by the methods given by
-\fB--gradients\fR.
-.PP
-Expected intensities at the CCD are then calculated using:
-.IP
-I(g) = I0 * r^2 * I_latt(g) * I_mol(g) * S
-.IP
-I0 = number of photons per unit area in the incident beam
- r = Thomson radius
- S = solid angle of corresponding pixel
-.PP
-Polarisation is not currently included in pattern_sim, although it is included
-in the analysis of Bragg peaks done by \fBindexamajig\fR.
-.PP
-Poisson counts are generated from the expected intensities using Knuth's
-algorithm. When the intensity is sufficiently high that Knuth's algorithm
-would result in machine precision problems, a normal distribution with
-standard deviation sqrt(I) is used instead.
-
-.SH GRADIENT METHODS
-
-The available options for \fB--gradients\fR as as follows:
-
-.IP \fBmosaic\fR
-.PD
-Take the intensity of the nearest Bragg position. This is the fastest method and the only one supported on the GPU, but the least accurate.
-
-.IP \fBinterpolate\fR
-.PD
-Interpolate trilinearly between six adjacent Bragg intensities. This method has intermediate accuracy.
-
-.IP \fBphased\fR
-.PD
-As 'interpolate', but take phase values into account. This is the most accurate method, but the slowest.
-
-.SH SPECTRUM TYPES
-
-The available options for \fB--spectrum\fR are:
-
-.IP \fBtophat\fR
-.PD
-The spectrum samples will be distributed equidistantly either side of the specified photon energy to give a uniform distribution.
-
-.IP \fBsase\fR
-.PD
-A self-amplified spontaneous emission (SASE) spectrum will be simulated, as follows. First, a central photon energy will be chosen using a Gaussian distribution centered on the specified photon energy with a standard deviation of 8 eV. A Gaussian spectrum will then be calculated using the specified bandwidth, and noise added to simulatie the SASE 'spikes'.
-
-.IP \fBtwocolour\fR
-.PD
-The spectrum will consist of two Gaussian peaks separated by the specified bandwidth, each with a standard deviation of one fifth the specified bandwidth.
-
-.SH AUTHOR
-This page was written by Thomas White.
-
-.SH REPORTING BUGS
-Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
-
-.SH COPYRIGHT AND DISCLAIMER
-Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
-.P
-pattern_sim, and this manual, are part of CrystFEL.
-.P
-CrystFEL is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-.P
-CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-.P
-You should have received a copy of the GNU General Public License along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
-
-.SH SEE ALSO
-.BR crystfel (7)
-and
-.BR crystfel_geometry (5).
diff --git a/meson.build b/meson.build
index 21f4534a..65a65823 100644
--- a/meson.build
+++ b/meson.build
@@ -48,14 +48,6 @@ if cairodep.found()
conf_data.set10('HAVE_CAIRO', true)
endif
-opencldep = dependency('OpenCL', required: get_option('opencl'))
-if opencldep.found()
- conf_data.set10('HAVE_OPENCL', true)
- if cc.has_header('CL/cl.h', dependencies: opencldep)
- conf_data.set10('HAVE_CL_CL_H', true)
- endif
-endif
-
zmqdep = dependency('libzmq', required: false)
if zmqdep.found()
conf_data.set10('HAVE_ZMQ', true)
@@ -88,13 +80,6 @@ executable('cell_tool',
install: true,
install_rpath: '$ORIGIN/../lib64/:$ORIGIN/../lib')
-# partial_sim
-executable('partial_sim',
- ['src/partial_sim.c', versionc],
- dependencies: [mdep, libcrystfeldep, gsldep, pthreaddep],
- install: true,
- install_rpath: '$ORIGIN/../lib64/:$ORIGIN/../lib')
-
# process_hkl
process_hkl = executable('process_hkl',
['src/process_hkl.c', versionc],
@@ -241,21 +226,6 @@ if gtkdep.found()
install_dir: get_option('datadir') / 'applications')
endif
-# pattern_sim
-if hdf5dep.found()
- pattern_sim_sources = ['src/pattern_sim.c', 'src/diffraction.c', versionc]
- if opencldep.found()
- pattern_sim_sources += ['src/diffraction-gpu.c', 'src/cl-utils.c']
- endif
- executable('pattern_sim', pattern_sim_sources,
- dependencies: [mdep, libcrystfeldep, gsldep, hdf5dep, opencldep],
- install: true,
- install_rpath: '$ORIGIN/../lib64/:$ORIGIN/../lib')
-endif
-# For testing (see tests/meson.build)
-simulation_bits = files(['src/diffraction.c',
- 'src/diffraction-gpu.c',
- 'src/cl-utils.c'])
# ************************ Misc resources ************************
@@ -289,8 +259,6 @@ install_man(['doc/man/ambigator.1',
'doc/man/list_events.1',
'doc/man/make_pixelmap.1',
'doc/man/partialator.1',
- 'doc/man/partial_sim.1',
- 'doc/man/pattern_sim.1',
'doc/man/process_hkl.1',
'doc/man/render_hkl.1',
'doc/man/whirligig.1'])
diff --git a/src/cl-utils.c b/src/cl-utils.c
deleted file mode 100644
index 1b6667c4..00000000
--- a/src/cl-utils.c
+++ /dev/null
@@ -1,354 +0,0 @@
-/*
- * cl-utils.c
- *
- * OpenCL utility functions
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2010-2020 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#define CL_TARGET_OPENCL_VERSION 220
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-#include "utils.h"
-
-
-/* Return 1 if a GPU device is present, 0 if not, 2 on error. */
-int have_gpu_device()
-{
- cl_uint nplat;
- cl_platform_id platforms[8];
- cl_context_properties prop[3];
- cl_int err;
- int i;
-
- err = clGetPlatformIDs(8, platforms, &nplat);
- if ( err != CL_SUCCESS ) return 2;
- if ( nplat == 0 ) return 0;
-
- /* Find a GPU platform in the list */
- for ( i=0; i<nplat; i++ ) {
-
- prop[0] = CL_CONTEXT_PLATFORM;
- prop[1] = (cl_context_properties)platforms[i];
- prop[2] = 0;
-
- clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
- NULL, NULL, &err);
-
- if ( err != CL_SUCCESS ) {
- if ( err != CL_DEVICE_NOT_FOUND ) return 2;
- } else {
- return 1;
- }
- }
-
- return 0;
-}
-
-
-const char *clError(cl_int err)
-{
- switch ( err ) {
-
- case CL_SUCCESS :
- return "no error";
-
- case CL_DEVICE_NOT_AVAILABLE :
- return "device not available";
-
- case CL_DEVICE_NOT_FOUND :
- return "device not found";
-
- case CL_INVALID_DEVICE_TYPE :
- return "invalid device type";
-
- case CL_INVALID_PLATFORM :
- return "invalid platform";
-
- case CL_INVALID_KERNEL :
- return "invalid kernel";
-
- case CL_INVALID_ARG_INDEX :
- return "invalid argument index";
-
- case CL_INVALID_ARG_VALUE :
- return "invalid argument value";
-
- case CL_INVALID_MEM_OBJECT :
- return "invalid memory object";
-
- case CL_INVALID_SAMPLER :
- return "invalid sampler";
-
- case CL_INVALID_ARG_SIZE :
- return "invalid argument size";
-
- case CL_INVALID_COMMAND_QUEUE :
- return "invalid command queue";
-
- case CL_INVALID_CONTEXT :
- return "invalid context";
-
- case CL_INVALID_VALUE :
- return "invalid value";
-
- case CL_INVALID_EVENT_WAIT_LIST :
- return "invalid wait list";
-
- case CL_MAP_FAILURE :
- return "map failure";
-
- case CL_MEM_OBJECT_ALLOCATION_FAILURE :
- return "object allocation failure";
-
- case CL_OUT_OF_HOST_MEMORY :
- return "out of host memory";
-
- case CL_OUT_OF_RESOURCES :
- return "out of resources";
-
- case CL_INVALID_KERNEL_NAME :
- return "invalid kernel name";
-
- case CL_INVALID_KERNEL_ARGS :
- return "invalid kernel arguments";
-
- case CL_INVALID_WORK_GROUP_SIZE :
- return "invalid work group size";
-
- case CL_IMAGE_FORMAT_NOT_SUPPORTED :
- return "image format not supported";
-
- case CL_INVALID_WORK_DIMENSION :
- return "invalid work dimension";
-
- default :
- return "unknown error";
- }
-}
-
-
-static char *get_device_string(cl_device_id dev, cl_device_info info)
-{
- int r;
- size_t size;
- char *val;
-
- r = clGetDeviceInfo(dev, info, 0, NULL, &size);
- if ( r != CL_SUCCESS ) {
- ERROR("Couldn't get device vendor size: %s\n",
- clError(r));
- return NULL;
- }
- val = malloc(size);
- r = clGetDeviceInfo(dev, info, size, val, NULL);
- if ( r != CL_SUCCESS ) {
- ERROR("Couldn't get dev vendor: %s\n", clError(r));
- return NULL;
- }
-
- return val;
-}
-
-
-cl_device_id get_cl_dev(cl_context ctx, int n)
-{
- cl_device_id *dev;
- cl_int r;
- size_t size;
- int num_devs;
-
- /* Get the required size of the array */
- r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, 0, NULL, &size);
- if ( r != CL_SUCCESS ) {
- ERROR("Couldn't get array size for devices: %s\n", clError(r));
- return 0;
- }
-
- dev = malloc(size);
- r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, size, dev, NULL);
- if ( r != CL_SUCCESS ) {
- ERROR("Couldn't get device: %s\n", clError(r));
- return 0;
- }
- num_devs = size/sizeof(cl_device_id);
-
- if ( n >= num_devs ) {
- ERROR("Device ID out of range\n");
- return 0;
- }
-
- if ( n < 0 ) {
-
- int i;
-
- STATUS("Available devices:\n");
- for ( i=0; i<num_devs; i++ ) {
-
- char *vendor;
- char *name;
-
- vendor = get_device_string(dev[i], CL_DEVICE_VENDOR);
- name = get_device_string(dev[i], CL_DEVICE_NAME);
-
- STATUS("Device %i: %s %s\n", i, vendor, name);
-
- }
- n = 0;
-
- STATUS("Using device 0. Use --gpu-dev to choose another.\n");
-
- } else {
-
- char *vendor;
- char *name;
-
- vendor = get_device_string(dev[n], CL_DEVICE_VENDOR);
- name = get_device_string(dev[n], CL_DEVICE_NAME);
-
- STATUS("Using device %i: %s %s\n", n, vendor, name);
-
- }
-
- return dev[n];
-}
-
-
-static void show_build_log(cl_program prog, cl_device_id dev)
-{
- cl_int r;
- char log[4096];
- size_t s;
-
- r = clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, 4096, log,
- &s);
-
- STATUS("Status: %i\n", r);
- STATUS("%s\n", log);
-}
-
-
-cl_program load_program_from_string(const char *source_in, size_t len,
- cl_context ctx, cl_device_id dev,
- cl_int *err, const char *extra_cflags,
- const char *insert_stuff)
-{
- cl_program prog;
- cl_int r;
- char cflags[1024] = "";
- char *insert_pos;
- char *source;
-
- /* Copy the code because we need to zero-terminate it */
- source = malloc(len+1);
- if ( source == NULL ) return 0;
- memcpy(source, source_in, len);
- source[len] = '\0';
-
- if ( insert_stuff != NULL ) {
- insert_pos = strstr(source, "INSERT_HERE");
-
- if ( insert_pos != NULL ) {
-
- char *source2;
- size_t il;
-
- source2 = malloc(strlen(source)+strlen(insert_stuff)+1);
- if ( source2 == NULL ) {
- free(source);
- return 0;
- }
-
- il = insert_pos - source;
- memcpy(source2, source, il);
- memcpy(source2+il, insert_stuff, strlen(insert_stuff)+1);
- memcpy(source2+il+strlen(insert_stuff),
- source+il+11, strlen(source+il+11)+1);
- free(source);
- source = source2;
-
- }
- }
-
- prog = clCreateProgramWithSource(ctx, 1, (const char **)&source,
- NULL, err);
- if ( *err != CL_SUCCESS ) {
- ERROR("Couldn't load source\n");
- return 0;
- }
-
- cflags[0] = '\0';
- strncat(cflags, "-cl-no-signed-zeros ", 1023-strlen(cflags));
- strncat(cflags, extra_cflags, 1023-strlen(cflags));
-
- r = clBuildProgram(prog, 0, NULL, cflags, NULL, NULL);
- if ( r != CL_SUCCESS ) {
- ERROR("Couldn't build program\n");
- show_build_log(prog, dev);
- *err = r;
- return 0;
- }
-
- free(source);
- *err = CL_SUCCESS;
- return prog;
-}
-
-
-cl_program load_program(const char *filename, cl_context ctx,
- cl_device_id dev, cl_int *err, const char *extra_cflags,
- const char *insert_stuff)
-{
- FILE *fh;
- char *source;
- size_t len;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) {
- ERROR("Couldn't open '%s'\n", filename);
- *err = CL_INVALID_PROGRAM;
- return 0;
- }
- source = malloc(16384);
- if ( source == NULL ) {
- fclose(fh);
- return 0;
- }
- len = fread(source, 1, 16383, fh);
- fclose(fh);
-
- return load_program_from_string(source, len, ctx, dev,err, extra_cflags,
- insert_stuff);
-}
diff --git a/src/cl-utils.h b/src/cl-utils.h
deleted file mode 100644
index d72ee6d5..00000000
--- a/src/cl-utils.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/*
- * cl-utils.h
- *
- * OpenCL utility functions
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2010-2020 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifndef CLUTILS_H
-#define CLUTILS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#define CL_TARGET_OPENCL_VERSION 220
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-extern int have_gpu_device(void);
-extern const char *clError(cl_int err);
-extern cl_device_id get_cl_dev(cl_context ctx, int n);
-extern cl_program load_program_from_string(const char *source_in, size_t len,
- cl_context ctx, cl_device_id dev,
- cl_int *err, const char *extra_cflags,
- const char *insert_stuff);
-extern cl_program load_program(const char *filename, cl_context ctx,
- cl_device_id dev, cl_int *err,
- const char *extra_cflags, const char *insert_stuff);
-
-
-#endif /* CLUTILS_H */
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
deleted file mode 100644
index 97badc36..00000000
--- a/src/diffraction-gpu.c
+++ /dev/null
@@ -1,596 +0,0 @@
-/*
- * diffraction-gpu.c
- *
- * Calculate diffraction patterns by Fourier methods (GPU version)
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2009-2020 Thomas White <taw@physics.org>
- * 2013 Alexandra Tolstikova
- * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <complex.h>
-
-#define CL_TARGET_OPENCL_VERSION 220
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "diffraction.h"
-#include "cl-utils.h"
-#include "pattern_sim.h"
-
-#include "diffraction.cl.h"
-
-#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,
- int no_fringes, int flat)
-{
- 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; i<n; i++ ) {
- gctx->sinc_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; i<SINC_LUT_ELEMENTS; i++ ) {
- gctx->sinc_lut_ptrs[n-1][i] = 1.0;
- }
- } else {
- for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
- double x, val;
- x = (double)i/SINC_LUT_ELEMENTS;
- if ( (flat || no_fringes) && (x > 1.0/n) && (1.0-x > 1.0/n) ) {
- val = 0.0;
- } else if ( flat ) {
- val = n;
- } else {
- val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
- }
- gctx->sinc_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;
-}
-
-
-static int do_panels(struct gpu_context *gctx, struct image *image,
- double k, double weight,
- int *n_inf, int *n_neg, int *n_nan)
-{
- int i;
- const int sampling = 4; /* This, squared, number of samples / pixel */
-
- if ( set_arg_float(gctx, 1, k) ) return 1;
- if ( set_arg_float(gctx, 2, weight) ) return 1;
-
- /* Iterate over panels */
- for ( i=0; i<image->detgeom->n_panels; i++ ) {
-
- size_t dims[2];
- size_t ldims[2];
- struct detgeom_panel *p;
- cl_mem diff;
- size_t diff_size;
- float *diff_ptr;
- int fs, ss;
- cl_int err;
-
- p = &image->detgeom->panels[i];
-
- /* Buffer for the results of this panel */
- diff_size = p->w * p->h * 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 1;
- }
-
- if ( set_arg_mem(gctx, 0, diff) ) return 1;
-
- if ( set_arg_int(gctx, 3, p->w) ) return 1;
- if ( set_arg_float(gctx, 4, p->cnx) ) return 1;
- if ( set_arg_float(gctx, 5, p->cny) ) return 1;
- if ( set_arg_float(gctx, 6, p->fsx) ) return 1;
- if ( set_arg_float(gctx, 7, p->fsy) ) return 1;
- if ( set_arg_float(gctx, 8, p->fsz) ) return 1;
- 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, 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;
-
- ldims[0] = sampling;
- ldims[1] = sampling;
-
- err = clSetKernelArg(gctx->kern, 20,
- sampling*sampling*sizeof(cl_float), NULL);
- if ( err != CL_SUCCESS ) {
- ERROR("Couldn't set local memory: %s\n", clError(err));
- return 1;
- }
-
- err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 2, NULL,
- dims, ldims, 0, NULL, NULL);
- if ( err != CL_SUCCESS ) {
- ERROR("Couldn't enqueue diffraction kernel: %s\n",
- clError(err));
- return 1;
- }
-
- 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 1;
- }
-
- for ( ss=0; ss<p->h; ss++ ) {
- for ( fs=0; fs<p->w; fs++ ) {
-
- float val;
-
- val = diff_ptr[fs + p->w*ss];
- if ( isinf(val) ) (*n_inf)++;
- if ( val < 0.0 ) (*n_neg)++;
- if ( isnan(val) ) (*n_nan)++;
-
- image->dp[i][fs + p->w*ss] += val;
-
- }
- }
-
- clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr,
- 0, NULL, NULL);
-
- clReleaseMemObject(diff);
-
- }
-
- return 0;
-}
-
-
-int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
- int na, int nb, int nc, UnitCell *ucell,
- int no_fringes, int flat, int n_samples)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- cl_float16 cell;
- cl_int err;
- int n_inf = 0;
- int n_neg = 0;
- int n_nan = 0;
- int i;
- double kmin, kmax, step, norm;
-
- if ( gctx == NULL ) {
- ERROR("GPU setup failed.\n");
- return 1;
- }
-
- /* Ensure all required LUTs are available */
- check_sinc_lut(gctx, na, no_fringes, flat);
- check_sinc_lut(gctx, nb, no_fringes, flat);
- check_sinc_lut(gctx, nc, no_fringes, flat);
-
- /* Unit cell */
- 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;
-
- err = clSetKernelArg(gctx->kern, 14, sizeof(cl_float16), &cell);
- if ( err != CL_SUCCESS ) {
- ERROR("Couldn't set unit cell: %s\n", clError(err));
- return 1;
- }
-
- if ( set_arg_mem(gctx, 15, gctx->intensities) ) return 1;
- if ( set_arg_mem(gctx, 16, gctx->flags) ) return 1;
- if ( set_arg_mem(gctx, 17, gctx->sinc_luts[na-1]) ) return 1;
- if ( set_arg_mem(gctx, 18, gctx->sinc_luts[nb-1]) ) return 1;
- if ( set_arg_mem(gctx, 19, gctx->sinc_luts[nc-1]) ) return 1;
-
- /* Allocate memory for the result */
- 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->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);
- return 1;
- }
- }
-
- spectrum_get_range(image->spectrum, &kmin, &kmax);
- step = (kmax-kmin)/(n_samples+1);
-
- /* Determine normalisation factor such that weights add up to 1 after
- * sampling (bins must have constant width) */
- norm = 0.0;
- for ( i=1; i<=n_samples; i++ ) {
- double k = kmin + i*step;
- norm += spectrum_get_density_at_k(image->spectrum, k);
- }
- for ( i=1; i<=n_samples; i++ ) {
-
- double k = kmin + i*step;
- double prob;
-
- /* Probability = p.d.f. times step width */
- prob = spectrum_get_density_at_k(image->spectrum, k)/norm;
-
- STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob);
-
- err = do_panels(gctx, image, k, prob, &n_inf, &n_neg, &n_nan);
- if ( err ) return 1;
-
- }
-
- 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);
- }
-
- return 0;
-}
-
-
-/* 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 iplat;
- int have_ctx = 0;
- char cflags[512] = "";
- char *insert_stuff = NULL;
-
- 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;
- }
-
- /* Find a GPU platform in the list */
- for ( iplat=0; iplat<nplat; iplat++ ) {
-
- prop[0] = CL_CONTEXT_PLATFORM;
- prop[1] = (cl_context_properties)platforms[iplat];
- prop[2] = 0;
-
- gctx = malloc(sizeof(*gctx));
- gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
- NULL, NULL, &err);
-
- if ( err != CL_SUCCESS ) {
- if ( err == CL_DEVICE_NOT_FOUND ) {
- /* No GPU device in this platform */
- continue; /* Try next platform */
- } else {
- ERROR("Couldn't create OpenCL context: %s (%i)\n",
- clError(err), err);
- free(gctx);
- return NULL;
- }
- } else {
- STATUS("Using OpenCL platform %i (%i total)\n",
- iplat, nplat);
- have_ctx = 1;
- break;
- }
- }
-
- if ( !have_ctx ) {
- ERROR("Couldn't find a GPU device in any platform.\n");
- 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 ) {
- int i;
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- intensities_ptr[i] = intensities[i];
- }
- } else {
- int i;
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- intensities_ptr[i] = 100.0; /* Does nothing */
- }
- strncat(cflags, "-DFLAT_INTENSITIES ", 511-strlen(cflags));
- }
- gctx->intensities = 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 ) {
-
- int i, n;
- SymOpList *pg;
- size_t islen = 0;
-
- insert_stuff = malloc(16384);
- if ( insert_stuff == NULL ) return NULL;
- insert_stuff[0] = '\0';
-
- pg = get_pointgroup(sym);
- n = num_equivs(pg, NULL);
- for ( i=0; i<n; i++ ) {
-
- IntegerMatrix *op = get_symop(pg, NULL, i);
- char line[1024];
-
- snprintf(line, 1023,
- "val += lookup_flagged_intensity(intensities, "
- "flags, %s, %s, %s);\n\t",
- get_matrix_name(op, 0),
- get_matrix_name(op, 1),
- get_matrix_name(op, 2));
-
- islen += strlen(line);
- if ( islen > 16383 ) {
- ERROR("Too many symmetry operators.\n");
- return NULL;
- }
- strcat(insert_stuff, line);
-
- }
-
- free_symoplist(pg);
-
- printf("Inserting --->%s<---\n", insert_stuff);
-
- } 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 ) {
- int i;
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- flags_ptr[i] = flags[i];
- }
- } else {
- int i;
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- flags_ptr[i] = 1.0;
- }
- }
- gctx->flags = 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_from_string((char *)data_diffraction_cl,
- data_diffraction_cl_len, gctx->ctx,
- dev, &err, cflags, insert_stuff);
- 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
deleted file mode 100644
index 1128c8b1..00000000
--- a/src/diffraction-gpu.h
+++ /dev/null
@@ -1,77 +0,0 @@
-/*
- * diffraction-gpu.h
- *
- * Calculate diffraction patterns by Fourier methods (GPU version)
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2010-2019 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DIFFRACTION_GPU_H
-#define DIFFRACTION_GPU_H
-
-#include "image.h"
-#include "cell.h"
-
-struct gpu_context;
-
-#ifdef HAVE_OPENCL
-
-extern int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
- int na, int nb, int nc, UnitCell *ucell,
- int no_fringes, int flat, int n_samples);
-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 int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
- int na, int nb, int nc, UnitCell *ucell,
- int no_fringes, int flat, int n_samples)
-{
- /* Do nothing */
- ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n");
- return 1;
-}
-
-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
deleted file mode 100644
index c0dba8aa..00000000
--- a/src/diffraction.c
+++ /dev/null
@@ -1,502 +0,0 @@
-/*
- * diffraction.c
- *
- * Calculate diffraction patterns by Fourier methods
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2009-2020 Thomas White <taw@physics.org>
- * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- * 2013 Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <complex.h>
-#include <assert.h>
-#include <fenv.h>
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "diffraction.h"
-#include "symmetry.h"
-#include "pattern_sim.h"
-
-
-#define SINC_LUT_ELEMENTS (4096)
-
-
-static double *get_sinc_lut(int n, int no_fringes, int flat)
-{
- int i;
- double *lut;
-
- lut = malloc(SINC_LUT_ELEMENTS*sizeof(double));
- lut[0] = n;
- if ( n == 1 ) {
- for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
- lut[i] = 1.0;
- }
- } else {
- for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
- double x, val;
- x = (double)i/SINC_LUT_ELEMENTS;
- if ( (flat || no_fringes) && (x > 1.0/n) && (1.0-x > 1.0/n) ) {
- val = 0.0;
- } else if ( flat ) {
- val = n;
- } else {
- val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
- }
- lut[i] = val;
- }
- }
-
- return lut;
-}
-
-
-static double interpolate_lut(double *lut, double val)
-{
- double i, pos, f;
- unsigned int low, high;
-
- pos = SINC_LUT_ELEMENTS * modf(fabs(val), &i);
- low = (int)pos; /* Discard fractional part */
- high = low + 1;
- f = modf(pos, &i); /* Fraction */
- if ( high == SINC_LUT_ELEMENTS ) high = 0;
-
- return (1.0-f)*lut[low] + f*lut[high];
-}
-
-
-static double lattice_factor(struct rvec q, double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz,
- double *lut_a, double *lut_b,
- double *lut_c)
-{
- struct rvec Udotq;
- double f1, f2, f3;
-
- Udotq.u = ax*q.u + ay*q.v + az*q.w;
- Udotq.v = bx*q.u + by*q.v + bz*q.w;
- Udotq.w = cx*q.u + cy*q.v + cz*q.w;
-
- f1 = interpolate_lut(lut_a, Udotq.u);
- f2 = interpolate_lut(lut_b, Udotq.v);
- f3 = interpolate_lut(lut_c, Udotq.w);
-
- return f1 * f2 * f3;
-}
-
-
-static double sym_lookup_intensity(const double *intensities,
- const unsigned char *flags,
- const SymOpList *sym,
- signed int h, signed int k, signed int l)
-{
- int i;
- double ret = 0.0;
-
- for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
- signed int he;
- signed int ke;
- signed int le;
- double f, val;
-
- get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-
- f = (double)lookup_arr_flag(flags, he, ke, le);
- val = lookup_arr_intensity(intensities, he, ke, le);
-
- ret += f*val;
-
- }
-
- return ret;
-}
-
-
-static double sym_lookup_phase(const double *phases,
- const unsigned char *flags, const SymOpList *sym,
- signed int h, signed int k, signed int l)
-{
- int i;
-
- for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
- signed int he;
- signed int ke;
- signed int le;
- int f;
-
- get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-
- f = lookup_arr_flag(flags, he, ke, le);
-
- if ( f ) return lookup_arr_phase(phases, he, ke, le);
-
- }
-
- return 0.0;
-}
-
-
-static double interpolate_linear(const double *ref, const unsigned char *flags,
- const SymOpList *sym, float hd,
- signed int k, signed int l)
-{
- signed int h;
- double val1, val2;
- float f;
-
- 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);
-
- 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);
-
- /* 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 100.0;
-
- 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;
-}
-
-
-static void diffraction_panel(struct image *image, const double *intensities,
- const double *phases, const unsigned char *flags,
- UnitCell *cell, GradientMethod m,
- const SymOpList *sym, double k,
- double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz,
- double *lut_a, double *lut_b, double *lut_c,
- int pn, double weight)
-{
- int fs, ss;
- const int nxs = 4;
- const int nys = 4;
- struct detgeom_panel *p = &image->detgeom->panels[pn];
-
- weight /= nxs*nys;
-
- for ( ss=0; ss<p->h; ss++ ) {
- for ( fs=0; fs<p->w; fs++ ) {
-
- int idx;
- double f_lattice, I_lattice;
- double I_molecule;
- 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;
-
- detgeom_transform_coords(p, fs+xo, ss+yo,
- 1.0/k, 0.0, 0.0, qv);
-
- q.u = qv[0]; q.v = qv[1]; q.w = qv[2];
-
- f_lattice = lattice_factor(q, ax, ay, az,
- bx, by, bz,
- cx, cy, cz,
- lut_a, lut_b, lut_c);
-
- I_molecule = molecule_factor(intensities,
- phases, flags, q,
- ax, ay, az,
- bx, by, bz,
- cx, cy, cz,
- m, sym);
-
- I_lattice = pow(f_lattice, 2.0);
-
- idx = fs + p->w*ss;
- image->dp[pn][idx] += I_lattice * I_molecule * weight;
-
- }
- }
- }
- progress_bar(ss, p->h-1, "Calculating diffraction");
- }
-}
-
-
-static void diffraction_at_k(struct image *image, const double *intensities,
- const double *phases, const unsigned char *flags,
- UnitCell *cell, GradientMethod m,
- const SymOpList *sym, double k,
- double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz,
- double *lut_a, double *lut_b, double *lut_c,
- double weight)
-{
- int 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);
- }
-}
-
-
-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,
- int no_fringes, int flat, int n_samples)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- double *lut_a;
- double *lut_b;
- double *lut_c;
- int i;
- double kmin, kmax, step;
- double norm = 0.0;
-
- cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
- lut_a = get_sinc_lut(na, no_fringes, flat);
- lut_b = get_sinc_lut(nb, no_fringes, flat);
- lut_c = get_sinc_lut(nc, no_fringes, flat);
-
- spectrum_get_range(image->spectrum, &kmin, &kmax);
- step = (kmax-kmin)/(n_samples+1);
-
- /* Determine normalisation factor such that weights add up to 1 after
- * sampling (bins must have constant width) */
- for ( i=1; i<=n_samples; i++ ) {
- double k = kmin + i*step;
- norm += spectrum_get_density_at_k(image->spectrum, k);
- }
- for ( i=1; i<=n_samples; i++ ) {
-
- double k = kmin + i*step;
- double prob;
-
- /* Probability = p.d.f. times step width */
- prob = spectrum_get_density_at_k(image->spectrum, k)/norm;
-
- STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob);
-
- diffraction_at_k(image, intensities, phases,
- flags, cell, m, sym, k,
- ax, ay, az, bx, by, bz, cx, cy, cz,
- lut_a, lut_b, lut_c, prob);
-
- }
-
- free(lut_a);
- free(lut_b);
- free(lut_c);
-}
diff --git a/src/diffraction.cl.h b/src/diffraction.cl.h
deleted file mode 100644
index 1f788b67..00000000
--- a/src/diffraction.cl.h
+++ /dev/null
@@ -1,513 +0,0 @@
-/*
- * This file was generated from data/diffraction.cl
- * using the following command:
- * xxd -i data/diffraction.cl src/diffraction.cl.h
- *
- * If you have 'xxd' installed, you can run the script
- * data/gen-resources to re-create this file.
- */
-
-unsigned char data_diffraction_cl[] = {
- 0x2f, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x64, 0x69, 0x66, 0x66, 0x72, 0x61,
- 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x2e, 0x63, 0x6c, 0x0a, 0x20, 0x2a, 0x0a,
- 0x20, 0x2a, 0x20, 0x47, 0x50, 0x55, 0x20, 0x63, 0x61, 0x6c, 0x63, 0x75,
- 0x6c, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x6b, 0x65, 0x72, 0x6e, 0x65,
- 0x6c, 0x20, 0x66, 0x6f, 0x72, 0x20, 0x74, 0x72, 0x75, 0x6e, 0x63, 0x61,
- 0x74, 0x65, 0x64, 0x20, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20,
- 0x64, 0x69, 0x66, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x0a,
- 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x6f, 0x70, 0x79, 0x72, 0x69,
- 0x67, 0x68, 0x74, 0x20, 0xc2, 0xa9, 0x20, 0x32, 0x30, 0x31, 0x32, 0x2d,
- 0x32, 0x30, 0x31, 0x34, 0x20, 0x44, 0x65, 0x75, 0x74, 0x73, 0x63, 0x68,
- 0x65, 0x73, 0x20, 0x45, 0x6c, 0x65, 0x6b, 0x74, 0x72, 0x6f, 0x6e, 0x65,
- 0x6e, 0x2d, 0x53, 0x79, 0x6e, 0x63, 0x68, 0x72, 0x6f, 0x74, 0x72, 0x6f,
- 0x6e, 0x20, 0x44, 0x45, 0x53, 0x59, 0x2c, 0x0a, 0x20, 0x2a, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x61, 0x20, 0x72,
- 0x65, 0x73, 0x65, 0x61, 0x72, 0x63, 0x68, 0x20, 0x63, 0x65, 0x6e, 0x74,
- 0x72, 0x65, 0x20, 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, 0x48, 0x65,
- 0x6c, 0x6d, 0x68, 0x6f, 0x6c, 0x74, 0x7a, 0x20, 0x41, 0x73, 0x73, 0x6f,
- 0x63, 0x69, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2e, 0x0a, 0x20, 0x2a, 0x0a,
- 0x20, 0x2a, 0x20, 0x41, 0x75, 0x74, 0x68, 0x6f, 0x72, 0x73, 0x3a, 0x0a,
- 0x20, 0x2a, 0x20, 0x20, 0x20, 0x32, 0x30, 0x30, 0x39, 0x2d, 0x32, 0x30,
- 0x31, 0x34, 0x20, 0x54, 0x68, 0x6f, 0x6d, 0x61, 0x73, 0x20, 0x57, 0x68,
- 0x69, 0x74, 0x65, 0x20, 0x3c, 0x74, 0x61, 0x77, 0x40, 0x70, 0x68, 0x79,
- 0x73, 0x69, 0x63, 0x73, 0x2e, 0x6f, 0x72, 0x67, 0x3e, 0x0a, 0x20, 0x2a,
- 0x20, 0x20, 0x20, 0x32, 0x30, 0x31, 0x33, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x41, 0x6c, 0x65, 0x78, 0x61, 0x6e, 0x64, 0x72, 0x61, 0x20, 0x54,
- 0x6f, 0x6c, 0x73, 0x74, 0x69, 0x6b, 0x6f, 0x76, 0x61, 0x0a, 0x20, 0x2a,
- 0x0a, 0x20, 0x2a, 0x20, 0x54, 0x68, 0x69, 0x73, 0x20, 0x66, 0x69, 0x6c,
- 0x65, 0x20, 0x69, 0x73, 0x20, 0x70, 0x61, 0x72, 0x74, 0x20, 0x6f, 0x66,
- 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, 0x4c, 0x2e, 0x0a, 0x20,
- 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45,
- 0x4c, 0x20, 0x69, 0x73, 0x20, 0x66, 0x72, 0x65, 0x65, 0x20, 0x73, 0x6f,
- 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x3a, 0x20, 0x79, 0x6f, 0x75, 0x20,
- 0x63, 0x61, 0x6e, 0x20, 0x72, 0x65, 0x64, 0x69, 0x73, 0x74, 0x72, 0x69,
- 0x62, 0x75, 0x74, 0x65, 0x20, 0x69, 0x74, 0x20, 0x61, 0x6e, 0x64, 0x2f,
- 0x6f, 0x72, 0x20, 0x6d, 0x6f, 0x64, 0x69, 0x66, 0x79, 0x0a, 0x20, 0x2a,
- 0x20, 0x69, 0x74, 0x20, 0x75, 0x6e, 0x64, 0x65, 0x72, 0x20, 0x74, 0x68,
- 0x65, 0x20, 0x74, 0x65, 0x72, 0x6d, 0x73, 0x20, 0x6f, 0x66, 0x20, 0x74,
- 0x68, 0x65, 0x20, 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, 0x6e, 0x65, 0x72,
- 0x61, 0x6c, 0x20, 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, 0x20, 0x4c, 0x69,
- 0x63, 0x65, 0x6e, 0x73, 0x65, 0x20, 0x61, 0x73, 0x20, 0x70, 0x75, 0x62,
- 0x6c, 0x69, 0x73, 0x68, 0x65, 0x64, 0x20, 0x62, 0x79, 0x0a, 0x20, 0x2a,
- 0x20, 0x74, 0x68, 0x65, 0x20, 0x46, 0x72, 0x65, 0x65, 0x20, 0x53, 0x6f,
- 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x20, 0x46, 0x6f, 0x75, 0x6e, 0x64,
- 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2c, 0x20, 0x65, 0x69, 0x74, 0x68, 0x65,
- 0x72, 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x20, 0x33, 0x20,
- 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e,
- 0x73, 0x65, 0x2c, 0x20, 0x6f, 0x72, 0x0a, 0x20, 0x2a, 0x20, 0x28, 0x61,
- 0x74, 0x20, 0x79, 0x6f, 0x75, 0x72, 0x20, 0x6f, 0x70, 0x74, 0x69, 0x6f,
- 0x6e, 0x29, 0x20, 0x61, 0x6e, 0x79, 0x20, 0x6c, 0x61, 0x74, 0x65, 0x72,
- 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x2e, 0x0a, 0x20, 0x2a,
- 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, 0x4c,
- 0x20, 0x69, 0x73, 0x20, 0x64, 0x69, 0x73, 0x74, 0x72, 0x69, 0x62, 0x75,
- 0x74, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x20, 0x74, 0x68, 0x65, 0x20, 0x68,
- 0x6f, 0x70, 0x65, 0x20, 0x74, 0x68, 0x61, 0x74, 0x20, 0x69, 0x74, 0x20,
- 0x77, 0x69, 0x6c, 0x6c, 0x20, 0x62, 0x65, 0x20, 0x75, 0x73, 0x65, 0x66,
- 0x75, 0x6c, 0x2c, 0x0a, 0x20, 0x2a, 0x20, 0x62, 0x75, 0x74, 0x20, 0x57,
- 0x49, 0x54, 0x48, 0x4f, 0x55, 0x54, 0x20, 0x41, 0x4e, 0x59, 0x20, 0x57,
- 0x41, 0x52, 0x52, 0x41, 0x4e, 0x54, 0x59, 0x3b, 0x20, 0x77, 0x69, 0x74,
- 0x68, 0x6f, 0x75, 0x74, 0x20, 0x65, 0x76, 0x65, 0x6e, 0x20, 0x74, 0x68,
- 0x65, 0x20, 0x69, 0x6d, 0x70, 0x6c, 0x69, 0x65, 0x64, 0x20, 0x77, 0x61,
- 0x72, 0x72, 0x61, 0x6e, 0x74, 0x79, 0x20, 0x6f, 0x66, 0x0a, 0x20, 0x2a,
- 0x20, 0x4d, 0x45, 0x52, 0x43, 0x48, 0x41, 0x4e, 0x54, 0x41, 0x42, 0x49,
- 0x4c, 0x49, 0x54, 0x59, 0x20, 0x6f, 0x72, 0x20, 0x46, 0x49, 0x54, 0x4e,
- 0x45, 0x53, 0x53, 0x20, 0x46, 0x4f, 0x52, 0x20, 0x41, 0x20, 0x50, 0x41,
- 0x52, 0x54, 0x49, 0x43, 0x55, 0x4c, 0x41, 0x52, 0x20, 0x50, 0x55, 0x52,
- 0x50, 0x4f, 0x53, 0x45, 0x2e, 0x20, 0x20, 0x53, 0x65, 0x65, 0x20, 0x74,
- 0x68, 0x65, 0x0a, 0x20, 0x2a, 0x20, 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65,
- 0x6e, 0x65, 0x72, 0x61, 0x6c, 0x20, 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63,
- 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e, 0x73, 0x65, 0x20, 0x66, 0x6f, 0x72,
- 0x20, 0x6d, 0x6f, 0x72, 0x65, 0x20, 0x64, 0x65, 0x74, 0x61, 0x69, 0x6c,
- 0x73, 0x2e, 0x0a, 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x59, 0x6f, 0x75,
- 0x20, 0x73, 0x68, 0x6f, 0x75, 0x6c, 0x64, 0x20, 0x68, 0x61, 0x76, 0x65,
- 0x20, 0x72, 0x65, 0x63, 0x65, 0x69, 0x76, 0x65, 0x64, 0x20, 0x61, 0x20,
- 0x63, 0x6f, 0x70, 0x79, 0x20, 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20,
- 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, 0x6e, 0x65, 0x72, 0x61, 0x6c, 0x20,
- 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e,
- 0x73, 0x65, 0x0a, 0x20, 0x2a, 0x20, 0x61, 0x6c, 0x6f, 0x6e, 0x67, 0x20,
- 0x77, 0x69, 0x74, 0x68, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45,
- 0x4c, 0x2e, 0x20, 0x20, 0x49, 0x66, 0x20, 0x6e, 0x6f, 0x74, 0x2c, 0x20,
- 0x73, 0x65, 0x65, 0x20, 0x3c, 0x68, 0x74, 0x74, 0x70, 0x3a, 0x2f, 0x2f,
- 0x77, 0x77, 0x77, 0x2e, 0x67, 0x6e, 0x75, 0x2e, 0x6f, 0x72, 0x67, 0x2f,
- 0x6c, 0x69, 0x63, 0x65, 0x6e, 0x73, 0x65, 0x73, 0x2f, 0x3e, 0x2e, 0x0a,
- 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x2f, 0x0a, 0x0a, 0x0a, 0x2f, 0x2a, 0x20,
- 0x4d, 0x61, 0x78, 0x6d, 0x69, 0x6d, 0x75, 0x6d, 0x20, 0x69, 0x6e, 0x64,
- 0x65, 0x78, 0x20, 0x74, 0x6f, 0x20, 0x68, 0x6f, 0x6c, 0x64, 0x20, 0x76,
- 0x61, 0x6c, 0x75, 0x65, 0x73, 0x20, 0x75, 0x70, 0x20, 0x74, 0x6f, 0x20,
- 0x28, 0x63, 0x61, 0x6e, 0x20, 0x62, 0x65, 0x20, 0x69, 0x6e, 0x63, 0x72,
- 0x65, 0x61, 0x73, 0x65, 0x64, 0x20, 0x69, 0x66, 0x20, 0x6e, 0x65, 0x63,
- 0x65, 0x73, 0x73, 0x61, 0x72, 0x79, 0x29, 0x0a, 0x20, 0x2a, 0x20, 0x57,
- 0x41, 0x52, 0x4e, 0x49, 0x4e, 0x47, 0x3a, 0x20, 0x41, 0x6c, 0x74, 0x65,
- 0x72, 0x69, 0x6e, 0x67, 0x20, 0x74, 0x68, 0x69, 0x73, 0x20, 0x76, 0x61,
- 0x6c, 0x75, 0x65, 0x20, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x69, 0x74, 0x75,
- 0x74, 0x65, 0x73, 0x20, 0x61, 0x6e, 0x20, 0x41, 0x42, 0x49, 0x20, 0x63,
- 0x68, 0x61, 0x6e, 0x67, 0x65, 0x2c, 0x20, 0x61, 0x6e, 0x64, 0x20, 0x6d,
- 0x65, 0x61, 0x6e, 0x73, 0x20, 0x79, 0x6f, 0x75, 0x20, 0x6d, 0x75, 0x73,
- 0x74, 0x0a, 0x20, 0x2a, 0x20, 0x75, 0x70, 0x64, 0x61, 0x74, 0x65, 0x20,
- 0x73, 0x72, 0x63, 0x2f, 0x70, 0x61, 0x74, 0x74, 0x65, 0x72, 0x6e, 0x5f,
- 0x73, 0x69, 0x6d, 0x2e, 0x68, 0x20, 0x74, 0x68, 0x65, 0x6e, 0x20, 0x72,
- 0x65, 0x63, 0x6f, 0x6d, 0x70, 0x69, 0x6c, 0x65, 0x20, 0x61, 0x6e, 0x64,
- 0x20, 0x72, 0x65, 0x69, 0x6e, 0x73, 0x74, 0x61, 0x6c, 0x6c, 0x20, 0x65,
- 0x76, 0x65, 0x72, 0x79, 0x74, 0x68, 0x69, 0x6e, 0x67, 0x2e, 0x20, 0x2a,
- 0x2f, 0x0a, 0x23, 0x64, 0x65, 0x66, 0x69, 0x6e, 0x65, 0x20, 0x49, 0x4e,
- 0x44, 0x4d, 0x41, 0x58, 0x20, 0x31, 0x33, 0x30, 0x0a, 0x23, 0x64, 0x65,
- 0x66, 0x69, 0x6e, 0x65, 0x20, 0x49, 0x44, 0x49, 0x4d, 0x20, 0x28, 0x49,
- 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x2a, 0x32, 0x20, 0x2b, 0x31, 0x29, 0x0a,
- 0x0a, 0x23, 0x69, 0x66, 0x6e, 0x64, 0x65, 0x66, 0x20, 0x4d, 0x5f, 0x50,
- 0x49, 0x0a, 0x23, 0x64, 0x65, 0x66, 0x69, 0x6e, 0x65, 0x20, 0x4d, 0x5f,
- 0x50, 0x49, 0x20, 0x28, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x29, 0x28,
- 0x33, 0x2e, 0x31, 0x34, 0x31, 0x35, 0x39, 0x32, 0x36, 0x35, 0x29, 0x29,
- 0x0a, 0x23, 0x65, 0x6e, 0x64, 0x69, 0x66, 0x0a, 0x0a, 0x63, 0x6f, 0x6e,
- 0x73, 0x74, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74,
- 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x61, 0x20, 0x3d,
- 0x20, 0x43, 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49,
- 0x5a, 0x45, 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54,
- 0x52, 0x55, 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43,
- 0x4c, 0x4b, 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52,
- 0x45, 0x50, 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c,
- 0x20, 0x43, 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f,
- 0x4c, 0x49, 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x63, 0x6f, 0x6e, 0x73,
- 0x74, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, 0x20,
- 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x62, 0x20, 0x3d, 0x20,
- 0x43, 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, 0x5a,
- 0x45, 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, 0x52,
- 0x55, 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, 0x4c,
- 0x4b, 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, 0x45,
- 0x50, 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20,
- 0x43, 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, 0x4c,
- 0x49, 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x63, 0x6f, 0x6e, 0x73, 0x74,
- 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, 0x20, 0x73,
- 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x63, 0x20, 0x3d, 0x20, 0x43,
- 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, 0x5a, 0x45,
- 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, 0x52, 0x55,
- 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, 0x4c, 0x4b,
- 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, 0x45, 0x50,
- 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43,
- 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, 0x4c, 0x49,
- 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61,
- 0x74, 0x34, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x71, 0x28, 0x66, 0x6c, 0x6f,
- 0x61, 0x74, 0x20, 0x66, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74,
- 0x20, 0x73, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x72,
- 0x65, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6c,
- 0x65, 0x6e, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x6b, 0x2c,
- 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, 0x72, 0x6e,
- 0x65, 0x72, 0x5f, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
- 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x2c, 0x0a, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x66,
- 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c,
- 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f,
- 0x61, 0x74, 0x20, 0x66, 0x73, 0x7a, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61,
- 0x74, 0x20, 0x73, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74,
- 0x20, 0x73, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
- 0x73, 0x73, 0x7a, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61,
- 0x74, 0x20, 0x72, 0x78, 0x2c, 0x20, 0x72, 0x79, 0x2c, 0x20, 0x72, 0x7a,
- 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x61, 0x7a, 0x2c,
- 0x20, 0x74, 0x74, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34,
- 0x20, 0x71, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x78,
- 0x73, 0x2c, 0x20, 0x79, 0x73, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61,
- 0x74, 0x20, 0x6b, 0x78, 0x2c, 0x20, 0x6b, 0x79, 0x2c, 0x20, 0x6b, 0x7a,
- 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x74, 0x74,
- 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75,
- 0x6c, 0x61, 0x74, 0x65, 0x20, 0x33, 0x44, 0x20, 0x70, 0x6f, 0x73, 0x69,
- 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x6f, 0x66, 0x20, 0x67, 0x69, 0x76, 0x65,
- 0x6e, 0x20, 0x70, 0x6f, 0x73, 0x69, 0x74, 0x69, 0x6f, 0x6e, 0x2c, 0x20,
- 0x69, 0x6e, 0x20, 0x6d, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x72, 0x78, 0x20,
- 0x3d, 0x20, 0x28, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x78, 0x20,
- 0x2b, 0x20, 0x66, 0x73, 0x2a, 0x66, 0x73, 0x78, 0x20, 0x2b, 0x20, 0x73,
- 0x73, 0x2a, 0x73, 0x73, 0x78, 0x29, 0x20, 0x2f, 0x20, 0x72, 0x65, 0x73,
- 0x3b, 0x0a, 0x09, 0x72, 0x79, 0x20, 0x3d, 0x20, 0x28, 0x63, 0x6f, 0x72,
- 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x20, 0x2b, 0x20, 0x66, 0x73, 0x2a, 0x66,
- 0x73, 0x79, 0x20, 0x2b, 0x20, 0x73, 0x73, 0x2a, 0x73, 0x73, 0x79, 0x29,
- 0x20, 0x2f, 0x20, 0x72, 0x65, 0x73, 0x3b, 0x0a, 0x09, 0x72, 0x7a, 0x20,
- 0x3d, 0x20, 0x63, 0x6c, 0x65, 0x6e, 0x20, 0x2b, 0x20, 0x28, 0x66, 0x73,
- 0x2a, 0x66, 0x73, 0x7a, 0x20, 0x2b, 0x20, 0x73, 0x73, 0x2a, 0x73, 0x73,
- 0x7a, 0x29, 0x2f, 0x72, 0x65, 0x73, 0x3b, 0x0a, 0x0a, 0x09, 0x63, 0x74,
- 0x74, 0x20, 0x3d, 0x20, 0x72, 0x7a, 0x20, 0x2f, 0x20, 0x73, 0x71, 0x72,
- 0x74, 0x28, 0x72, 0x78, 0x2a, 0x72, 0x78, 0x20, 0x2b, 0x20, 0x72, 0x79,
- 0x2a, 0x72, 0x79, 0x20, 0x2b, 0x20, 0x72, 0x7a, 0x2a, 0x72, 0x7a, 0x29,
- 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, 0x63, 0x6f, 0x73, 0x28, 0x32, 0x74,
- 0x68, 0x65, 0x74, 0x61, 0x29, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x74, 0x74,
- 0x20, 0x3d, 0x20, 0x61, 0x63, 0x6f, 0x73, 0x28, 0x63, 0x74, 0x74, 0x29,
- 0x3b, 0x0a, 0x09, 0x61, 0x7a, 0x20, 0x3d, 0x20, 0x61, 0x74, 0x61, 0x6e,
- 0x32, 0x28, 0x72, 0x79, 0x2c, 0x20, 0x72, 0x78, 0x29, 0x3b, 0x0a, 0x0a,
- 0x09, 0x6b, 0x78, 0x20, 0x3d, 0x20, 0x6b, 0x2a, 0x6e, 0x61, 0x74, 0x69,
- 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, 0x28, 0x74, 0x74, 0x29, 0x2a, 0x6e,
- 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x63, 0x6f, 0x73, 0x28, 0x61, 0x7a,
- 0x29, 0x3b, 0x0a, 0x09, 0x6b, 0x79, 0x20, 0x3d, 0x20, 0x6b, 0x2a, 0x6e,
- 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, 0x28, 0x74, 0x74,
- 0x29, 0x2a, 0x6e, 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e,
- 0x28, 0x61, 0x7a, 0x29, 0x3b, 0x0a, 0x09, 0x6b, 0x7a, 0x20, 0x3d, 0x20,
- 0x6b, 0x2a, 0x28, 0x63, 0x74, 0x74, 0x20, 0x2d, 0x20, 0x31, 0x2e, 0x30,
- 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x71, 0x20, 0x3d, 0x20, 0x28, 0x66, 0x6c,
- 0x6f, 0x61, 0x74, 0x34, 0x29, 0x28, 0x6b, 0x78, 0x2c, 0x20, 0x6b, 0x79,
- 0x2c, 0x20, 0x6b, 0x7a, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, 0x3b, 0x0a,
- 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x71, 0x3b, 0x0a,
- 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x6c, 0x61,
- 0x74, 0x74, 0x69, 0x63, 0x65, 0x5f, 0x66, 0x61, 0x63, 0x74, 0x6f, 0x72,
- 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, 0x65, 0x6c,
- 0x6c, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, 0x71, 0x2c,
- 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, 0x65,
- 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, 0x67,
- 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x61,
- 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72,
- 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61,
- 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f,
- 0x62, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x72, 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d,
- 0x61, 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63,
- 0x5f, 0x63, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74,
- 0x20, 0x66, 0x31, 0x2c, 0x20, 0x66, 0x32, 0x2c, 0x20, 0x66, 0x33, 0x2c,
- 0x20, 0x76, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20,
- 0x55, 0x64, 0x6f, 0x74, 0x71, 0x3b, 0x0a, 0x0a, 0x09, 0x55, 0x64, 0x6f,
- 0x74, 0x71, 0x2e, 0x78, 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
- 0x73, 0x30, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c,
- 0x6c, 0x2e, 0x73, 0x31, 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63,
- 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x32, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x0a,
- 0x09, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x79, 0x20, 0x3d, 0x20, 0x63,
- 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x33, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b,
- 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x34, 0x2a, 0x71, 0x2e, 0x79,
- 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x35, 0x2a, 0x71,
- 0x2e, 0x7a, 0x3b, 0x0a, 0x09, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x7a,
- 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x36, 0x2a, 0x71,
- 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x37,
- 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
- 0x73, 0x38, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a,
- 0x20, 0x4c, 0x6f, 0x6f, 0x6b, 0x20, 0x75, 0x70, 0x20, 0x76, 0x61, 0x6c,
- 0x75, 0x65, 0x73, 0x20, 0x66, 0x72, 0x6f, 0x6d, 0x20, 0x70, 0x72, 0x65,
- 0x63, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, 0x74, 0x65, 0x64, 0x20, 0x73,
- 0x69, 0x6e, 0x63, 0x28, 0x29, 0x20, 0x74, 0x61, 0x62, 0x6c, 0x65, 0x20,
- 0x2a, 0x2f, 0x0a, 0x09, 0x66, 0x31, 0x20, 0x3d, 0x20, 0x72, 0x65, 0x61,
- 0x64, 0x5f, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x66, 0x28, 0x66, 0x75, 0x6e,
- 0x63, 0x5f, 0x61, 0x2c, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72,
- 0x5f, 0x61, 0x2c, 0x20, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x32, 0x29,
- 0x28, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x78, 0x2c, 0x20, 0x30, 0x2e,
- 0x30, 0x29, 0x29, 0x2e, 0x73, 0x30, 0x3b, 0x0a, 0x09, 0x66, 0x32, 0x20,
- 0x3d, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, 0x69, 0x6d, 0x61, 0x67, 0x65,
- 0x66, 0x28, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, 0x2c, 0x20, 0x73, 0x61,
- 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x62, 0x2c, 0x20, 0x28, 0x66, 0x6c,
- 0x6f, 0x61, 0x74, 0x32, 0x29, 0x28, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e,
- 0x79, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, 0x29, 0x2e, 0x73, 0x30, 0x3b,
- 0x0a, 0x09, 0x66, 0x33, 0x20, 0x3d, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f,
- 0x69, 0x6d, 0x61, 0x67, 0x65, 0x66, 0x28, 0x66, 0x75, 0x6e, 0x63, 0x5f,
- 0x63, 0x2c, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x63,
- 0x2c, 0x20, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x32, 0x29, 0x28, 0x55,
- 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x7a, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29,
- 0x29, 0x2e, 0x73, 0x30, 0x3b, 0x0a, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75,
- 0x72, 0x6e, 0x20, 0x66, 0x31, 0x20, 0x2a, 0x20, 0x66, 0x32, 0x20, 0x2a,
- 0x20, 0x66, 0x33, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f,
- 0x61, 0x74, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x69, 0x6e,
- 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x67, 0x6c, 0x6f, 0x62,
- 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e,
- 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x0a, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x73, 0x69,
- 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, 0x20,
- 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6b,
- 0x2c, 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74,
- 0x20, 0x6c, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x69, 0x6e, 0x74, 0x20, 0x69,
- 0x64, 0x78, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x4f, 0x75, 0x74,
- 0x20, 0x6f, 0x66, 0x20, 0x72, 0x61, 0x6e, 0x67, 0x65, 0x3f, 0x20, 0x2a,
- 0x2f, 0x0a, 0x09, 0x69, 0x66, 0x20, 0x28, 0x20, 0x28, 0x61, 0x62, 0x73,
- 0x28, 0x68, 0x29, 0x20, 0x3e, 0x20, 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58,
- 0x29, 0x20, 0x7c, 0x7c, 0x20, 0x28, 0x61, 0x62, 0x73, 0x28, 0x6b, 0x29,
- 0x20, 0x3e, 0x20, 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x29, 0x20, 0x7c,
- 0x7c, 0x20, 0x28, 0x61, 0x62, 0x73, 0x28, 0x6c, 0x29, 0x20, 0x3e, 0x20,
- 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x29, 0x20, 0x29, 0x20, 0x7b, 0x0a,
- 0x09, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x30, 0x2e, 0x30,
- 0x3b, 0x0a, 0x09, 0x7d, 0x0a, 0x0a, 0x09, 0x68, 0x20, 0x3d, 0x20, 0x28,
- 0x68, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, 0x68, 0x20, 0x3a, 0x20,
- 0x68, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, 0x09, 0x6b, 0x20, 0x3d,
- 0x20, 0x28, 0x6b, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, 0x6b, 0x20,
- 0x3a, 0x20, 0x6b, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, 0x09, 0x6c,
- 0x20, 0x3d, 0x20, 0x28, 0x6c, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20,
- 0x6c, 0x20, 0x3a, 0x20, 0x6c, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a,
- 0x0a, 0x09, 0x69, 0x64, 0x78, 0x20, 0x3d, 0x20, 0x68, 0x20, 0x2b, 0x20,
- 0x28, 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x6b, 0x29, 0x20, 0x2b, 0x20, 0x28,
- 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x6c, 0x29,
- 0x3b, 0x0a, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x69,
- 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x5b, 0x69,
- 0x64, 0x78, 0x5d, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f,
- 0x61, 0x74, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x66, 0x6c,
- 0x61, 0x67, 0x67, 0x65, 0x64, 0x5f, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73,
- 0x69, 0x74, 0x79, 0x28, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66,
- 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73,
- 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61,
- 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, 0x61,
- 0x67, 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x73,
- 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c,
- 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20,
- 0x6b, 0x2c, 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e,
- 0x74, 0x20, 0x6c, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75,
- 0x72, 0x6e, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x69, 0x6e,
- 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x69, 0x6e, 0x74, 0x65,
- 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x68, 0x2c, 0x20,
- 0x6b, 0x2c, 0x20, 0x6c, 0x29, 0x0a, 0x09, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x2a, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f,
- 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x66, 0x6c,
- 0x61, 0x67, 0x73, 0x2c, 0x20, 0x68, 0x2c, 0x20, 0x6b, 0x2c, 0x20, 0x6c,
- 0x29, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, 0x74,
- 0x20, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x5f, 0x66, 0x61,
- 0x63, 0x74, 0x6f, 0x72, 0x28, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20,
- 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e,
- 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62,
- 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c,
- 0x61, 0x67, 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63,
- 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20,
- 0x71, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
- 0x68, 0x66, 0x2c, 0x20, 0x6b, 0x66, 0x2c, 0x20, 0x6c, 0x66, 0x3b, 0x0a,
- 0x09, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, 0x20, 0x6b, 0x2c, 0x20, 0x6c,
- 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x76, 0x61, 0x6c,
- 0x20, 0x3d, 0x20, 0x30, 0x2e, 0x30, 0x3b, 0x0a, 0x0a, 0x09, 0x23, 0x69,
- 0x66, 0x64, 0x65, 0x66, 0x20, 0x46, 0x4c, 0x41, 0x54, 0x5f, 0x49, 0x4e,
- 0x54, 0x45, 0x4e, 0x53, 0x49, 0x54, 0x49, 0x45, 0x53, 0x0a, 0x09, 0x72,
- 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x31, 0x30, 0x30, 0x2e, 0x30, 0x3b,
- 0x0a, 0x09, 0x23, 0x65, 0x6c, 0x73, 0x65, 0x0a, 0x0a, 0x09, 0x68, 0x66,
- 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x30, 0x2a, 0x71,
- 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x31,
- 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
- 0x73, 0x32, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20,
- 0x68, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x6b, 0x66, 0x20, 0x3d, 0x20, 0x63,
- 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x33, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b,
- 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x34, 0x2a, 0x71, 0x2e, 0x79,
- 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x35, 0x2a, 0x71,
- 0x2e, 0x7a, 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, 0x6b, 0x20, 0x2a, 0x2f,
- 0x0a, 0x09, 0x6c, 0x66, 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
- 0x73, 0x36, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c,
- 0x6c, 0x2e, 0x73, 0x37, 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63,
- 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x38, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x20,
- 0x20, 0x2f, 0x2a, 0x20, 0x6c, 0x20, 0x2a, 0x2f, 0x0a, 0x0a, 0x09, 0x68,
- 0x20, 0x3d, 0x20, 0x72, 0x6f, 0x75, 0x6e, 0x64, 0x28, 0x68, 0x66, 0x29,
- 0x3b, 0x0a, 0x09, 0x6b, 0x20, 0x3d, 0x20, 0x72, 0x6f, 0x75, 0x6e, 0x64,
- 0x28, 0x6b, 0x66, 0x29, 0x3b, 0x0a, 0x09, 0x6c, 0x20, 0x3d, 0x20, 0x72,
- 0x6f, 0x75, 0x6e, 0x64, 0x28, 0x6c, 0x66, 0x29, 0x3b, 0x0a, 0x0a, 0x09,
- 0x2f, 0x2a, 0x20, 0x53, 0x79, 0x6d, 0x6d, 0x65, 0x74, 0x72, 0x79, 0x20,
- 0x73, 0x74, 0x75, 0x66, 0x66, 0x20, 0x67, 0x6f, 0x65, 0x73, 0x20, 0x68,
- 0x65, 0x72, 0x65, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x49, 0x4e, 0x53, 0x45,
- 0x52, 0x54, 0x5f, 0x48, 0x45, 0x52, 0x45, 0x0a, 0x0a, 0x09, 0x72, 0x65,
- 0x74, 0x75, 0x72, 0x6e, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x09, 0x23,
- 0x65, 0x6e, 0x64, 0x69, 0x66, 0x20, 0x2f, 0x2a, 0x20, 0x46, 0x4c, 0x41,
- 0x54, 0x5f, 0x49, 0x4e, 0x54, 0x45, 0x4e, 0x53, 0x49, 0x54, 0x49, 0x49,
- 0x45, 0x53, 0x20, 0x2a, 0x2f, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x6b, 0x65,
- 0x72, 0x6e, 0x65, 0x6c, 0x20, 0x76, 0x6f, 0x69, 0x64, 0x20, 0x64, 0x69,
- 0x66, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x28, 0x67, 0x6c,
- 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a,
- 0x64, 0x69, 0x66, 0x66, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
- 0x6b, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x77, 0x65, 0x69,
- 0x67, 0x68, 0x74, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x77, 0x2c, 0x20,
- 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72,
- 0x5f, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f,
- 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61,
- 0x74, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74,
- 0x20, 0x66, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
- 0x66, 0x73, 0x7a, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73,
- 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, 0x73,
- 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, 0x73, 0x7a,
- 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x72, 0x65, 0x73, 0x2c,
- 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6c, 0x65, 0x6e, 0x2c,
- 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, 0x65, 0x6c,
- 0x6c, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, 0x6c,
- 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69,
- 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c,
- 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, 0x61, 0x67,
- 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79,
- 0x20, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66,
- 0x75, 0x6e, 0x63, 0x5f, 0x61, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f,
- 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x32, 0x64,
- 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, 0x2c, 0x0a, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72,
- 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61,
- 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f,
- 0x63, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x6c, 0x6f, 0x63, 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f,
- 0x61, 0x74, 0x20, 0x2a, 0x74, 0x6d, 0x70, 0x29, 0x0a, 0x7b, 0x0a, 0x09,
- 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x2c, 0x20, 0x73, 0x73,
- 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x5f, 0x6c,
- 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x2c, 0x20, 0x49, 0x5f, 0x6c, 0x61,
- 0x74, 0x74, 0x69, 0x63, 0x65, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61,
- 0x74, 0x20, 0x49, 0x5f, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65,
- 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, 0x71, 0x3b,
- 0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20,
- 0x6c, 0x73, 0x30, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f,
- 0x63, 0x61, 0x6c, 0x5f, 0x73, 0x69, 0x7a, 0x65, 0x28, 0x30, 0x29, 0x3b,
- 0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20,
- 0x6c, 0x73, 0x31, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f,
- 0x63, 0x61, 0x6c, 0x5f, 0x73, 0x69, 0x7a, 0x65, 0x28, 0x31, 0x29, 0x3b,
- 0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20,
- 0x6c, 0x69, 0x30, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f,
- 0x63, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28, 0x30, 0x29, 0x3b, 0x0a, 0x09,
- 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6c, 0x69,
- 0x31, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f, 0x63, 0x61,
- 0x6c, 0x5f, 0x69, 0x64, 0x28, 0x31, 0x29, 0x3b, 0x0a, 0x09, 0x63, 0x6f,
- 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6c, 0x73, 0x20, 0x3d,
- 0x20, 0x6c, 0x73, 0x30, 0x20, 0x2a, 0x20, 0x6c, 0x73, 0x31, 0x3b, 0x0a,
- 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61,
- 0x74, 0x65, 0x20, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x61,
- 0x6c, 0x20, 0x63, 0x6f, 0x6f, 0x72, 0x64, 0x69, 0x6e, 0x61, 0x74, 0x65,
- 0x73, 0x20, 0x69, 0x6e, 0x20, 0x66, 0x73, 0x2f, 0x73, 0x73, 0x20, 0x2a,
- 0x2f, 0x0a, 0x09, 0x66, 0x73, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76,
- 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x67, 0x65,
- 0x74, 0x5f, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28,
- 0x30, 0x29, 0x29, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, 0x72,
- 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, 0x73, 0x30, 0x29,
- 0x3b, 0x0a, 0x09, 0x73, 0x73, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76,
- 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x67, 0x65,
- 0x74, 0x5f, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28,
- 0x31, 0x29, 0x29, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, 0x72,
- 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, 0x73, 0x31, 0x29,
- 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x47, 0x65, 0x74, 0x20, 0x74,
- 0x68, 0x65, 0x20, 0x73, 0x63, 0x61, 0x74, 0x74, 0x65, 0x72, 0x69, 0x6e,
- 0x67, 0x20, 0x76, 0x65, 0x63, 0x74, 0x6f, 0x72, 0x20, 0x2a, 0x2f, 0x0a,
- 0x09, 0x71, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x71, 0x28, 0x66,
- 0x73, 0x2c, 0x20, 0x73, 0x73, 0x2c, 0x20, 0x72, 0x65, 0x73, 0x2c, 0x20,
- 0x63, 0x6c, 0x65, 0x6e, 0x2c, 0x20, 0x6b, 0x2c, 0x0a, 0x09, 0x20, 0x20,
- 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x63, 0x6f, 0x72, 0x6e,
- 0x65, 0x72, 0x5f, 0x78, 0x2c, 0x20, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72,
- 0x5f, 0x79, 0x2c, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x73, 0x79,
- 0x2c, 0x20, 0x66, 0x73, 0x7a, 0x2c, 0x20, 0x73, 0x73, 0x78, 0x2c, 0x20,
- 0x73, 0x73, 0x79, 0x2c, 0x20, 0x73, 0x73, 0x7a, 0x29, 0x3b, 0x0a, 0x0a,
- 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, 0x74,
- 0x65, 0x20, 0x74, 0x68, 0x65, 0x20, 0x64, 0x69, 0x66, 0x66, 0x72, 0x61,
- 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x66, 0x5f,
- 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, 0x3d, 0x20, 0x6c, 0x61,
- 0x74, 0x74, 0x69, 0x63, 0x65, 0x5f, 0x66, 0x61, 0x63, 0x74, 0x6f, 0x72,
- 0x28, 0x63, 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x71, 0x2c, 0x20, 0x66, 0x75,
- 0x6e, 0x63, 0x5f, 0x61, 0x2c, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62,
- 0x2c, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x63, 0x29, 0x3b, 0x0a, 0x09,
- 0x49, 0x5f, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x20, 0x3d,
- 0x20, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x5f, 0x66, 0x61,
- 0x63, 0x74, 0x6f, 0x72, 0x28, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69,
- 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x61, 0x67, 0x73, 0x2c,
- 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x71, 0x29, 0x3b, 0x0a, 0x09,
- 0x49, 0x5f, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, 0x3d, 0x20,
- 0x70, 0x6f, 0x77, 0x28, 0x66, 0x5f, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63,
- 0x65, 0x2c, 0x20, 0x32, 0x2e, 0x30, 0x66, 0x29, 0x3b, 0x0a, 0x0a, 0x09,
- 0x74, 0x6d, 0x70, 0x5b, 0x6c, 0x69, 0x30, 0x20, 0x2b, 0x20, 0x6c, 0x73,
- 0x30, 0x2a, 0x6c, 0x69, 0x31, 0x5d, 0x20, 0x3d, 0x20, 0x49, 0x5f, 0x6d,
- 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x20, 0x2a, 0x20, 0x49, 0x5f,
- 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x3b, 0x0a, 0x0a, 0x09, 0x62,
- 0x61, 0x72, 0x72, 0x69, 0x65, 0x72, 0x28, 0x43, 0x4c, 0x4b, 0x5f, 0x4c,
- 0x4f, 0x43, 0x41, 0x4c, 0x5f, 0x4d, 0x45, 0x4d, 0x5f, 0x46, 0x45, 0x4e,
- 0x43, 0x45, 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x46, 0x69,
- 0x72, 0x73, 0x74, 0x20, 0x74, 0x68, 0x72, 0x65, 0x61, 0x64, 0x20, 0x69,
- 0x6e, 0x20, 0x67, 0x72, 0x6f, 0x75, 0x70, 0x20, 0x73, 0x75, 0x6d, 0x73,
- 0x20, 0x74, 0x68, 0x65, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x73,
- 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x69, 0x66, 0x20, 0x28, 0x20, 0x6c, 0x69,
- 0x30, 0x20, 0x2b, 0x20, 0x6c, 0x69, 0x31, 0x20, 0x3d, 0x3d, 0x20, 0x30,
- 0x20, 0x29, 0x20, 0x7b, 0x0a, 0x0a, 0x09, 0x09, 0x69, 0x6e, 0x74, 0x20,
- 0x69, 0x3b, 0x0a, 0x09, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73,
- 0x75, 0x6d, 0x20, 0x3d, 0x20, 0x30, 0x2e, 0x30, 0x3b, 0x0a, 0x09, 0x09,
- 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x09,
- 0x09, 0x69, 0x6e, 0x74, 0x20, 0x69, 0x64, 0x78, 0x3b, 0x0a, 0x0a, 0x09,
- 0x09, 0x69, 0x64, 0x78, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65,
- 0x72, 0x74, 0x5f, 0x69, 0x6e, 0x74, 0x5f, 0x72, 0x74, 0x7a, 0x28, 0x66,
- 0x73, 0x29, 0x20, 0x2b, 0x20, 0x77, 0x2a, 0x63, 0x6f, 0x6e, 0x76, 0x65,
- 0x72, 0x74, 0x5f, 0x69, 0x6e, 0x74, 0x5f, 0x72, 0x74, 0x7a, 0x28, 0x73,
- 0x73, 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x09, 0x66, 0x6f, 0x72, 0x20, 0x28,
- 0x20, 0x69, 0x3d, 0x30, 0x3b, 0x20, 0x69, 0x3c, 0x6c, 0x73, 0x3b, 0x20,
- 0x69, 0x2b, 0x2b, 0x20, 0x29, 0x20, 0x73, 0x75, 0x6d, 0x20, 0x2b, 0x3d,
- 0x20, 0x74, 0x6d, 0x70, 0x5b, 0x69, 0x5d, 0x3b, 0x0a, 0x0a, 0x09, 0x09,
- 0x76, 0x61, 0x6c, 0x20, 0x3d, 0x20, 0x77, 0x65, 0x69, 0x67, 0x68, 0x74,
- 0x20, 0x2a, 0x20, 0x73, 0x75, 0x6d, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e,
- 0x76, 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c,
- 0x73, 0x29, 0x3b, 0x0a, 0x09, 0x09, 0x64, 0x69, 0x66, 0x66, 0x5b, 0x69,
- 0x64, 0x78, 0x5d, 0x20, 0x3d, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x0a,
- 0x09, 0x7d, 0x0a, 0x0a, 0x7d, 0x0a
-};
-unsigned int data_diffraction_cl_len = 6006;
diff --git a/src/diffraction.h b/src/diffraction.h
deleted file mode 100644
index 63b022c9..00000000
--- a/src/diffraction.h
+++ /dev/null
@@ -1,66 +0,0 @@
-/*
- * diffraction.h
- *
- * Calculate diffraction patterns by Fourier methods
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2009-2019 Thomas White <taw@physics.org>
- * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- * 2013 Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DIFFRACTION_H
-#define DIFFRACTION_H
-
-#include <gsl/gsl_rng.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,
- int no_fringes, int flat, int n_samples);
-
-extern struct sample *generate_tophat(struct image *image);
-
-extern struct sample *generate_SASE(struct image *image, gsl_rng *rng);
-
-extern struct sample *generate_twocolour(struct image *image);
-
-extern struct sample *generate_spectrum_fromfile(struct image *image,
- char *spectrum_fn);
-
-#endif /* DIFFRACTION_H */
diff --git a/src/list_tmp.h b/src/list_tmp.h
deleted file mode 100644
index 0c1a84e4..00000000
--- a/src/list_tmp.h
+++ /dev/null
@@ -1,104 +0,0 @@
-/*
- * 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_<LABEL>(), for creating a new list of the given type,
- * - set_<LABEL>(), for setting a value in a list,
- * - lookup_<LABEL>(), for retrieving values from a list,
- * ..and so on.
- *
- * See src/utils.h for more information.
- *
- */
-
-/*
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2009-2012 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#include <stdio.h>
-
-#define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
-
-
-static inline void LABEL(set_arr)(TYPE *ref, signed int h,
- signed int k, signed int l,
- TYPE i)
-{
- int idx;
-
- if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
- ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
- ERROR_T("You need to re-configure INDMAX and re-run.\n");
- exit(1);
- }
-
- if ( h < 0 ) h += IDIM;
- if ( k < 0 ) k += IDIM;
- if ( l < 0 ) l += IDIM;
-
- idx = h + (IDIM*k) + (IDIM*IDIM*l);
- ref[idx] = i;
-}
-
-
-static inline TYPE LABEL(lookup_arr)(const TYPE *ref, signed int h,
- signed int k, signed int l)
-{
- int idx;
-
- if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
- ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
- ERROR_T("You need to re-configure INDMAX and re-run.\n");
- exit(1);
- }
-
- if ( h < 0 ) h += IDIM;
- if ( k < 0 ) k += IDIM;
- if ( l < 0 ) l += IDIM;
-
- idx = h + (IDIM*k) + (IDIM*IDIM*l);
- return ref[idx];
-}
-
-
-static inline TYPE *LABEL(new_arr)(void)
-{
- TYPE *r;
- size_t r_size;
- r_size = IDIM*IDIM*IDIM*sizeof(TYPE);
- r = malloc(r_size);
- memset(r, 0, r_size);
- return r;
-}
-
-
-#undef LABEL
-#undef TYPE
-#undef ERROR_T
diff --git a/src/partial_sim.c b/src/partial_sim.c
deleted file mode 100644
index f5588ceb..00000000
--- a/src/partial_sim.c
+++ /dev/null
@@ -1,1034 +0,0 @@
-/*
- * partial_sim.c
- *
- * Generate partials for testing scaling
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2011-2020 Thomas White <taw@physics.org>
- * 2014 Valerio Mariani
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <getopt.h>
-#include <assert.h>
-#include <pthread.h>
-#include <gsl/gsl_rng.h>
-
-#include <image.h>
-#include <utils.h>
-#include <reflist-utils.h>
-#include <symmetry.h>
-#include <geometry.h>
-#include <stream.h>
-#include <thread-pool.h>
-#include <cell-utils.h>
-
-#include "version.h"
-
-
-/* Number of bins for partiality graph */
-#define NBINS 50
-
-
-static void mess_up_cell(Crystal *cr, double cnoise, gsl_rng *rng)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- UnitCell *cell = crystal_get_cell(cr);
-
- //STATUS("Real:\n");
- //cell_print(cell);
-
- cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- ax = flat_noise(rng, ax, cnoise*fabs(ax)/100.0);
- ay = flat_noise(rng, ay, cnoise*fabs(ay)/100.0);
- az = flat_noise(rng, az, cnoise*fabs(az)/100.0);
- bx = flat_noise(rng, bx, cnoise*fabs(bx)/100.0);
- by = flat_noise(rng, by, cnoise*fabs(by)/100.0);
- bz = flat_noise(rng, bz, cnoise*fabs(bz)/100.0);
- cx = flat_noise(rng, cx, cnoise*fabs(cx)/100.0);
- cy = flat_noise(rng, cy, cnoise*fabs(cy)/100.0);
- cz = flat_noise(rng, cz, cnoise*fabs(cz)/100.0);
- cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
-
- //STATUS("Changed:\n");
- //cell_print(cell);
-}
-
-
-/* For each reflection in "partial", fill in what the intensity would be
- * according to "full" */
-static void calculate_partials(Crystal *cr,
- RefList *full, const SymOpList *sym,
- int random_intensities,
- pthread_rwlock_t *full_lock,
- unsigned long int *n_ref, double *p_hist,
- double *p_max, double max_q, double full_stddev,
- double noise_stddev, gsl_rng *rng,
- UnitCell *template_cell, RefList *template_reflist)
-{
- Reflection *refl;
- RefListIterator *iter;
-
- for ( refl = first_refl(crystal_get_reflections(cr), &iter);
- refl != NULL;
- refl = next_refl(refl, iter) )
- {
- signed int h, k, l;
- Reflection *rfull;
- double L, p, Ip, If;
- int bin;
- double res;
-
- get_indices(refl, &h, &k, &l);
- get_asymm(sym, h, k, l, &h, &k, &l);
- p = get_partiality(refl);
- L = get_lorentz(refl);
-
- pthread_rwlock_rdlock(full_lock);
- rfull = find_refl(full, h, k, l);
- pthread_rwlock_unlock(full_lock);
-
- if ( rfull == NULL ) {
- if ( random_intensities ) {
-
- pthread_rwlock_wrlock(full_lock);
-
- /* In the gap between the unlock and the wrlock,
- * the reflection might have been created by
- * another thread. So, we must check again */
- rfull = find_refl(full, h, k, l);
- if ( rfull == NULL ) {
- rfull = add_refl(full, h, k, l);
- If = fabs(gaussian_noise(rng, 0.0,
- full_stddev));
- set_intensity(rfull, If);
- set_redundancy(rfull, 1);
- } else {
- If = get_intensity(rfull);
- }
- pthread_rwlock_unlock(full_lock);
-
- } else {
- set_redundancy(refl, 0);
- If = 0.0;
- }
- } else {
- If = get_intensity(rfull);
- if ( random_intensities ) {
- lock_reflection(rfull);
- int red = get_redundancy(rfull);
- set_redundancy(rfull, red+1);
- unlock_reflection(rfull);
- }
- }
-
- Ip = crystal_get_osf(cr) * L * p * If;
-
- res = resolution(crystal_get_cell(cr), h, k, l);
- bin = NBINS*2.0*res/max_q;
- if ( (bin < NBINS) && (bin>=0) ) {
- p_hist[bin] += p;
- n_ref[bin]++;
- if ( p > p_max[bin] ) p_max[bin] = p;
- } else {
- STATUS("Reflection out of histogram range: %e %i %f\n",
- res, bin, p);
- }
-
- Ip = gaussian_noise(rng, Ip, noise_stddev);
-
- set_intensity(refl, Ip);
- set_esd_intensity(refl, noise_stddev);
- }
-}
-
-
-static void draw_and_write_image(struct image *image,
- const DataTemplate *dtempl,
- RefList *reflections,
- gsl_rng *rng, double background)
-{
- Reflection *refl;
- RefListIterator *iter;
- int i;
-
- image->dp = malloc(image->detgeom->n_panels*sizeof(float *));
- if ( image->dp == NULL ) {
- ERROR("Failed to allocate data\n");
- return;
- }
- for ( i=0; i<image->detgeom->n_panels; i++ ) {
-
- int j;
- struct detgeom_panel *p = &image->detgeom->panels[i];
-
- image->dp[i] = calloc(p->w * p->h, sizeof(float));
- if ( image->dp[i] == NULL ) {
- ERROR("Failed to allocate data\n");
- return;
- }
- for ( j=0; j<p->w*p->h; j++ ) {
- image->dp[i][j] = poisson_noise(rng, background);
- }
-
- }
-
- for ( refl = first_refl(reflections, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) )
- {
- double Ip;
- double dfs, dss;
- int fs, ss;
- struct detgeom_panel *p;
- signed int pn;
-
- Ip = get_intensity(refl);
-
- get_detector_pos(refl, &dfs, &dss);
- pn = get_panel_number(refl);
- assert(pn < image->detgeom->n_panels);
- p = &image->detgeom->panels[pn];
-
- /* Explicit rounding, downwards */
- fs = dfs; ss = dss;
- assert(fs >= 0);
- assert(ss >= 0);
- assert(fs < p->w);
- assert(ss < p->h);
-
- image->dp[pn][fs + p->w*ss] += Ip;
- }
-
- image_write(image, dtempl, image->filename);
-
- for ( i=0; i<image->detgeom->n_panels; i++ ) {
- free(image->dp[i]);
- }
- free(image->dp);
-}
-
-
-static void show_help(const char *s)
-{
- printf("Syntax: %s [options]\n\n", s);
- printf(
-"Generate a stream containing partials from a reflection list.\n"
-"\n"
-" -h, --help Display this help message.\n"
-" --version Print CrystFEL version number and exit.\n"
-"\n"
-"You need to provide the following basic options:\n"
-" -i, --input=<file> Read reflections from <file>.\n"
-" Default: generate random ones instead (see -r).\n"
-" -o, --output=<file> Write partials in stream format to <file>.\n"
-" --images=<prefix> Write images to <prefix>NNN.h5.\n"
-" -g. --geometry=<file> Get detector geometry from file.\n"
-" -p, --pdb=<file> PDB file from which to get the unit cell.\n"
-"\n"
-" -y, --symmetry=<sym> Symmetry of the input reflection list.\n"
-" -n <n> Simulate <n> patterns. Default: 2.\n"
-" -r, --save-random=<file> Save randomly generated intensities to file.\n"
-" --pgraph=<file> Save a histogram of partiality values to file.\n"
-" -c, --cnoise=<val> Amount of reciprocal space cell noise, in percent.\n"
-" --osf-stddev=<val> Standard deviation of the scaling factors.\n"
-" --full-stddev=<val> Standard deviation of the randomly\n"
-" generated full intensities, if not using -i.\n"
-" --noise-stddev=<val> Set the standard deviation of the noise.\n"
-" --background=<val> Background level in photons. Default 3000.\n"
-" --profile-radius Reciprocal space reflection profile radius in m^-1.\n"
-" Default 0.001e9 m^-1\n"
-" --really-random Be non-deterministic.\n"
-"\n"
-);
-}
-
-
-struct partial_sim_queue_args
-{
- RefList *full;
- pthread_rwlock_t full_lock;
- const DataTemplate *dtempl;
-
- int n_done;
- int n_started;
- int n_to_do;
-
- SymOpList *sym;
- int random_intensities;
- UnitCell *cell;
- double cnoise;
- double osf_stddev;
- double full_stddev;
- double noise_stddev;
- double background;
- double profile_radius;
-
- double max_q;
-
- char *image_prefix;
-
- /* The overall histogram */
- double p_hist[NBINS];
- unsigned long int n_ref[NBINS];
- double p_max[NBINS];
-
- Stream *stream;
- gsl_rng **rngs;
- Stream *template_stream;
-};
-
-
-struct partial_sim_worker_args
-{
- struct partial_sim_queue_args *qargs;
- struct image *image;
-
- UnitCell *template_cell;
- RefList *template_reflist;
-
- /* Histogram for this image */
- double p_hist[NBINS];
- unsigned long int n_ref[NBINS];
- double p_max[NBINS];
-
- int n;
-};
-
-
-static void *create_job(void *vqargs)
-{
- struct partial_sim_worker_args *wargs;
- struct partial_sim_queue_args *qargs = vqargs;
-
- /* All done already? */
- if ( qargs->n_started == qargs->n_to_do ) return NULL;
-
- wargs = malloc(sizeof(struct partial_sim_worker_args));
-
- wargs->qargs = qargs;
-
- if ( qargs->template_stream != NULL ) {
-
- struct image *image;
-
- image = stream_read_chunk(qargs->template_stream,
- STREAM_REFLECTIONS);
- if ( image == NULL ) {
- ERROR("Failed to read template chunk!\n");
- free(wargs);
- return NULL;
- }
- if ( image->n_crystals != 1 ) {
- ERROR("Template stream must have exactly one crystal "
- "per frame.\n");
- free(wargs);
- return NULL;
- }
-
- wargs->template_cell = crystal_get_cell(image->crystals[0]);
- wargs->template_reflist = crystal_get_reflections(image->crystals[0]);
-
- image_free(image);
-
- } else {
- wargs->template_cell = NULL;
- wargs->template_reflist = NULL;
- }
-
- qargs->n_started++;
- wargs->n = qargs->n_started;
-
- return wargs;
-}
-
-
-static void run_job(void *vwargs, int cookie)
-{
- struct partial_sim_worker_args *wargs = vwargs;
- struct partial_sim_queue_args *qargs = wargs->qargs;
- int i;
- Crystal *cr;
- double osf;
- struct image *image;
-
- image = image_create_for_simulation(qargs->dtempl);
- if ( image == NULL ) {
- ERROR("Failed to create image.\n");
- return;
- }
-
- cr = crystal_new();
- if ( cr == NULL ) {
- ERROR("Failed to create crystal.\n");
- return;
- }
- crystal_set_image(cr, image);
- image_add_crystal(image, cr);
-
- do {
- osf = gaussian_noise(qargs->rngs[cookie], 1.0,
- qargs->osf_stddev);
- } while ( osf <= 0.0 );
- crystal_set_osf(cr, osf);
- crystal_set_mosaicity(cr, 0.0);
- crystal_set_profile_radius(cr, qargs->profile_radius);
-
- if ( wargs->template_cell == NULL ) {
- /* Set up a random orientation */
- struct quaternion orientation;
- char tmp[128];
- orientation = random_quaternion(qargs->rngs[cookie]);
- crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
- snprintf(tmp, 127, "quaternion = %f %f %f %f",
- orientation.w, orientation.x, orientation.y, orientation.z);
- crystal_add_notes(cr, tmp);
- } else {
- crystal_set_cell(cr, wargs->template_cell);
- }
-
- image->filename = malloc(256);
- if ( image->filename == NULL ) {
- ERROR("Failed to allocate filename.\n");
- return;
- }
- if ( qargs->image_prefix != NULL ) {
- snprintf(image->filename, 255, "%s%i.h5",
- qargs->image_prefix, wargs->n);
- } else {
- snprintf(image->filename, 255, "dummy-%i.h5", wargs->n);
- }
-
- if ( wargs->template_reflist == NULL ) {
- RefList *reflections;
- reflections = predict_to_res(cr, qargs->max_q);
- crystal_set_reflections(cr, reflections);
- calculate_partialities(cr, PMODEL_XSPHERE);
- } else {
- crystal_set_reflections(cr, wargs->template_reflist);
- update_predictions(cr);
- calculate_partialities(cr, PMODEL_XSPHERE);
- }
-
- for ( i=0; i<NBINS; i++ ) {
- wargs->n_ref[i] = 0;
- wargs->p_hist[i] = 0.0;
- wargs->p_max[i] = 0.0;
- }
-
- calculate_partials(cr, qargs->full,
- qargs->sym, qargs->random_intensities,
- &qargs->full_lock,
- wargs->n_ref, wargs->p_hist, wargs->p_max,
- qargs->max_q, qargs->full_stddev,
- qargs->noise_stddev, qargs->rngs[cookie],
- wargs->template_cell, wargs->template_reflist);
-
- if ( qargs->image_prefix != NULL ) {
- draw_and_write_image(image, qargs->dtempl,
- crystal_get_reflections(cr),
- qargs->rngs[cookie],
- qargs->background);
- }
-
- /* Give a slightly incorrect cell in the stream */
- mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]);
-
- wargs->image = image;
-}
-
-
-static void finalise_job(void *vqargs, void *vwargs)
-{
- struct partial_sim_worker_args *wargs = vwargs;
- struct partial_sim_queue_args *qargs = vqargs;
- int i;
- int ret;
-
- ret = stream_write_chunk(qargs->stream, wargs->image,
- STREAM_REFLECTIONS);
- if ( ret != 0 ) {
- ERROR("WARNING: error writing stream file.\n");
- }
-
- for ( i=0; i<NBINS; i++ ) {
- qargs->n_ref[i] += wargs->n_ref[i];
- qargs->p_hist[i] += wargs->p_hist[i];
- if ( wargs->p_max[i] > qargs->p_max[i] ) {
- qargs->p_max[i] = wargs->p_max[i];
- }
- }
-
- qargs->n_done++;
- progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
-
- image_free(wargs->image);
-}
-
-
-int main(int argc, char *argv[])
-{
- int c;
- char *input_file = NULL;
- char *output_file = NULL;
- char *geomfile = NULL;
- char *cellfile = NULL;
- DataTemplate *dtempl;
- RefList *full = NULL;
- char *sym_str = NULL;
- SymOpList *sym;
- UnitCell *cell = NULL;
- Stream *stream;
- int n = 2;
- int random_intensities = 0;
- char *save_file = NULL;
- struct partial_sim_queue_args qargs;
- int n_threads = 1;
- char *rval;
- int i;
- char *phist_file = NULL;
- int config_random = 0;
- char *image_prefix = NULL;
- Stream *template_stream = NULL;
- char *template = NULL;
- struct image *test_image;
-
- /* Default simulation parameters */
- double profile_radius = 0.001e9;
- double osf_stddev = 2.0;
- double full_stddev = 1000.0;
- double noise_stddev = 20.0;
- double background = 3000.0;
- double cnoise = 0.0;
-
- /* Long options */
- const struct option longopts[] = {
- {"help", 0, NULL, 'h'},
- {"version", 0, NULL, 'v'},
- {"beam", 1, NULL, 'b'},
- {"output", 1, NULL, 'o'},
- {"input", 1, NULL, 'i'},
- {"pdb", 1, NULL, 'p'},
- {"geometry", 1, NULL, 'g'},
- {"symmetry", 1, NULL, 'y'},
- {"save-random", 1, NULL, 'r'},
- {"cnoise", 1, NULL, 'c'},
-
- {"pgraph", 1, NULL, 2},
- {"osf-stddev", 1, NULL, 3},
- {"full-stddev", 1, NULL, 4},
- {"noise-stddev", 1, NULL, 5},
- {"images", 1, NULL, 6},
- {"background", 1, NULL, 7},
- {"beam-divergence", 1, NULL, 8},
- {"beam-bandwidth", 1, NULL, 9},
- {"profile-radius", 1, NULL, 10},
- {"photon-energy", 1, NULL, 11},
- {"template-stream", 1, NULL, 12},
-
- {"really-random", 0, &config_random, 1},
-
- {0, 0, NULL, 0}
- };
-
- /* Short options */
- while ((c = getopt_long(argc, argv, "hi:o:p:g:y:n:r:j:c:vb:",
- longopts, NULL)) != -1)
- {
- switch (c) {
-
- case 'h' :
- show_help(argv[0]);
- return 0;
-
- case 'v' :
- printf("CrystFEL: %s\n",
- crystfel_version_string());
- printf("%s\n",
- crystfel_licence_string());
- return 0;
-
- case 'b' :
- ERROR("WARNING: This version of CrystFEL no longer "
- "uses beam files. Please remove the beam file "
- "from your partial_sim command line.\n");
- return 1;
-
- case 'i' :
- input_file = strdup(optarg);
- break;
-
- case 'o' :
- output_file = strdup(optarg);
- break;
-
- case 'p' :
- cellfile = strdup(optarg);
- break;
-
- case 'g' :
- geomfile = strdup(optarg);
- break;
-
- case 'y' :
- sym_str = strdup(optarg);
- break;
-
- case 'n' :
- n = atoi(optarg);
- break;
-
- case 'r' :
- save_file = strdup(optarg);
- break;
-
- case 'j' :
- n_threads = atoi(optarg);
- break;
-
- case 'c' :
- cnoise = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid cell noise value.\n");
- return 1;
- }
- break;
-
- case 2 :
- phist_file = strdup(optarg);
- break;
-
- case 3 :
- osf_stddev = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid OSF standard deviation.\n");
- return 1;
- }
- if ( osf_stddev < 0.0 ) {
- ERROR("Invalid OSF standard deviation.");
- ERROR(" (must be positive).\n");
- return 1;
- }
- break;
-
- case 4 :
- full_stddev = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid full standard deviation.\n");
- return 1;
- }
- if ( full_stddev < 0.0 ) {
- ERROR("Invalid full standard deviation.");
- ERROR(" (must be positive).\n");
- return 1;
- }
- break;
-
- case 5 :
- noise_stddev = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid noise standard deviation.\n");
- return 1;
- }
- if ( noise_stddev < 0.0 ) {
- ERROR("Invalid noise standard deviation.");
- ERROR(" (must be positive).\n");
- return 1;
- }
- break;
-
- case 6 :
- image_prefix = strdup(optarg);
- break;
-
- case 7 :
- background = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid background level.\n");
- return 1;
- }
- if ( background < 0.0 ) {
- ERROR("Background level must be positive.\n");
- return 1;
- }
- break;
-
- case 8 :
- ERROR("--beam-divergence is no longer used.\n");
- ERROR("The 'xsphere' partiality model does not take divergence into account.\n");
- return 1;
-
- case 9 :
- ERROR("--beam-bandwidth is no longer used.\n");
- ERROR("Set the bandwidth in the geometry file instead.\n");
- return 1;
-
- case 10 :
- profile_radius = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid profile radius.\n");
- return 1;
- }
- if ( profile_radius < 0.0 ) {
- ERROR("Profile radius must be positive.\n");
- return 1;
- }
- break;
-
- case 11 :
- ERROR("--photon-energy is no longer used.\n");
- ERROR("Set the photon energy in the geometry file instead.\n");
- return 1;
-
- case 12 :
- template = strdup(optarg);
- break;
-
- case 0 :
- break;
-
- case '?' :
- return 1;
-
- default :
- ERROR("Unhandled option '%c'\n", c);
- break;
-
- }
- }
-
- if ( n_threads < 1 ) {
- ERROR("Invalid number of threads.\n");
- return 1;
- }
-
- if ( (n_threads > 1) && (image_prefix != NULL) ) {
- ERROR("Option \"--images\" is incompatible with \"-j\".\n");
- return 1;
- }
-
- /* Load cell */
- if ( cellfile == NULL ) {
- ERROR("You need to give a PDB file with the unit cell.\n");
- return 1;
- }
- cell = load_cell_from_file(cellfile);
- if ( cell == NULL ) {
- ERROR("Failed to get cell from '%s'\n", cellfile);
- return 1;
- }
- free(cellfile);
-
- if ( !cell_is_sensible(cell) ) {
- ERROR("Invalid unit cell parameters:\n");
- cell_print(cell);
- return 1;
- }
-
- /* Load geometry */
- if ( geomfile == NULL ) {
- ERROR("You need to give a geometry file.\n");
- return 1;
- }
- dtempl = data_template_new_from_file(geomfile);
- if ( dtempl == NULL ) {
- ERROR("Failed to read geometry from '%s'\n", geomfile);
- return 1;
- }
-
- if ( save_file == NULL ) save_file = strdup("partial_sim.hkl");
-
- /* Load (full) reflections */
- if ( input_file != NULL ) {
-
- RefList *as;
- char *sym_str_fromfile = NULL;
-
- full = read_reflections_2(input_file, &sym_str_fromfile);
- if ( full == NULL ) {
- ERROR("Failed to read reflections from '%s'\n",
- input_file);
- return 1;
- }
-
- /* If we don't have a point group yet, and if the file provides
- * one, use the one from the file */
- if ( (sym_str == NULL) && (sym_str_fromfile != NULL) ) {
- sym_str = sym_str_fromfile;
- STATUS("Using symmetry from reflection file: %s\n",
- sym_str);
- }
-
- /* If we still don't have a point group, use "1" */
- if ( sym_str == NULL ) sym_str = strdup("1");
-
- pointgroup_warning(sym_str);
- sym = get_pointgroup(sym_str);
-
- if ( check_list_symmetry(full, sym) ) {
- ERROR("The input reflection list does not appear to"
- " have symmetry %s\n", symmetry_name(sym));
- if ( cell_get_lattice_type(cell) == L_MONOCLINIC ) {
- ERROR("You may need to specify the unique axis "
- "in your point group. The default is "
- "unique axis c.\n");
- ERROR("See 'man crystfel' for more details.\n");
- }
- return 1;
- }
-
- as = asymmetric_indices(full, sym);
- reflist_free(full);
- full = as;
-
- } else {
- random_intensities = 1;
- if ( sym_str == NULL ) sym_str = strdup("1");
- sym = get_pointgroup(sym_str);
- }
-
- if ( n < 1 ) {
- ERROR("Number of patterns must be at least 1.\n");
- return 1;
- }
-
- if ( output_file == NULL ) {
- ERROR("You must give a filename for the output.\n");
- return 1;
- }
- stream = stream_open_for_write(output_file, dtempl);
- if ( stream == NULL ) {
- ERROR("Couldn't open output file '%s'\n", output_file);
- return 1;
- }
- free(output_file);
-
- stream_write_geometry_file(stream, geomfile);
- stream_write_commandline_args(stream, argc, argv);
-
- if ( template != NULL ) {
- template_stream = stream_open_for_read(template);
- if ( template_stream == NULL ) {
- ERROR("Couldn't open template stream '%s'\n", template);
- return 1;
- }
- }
-
- test_image = image_create_for_simulation(dtempl);
- if ( test_image == NULL ) {
- ERROR("Could not create image structure.\n");
- ERROR("Does the geometry file contain references?\n");
- return 1;
- }
-
- STATUS("Simulation parameters:\n");
- STATUS(" Wavelength: %.5f A (photon energy %.2f eV)\n",
- test_image->lambda*1e10,
- ph_lambda_to_eV(test_image->lambda));
- STATUS(" Beam divergence: 0 (not modelled)\n");
- STATUS(" Beam bandwidth: %.5f %%\n",
- test_image->bw*100.0);
- STATUS("Reciprocal space profile radius: %e m^-1\n",
- profile_radius);
- if ( image_prefix != NULL ) {
- STATUS(" Background: %.2f detector units\n",
- background);
- } else {
- STATUS(" Background: none (no image "
- "output)\n");
- }
- STATUS(" Partiality model: xsphere (hardcoded)\n");
- STATUS(" Noise standard deviation: %.2f detector units\n",
- noise_stddev);
- if ( random_intensities ) {
- STATUS(" Full intensities: randomly generated: "
- "abs(Gaussian(sigma=%.2f)), symmetry %s\n",
- full_stddev, sym_str);
- } else {
- STATUS(" Full intensities: from %s (symmetry %s)\n",
- input_file, sym_str);
- }
- STATUS(" Max error in cell components: %.2f %%\n", cnoise);
- STATUS("Scale factor standard deviation: %.2f\n", osf_stddev);
- if ( template_stream != NULL ) {
- STATUS("Crystal orientations and reflections to use from %s\n",
- template);
- }
-
- if ( random_intensities ) {
- full = reflist_new();
- }
-
- qargs.full = full;
- pthread_rwlock_init(&qargs.full_lock, NULL);
- qargs.n_to_do = n;
- qargs.n_done = 0;
- qargs.n_started = 0;
- qargs.sym = sym;
- qargs.dtempl = dtempl;
- qargs.random_intensities = random_intensities;
- qargs.cell = cell;
- qargs.stream = stream;
- qargs.cnoise = cnoise;
- qargs.osf_stddev = osf_stddev;
- qargs.full_stddev = full_stddev;
- qargs.noise_stddev = noise_stddev;
- qargs.background = background;
- qargs.max_q = detgeom_max_resolution(test_image->detgeom,
- test_image->lambda);
- qargs.image_prefix = image_prefix;
- qargs.profile_radius = profile_radius;
- qargs.template_stream = template_stream;
-
- qargs.rngs = malloc(n_threads * sizeof(gsl_rng *));
- if ( qargs.rngs == NULL ) {
- ERROR("Failed to allocate RNGs\n");
- return 1;
- }
-
- if ( config_random ) {
-
- FILE *fh;
-
- fh = fopen("/dev/urandom", "r");
- if ( fh == NULL ) {
- ERROR("Failed to open /dev/urandom. Try again without"
- " --really-random.\n");
- free(qargs.rngs);
- return 1;
- }
-
- for ( i=0; i<n_threads; i++ ) {
-
- unsigned long int seed;
- qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
-
- if ( fread(&seed, sizeof(seed), 1, fh) == 1 ) {
- gsl_rng_set(qargs.rngs[i], seed);
- } else {
- ERROR("Failed to seed RNG %i\n", i);
- }
-
- }
-
- fclose(fh);
-
- } else {
- gsl_rng *rng_for_seeds;
- rng_for_seeds = gsl_rng_alloc(gsl_rng_mt19937);
- for ( i=0; i<n_threads; i++ ) {
- qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
- gsl_rng_set(qargs.rngs[i], gsl_rng_get(rng_for_seeds));
- }
- gsl_rng_free(rng_for_seeds);
- }
-
- for ( i=0; i<NBINS; i++ ) {
- qargs.n_ref[i] = 0;
- qargs.p_hist[i] = 0.0;
- qargs.p_max[i] = 0.0;
- }
-
- run_threads(n_threads, run_job, create_job, finalise_job,
- &qargs, n, 0, 0, 0);
-
- if ( random_intensities ) {
- STATUS("Writing full intensities to %s\n", save_file);
- write_reflist_2(save_file, full, sym);
- }
-
- if ( phist_file != NULL ) {
-
- FILE *fh;
-
- fh = fopen(phist_file, "w");
-
- if ( fh != NULL ) {
-
- double overall_max = 0.0;
- double overall_mean = 0.0;
- long long int overall_total = 0;
-
- for ( i=0; i<NBINS; i++ ) {
-
- double rcen;
-
- if ( qargs.p_max[i] > overall_max ) {
- overall_max = qargs.p_max[i];
- }
-
- overall_mean += qargs.p_hist[i];
- overall_total += qargs.n_ref[i];
-
- rcen = i/(double)NBINS*qargs.max_q
- + qargs.max_q/(2.0*NBINS);
- fprintf(fh, "%.2f %7lu %.3f %.3f\n", rcen/1.0e9,
- qargs.n_ref[i],
- qargs.p_hist[i]/qargs.n_ref[i],
- qargs.p_max[i]);
-
- }
-
- fclose(fh);
-
- overall_mean /= overall_total;
-
- STATUS("Overall max partiality = %.2f\n",
- overall_max);
- STATUS("Overall mean partiality = %.2f\n", overall_mean);
- STATUS("Total number of reflections = %lli\n",
- overall_total);
-
- } else {
- ERROR("Failed to open file '%s' for writing.\n",
- phist_file);
- }
-
- }
-
- for ( i=0; i<n_threads; i++ ) {
- gsl_rng_free(qargs.rngs[i]);
- }
- free(qargs.rngs);
- pthread_rwlock_destroy(&qargs.full_lock);
- stream_close(stream);
- cell_free(cell);
- data_template_free(dtempl);
- free_symoplist(sym);
- reflist_free(full);
- free(save_file);
- free(geomfile);
- free(input_file);
-
- return 0;
-}
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
deleted file mode 100644
index 9769a0d4..00000000
--- a/src/pattern_sim.c
+++ /dev/null
@@ -1,1183 +0,0 @@
-/*
- * pattern_sim.c
- *
- * Simulate diffraction patterns from small crystals
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2009-2020 Thomas White <taw@physics.org>
- * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- * 2014 Valerio Mariani
- * 2013 Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <getopt.h>
-#include <hdf5.h>
-
-#include <image.h>
-#include <cell.h>
-#include <cell-utils.h>
-#include <utils.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"
-#include "pattern_sim.h"
-#include "version.h"
-
-
-enum spectrum_type {
- SPECTRUM_TOPHAT, /**< A top hat distribution */
- SPECTRUM_GAUSSIAN, /**< A Gaussian distribution */
- SPECTRUM_SASE, /**< A simulated SASE spectrum */
- SPECTRUM_TWOCOLOUR, /**< A spectrum containing two peaks */
- SPECTRUM_FROMFILE /**< An arbitrary spectrum read from a file */
-};
-
-static void show_help(const char *s)
-{
- printf("Syntax: %s [options]\n\n", s);
- printf(
-"Simulate diffraction patterns from small crystals probed with femtosecond\n"
-"pulses of X-rays from a free electron laser.\n"
-"\n"
-" -h, --help Display this help message.\n"
-" --version Print CrystFEL version number and exit.\n"
-"\n"
-" -p, --pdb=<file> File from which to get the unit cell.\n"
-" (The actual Bragg intensities come from the\n"
-" intensities file)\n"
-" --gpu Use the GPU to speed up the calculation.\n"
-" --gpu-dev=<n> Use GPU device <n>. Omit this option to see the\n"
-" available devices.\n"
-" -g, --geometry=<file> Get detector geometry from file.\n"
-" -n, --number=<N> Generate N images. Default 1.\n"
-" --no-images Do not output any HDF5 files.\n"
-" -o, --output=<filename> Output HDF5 filename. Default: sim.h5.\n"
-" -r, --random-orientation Use randomly generated orientations.\n"
-" --powder=<file> Write a summed pattern of all images simulated by\n"
-" this invocation as the given filename.\n"
-" -i, --intensities=<file> Specify file containing reflection intensities\n"
-" (and phases) to use.\n"
-" -y, --symmetry=<sym> The symmetry of the intensities file.\n"
-" -t, --gradients=<method> Select method for calculation of shape transforms\n"
-" --really-random Seed the random number generator with /dev/urandom.\n"
-" --min-size=<s> Minimum crystal size in nm.\n"
-" --max-size=<s> Maximum crystal size in nm.\n"
-" --no-noise Do not calculate Poisson noise.\n"
-" -s, --sample-spectrum=<N> Use N samples from spectrum. Default 3.\n"
-" -x, --spectrum=<type> Type of spectrum to simulate.\n"
-" --background=<N> Add N photons of Poisson background (default 0).\n"
-" --template=<file> Take orientations from stream <file>.\n"
-" --no-fringes Exclude the side maxima of Bragg peaks.\n"
-" --flat Make Bragg peaks flat.\n"
-" --sase-spike-width SASE spike width (standard deviation in m^-1)\n"
-" Default 2e6 m^-1\n"
-" --twocol-separation Separation between peaks in two-colour mode in m^-1\n"
-" Default 8e6 m^-1\n"
-" --nphotons Number of photons per X-ray pulse. Default 1e12.\n"
-" --beam-radius Radius of beam in metres (default 1e-6).\n"
-);
-}
-
-
-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,
- int *n_neg2, int *n_inf2, int *n_nan2,
- double *max_tt)
-{
- int fs, ss;
-
- for ( ss=0; ss<p->h; ss++ ) {
- for ( fs=0; fs<p->w; fs++ ) {
-
- double counts;
- double intensity, sa;
- double pix_area, Lsq;
- double xs, ys, rx, ry;
- double dsq, proj_area;
- float dval;
- double twotheta;
-
- intensity = (double)dp[fs + p->w*ss];
- if ( isinf(intensity) ) (*n_inf1)++;
- if ( intensity < 0.0 ) (*n_neg1)++;
- if ( isnan(intensity) ) (*n_nan1)++;
-
- /* Area of one pixel */
- 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->pixel_pitch;
- ry = (ys + p->cny) * p->pixel_pitch;
- dsq = pow(rx, 2.0) + pow(ry, 2.0);
- twotheta = atan2(sqrt(dsq), p->cnz*p->pixel_pitch);
-
- /* Area of pixel as seen from crystal (approximate) */
- proj_area = pix_area * cos(twotheta);
-
- /* Projected area of pixel divided by distance squared */
- sa = proj_area / (dsq + Lsq);
-
- if ( do_poisson ) {
- counts = poisson_noise(rng, intensity * ph_per_e * sa);
- } else {
- counts = intensity * ph_per_e * sa;
- }
-
- /* Number of photons in pixel */
- dval = counts + poisson_noise(rng, background);
-
- /* Convert to ADU */
- dval *= p->adu_per_photon;
-
- /* Saturation */
- if ( dval > p->max_adu ) dval = p->max_adu;
-
- dp[fs + p->w*ss] = dval;
-
- /* Sanity checks */
- if ( isinf(dp[fs + p->w*ss]) ) n_inf2++;
- if ( isnan(dp[fs + p->w*ss]) ) n_nan2++;
- if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++;
-
- if ( twotheta > *max_tt ) *max_tt = twotheta;
-
- }
- }
-}
-
-
-void record_image(struct image *image, int do_poisson, double background,
- gsl_rng *rng, double beam_radius, double nphotons)
-{
- double total_energy, energy_density;
- double ph_per_e;
- double area;
- double max_tt = 0.0;
- int n_inf1 = 0;
- int n_neg1 = 0;
- int n_nan1 = 0;
- int n_inf2 = 0;
- int n_neg2 = 0;
- int n_nan2 = 0;
- int pn;
-
- /* How many photons are scattered per electron? */
- area = M_PI*pow(beam_radius, 2.0);
- total_energy = nphotons * ph_lambda_to_en(image->lambda);
- energy_density = total_energy / area;
- ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0);
- STATUS("Fluence = %8.2e photons, "
- "Energy density = %5.3f kJ/cm^2, "
- "Total energy = %5.3f microJ\n",
- nphotons, energy_density/1e7, total_energy*1e6);
-
- for ( pn=0; pn<image->detgeom->n_panels; 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,
- &n_neg2, &n_inf2, &n_nan2, &max_tt);
- }
-
- STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
- rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
-
- STATUS("Halve the d values to get the voxel size for a synthesis.\n");
-
- if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) {
- ERROR("WARNING: The raw calculation produced %i negative"
- " values, %i infinities and %i NaNs.\n",
- n_neg1, n_inf1, n_nan1);
- ERROR("WARNING: After processing, there were %i negative"
- " values, %i infinities and %i NaNs.\n",
- n_neg2, n_inf2, n_nan2);
- }
-}
-
-
-static double *intensities_from_list(RefList *list, SymOpList *sym)
-{
- Reflection *refl;
- RefListIterator *iter;
- double *out = new_arr_intensity();
- SymOpMask *m = new_symopmask(sym);
- int neq = num_equivs(sym, NULL);
-
- for ( refl = first_refl(list, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) ) {
-
- signed int h, k, l;
- int eps;
- double intensity = get_intensity(refl);
-
- get_indices(refl, &h, &k, &l);
- special_position(sym, m, h, k, l);
- eps = neq / num_equivs(sym, m);
-
- set_arr_intensity(out, h, k, l, intensity / eps);
-
- }
-
- return out;
-}
-
-
-static double *phases_from_list(RefList *list)
-{
- Reflection *refl;
- RefListIterator *iter;
- double *out = new_arr_phase();
-
- for ( refl = first_refl(list, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) ) {
-
- signed int h, k, l;
- double phase = get_phase(refl, NULL);
-
- get_indices(refl, &h, &k, &l);
-
- set_arr_phase(out, h, k, l, phase);
-
- }
-
- return out;
-
-}
-
-
-static unsigned char *flags_from_list(RefList *list)
-{
- Reflection *refl;
- RefListIterator *iter;
- unsigned char *out = new_arr_flag();
-
- for ( refl = first_refl(list, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) ) {
-
- signed int h, k, l;
-
- get_indices(refl, &h, &k, &l);
-
- set_arr_flag(out, h, k, l, 1);
-
- }
-
- return out;
-
-}
-
-
-static struct quaternion read_quaternion()
-{
- do {
-
- int r;
- float w, x, y, z;
- char line[1024];
- char *rval;
-
- printf("Please input quaternion: w x y z\n");
- rval = fgets(line, 1023, stdin);
- if ( rval == NULL ) return invalid_quaternion();
- chomp(line);
-
- r = sscanf(line, "%f %f %f %f", &w, &x, &y, &z);
- if ( r == 4 ) {
-
- struct quaternion quat;
-
- quat.w = w;
- quat.x = x;
- quat.y = y;
- quat.z = z;
-
- return quat;
-
- } else {
- ERROR("Invalid rotation '%s'\n", line);
- }
-
- } while ( 1 );
-}
-
-
-static int random_ncells(double len, double min_m, double max_m)
-{
- int min_cells, max_cells, wid;
-
- min_cells = min_m / len;
- max_cells = max_m / len;
- wid = max_cells - min_cells;
-
- return wid*random()/RAND_MAX + min_cells;
-}
-
-
-static void add_metadata(const char *filename,
- struct quaternion q,
- UnitCell *cell)
-{
- hid_t fh, gh, sh, dh;
- herr_t r;
- hsize_t size[1];
- float data[4];
- double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
-
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
-
- fh = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
- if ( fh < 0 ) {
- ERROR("Failed to open file to add metadata.\n");
- return;
- }
-
- gh = H5Gcreate2(fh, "pattern_sim", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if ( gh < 0 ) {
- ERROR("Failed to create metadata group.\n");
- H5Fclose(fh);
- return;
- }
-
- size[0] = 4;
- sh = H5Screate_simple(1, size, NULL);
-
- dh = H5Dcreate2(gh, "quaternion", H5T_NATIVE_FLOAT, sh,
- H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't create dataset\n");
- H5Fclose(fh);
- return;
- }
-
- data[0] = q.w;
- data[1] = q.x;
- data[2] = q.y;
- data[3] = q.z;
- r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
- if ( r < 0 ) {
- ERROR("Failed to write quaternion to file\n");
- }
- H5Dclose(dh);
-
- size[0] = 3;
- sh = H5Screate_simple(1, size, NULL);
-
- dh = H5Dcreate2(gh, "astar", H5T_NATIVE_FLOAT, sh,
- H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't create dataset\n");
- H5Fclose(fh);
- return;
- }
- data[0] = asx;
- data[1] = asy;
- data[2] = asz;
- r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
- if ( r < 0 ) {
- ERROR("Failed to write astar to file.\n");
- }
- H5Dclose(dh);
-
- dh = H5Dcreate2(gh, "bstar", H5T_NATIVE_FLOAT, sh,
- H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't create dataset\n");
- H5Fclose(fh);
- return;
- }
- data[0] = bsx;
- data[1] = bsy;
- data[2] = bsz;
- r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
- if ( r < 0 ) {
- ERROR("Failed to write bstar to file.\n");
- }
- H5Dclose(dh);
-
- dh = H5Dcreate2(gh, "cstar", H5T_NATIVE_FLOAT, sh,
- H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't create dataset\n");
- H5Fclose(fh);
- return;
- }
- data[0] = csx;
- data[1] = csy;
- data[2] = csz;
- r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
- if ( r < 0 ) {
- ERROR("Failed to write cstar to file.\n");
- }
- H5Dclose(dh);
-
- H5Gclose(gh);
- H5Fclose(fh);
-}
-
-
-int main(int argc, char *argv[])
-{
- int argn;
- struct image *image;
- DataTemplate *dtempl;
- struct gpu_context *gctx = NULL;
- struct image *powder;
- char *intfile = NULL;
- double *intensities;
- char *rval;
- double *phases;
- unsigned char *flags;
- int config_randomquat = 0;
- int config_noimages = 0;
- int config_nonoise = 0;
- int config_nosfac = 0;
- int config_gpu = 0;
- int config_random = 0;
- char *powder_fn = NULL;
- char *cell_filename = NULL;
- char *grad_str = NULL;
- char *outfile = NULL;
- char *geometry = NULL;
- char *spectrum_str = NULL;
- char *spectrum_fn = NULL;
- GradientMethod grad;
- enum spectrum_type spectrum_type;
- int ndone = 0; /* Number of simulations done (images or not) */
- int number = 1; /* Number used for filename of image */
- int n_images = 1; /* Generate one image by default */
- int done = 0;
- UnitCell *input_cell;
- struct quaternion orientation;
- int gpu_dev = -1;
- int random_size = 0;
- double min_size = 0.0;
- double max_size = 0.0;
- char *sym_str = NULL;
- SymOpList *sym;
- int nsamples = 3;
- gsl_rng *rng;
- double background = 0.0;
- char *template_file = NULL;
- Stream *st = NULL;
- int no_fringes = 0;
- int flat = 0;
- double nphotons = 1e12;
- double beam_radius = 1e-6; /* metres */
- double sase_spike_width = 2e6; /* m^-1 */
- double twocol_sep = 8e6; /* m^-1 */
-
- /* Long options */
- const struct option longopts[] = {
- {"help", 0, NULL, 'h'},
- {"version", 0, NULL, 'v'},
- {"gpu", 0, &config_gpu, 1},
- {"beam", 1, NULL, 'b'},
- {"random-orientation", 0, NULL, 'r'},
- {"number", 1, NULL, 'n'},
- {"no-images", 0, &config_noimages, 1},
- {"no-noise", 0, &config_nonoise, 1},
- {"intensities", 1, NULL, 'i'},
- {"symmetry", 1, NULL, 'y'},
- {"powder", 1, NULL, 'w'},
- {"gradients", 1, NULL, 't'},
- {"pdb", 1, NULL, 'p'},
- {"output", 1, NULL, 'o'},
- {"geometry", 1, NULL, 'g'},
- {"sample-spectrum", 1, NULL, 's'},
- {"type-spectrum", 1, NULL, 'x'},
- {"spectrum", 1, NULL, 'x'},
- {"really-random", 0, &config_random, 1},
- {"no-fringes", 0, &no_fringes, 1},
- {"flat", 0, &flat, 1},
-
- {"gpu-dev", 1, NULL, 2},
- {"min-size", 1, NULL, 3},
- {"max-size", 1, NULL, 4},
- {"background", 1, NULL, 5},
- {"template", 1, NULL, 6},
- {"beam-bandwidth", 1, NULL, 7},
- {"photon-energy", 1, NULL, 9},
- {"nphotons", 1, NULL, 10},
- {"beam-radius", 1, NULL, 11},
- {"spectrum-file", 1, NULL, 12},
- {"sase-spike-width", 1, NULL, 13},
- {"twocol-separation", 1, NULL, 14},
-
- {0, 0, NULL, 0}
- };
-
- /* Short options */
- while ((argn = getopt_long(argc, argv, "hrn:i:t:p:o:g:y:s:x:vb:",
- longopts, NULL)) != -1) {
-
- switch (argn) {
-
- case 'h' :
- show_help(argv[0]);
- return 0;
-
- case 'v' :
- printf("CrystFEL: %s\n",
- crystfel_version_string());
- printf("%s\n",
- crystfel_licence_string());
- return 0;
-
- case 'b' :
- ERROR("WARNING: This version of CrystFEL no longer "
- "uses beam files. Please remove the beam file "
- "from your pattern_sim command line.\n");
- return 1;
-
- case 'r' :
- config_randomquat = 1;
- break;
-
- case 'n' :
- n_images = strtol(optarg, &rval, 10);
- if ( *rval != '\0' ) {
- ERROR("Invalid number of images.\n");
- return 1;
- }
- break;
-
- case 'i' :
- intfile = strdup(optarg);
- break;
-
- case 't' :
- grad_str = strdup(optarg);
- break;
-
- case 'p' :
- cell_filename = strdup(optarg);
- break;
-
- case 'o' :
- outfile = strdup(optarg);
- break;
-
- case 'w' :
- powder_fn = strdup(optarg);
- break;
-
- case 'g' :
- geometry = strdup(optarg);
- break;
-
- case 'y' :
- sym_str = strdup(optarg);
- break;
-
- case 's' :
- nsamples = strtol(optarg, &rval, 10);
- if ( *rval != '\0' ) {
- ERROR("Invalid number of spectrum samples.\n");
- return 1;
- }
- break;
-
- case 'x' :
- spectrum_str = strdup(optarg);
- break;
-
- case 2 :
- gpu_dev = atoi(optarg);
- break;
-
- case 3 :
- min_size = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid minimum size.\n");
- return 1;
- }
- min_size /= 1e9;
- random_size++;
- break;
-
- case 4 :
- max_size = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid maximum size.\n");
- return 1;
- }
- max_size /= 1e9;
- random_size++;
- break;
-
- case 5 :
- background = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid background level.\n");
- return 1;
- }
- break;
-
- case 6 :
- template_file = strdup(optarg);
- break;
-
- case 7 :
- ERROR("--beam-bandwidth is no longer used.\n");
- ERROR("Set the bandwidth in the geometry file instead.\n");
- return 1;
-
- case 9 :
- ERROR("--photon-energy is no longer used.\n");
- ERROR("Set the photon energy in the geometry file instead.\n");
- return 1;
-
- case 10 :
- nphotons = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid number of photons.\n");
- return 1;
- }
- if ( nphotons < 0.0 ) {
- ERROR("Number of photons must be positive.\n");
- return 1;
- }
- break;
-
- case 11 :
- beam_radius = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid beam radius.\n");
- return 1;
- }
- if ( beam_radius < 0.0 ) {
- ERROR("Beam radius must be positive.\n");
- return 1;
- }
- break;
-
- case 12 :
- spectrum_fn = strdup(optarg);
- break;
-
- case 13 :
- sase_spike_width = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid SASE spike width.\n");
- return 1;
- }
- if ( beam_radius < 0.0 ) {
- ERROR("SASE spike width must be positive.\n");
- return 1;
- }
- break;
-
- case 14 :
- twocol_sep = strtod(optarg, &rval);
- if ( *rval != '\0' ) {
- ERROR("Invalid two colour separation.\n");
- return 1;
- }
- if ( beam_radius < 0.0 ) {
- ERROR("Two-colour separation must be positive.\n");
- return 1;
- }
- break;
-
- case 0 :
- break;
-
- case '?' :
- break;
-
- default :
- ERROR("Unhandled option '%c'\n", argn);
- break;
-
- }
-
- }
-
- if ( random_size == 1 ) {
- ERROR("You must specify both --min-size and --max-size.\n");
- return 1;
- }
-
- if ( cell_filename == NULL ) {
- cell_filename = strdup("molecule.pdb");
- }
-
- if ( outfile == NULL ) {
- if ( n_images == 1 ) {
- outfile = strdup("sim.h5");
- } else {
- outfile = strdup("sim");
- }
- }
-
- if ( template_file != NULL ) {
- if ( config_randomquat ) {
- ERROR("You cannot use -r and --template together.\n");
- return 1;
- }
- st = stream_open_for_read(template_file);
- if ( st == NULL ) {
- ERROR("Failed to open stream.\n");
- return 1;
- }
- free(template_file);
- }
-
- if ( grad_str == NULL ) {
- STATUS("You didn't specify a gradient calculation method, so"
- " I'm using the 'mosaic' method, which is fastest.\n");
- grad = GRADIENT_MOSAIC;
- } else if ( strcmp(grad_str, "mosaic") == 0 ) {
- grad = GRADIENT_MOSAIC;
- } else if ( strcmp(grad_str, "interpolate") == 0) {
- grad = GRADIENT_INTERPOLATE;
- } else if ( strcmp(grad_str, "phased") == 0) {
- grad = GRADIENT_PHASED;
- } else {
- ERROR("Unrecognised gradient method '%s'\n", grad_str);
- return 1;
- }
- free(grad_str);
-
- if ( config_gpu && (grad != GRADIENT_MOSAIC) ) {
- ERROR("Only the mosaic method can be used for gradients when"
- "calculating on the GPU.\n");
- return 1;
- }
-
- if ( geometry == NULL ) {
- ERROR("You need to specify a geometry file with --geometry\n");
- return 1;
- }
- 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 ( spectrum_str == NULL ) {
- STATUS("You didn't specify a spectrum type, so"
- " I'm using a 'tophat' spectrum.\n");
- spectrum_type = SPECTRUM_TOPHAT;
- } else if ( strcasecmp(spectrum_str, "tophat") == 0) {
- spectrum_type = SPECTRUM_TOPHAT;
- } else if ( strcasecmp(spectrum_str, "gaussian") == 0) {
- spectrum_type = SPECTRUM_GAUSSIAN;
- } else if ( strcasecmp(spectrum_str, "sase") == 0) {
- spectrum_type = SPECTRUM_SASE;
- } else if ( strcasecmp(spectrum_str, "twocolour") == 0 ||
- strcasecmp(spectrum_str, "twocolor") == 0 ||
- strcasecmp(spectrum_str, "twocolours") == 0 ||
- strcasecmp(spectrum_str, "twocolors") == 0) {
- spectrum_type = SPECTRUM_TWOCOLOUR;
- } else if ( strcasecmp(spectrum_str, "fromfile") == 0) {
- spectrum_type = SPECTRUM_FROMFILE;
- } else {
- ERROR("Unrecognised spectrum type '%s'\n", spectrum_str);
- return 1;
- }
- free(spectrum_str);
-
- /* Load unit cell */
- input_cell = load_cell_from_file(cell_filename);
- if ( input_cell == NULL ) {
- exit(1);
- }
-
- if ( intfile == NULL ) {
-
- /* Gentle reminder */
- STATUS("You didn't specify the file containing the ");
- STATUS("reflection intensities (with --intensities).\n");
- STATUS("I'll simulate a flat intensity distribution.\n");
- intensities = NULL;
- phases = NULL;
- flags = NULL;
-
- if ( sym_str == NULL ) sym_str = strdup("1");
- pointgroup_warning(sym_str);
- sym = get_pointgroup(sym_str);
-
- } else {
-
- RefList *reflections;
- char *sym_str_fromfile = NULL;
-
- reflections = read_reflections_2(intfile, &sym_str_fromfile);
- if ( reflections == NULL ) {
- ERROR("Problem reading input file %s\n", intfile);
- return 1;
- }
-
- /* If we don't have a point group yet, and if the file provides
- * one, use the one from the file */
- if ( (sym_str == NULL) && (sym_str_fromfile != NULL) ) {
- sym_str = sym_str_fromfile;
- STATUS("Using symmetry from reflection file: %s\n",
- sym_str);
- }
-
- /* If we still don't have a point group, use "1" */
- if ( sym_str == NULL ) sym_str = strdup("1");
-
- pointgroup_warning(sym_str);
- sym = get_pointgroup(sym_str);
-
- if ( grad == GRADIENT_PHASED ) {
- phases = phases_from_list(reflections);
- } else {
- phases = NULL;
- }
- intensities = intensities_from_list(reflections, sym);
- flags = flags_from_list(reflections);
-
- /* Check that the intensities have the correct symmetry */
- if ( check_list_symmetry(reflections, sym) ) {
- ERROR("The input reflection list does not appear to"
- " have symmetry %s\n", symmetry_name(sym));
- if ( cell_get_lattice_type(input_cell) == L_MONOCLINIC )
- {
- ERROR("You may need to specify the unique axis "
- "in your point group. The default is "
- "unique axis c.\n");
- ERROR("See 'man crystfel' for more details.\n");
- }
- return 1;
- }
-
- reflist_free(reflections);
-
- }
-
- rng = gsl_rng_alloc(gsl_rng_mt19937);
- if ( config_random ) {
- FILE *fh;
- unsigned long int seed;
- fh = fopen("/dev/urandom", "r");
- if ( fread(&seed, sizeof(seed), 1, fh) == 1 ) {
- gsl_rng_set(rng, seed);
- } else {
- ERROR("Failed to seed random number generator\n");
- }
- fclose(fh);
- }
-
- image = image_create_for_simulation(dtempl);
- powder = image_create_for_simulation(dtempl);
-
- /* Splurge a few useful numbers */
- STATUS("Simulation parameters:\n");
- 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(ph_lambda_to_eV(image->lambda))*nphotons*1e3);
- STATUS(" Beam divergence: not simulated\n");
- STATUS(" Beam radius: %.2f microns\n",
- beam_radius*1e6);
- STATUS(" Background: %.2f photons\n", background);
-
- switch ( spectrum_type ) {
-
- case SPECTRUM_TOPHAT:
- STATUS(" X-ray spectrum: top hat, "
- "width %.5f %%\n", image->bw*100.0);
- break;
-
- case SPECTRUM_GAUSSIAN:
- STATUS(" X-ray spectrum: Gaussian, "
- "bandwidth %.5f %%\n", image->bw*100.0);
- break;
-
- case SPECTRUM_SASE:
- STATUS(" X-ray spectrum: SASE, "
- "bandwidth %.5f %%\n", image->bw*100.0);
- break;
-
- case SPECTRUM_TWOCOLOUR:
- STATUS(" X-ray spectrum: two colour, "
- "bandwidth %.5f %%, separation %.5e m^-1\n",
- image->bw*100.0, twocol_sep);
- break;
-
- case SPECTRUM_FROMFILE:
- STATUS(" X-ray spectrum: from %s\n",
- spectrum_fn);
- break;
- }
- if ( random_size ) {
- STATUS(" Crystal size: random, between "
- "%.2f and %.2f nm along each of a, b and c\n",
- min_size*1e9, max_size*1e9);
- } else {
- STATUS(" Crystal size: 8 unit cells along "
- "each of a, b and c\n");
- }
- if ( intfile == NULL ) {
- STATUS(" Full intensities: all equal\n");
- } else {
- STATUS(" Full intensities: from %s (symmetry %s)\n",
- intfile, sym_str);
- }
- do {
-
- int na, nb, nc;
- double a, b, c, dis;
- UnitCell *cell;
- int err = 0;
- int pi;
-
- /* Reset image data to zero for new pattern */
- for ( pi=0; pi<image->detgeom->n_panels; pi++ ) {
- long j;
- long np = image->detgeom->panels[pi].w
- * image->detgeom->panels[pi].h;
- for ( j=0; j<np; j++ ) {
- image->dp[pi][j] = 0.0;
- }
- }
-
- if ( random_size ) {
-
- double alen, blen, clen, discard;
-
- cell_get_parameters(input_cell, &alen, &blen, &clen,
- &discard, &discard, &discard);
-
- na = random_ncells(alen, min_size, max_size);
- nb = random_ncells(blen, min_size, max_size);
- nc = random_ncells(clen, min_size, max_size);
-
- } else {
-
- na = 8;
- nb = 8;
- nc = 8;
-
- }
-
- if ( st == NULL ) {
-
- if ( config_randomquat ) {
- orientation = random_quaternion(rng);
- } else {
- orientation = read_quaternion();
- }
-
- STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
- orientation.w, orientation.x,
- orientation.y, orientation.z);
-
- if ( !quaternion_valid(orientation) ) {
- ERROR("Orientation modulus is not zero!\n");
- return 1;
- }
-
- cell = cell_rotate(input_cell, orientation);
-
- } else {
-
- struct image *templ_image;
- Crystal *cr;
-
- /* Get data from next chunk */
- templ_image = stream_read_chunk(st, 0);
- if ( templ_image == NULL ) break;
- if ( templ_image->n_crystals == 0 ) continue;
-
- cr = templ_image->crystals[0];
- cell = cell_new_from_cell(crystal_get_cell(cr));
-
- if ( templ_image->n_crystals > 1 ) {
- ERROR("Using the first crystal only.\n");
- }
-
- image_free(templ_image);
-
- }
-
- switch ( spectrum_type ) {
-
- case SPECTRUM_TOPHAT :
- image->spectrum = spectrum_generate_tophat(image->lambda,
- image->bw);
- break;
-
- case SPECTRUM_GAUSSIAN :
- 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);
- break;
-
- case SPECTRUM_TWOCOLOUR :
- image->spectrum = spectrum_generate_twocolour(image->lambda,
- image->bw,
- twocol_sep);
- break;
-
- case SPECTRUM_FROMFILE :
- image->spectrum = spectrum_load(spectrum_fn);
- break;
-
- }
-
- cell_get_parameters(cell, &a, &b, &c, &dis, &dis, &dis);
- STATUS("Particle size = %i x %i x %i"
- " ( = %5.2f x %5.2f x %5.2f nm)\n",
- na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);
-
- if ( config_gpu ) {
-
- if ( gctx == NULL ) {
- gctx = setup_gpu(config_nosfac,
- intensities, flags, sym_str,
- gpu_dev);
- }
- err = get_diffraction_gpu(gctx, image, na, nb, nc,
- cell, no_fringes, flat,
- nsamples);
-
- } else {
- get_diffraction(image, na, nb, nc, intensities, phases,
- flags, cell, grad, sym, no_fringes, flat,
- nsamples);
- }
- if ( err ) {
- ERROR("Diffraction calculation failed.\n");
- done = 1;
- goto skip;
- }
-
- record_image(image, !config_nonoise, background, rng,
- beam_radius, nphotons);
-
- if ( powder_fn != NULL ) {
-
- int pn;
-
- for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
-
- size_t w, i;
-
- 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];
- }
-
- }
-
- if ( !(ndone % 10) ) {
- image_write(powder, dtempl, powder_fn);
- }
- }
-
- if ( !config_noimages ) {
-
- char h5filename[1024];
-
- if ( n_images != 1 ) {
- snprintf(h5filename, 1023, "%s-%i.h5",
- outfile, number);
- } else {
- strncpy(h5filename, outfile, 1023);
- }
-
- number++;
-
- /* Write the output file */
- image_write(image, dtempl, h5filename);
-
- /* Add some pattern_sim-specific metadata */
- add_metadata(h5filename, orientation, cell);
-
- }
-
- cell_free(cell);
-
-skip:
- ndone++;
-
- if ( n_images && (ndone >= n_images) ) done = 1;
-
- } while ( !done );
-
- if ( powder_fn != NULL ) {
- image_write(powder, dtempl, powder_fn);
- }
-
- if ( gctx != NULL ) {
- cleanup_gpu(gctx);
- }
-
- image_free(image);
- image_free(powder);
- free(intfile);
- cell_free(input_cell);
- free(intensities);
- free(outfile);
- free(cell_filename);
- free(sym_str);
- free_symoplist(sym);
-
- gsl_rng_free(rng);
-
- return 0;
-}
diff --git a/src/pattern_sim.h b/src/pattern_sim.h
deleted file mode 100644
index f013c0fb..00000000
--- a/src/pattern_sim.h
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * pattern_sim.h
- *
- * Simulate diffraction patterns from small crystals
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2011-2021 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef PATTERN_SIM_H
-#define PATTERN_SIM_H
-
-
-/* Maxmimum index to hold values up to (can be increased if necessary)
- * WARNING: Altering this value constitutes an ABI change, and means you must
- * update data/diffraction.cl then recompile and reinstall everything. */
-#define INDMAX 130
-
-/* Array size */
-#define IDIM (INDMAX*2 +1)
-#define LIST_SIZE (IDIM*IDIM*IDIM)
-
-/* Create functions for storing reflection intensities indexed as h,k,l */
-#define LABEL(x) x##_intensity
-#define TYPE double
-#include "list_tmp.h"
-
-/* As above, but for phase values */
-#define LABEL(x) x##_phase
-#define TYPE double
-#include "list_tmp.h"
-
-/* As above, but for (unsigned) integer counts */
-#define LABEL(x) x##_count
-#define TYPE unsigned int
-#include "list_tmp.h"
-
-/* As above, but for simple flags */
-#define LABEL(x) x##_flag
-#define TYPE unsigned char
-#include "list_tmp.h"
-
-
-#endif /* PATTERN_SIM_H */
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
deleted file mode 100644
index f0b3193e..00000000
--- a/tests/gpu_sim_check.c
+++ /dev/null
@@ -1,198 +0,0 @@
-/*
- * gpu_sim_check.c
- *
- * Check that GPU simulation agrees with CPU version
- *
- * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2012-2019 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#ifdef HAVE_CLOCK_GETTIME
-#include <time.h>
-#else
-#include <sys/time.h>
-#endif
-
-#define CL_TARGET_OPENCL_VERSION 220
-
-#include "../src/diffraction.h"
-#include "../src/diffraction-gpu.h"
-#include "../src/cl-utils.h"
-
-#include <datatemplate.h>
-#include <utils.h>
-#include <image.h>
-#include <symmetry.h>
-#include <cell-utils.h>
-
-
-#ifdef HAVE_CLOCK_GETTIME
-
-static double get_hires_seconds()
-{
- struct timespec tp;
- clock_gettime(CLOCK_MONOTONIC, &tp);
- return (double)tp.tv_sec + ((double)tp.tv_nsec/1e9);
-}
-
-#else
-
-/* Fallback version of the above. The time according to gettimeofday() is not
- * monotonic, so measuring intervals based on it will screw up if there's a
- * timezone change (e.g. daylight savings) while the program is running. */
-static double get_hires_seconds()
-{
- struct timeval tp;
- gettimeofday(&tp, NULL);
- return (double)tp.tv_sec + ((double)tp.tv_usec/1e6);
-}
-
-#endif
-
-
-int main(int argc, char *argv[])
-{
- struct gpu_context *gctx;
- struct image *gpu_image;
- struct image *cpu_image;
- DataTemplate *dtempl;
- UnitCell *cell;
- UnitCell *cell_raw;
- int i;
- double gpu_min, gpu_max, gpu_tot;
- double cpu_min, cpu_max, cpu_tot;
- double dev, perc;
- double start, end;
- double gpu_time, cpu_time;
- SymOpList *sym;
- gsl_rng *rng;
-
- if ( have_gpu_device() == 0 ) {
- ERROR("No GPU device found - skipping test.\n");
- return 0;
- }
-
- rng = gsl_rng_alloc(gsl_rng_mt19937);
-
- gctx = setup_gpu(1, NULL, NULL, NULL, 0);
- if ( gctx == NULL ) {
- ERROR("Couldn't set up GPU.\n");
- return 1;
- }
-
- cell_raw = cell_new_from_parameters(28.1e-9, 28.1e-9, 16.5e-9,
- deg2rad(90.0), deg2rad(90.0), deg2rad(120.0));
-
- cell = cell_rotate(cell_raw, random_quaternion(rng));
-
- dtempl = data_template_new_from_file(argv[1]);
- if ( dtempl == NULL ) return 1;
-
- cpu_image = image_create_for_simulation(dtempl);
- gpu_image = image_create_for_simulation(dtempl);
-
- start = get_hires_seconds();
- if ( get_diffraction_gpu(gctx, gpu_image, 8, 8, 8, cell, 0, 0, 10) ) {
- return 1;
- }
- end = get_hires_seconds();
- gpu_time = end - start;
-
- sym = get_pointgroup("1");
-
- cpu_image->dp = malloc(cpu_image->detgeom->n_panels * sizeof(float *));
- if ( cpu_image->dp == NULL ) {
- ERROR("Couldn't allocate memory for result.\n");
- return 1;
- }
- for ( i=0; i<cpu_image->detgeom->n_panels; i++ ) {
- struct detgeom_panel *p = &cpu_image->detgeom->panels[i];
- cpu_image->dp[i] = calloc(p->w * p->h, sizeof(float));
- if ( cpu_image->dp[i] == NULL ) {
- ERROR("Couldn't allocate memory for panel %i\n", i);
- return 1;
- }
- }
-
- start = get_hires_seconds();
- get_diffraction(cpu_image, 8, 8, 8, NULL, NULL, NULL, cell,
- GRADIENT_MOSAIC, sym, 0, 0, 10);
- end = get_hires_seconds();
- cpu_time = end - start;
-
- free_symoplist(sym);
-
- STATUS("The GPU version was %5.2f times faster.\n", cpu_time/gpu_time);
-
- gpu_min = +INFINITY; gpu_max = -INFINITY; gpu_tot = 0.0;
- cpu_min = +INFINITY; cpu_max = -INFINITY; cpu_tot = 0.0;
- dev = 0.0;
- for ( i=0; i<cpu_image->detgeom->n_panels; i++ ) {
-
- int j;
- struct detgeom_panel *p = &cpu_image->detgeom->panels[i];
-
- for ( j=0; j<p->w*p->h; j++ ) {
-
- const double cpu = cpu_image->dp[i][j];
- const double gpu = gpu_image->dp[i][j];
-
- if ( cpu > cpu_max ) cpu_max = cpu;
- if ( cpu < cpu_min ) cpu_min = cpu;
- if ( gpu > gpu_max ) gpu_max = gpu;
- if ( gpu < gpu_min ) gpu_min = gpu;
- gpu_tot += gpu;
- cpu_tot += cpu;
- dev += fabs(gpu - cpu);
-
- }
-
- }
- perc = 100.0*dev/cpu_tot;
-
- STATUS("GPU: min=%8e, max=%8e, total=%8e\n", gpu_min, gpu_max, gpu_tot);
- STATUS("CPU: min=%8e, max=%8e, total=%8e\n", cpu_min, cpu_max, cpu_tot);
- STATUS("dev = %8e (%5.2f%% of CPU total)\n", dev, perc);
-
- if ( perc > 1.2 ) {
-
- STATUS("Test failed! I'm writing cpu-sim.h5 and gpu-sim.h5"
- " for you to inspect.\n");
-
- image_write(cpu_image, dtempl, "cpu-sim.h5");
- image_write(gpu_image, dtempl, "gpu-sim.h5");
-
- return 1;
-
- }
-
- gsl_rng_free(rng);
- cell_free(cell);
- image_free(cpu_image);
- image_free(gpu_image);
- data_template_free(dtempl);
-
- return 0;
-}
diff --git a/tests/gpu_sim_check.geom b/tests/gpu_sim_check.geom
deleted file mode 100644
index 57064e58..00000000
--- a/tests/gpu_sim_check.geom
+++ /dev/null
@@ -1,28 +0,0 @@
-photon_energy = 6 keV
-bandwidth = 0.01
-
-panel0/min_fs = 0
-panel0/max_fs = 1023
-panel0/min_ss = 0
-panel0/max_ss = 511
-panel0/corner_x = -512.0
-panel0/corner_y = -532.0
-panel0/clen = 100 mm
-panel0/res = 9090.91
-panel0/adu_per_eV = 1
-panel0/data = /data/panel0
-panel0/fs = x
-panel0/ss = y
-
-panel1/min_fs = 0
-panel1/max_fs = 1023
-panel1/min_ss = 0
-panel1/max_ss = 511
-panel1/corner_x = -512.0
-panel1/corner_y = 20.0
-panel1/clen = 100 mm
-panel1/res = 9090.91
-panel1/adu_per_eV = 1
-panel1/data = /data/panel1
-panel1/fs = x
-panel1/ss = y \ No newline at end of file
diff --git a/tests/meson.build b/tests/meson.build
index b3b3cc79..917c59db 100644
--- a/tests/meson.build
+++ b/tests/meson.build
@@ -96,16 +96,6 @@ exe = executable('prof2d_check',
dependencies : [libcrystfeldep, mdep, gsldep])
test('prof2d_check', exe)
-if opencldep.found()
- exe = executable('gpu_sim_check',
- ['gpu_sim_check.c',
- simulation_bits],
- dependencies : [libcrystfeldep, mdep, gsldep, opencldep],
- include_directories: conf_inc)
- geom = files('gpu_sim_check.geom')
- test('gpu_sim_check', exe, args: [geom])
-endif
-
# Event enumeration tests
if hdf5dep.found()
ev_enum_tests = ['ev_enum1',
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-20 20:49:04 +0000