pattern_sim.c 17.2 KB
Newer Older
1
/*
Thomas White's avatar
Thomas White committed
2
3
4
 * pattern_sim.c
 *
 * Simulate diffraction patterns from small crystals
5
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2011 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
 * Part of CrystFEL - crystallography with a FEL
9
10
11
12
13
14
15
16
17
18
19
20
21
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
Thomas White's avatar
Thomas White committed
22
#include <getopt.h>
23

24
#include "image.h"
Thomas White's avatar
Thomas White committed
25
#include "diffraction.h"
26
#include "diffraction-gpu.h"
27
#include "cell.h"
28
29
#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
30
#include "detector.h"
31
#include "peaks.h"
32
#include "beam-parameters.h"
Thomas White's avatar
Thomas White committed
33
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
34
35
#include "reflist.h"
#include "reflist-utils.h"
36
#include "pattern_sim.h"
37
38


Thomas White's avatar
Thomas White committed
39
static void show_help(const char *s)
40
{
41
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
42
	printf(
Thomas White's avatar
Thomas White committed
43
"Simulate diffraction patterns from small crystals probed with femtosecond\n"
Thomas White's avatar
Thomas White committed
44
45
"pulses of X-rays from a free electron laser.\n"
"\n"
46
" -h, --help                Display this help message.\n"
Thomas White's avatar
Thomas White committed
47
"\n"
48
49
50
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
51
"     --simulation-details  Show technical details of the simulation.\n"
52
"     --gpu                 Use the GPU to speed up the calculation.\n"
Thomas White's avatar
Thomas White committed
53
54
"     --gpu-dev=<n>         Use GPU device <n>.  Omit this option to see the\n"
"                            available devices.\n"
55
56
" -g, --geometry=<file>     Get detector geometry from file.\n"
" -b, --beam=<file>         Get beam parameters from file.\n"
57
58
59
"\n"
" -n, --number=<N>          Generate N images.  Default 1.\n"
"     --no-images           Do not output any HDF5 files.\n"
60
61
62
63
64
65
" -o, --output=<filename>   Output HDF5 filename.  Default: sim.h5.\n"
"                            If you choose to simulate more than one pattern,\n"
"                            the filename given will be postfixed with a\n"
"                            hyphen, the image number and '.h5'.  In this\n"
"                            case, the default value is 'sim', such that the\n"
"                            files produced are sim-1.h5, sim-2.h5 and so on.\n"
Thomas White's avatar
Thomas White committed
66
" -r, --random-orientation  Use a randomly generated orientation\n"
67
"                            (a new orientation will be used for each image).\n"
68
69
"     --powder=<file>       Write a summed pattern of all images simulated by\n"
"                            this invocation as the given filename.\n"
70
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
Thomas White's avatar
Thomas White committed
71
"                            (and phases) to use.\n"
72
" -y, --symmetry=<sym>      The symmetry of the intensities file.\n"
Thomas White's avatar
Thomas White committed
73
" -t, --gradients=<method>  Use <method> for the calculation of shape\n"
74
75
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
Thomas White's avatar
Thomas White committed
76
77
78
"                                           Bragg position.  This is the\n"
"                                           fastest method and the only one\n"
"                                           supported on the GPU.\n"
79
80
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
Thomas White's avatar
Thomas White committed
81
82
"                                           This method has intermediate\n"
"                                           accuracy.\n"
83
"                             phased      : As 'interpolate', but take phase\n"
Thomas White's avatar
Thomas White committed
84
85
86
"                                           values into account.  This is the\n"
"                                           most accurate method, but the\n"
"                                           slowest.\n"
87
88
89
"     --really-random       Use really random numbers for the orientation and\n"
"                            crystal size.  By default, the same sequence\n"
"                            will be used for each new run.\n"
90
91
92
93
94
"     --min-size=<s>        Generate random crystal sizes using <s> as the\n"
"                            minimum crystal size in nm. --max-size is also\n"
"                            required.\n"
"     --max-size=<s>        Use <s> as the maximum crystal size in nm.\n"
"                            --min-size is also required.\n"
Thomas White's avatar
Thomas White committed
95
96
97
98
99
"\n"
"By default, the simulation aims to be as accurate as possible.  For greater\n"
"speed, or for testing, you can choose to disable certain things using the\n"
"following options.\n"
"\n"
Thomas White's avatar
Thomas White committed
100
"     --no-noise            Do not calculate Poisson noise.\n"
101
);
Thomas White's avatar
Thomas White committed
102
103
104
}


105
static void show_details()
Thomas White's avatar
Thomas White committed
106
{
Thomas White's avatar
Thomas White committed
107
108
109
110
111
	printf(
"This program simulates diffraction patterns from small crystals illuminated\n"
"with femtosecond X-ray pulses from a free electron laser.\n"
"\n"
"The lattice transform from the specified number of unit cells is calculated\n"
112
113
"using the closed-form solution for a truncated lattice faceted on the\n"
"(001), (010) and (100) planes:\n"
Thomas White's avatar
Thomas White committed
114
"\n"
115
116
117
"I_latt(q) =  sin^2(pi*na*q.a)/sin^2(pi*q.a)\n"
"           * sin^2(pi*nb*q.b)/sin^2(pi*q.b)\n"
"           * sin^2(pi*nc*q.c)/sin^2(pi*q.c)\n"
Thomas White's avatar
Thomas White committed
118
119
120
121
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
122
123
124
"This is multiplied by a model of the underlying molecular transform, I_mol(q).\n"
"This can be approximated to varying levels of accuracy by the methods given by\n"
"the '--gradients' option.\n"
Thomas White's avatar
Thomas White committed
125
126
127
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
128
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
129
130
131
132
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
133
134
135
"\n"
"Polarisation is not currently included in pattern_sim, although it is included\n"
"in the analysis of Bragg peaks by 'indexamajig'.\n"
Thomas White's avatar
Thomas White committed
136
137
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
138
139
140
"algorithm.  When the intensity is sufficiently high that Knuth's algorithm\n"
"would result in machine precision problems, a normal distribution with\n"
"standard deviation sqrt(I) is used instead.\n"
141
);
142
143
144
}


145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
static double *intensities_from_list(RefList *list)
{
	Reflection *refl;
	RefListIterator *iter;
	double *out = new_list_intensity();

	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {

		signed int h, k, l;
		double intensity = get_intensity(refl);

		get_indices(refl, &h, &k, &l);

		set_intensity(out, h, k, l, intensity);

	}

	return out;
}


static double *phases_from_list(RefList *list)
{
	Reflection *refl;
	RefListIterator *iter;
	double *out = new_list_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_phase(out, h, k, l, phase);

	}

	return out;

}


static unsigned char *flags_from_list(RefList *list)
{
	Reflection *refl;
	RefListIterator *iter;
	unsigned char *out = new_list_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_flag(out, h, k, l, 1);

	}

	return out;

}


215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 {
242
			ERROR("Invalid rotation '%s'\n", line);
243
244
245
		}

	} while ( 1 );
246
247
248
}


249
250
251
252
253
254
255
256
257
258
259
260
static int random_ncells(double len, double min_nm, double max_nm)
{
	int min_cells, max_cells, wid;

	min_cells = min_nm*1e-9 / len;
	max_cells = max_nm*1e-9 / len;
	wid = max_cells - min_cells;

	return wid*random()/RAND_MAX + min_cells;
}


261
262
int main(int argc, char *argv[])
{
263
	int c;
264
	struct image image;
265
	struct gpu_context *gctx = NULL;
266
	double *powder;
267
268
	char *intfile = NULL;
	double *intensities;
269
	char *rval;
270
	double *phases;
271
	unsigned char *flags;
Thomas White's avatar
Thomas White committed
272
	int config_simdetails = 0;
273
	int config_randomquat = 0;
274
	int config_noimages = 0;
Thomas White's avatar
Thomas White committed
275
	int config_nonoise = 0;
276
	int config_nosfac = 0;
277
	int config_gpu = 0;
278
	int config_random = 0;
279
	char *powder_fn = NULL;
280
	char *filename = NULL;
281
	char *grad_str = NULL;
282
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
283
	char *geometry = NULL;
284
	char *beamfile = NULL;
285
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
286
287
288
	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 */
