pattern_sim.c 16.8 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-2010 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"
Thomas White's avatar
Thomas White committed
32
#include "sfac.h"
33
#include "reflections.h"
34
#include "beam-parameters.h"
Thomas White's avatar
Thomas White committed
35
#include "symmetry.h"
36
37


Thomas White's avatar
Thomas White committed
38
static void show_help(const char *s)
39
{
40
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
41
	printf(
Thomas White's avatar
Thomas White committed
42
"Simulate diffraction patterns from small crystals probed with femtosecond\n"
Thomas White's avatar
Thomas White committed
43
44
"pulses of X-rays from a free electron laser.\n"
"\n"
45
" -h, --help                Display this help message.\n"
Thomas White's avatar
Thomas White committed
46
"\n"
47
48
49
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
50
"     --simulation-details  Show technical details of the simulation.\n"
51
"     --gpu                 Use the GPU to speed up the calculation.\n"
Thomas White's avatar
Thomas White committed
52
53
"     --gpu-dev=<n>         Use GPU device <n>.  Omit this option to see the\n"
"                            available devices.\n"
54
55
" -g, --geometry=<file>     Get detector geometry from file.\n"
" -b, --beam=<file>         Get beam parameters from file.\n"
56
57
58
59
"\n"
"     --near-bragg          Output h,k,l,I near Bragg conditions.\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"
100
"     --no-water            Do not simulate water background.\n"
Thomas White's avatar
Thomas White committed
101
"     --no-noise            Do not calculate Poisson noise.\n"
102
);
Thomas White's avatar
Thomas White committed
103
104
105
}


106
static void show_details()
Thomas White's avatar
Thomas White committed
107
{
Thomas White's avatar
Thomas White committed
108
109
110
111
112
	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"
113
114
"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
115
"\n"
116
117
118
"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
119
120
121
122
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
123
124
125
"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
126
"\n"
127
128
"Intensity from water is added according to the first term of equation 5\n"
"from Phys Chem Chem Phys 2003 (5) 1981--1991.  This simulates the\n"
129
"coherent, elastic part of the diffuse scattering from the water jet only.\n"
Thomas White's avatar
Thomas White committed
130
131
132
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
133
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
134
135
136
137
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
138
139
140
"\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
141
142
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
143
144
145
"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"
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
}


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 {
177
			ERROR("Invalid rotation '%s'\n", line);
178
179
180
		}

	} while ( 1 );
181
182
183
}


184
185
186
187
188
189
190
191
192
193
194
195
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;
}


