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


Thomas White's avatar
Thomas White committed
37
static void show_help(const char *s)
38
{
39
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
40
	printf(
Thomas White's avatar
Thomas White committed
41
"Simulate diffraction patterns from small crystals probed with femtosecond\n"
Thomas White's avatar
Thomas White committed
42
43
"pulses of X-rays from a free electron laser.\n"
"\n"
44
" -h, --help                Display this help message.\n"
Thomas White's avatar
Thomas White committed
45
"\n"
46
47
48
" -p, --pdb=<file>          PDB file from which to get the unit cell.\n"
"                            (The actual Bragg intensities come from the\n"
"                            intensities file)\n"
49
"     --simulation-details  Show technical details of the simulation.\n"
50
"     --gpu                 Use the GPU to speed up the calculation.\n"
Thomas White's avatar
Thomas White committed
51
52
"     --gpu-dev=<n>         Use GPU device <n>.  Omit this option to see the\n"
"                            available devices.\n"
53
54
" -g, --geometry=<file>     Get detector geometry from file.\n"
" -b, --beam=<file>         Get beam parameters from file.\n"
55
56
57
58
"\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"
59
60
61
62
63
64
" -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
65
" -r, --random-orientation  Use a randomly generated orientation\n"
66
"                            (a new orientation will be used for each image).\n"
67
68
"     --powder=<file>       Write a summed pattern of all images simulated by\n"
"                            this invocation as the given filename.\n"
69
" -i, --intensities=<file>  Specify file containing reflection intensities\n"
Thomas White's avatar
Thomas White committed
70
"                            (and phases) to use.\n"
71
" -y, --symmetry=<sym>      The symmetry of the intensities file.\n"
Thomas White's avatar
Thomas White committed
72
" -t, --gradients=<method>  Use <method> for the calculation of shape\n"
73
74
"                            transform intensities.  Choose from:\n"
"                             mosaic      : Take the intensity of the nearest\n"
Thomas White's avatar
Thomas White committed
75
76
77
"                                           Bragg position.  This is the\n"
"                                           fastest method and the only one\n"
"                                           supported on the GPU.\n"
78
79
"                             interpolate : Interpolate trilinearly between\n"
"                                           six adjacent Bragg intensities.\n"
Thomas White's avatar
Thomas White committed
80
81
"                                           This method has intermediate\n"
"                                           accuracy.\n"
82
"                             phased      : As 'interpolate', but take phase\n"
Thomas White's avatar
Thomas White committed
83
84
85
"                                           values into account.  This is the\n"
"                                           most accurate method, but the\n"
"                                           slowest.\n"
86
87
88
"     --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"
89
90
91
92
93
"     --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
94
95
96
97
98
"\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
99
"     --no-noise            Do not calculate Poisson noise.\n"
100
);
Thomas White's avatar
Thomas White committed
101
102
103
}


104
static void show_details()
Thomas White's avatar
Thomas White committed
105
{
Thomas White's avatar
Thomas White committed
106
107
108
109
110
	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"
111
112
"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
113
"\n"
114
115
116
"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
117
118
119
120
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
121
122
123
"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
124
125
126
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
127
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
128
129
130
131
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
132
133
134
"\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
135
136
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
137
138
139
"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"
140
);
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
}


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 {
171
			ERROR("Invalid rotation '%s'\n", line);
172
173
174
		}

	} while ( 1 );
175
176
177
}


178
179
180
181
182
183
184
185
186
187
188
189
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;
}


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

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

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

Thomas White's avatar
Thomas White committed
257
		switch (c) {
Thomas White's avatar
Thomas White committed
258
		case 'h' :
Thomas White's avatar
Thomas White committed
259
260
			show_help(argv[0]);
			return 0;
261

Thomas White's avatar
Thomas White committed
262
		case 'r' :
263
264
265
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
266
		case 'n' :
267
268
269
270
271
			n_images = strtol(optarg, &rval, 10);
			if ( *rval != '\0' ) {
				ERROR("Invalid number of images.\n");
				return 1;
			}
272
273
			break;

Thomas White's avatar
Thomas White committed
274
		case 'i' :
275
276
277
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
278
		case 't' :
279
280
281
			grad_str = strdup(optarg);
			break;

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

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

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

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

298
299
300
301
		case 'b' :
			beamfile = strdup(optarg);
			break;

302
303
304
305
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
306
307
308
309
		case 2 :
			gpu_dev = atoi(optarg);
			break;

310
		case 3 :
311
312
313
314
315
			min_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid minimum size.\n");
				return 1;
			}
316
317
318
319
			random_size++;
			break;

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

Thomas White's avatar
Thomas White committed
328
		case 0 :
Thomas White's avatar
Thomas White committed
329
			break;
330

Thomas White's avatar
Thomas White committed
331
		default :
Thomas White's avatar
Thomas White committed
332
333
			return 1;
		}
334
335
336

	}