289
	int done = 0;
Thomas White's avatar
Thomas White committed
290
291
	UnitCell *input_cell;
	struct quaternion orientation;
Thomas White's avatar
Thomas White committed
292
	int gpu_dev = -1;
293
294
295
	int random_size = 0;
	double min_size = 0.0;
	double max_size = 0.0;
Thomas White's avatar
Thomas White committed
296
297
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
298

299
	/* Long options */
Thomas White's avatar
Thomas White committed
300
	const struct option longopts[] = {
301
302
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
303
		{"gpu",                0, &config_gpu,         1},
304
305
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
306
		{"no-images",          0, &config_noimages,    1},
Thomas White's avatar
Thomas White committed
307
		{"no-noise",           0, &config_nonoise,     1},
308
		{"intensities",        1, NULL,               'i'},
309
		{"symmetry",           1, NULL,               'y'},
310
		{"powder",             1, NULL,               'w'},
Thomas White's avatar
Thomas White committed
311
		{"gradients",          1, NULL,               't'},
312
		{"pdb",                1, NULL,               'p'},
313
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
314
		{"geometry",           1, NULL,               'g'},
315
		{"beam",               1, NULL,               'b'},
316
		{"really-random",      0, &config_random,      1},
Thomas White's avatar
Thomas White committed
317
		{"gpu-dev",            1, NULL,                2},
318
319
		{"min-size",           1, NULL,                3},
		{"max-size",           1, NULL,                4},
320
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
321
	};