196
197
int main(int argc, char *argv[])
{
198
	int c;
199
	struct image image;
200
	struct gpu_context *gctx = NULL;
201
	double *powder;
202
203
	char *intfile = NULL;
	double *intensities;
204
	char *rval;
205
	double *phases;
206
	unsigned char *flags;
Thomas White's avatar
Thomas White committed
207
	int config_simdetails = 0;
208
209
	int config_nearbragg = 0;
	int config_randomquat = 0;
210
	int config_noimages = 0;
211
	int config_nowater = 0;
Thomas White's avatar
Thomas White committed
212
	int config_nonoise = 0;
213
	int config_nosfac = 0;
214
	int config_gpu = 0;
215
	int config_random = 0;
216
	char *powder_fn = NULL;
217
	char *filename = NULL;
218
	char *grad_str = NULL;
219
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
220
	char *geometry = NULL;
221
	char *beamfile = NULL;
222
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
223
224
225
	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 */
226
	int done = 0;
Thomas White's avatar
Thomas White committed
227
228
	UnitCell *input_cell;
	struct quaternion orientation;
Thomas White's avatar
Thomas White committed
229
	int gpu_dev = -1;
230
231
232
	int random_size = 0;
	double min_size = 0.0;
	double max_size = 0.0;
233
	char *sym = NULL;
Thomas White's avatar
Thomas White committed
234

235
	/* Long options */
Thomas White's avatar
Thomas White committed
236
	const struct option longopts[] = {
237
238
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
239
		{"gpu",                0, &config_gpu,         1},
240
241
242
		{"near-bragg",         0, &config_nearbragg,   1},
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
243
		{"no-images",          0, &config_noimages,    1},
244
		{"no-water",           0, &config_nowater,     1},
Thomas White's avatar
Thomas White committed
245
		{"no-noise",           0, &config_nonoise,     1},
246
		{"intensities",        1, NULL,               'i'},
247
		{"symmetry",           1, NULL,               'y'},
248
		{"powder",             1, NULL,               'w'},
Thomas White's avatar
Thomas White committed
249
		{"gradients",          1, NULL,               't'},
250
		{"pdb",                1, NULL,               'p'},
251
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
252
		{"geometry",           1, NULL,               'g'},
253
		{"beam",               1, NULL,               'b'},
254
		{"really-random",      0, &config_random,      1},
Thomas White's avatar
Thomas White committed
255
		{"gpu-dev",            1, NULL,                2},
256
257
		{"min-size",           1, NULL,                3},
		{"max-size",           1, NULL,                4},
258
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
259
	};
260

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

Thomas White's avatar
Thomas White committed
265
		switch (c) {
Thomas White's avatar
Thomas White committed
266
		case 'h' :
Thomas White's avatar
Thomas White committed
267
268
			show_help(argv[0]);
			return 0;
269

Thomas White's avatar
Thomas White committed
270
		case 'r' :
271
272
273
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
274
		case 'n' :
275
276
277
278
279
			n_images = strtol(optarg, &rval, 10);
			if ( *rval != '\0' ) {
				ERROR("Invalid number of images.\n");
				return 1;
			}
280
281
			break;

Thomas White's avatar
Thomas White committed
282
		case 'i' :
283
284
285
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
286
		case 't' :
287
288
289
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
290
		case 'p' :
291
292
293
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
294
		case 'o' :
295
296
297
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
298
		case 'w' :
299
300
301
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
302
		case 'g' :
Thomas White's avatar
Thomas White committed
303
304
305
			geometry = strdup(optarg);
			break;

306
307
308
309
		case 'b' :
			beamfile = strdup(optarg);
			break;

310
311
312
313
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
314
315
316
317
		case 2 :
			gpu_dev = atoi(optarg);
			break;

318
		case 3 :
319
320
321
322
323
			min_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid minimum size.\n");
				return 1;
			}
324
325
326
327
			random_size++;
			break;

		case 4 :
328
329
330
331
332
			max_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid maximum size.\n");
				return 1;
			}
333
334
335
			random_size++;
			break;

Thomas White's avatar
Thomas White committed
336
		case 0 :
Thomas White's avatar
Thomas White committed
337
			break;
338

Thomas White's avatar
Thomas White committed
339
		default :
Thomas White's avatar
Thomas White committed
340
341
			return 1;
		}
342
343
344

	}

345
	if ( config_random ) {
346
347
348
349
350
351
352
353
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

354
355
356
357
358
	if ( random_size == 1 ) {
		ERROR("You must specify both --min-size and --max-size.\n");
		return 1;
	}

359
360
361
362
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

363
364
365
366
367
368
369
370
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

371
372
	if ( sym == NULL ) sym = strdup("1");

Thomas White's avatar
Thomas White committed
373
	if ( config_simdetails ) {
374
		show_details();
Thomas White's avatar
Thomas White committed
375
376
377
		return 0;
	}

378
379
380
381
382
383
	if ( (!config_nowater) && config_gpu ) {
		ERROR("Cannot simulate water scattering on the GPU.\n");
		ERROR("Please try again with the --no-water option.\n");
		return 1;
	}

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
	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
406
407
408
409
410
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

411
412
413
414
415
416
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

417
	if ( intfile == NULL ) {
418

419
420
421
422
423
424
		/* 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;
425
426
		flags = NULL;

427
	} else {
428
429

		int i;
Thomas White's avatar
Thomas White committed
430
		ReflItemList *items;
431

432
433
434
435
436
		if ( grad == GRADIENT_PHASED ) {
			phases = new_list_phase();
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
437
438
		intensities = new_list_intensity();
		phases = new_list_phase();
439
		flags = new_list_flag();
440
441
		items = read_reflections(intfile, intensities, phases,
		                         NULL, NULL);
442
		free(intfile);
443
444
445
446
447
448

		for ( i=0; i<num_items(items); i++ ) {
			struct refl_item *it = get_item(items, i);
			set_flag(flags, it->h, it->k, it->l, 1);
		}

Thomas White's avatar
Thomas White committed
449
450
451
452
453
454
455
		/* Check that the intensities have the correct symmetry */
		if ( check_symmetry(items, sym) ) {
			ERROR("The input reflection list does not appear to"
			      " have symmetry %s\n", sym);
			return 1;
		}

Thomas White's avatar
Thomas White committed
456
		delete_items(items);
457
458
	}

459
460
461
462
463
464
465
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

466
467
468
469
470
471
472
	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
473
	/* Define image parameters */
474
475
	image.width = image.det->max_x + 1;
	image.height = image.det->max_y + 1;
476
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
477
478
479
480

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
481
482
		exit(1);
	}