337
	if ( config_random ) {
338
339
340
341
342
343
344
345
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

346
347
348
349
350
	if ( random_size == 1 ) {
		ERROR("You must specify both --min-size and --max-size.\n");
		return 1;
	}

351
352
353
354
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

355
356
357
358
359
360
361
362
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

363
364
	if ( sym == NULL ) sym = strdup("1");

Thomas White's avatar
Thomas White committed
365
	if ( config_simdetails ) {
366
		show_details();
Thomas White's avatar
Thomas White committed
367
368
369
		return 0;
	}

370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
	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
392
393
394
395
396
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

397
398
399
400
401
402
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

403
	if ( intfile == NULL ) {
404

405
406
407
408
409
410
		/* 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;
411
412
		flags = NULL;

413
	} else {
414
415

		int i;
Thomas White's avatar
Thomas White committed
416
		ReflItemList *items;
417

418
419
420
421
422
		if ( grad == GRADIENT_PHASED ) {
			phases = new_list_phase();
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
423
424
		intensities = new_list_intensity();
		phases = new_list_phase();
425
		flags = new_list_flag();
426
427
		items = read_reflections(intfile, intensities, phases,
		                         NULL, NULL);
428
		free(intfile);
429
430
431
432
433
434

		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
435
436
437
438
439
440
441
		/* 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
442
		delete_items(items);
443
444
	}

445
446
447
448
449
450
451
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

452
453
454
455
456
457
458
	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
459
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
460
461
	image.width = image.det->max_fs + 1;
	image.height = image.det->max_ss + 1;
462
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
463
464
465
466

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
467
468
		exit(1);
	}
Thomas White's avatar
Thomas White committed
469
470

	/* Initialise stuff */
471
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
472
473
	image.features = NULL;
	image.flags = NULL;
Thomas White's avatar
Thomas White committed
474
475
	image.i0 = 1.0;
	image.i0_available = 1;
Thomas White's avatar
Thomas White committed
476

477
478
	powder = calloc(image.width*image.height, sizeof(*powder));

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

482
483
	do {

Thomas White's avatar
Thomas White committed
484
		int na, nb, nc;
485
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
486
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
487

488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
		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
506

507
508
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
509
			orientation = random_quaternion();
510
		} else {
Thomas White's avatar
Thomas White committed
511
			orientation = read_quaternion();
512
513
		}

Thomas White's avatar
Thomas White committed
514
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
515
516
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
517

Thomas White's avatar
Thomas White committed
518
		if ( !quaternion_valid(orientation) ) {
519
			ERROR("Orientation modulus is not zero!\n");
520
521
522
			return 1;
		}

Thomas White's avatar
Thomas White committed
523
524
		cell = cell_rotate(input_cell, orientation);

525
526
527
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
528

529
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
Thomas White's avatar
Thomas White committed
530
531
		STATUS("Particle size = %i x %i x %i"
		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
532
533
	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);

534
		if ( config_gpu ) {
535
			if ( gctx == NULL ) {
536
				gctx = setup_gpu(config_nosfac,
537
538
				                 intensities, flags, sym,
				                 gpu_dev);
539
			}
540
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
541
		} else {
Thomas White's avatar
Thomas White committed
542
			get_diffraction(&image, na, nb, nc, intensities, phases,
Thomas White's avatar
Thomas White committed
543
			                flags, cell, grad, sym);
544
		}
545
		if ( image.data == NULL ) {
546
547
548
549
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

550
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
551

Thomas White's avatar
Thomas White committed
552
		if ( config_nearbragg ) {
553

Thomas White's avatar
Thomas White committed
554
			/* Do EITHER: */
555

Thomas White's avatar
Thomas White committed
556
557
558
559
560
			//image.div = beam->divergence;
			//image.bw = beam->bandwidth;
			//image.profile_radius = 0.0001e9;
			//image.reflections = find_intersections(&image,
			//                               image.indexed_cell, 0);
561

Thomas White's avatar
Thomas White committed
562
563
			image.reflections = find_projected_peaks(&image,
			                            image.indexed_cell, 0, 0.1);
564

Thomas White's avatar
Thomas White committed
565
566
567
568
569
570
			integrate_reflections(&image, 0, 0);

			/* OR */

			//image.reflections = integrate_pixels(&image, 0, 0.1,
			//                                     config_polar);
571

Thomas White's avatar
Thomas White committed
572
573
		}

574
		if ( powder_fn != NULL ) {
575
576
577
578
579
580
581

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
582
				powder[x+w*y] += (double)image.data[x+w*y];
583
584
585
586
			}
			}

			if ( !(ndone % 10) ) {
587
				hdf5_write(powder_fn, powder,
588
				           image.width, image.height,
589
				           H5T_NATIVE_DOUBLE);
590
591
592
			}
		}

593
594
595
596
		if ( !config_noimages ) {

			char filename[1024];

597
598
599
600
601
602
603
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

604
			number++;
Thomas White's avatar
Thomas White committed
605

606
607
			/* Write the output file */
			hdf5_write(filename, image.data,
608
			           image.width, image.height, H5T_NATIVE_FLOAT);
609
610

		}
611

612
613
614
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
615
		cell_free(cell);
Thomas White's avatar
Thomas White committed
616

617
skip:
Thomas White's avatar
Thomas White committed
618
		ndone++;
619

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

622
	} while ( !done );
623

624
625
626
627
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
628
629
	free(image.det->panels);
	free(image.det);
630
	free(image.beam);
631
	free(powder);
Thomas White's avatar
Thomas White committed
632
	cell_free(input_cell);
633
	free(intensities);
Thomas White's avatar
Thomas White committed
634
635
	free(outfile);
	free(filename);
636
	free(sym);
637

638
	return 0;
639
}