322

323
	/* Short options */
324
	while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:y:",
Thomas White's avatar
Thomas White committed
325
	                        longopts, NULL)) != -1) {
326

Thomas White's avatar
Thomas White committed
327
		switch (c) {
Thomas White's avatar
Thomas White committed
328
		case 'h' :
Thomas White's avatar
Thomas White committed
329
330
			show_help(argv[0]);
			return 0;
331

Thomas White's avatar
Thomas White committed
332
		case 'r' :
333
334
335
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
336
		case 'n' :
337
338
339
340
341
			n_images = strtol(optarg, &rval, 10);
			if ( *rval != '\0' ) {
				ERROR("Invalid number of images.\n");
				return 1;
			}
342
343
			break;

Thomas White's avatar
Thomas White committed
344
		case 'i' :
345
346
347
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
348
		case 't' :
349
350
351
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
352
		case 'p' :
353
354
355
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
356
		case 'o' :
357
358
359
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
360
		case 'w' :
361
362
363
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
364
		case 'g' :
Thomas White's avatar
Thomas White committed
365
366
367
			geometry = strdup(optarg);
			break;

368
369
370
371
		case 'b' :
			beamfile = strdup(optarg);
			break;

372
		case 'y' :
Thomas White's avatar
Thomas White committed
373
			sym_str = strdup(optarg);
374
375
			break;

Thomas White's avatar
Thomas White committed
376
377
378
379
		case 2 :
			gpu_dev = atoi(optarg);
			break;

380
		case 3 :
381
382
383
384
385
			min_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid minimum size.\n");
				return 1;
			}
386
387
388
389
			random_size++;
			break;

		case 4 :
390
391
392
393
394
			max_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid maximum size.\n");
				return 1;
			}
395
396
397
			random_size++;
			break;

Thomas White's avatar
Thomas White committed
398
		case 0 :
Thomas White's avatar
Thomas White committed
399
			break;
400

Thomas White's avatar
Thomas White committed
401
		default :
Thomas White's avatar
Thomas White committed
402
403
			return 1;
		}
404
405
406

	}

407
	if ( config_random ) {
408
409
410
411
412
413
414
415
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

416
417
418
419
420
	if ( random_size == 1 ) {
		ERROR("You must specify both --min-size and --max-size.\n");
		return 1;
	}

421
422
423
424
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

425
426
427
428
429
430
431
432
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
433
434
435
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	/* sym_str is used below */
436

Thomas White's avatar
Thomas White committed
437
	if ( config_simdetails ) {
438
		show_details();
Thomas White's avatar
Thomas White committed
439
440
441
		return 0;
	}

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
	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;
	}

Thomas White's avatar
Thomas White committed
464
465
466
467
468
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

469
470
471
472
473
474
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