Thomas White's avatar
Thomas White committed
483
484

	/* Initialise stuff */
485
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
486
487
	image.features = NULL;
	image.flags = NULL;
488
489
	image.f0 = 1.0;
	image.f0_available = 1;
Thomas White's avatar
Thomas White committed
490

491
492
	powder = calloc(image.width*image.height, sizeof(*powder));

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

496
497
	do {

Thomas White's avatar
Thomas White committed
498
		int na, nb, nc;
499
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
500
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
501

502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
		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
520

521
522
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
523
			orientation = random_quaternion();
524
		} else {
Thomas White's avatar
Thomas White committed
525
			orientation = read_quaternion();
526
527
		}

Thomas White's avatar
Thomas White committed
528
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
529
530
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
531

Thomas White's avatar
Thomas White committed
532
		if ( !quaternion_valid(orientation) ) {
533
			ERROR("Orientation modulus is not zero!\n");
534
535
536
			return 1;
		}

Thomas White's avatar
Thomas White committed
537
538
		cell = cell_rotate(input_cell, orientation);

539
540
541
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
542

543
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
544
545
546
		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);

547
		if ( config_gpu ) {
548
549
			if ( gctx == NULL ) {
				gctx = setup_gpu(config_nosfac, &image,
550
551
				                 intensities, flags, sym,
				                 gpu_dev);
552
			}
553
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
554
		} else {
Thomas White's avatar
Thomas White committed
555
			get_diffraction(&image, na, nb, nc, intensities, phases,
556
557
			                flags, cell, !config_nowater, grad,
			                sym);
558
		}
559
		if ( image.data == NULL ) {
560
561
562
563
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

564
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
565

Thomas White's avatar
Thomas White committed
566
		if ( config_nearbragg ) {
567
568
569
570
571
572
573
574
575
576
577

			RefList *reflections;

			reflections = find_projected_peaks(&image, cell,
			                                   0, 0.1);

			output_intensities(&image, cell, reflections, NULL,
			                   0, 0, stdout);

			reflist_free(reflections);

Thomas White's avatar
Thomas White committed
578
579
		}

580
		if ( powder_fn != NULL ) {
581
582
583
584
585
586
587

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
588
				powder[x+w*y] += (double)image.data[x+w*y];
589
590
591
592
			}
			}

			if ( !(ndone % 10) ) {
593
				hdf5_write(powder_fn, powder,
594
				           image.width, image.height,
595
				           H5T_NATIVE_DOUBLE);
596
597
598
			}
		}

599
600
601
602
		if ( !config_noimages ) {

			char filename[1024];

603
604
605
606
607
608
609
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

610
			number++;
Thomas White's avatar
Thomas White committed
611

612
613
			/* Write the output file */
			hdf5_write(filename, image.data,
614
			           image.width, image.height, H5T_NATIVE_FLOAT);
615
616

		}
617

618
619
620
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
621
		cell_free(cell);
Thomas White's avatar
Thomas White committed
622

623
skip:
Thomas White's avatar
Thomas White committed
624
		ndone++;
625

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

628
	} while ( !done );
629

630
631
632
633
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
634
635
	free(image.det->panels);
	free(image.det);
636
	free(image.beam);
637
	free(powder);
Thomas White's avatar
Thomas White committed
638
	cell_free(input_cell);
639
	free(intensities);
Thomas White's avatar
Thomas White committed
640
641
	free(outfile);
	free(filename);
642
	free(sym);
643

644
	return 0;
645
}