Add partial_sim, for generating test streams by geometrical methods

1 files changed, 263 insertions, 0 deletions

diff --git a/src/partial_sim.c b/src/partial_sim.c
new file mode 100644
index 00000000..a023737f
--- /dev/null
+++ b/src/partial_sim.c
@@ -0,0 +1,263 @@
+/*
+ * partial_sim.c
+ *
+ * Generate partials for testing scaling
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#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 "utils.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+#include "beam-parameters.h"
+#include "detector.h"
+#include "geometry.h"
+#include "stream.h"
+
+
+/* For each reflection in "partial", fill in what the intensity would be
+ * according to "full" */
+static void calculate_partials(RefList *partial, RefList *full, const char *sym)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+
+	for ( refl = first_refl(partial, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+		Reflection *rfull;
+		double p;
+		double Ip;
+
+		get_indices(refl, &h, &k, &l);
+		get_asymm(h, k, l, &h, &k, &l, sym);
+		p = get_partiality(refl);
+
+		rfull = find_refl(full, h, k, l);
+		if ( rfull == NULL ) {
+			set_redundancy(refl, 0);
+		} else {
+			Ip = p * get_intensity(rfull);
+			set_int(refl, Ip);
+		}
+
+	}
+}
+
+
+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"
+"\n"
+"You need to provide the following basic options:\n"
+" -i, --input=<file>      Read reflections from <file>.\n"
+" -o, --output=<file>     Write partials in stream format to <file>.\n"
+" -g. --geometry=<file>   Get detector geometry from file.\n"
+" -b, --beam=<file>       Get beam parameters 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"
+);
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	char *input_file = NULL;
+	char *output_file = NULL;
+	char *beamfile = NULL;
+	char *geomfile = NULL;
+	char *cellfile = NULL;
+	struct detector *det = NULL;
+	struct beam_params *beam = NULL;
+	RefList *full;
+	char *sym = NULL;
+	UnitCell *cell = NULL;
+	UnitCell *rcell;
+	struct quaternion orientation;
+	struct image image;
+	FILE *ofh;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"output",             1, NULL,               'o'},
+		{"input",              1, NULL,               'i'},
+		{"beam",               1, NULL,               'b'},
+		{"pdb",                1, NULL,               'p'},
+		{"geometry",           1, NULL,               'g'},
+		{"symmetry",           1, NULL,               'y'},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hi:o:b:p:g:y:",
+	                        longopts, NULL)) != -1) {
+
+		switch (c) {
+		case 'h' :
+			show_help(argv[0]);
+			return 0;
+
+		case 'o' :
+			output_file = strdup(optarg);
+			break;
+
+		case 'i' :
+			input_file = strdup(optarg);
+			break;
+
+		case 'b' :
+			beamfile = strdup(optarg);
+			break;
+
+		case 'p' :
+			cellfile = strdup(optarg);
+			break;
+
+		case 'g' :
+			geomfile = strdup(optarg);
+			break;
+
+		case 'y' :
+			sym = strdup(optarg);
+			break;
+
+		case 0 :
+			break;
+
+		default :
+			return 1;
+		}
+
+	}
+
+	/* Load beam */
+	if ( beamfile == NULL ) {
+		ERROR("You need to provide a beam parameters file.\n");
+		return 1;
+	}
+	beam = get_beam_parameters(beamfile);
+	if ( beam == NULL ) {
+		ERROR("Failed to load beam parameters from '%s'\n", beamfile);
+		return 1;
+	}
+	free(beamfile);
+
+	/* Load cell */
+	if ( cellfile == NULL ) {
+		ERROR("You need to give a PDB file with the unit cell.\n");
+		return 1;
+	}
+	cell = load_cell_from_pdb(cellfile);
+	if ( cell == NULL ) {
+		ERROR("Failed to get cell from '%s'\n", cellfile);
+		return 1;
+	}
+	free(cellfile);
+
+	/* Load geometry */
+	if ( geomfile == NULL ) {
+		ERROR("You need to give a geometry file.\n");
+		return 1;
+	}
+	det = get_detector_geometry(geomfile);
+	if ( cell == NULL ) {
+		ERROR("Failed to read geometry from '%s'\n", geomfile);
+		return 1;
+	}
+	free(geomfile);
+
+	/* Load (full) reflections */
+	if ( input_file == NULL ) {
+		ERROR("You must provide a reflection list.\n");
+		return 1;
+	}
+	full = read_reflections(input_file);
+	if ( full == NULL ) {
+		ERROR("Failed to read reflections from '%s'\n", input_file);
+		return 1;
+	}
+	free(input_file);
+	if ( check_list_symmetry(full, sym) ) {
+		ERROR("The input reflection list does not appear to"
+		      " have symmetry %s\n", sym);
+		return 1;
+	}
+
+	if ( output_file == NULL ) {
+		ERROR("You must pgive a filename for the output.\n");
+		return 1;
+	}
+	ofh = fopen(output_file, "w");
+	if ( ofh == NULL ) {
+		ERROR("Couldn't open output file '%s'\n", output_file);
+		return 1;
+	}
+	free(output_file);
+	write_stream_header(ofh, argc, argv);
+
+	/* Set up a random orientation */
+	orientation = random_quaternion();
+	image.indexed_cell = cell_rotate(cell, orientation);
+
+	image.det = det;
+	image.width = det->max_fs;
+	image.height = det->max_ss;
+
+	image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
+	image.div = beam->divergence;
+	image.bw = beam->bandwidth;
+	image.profile_radius = 0.005e9;
+	image.i0_available = 0;
+	image.filename = "(simulated 1)";
+	image.reflections = find_intersections(&image, image.indexed_cell, 0);
+	calculate_partials(image.reflections, full, sym);
+	write_chunk(ofh, &image, STREAM_INTEGRATED);
+	reflist_free(image.reflections);
+
+	/* Rotate the cell by a tiny amount */
+	rcell = rotate_cell(image.indexed_cell,
+	                    deg2rad(0.2), deg2rad(0.0), deg2rad(0.0));
+	cell_free(image.indexed_cell);
+	image.indexed_cell = rcell;
+	image.filename = "(simulated 2)";
+
+	/* Write another chunk */
+	image.reflections = find_intersections(&image, image.indexed_cell, 0);
+	calculate_partials(image.reflections, full, sym);
+	write_chunk(ofh, &image, STREAM_INTEGRATED);
+	reflist_free(image.reflections);
+
+	fclose(ofh);
+	cell_free(cell);
+	free_detector_geometry(det);
+	free(beam);
+	free(sym);
+	reflist_free(full);
+
+	return 0;
+}