475
	if ( intfile == NULL ) {
476

477
478
479
480
481
482
		/* 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;
483
484
		flags = NULL;

485
	} else {
486

Thomas White's avatar
Thomas White committed
487
488
489
		RefList *reflections;

		reflections = read_reflections(intfile);
Thomas White's avatar
Thomas White committed
490
491
492
493
		if ( reflections == NULL ) {
			ERROR("Problem reading input file %s\n", intfile);
			return 1;
		}
Thomas White's avatar
Thomas White committed
494
		free(intfile);
495

496
		if ( grad == GRADIENT_PHASED ) {
Thomas White's avatar
Thomas White committed
497
			phases = phases_from_list(reflections);
498
499
500
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
501
502
503
		intensities = intensities_from_list(reflections);
		phases = phases_from_list(reflections);
		flags = flags_from_list(reflections);
504

Thomas White's avatar
Thomas White committed
505
		/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
506
		if ( check_list_symmetry(reflections, sym) ) {
Thomas White's avatar
Thomas White committed
507
			ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
508
			      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
509
510
511
			return 1;
		}

Thomas White's avatar
Thomas White committed
512
513
		reflist_free(reflections);

514
515
	}

516
517
518
519
520
521
522
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

523
524
525
526
527
528
529
	image.beam = get_beam_parameters(beamfile);
	if ( image.beam == NULL ) {
		ERROR("Failed to read beam parameters from '%s'\n", beamfile);
		return 1;
	}
	free(beamfile);

Thomas White's avatar
Thomas White committed
530
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
531
532
	image.width = image.det->max_fs + 1;
	image.height = image.det->max_ss + 1;
533
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
534
535
536
537

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
538
539
		exit(1);
	}
Thomas White's avatar
Thomas White committed
540
541

	/* Initialise stuff */
542
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
543
544
	image.features = NULL;
	image.flags = NULL;
Thomas White's avatar
Thomas White committed
545
546
	image.i0 = 1.0;
	image.i0_available = 1;
Thomas White's avatar
Thomas White committed
547

548
549
	powder = calloc(image.width*image.height, sizeof(*powder));

Thomas White's avatar
Thomas White committed
550
	/* Splurge a few useful numbers */
551
	STATUS("Wavelength is %f nm\n", image.lambda/1.0e-9);
Thomas White's avatar
Thomas White committed
552

553
554
	do {

Thomas White's avatar
Thomas White committed
555
		int na, nb, nc;
556
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
557
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
558

559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
		if ( random_size ) {

			double alen, blen, clen, dis;

			cell_get_parameters(input_cell, &alen, &blen, &clen,
			                    &dis, &dis, &dis);

			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;

		}
Thomas White's avatar
Thomas White committed
577

578
579
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
580
			orientation = random_quaternion();
581
		} else {
Thomas White's avatar
Thomas White committed
582
			orientation = read_quaternion();
583
584
		}

Thomas White's avatar
Thomas White committed
585
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
586
587
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
588

Thomas White's avatar
Thomas White committed
589
		if ( !quaternion_valid(orientation) ) {
590
			ERROR("Orientation modulus is not zero!\n");
591
592
593
			return 1;
		}

Thomas White's avatar
Thomas White committed
594
595
		cell = cell_rotate(input_cell, orientation);

596
597
598
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
599

600
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
Thomas White's avatar
Thomas White committed
601
602
		STATUS("Particle size = %i x %i x %i"
		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
603
604
	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);

605
		if ( config_gpu ) {
606
			if ( gctx == NULL ) {
607
				gctx = setup_gpu(config_nosfac,
Thomas White's avatar
Thomas White committed
608
				                 intensities, flags, sym_str,
609
				                 gpu_dev);
610
			}
611
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
612
		} else {
Thomas White's avatar
Thomas White committed
613
			get_diffraction(&image, na, nb, nc, intensities, phases,
Thomas White's avatar
Thomas White committed
614
			                flags, cell, grad, sym);
615
		}
616
		if ( image.data == NULL ) {
617
618
619
620
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

621
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
622

623
		if ( powder_fn != NULL ) {
624
625
626
627
628
629
630

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
631
				powder[x+w*y] += (double)image.data[x+w*y];
632
633
634
635
			}
			}

			if ( !(ndone % 10) ) {
636
				hdf5_write(powder_fn, powder,
637
				           image.width, image.height,
638
				           H5T_NATIVE_DOUBLE);
639
640
641
			}
		}

642
643
644
645
		if ( !config_noimages ) {

			char filename[1024];

646
647
648
649
650
651
652
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

653
			number++;
Thomas White's avatar
Thomas White committed
654

655
656
			/* Write the output file */
			hdf5_write(filename, image.data,
657
			           image.width, image.height, H5T_NATIVE_FLOAT);
658
659

		}
660

661
662
663
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
664
		cell_free(cell);
Thomas White's avatar
Thomas White committed
665

666
skip:
Thomas White's avatar
Thomas White committed
667
		ndone++;
668

669
		if ( n_images && (ndone >= n_images) ) done = 1;
Thomas White's avatar
Thomas White committed
670

671
	} while ( !done );
672

673
674
675
676
677
678
	if ( powder_fn != NULL ) {
		hdf5_write(powder_fn, powder,
		           image.width, image.height,
		           H5T_NATIVE_DOUBLE);
	}

679
680
681
682
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
683
684
	free(image.det->panels);
	free(image.det);
685
	free(image.beam);
686
	free(powder);
Thomas White's avatar
Thomas White committed
687
	cell_free(input_cell);
688
	free(intensities);
Thomas White's avatar
Thomas White committed
689
690
	free(outfile);
	free(filename);
Thomas White's avatar
Thomas White committed
691
692
	free(sym_str);
	free_symoplist(sym);
693

694
	return 0;
695
}