pattern_sim.c 16.3 KB
Newer Older
1
/*
Thomas White's avatar
Thomas White committed
2
3
4
 * pattern_sim.c
 *
 * Simulate diffraction patterns from small crystals
5
 *
6
7
8
9
10
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
 *   2009-2012 Thomas White <taw@physics.org>
11
 *
Thomas White's avatar
Thomas White committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 * 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/>.
26
27
28
29
30
31
32
33
34
35
36
37
38
 *
 */


#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
39
#include <getopt.h>
40

41
#include "image.h"
Thomas White's avatar
Thomas White committed
42
#include "diffraction.h"
43
#include "diffraction-gpu.h"
44
#include "cell.h"
45
#include "cell-utils.h"
46
47
#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
48
#include "detector.h"
49
#include "peaks.h"
50
#include "beam-parameters.h"
Thomas White's avatar
Thomas White committed
51
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
52
53
#include "reflist.h"
#include "reflist-utils.h"
54
#include "pattern_sim.h"
55
56


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


117
118
119
120
static double *intensities_from_list(RefList *list)
{
	Reflection *refl;
	RefListIterator *iter;
121
	double *out = new_arr_intensity();
122
123
124
125
126
127
128
129
130
131

	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);

132
		set_arr_intensity(out, h, k, l, intensity);
133
134
135
136
137
138
139
140
141
142
143

	}

	return out;
}


static double *phases_from_list(RefList *list)
{
	Reflection *refl;
	RefListIterator *iter;
144
	double *out = new_arr_phase();
145
146
147
148
149
150
151
152
153
154

	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);

155
		set_arr_phase(out, h, k, l, phase);
156
157
158
159
160
161
162
163
164
165
166
167

	}

	return out;

}


static unsigned char *flags_from_list(RefList *list)
{
	Reflection *refl;
	RefListIterator *iter;
168
	unsigned char *out = new_arr_flag();
169
170
171
172
173
174
175
176
177

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

		signed int h, k, l;

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

178
		set_arr_flag(out, h, k, l, 1);
179
180
181
182
183
184
185
186

	}

	return out;

}


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
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 {
214
			ERROR("Invalid rotation '%s'\n", line);
215
216
217
		}

	} while ( 1 );
218
219
220
}


221
222
223
224
225
226
227
228
229
230
231
232
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;
}


233
234
int main(int argc, char *argv[])
{
235
	int c;
236
	struct image image;
237
	struct gpu_context *gctx = NULL;
238
	double *powder;
239
240
	char *intfile = NULL;
	double *intensities;
241
	char *rval;
242
	double *phases;
243
	unsigned char *flags;
244
	int config_randomquat = 0;
245
	int config_noimages = 0;
Thomas White's avatar
Thomas White committed
246
	int config_nonoise = 0;
247
	int config_nosfac = 0;
248
	int config_gpu = 0;
249
	int config_random = 0;
250
	char *powder_fn = NULL;
251
	char *filename = NULL;
252
	char *grad_str = NULL;
253
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
254
	char *geometry = NULL;
255
	char *beamfile = NULL;
256
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
257
258
259
	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 */
260
	int done = 0;
Thomas White's avatar
Thomas White committed
261
262
	UnitCell *input_cell;
	struct quaternion orientation;
Thomas White's avatar
Thomas White committed
263
	int gpu_dev = -1;
264
265
266
	int random_size = 0;
	double min_size = 0.0;
	double max_size = 0.0;
Thomas White's avatar
Thomas White committed
267
268
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
269

270
	/* Long options */
Thomas White's avatar
Thomas White committed
271
	const struct option longopts[] = {
272
		{"help",               0, NULL,               'h'},
273
		{"gpu",                0, &config_gpu,         1},
274
275
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
276
		{"no-images",          0, &config_noimages,    1},
Thomas White's avatar
Thomas White committed
277
		{"no-noise",           0, &config_nonoise,     1},
278
		{"intensities",        1, NULL,               'i'},
279
		{"symmetry",           1, NULL,               'y'},
280
		{"powder",             1, NULL,               'w'},
Thomas White's avatar
Thomas White committed
281
		{"gradients",          1, NULL,               't'},
282
		{"pdb",                1, NULL,               'p'},
283
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
284
		{"geometry",           1, NULL,               'g'},
285
		{"beam",               1, NULL,               'b'},
286
		{"really-random",      0, &config_random,      1},
Thomas White's avatar
Thomas White committed
287
		{"gpu-dev",            1, NULL,                2},
288
289
		{"min-size",           1, NULL,                3},
		{"max-size",           1, NULL,                4},
290
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
291
	};
292

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

Thomas White's avatar
Thomas White committed
297
		switch (c) {
Thomas White's avatar
Thomas White committed
298
299

			case 'h' :
Thomas White's avatar
Thomas White committed
300
301
			show_help(argv[0]);
			return 0;
302

Thomas White's avatar
Thomas White committed
303
			case 'r' :
304
305
306
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
307
			case 'n' :
308
309
310
311
312
			n_images = strtol(optarg, &rval, 10);
			if ( *rval != '\0' ) {
				ERROR("Invalid number of images.\n");
				return 1;
			}
313
314
			break;

Thomas White's avatar
Thomas White committed
315
			case 'i' :
316
317
318
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
319
			case 't' :
320
321
322
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
323
			case 'p' :
324
325
326
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
327
			case 'o' :
328
329
330
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
331
			case 'w' :
332
333
334
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
335
			case 'g' :
Thomas White's avatar
Thomas White committed
336
337
338
			geometry = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
339
			case 'b' :
340
341
342
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
343
			case 'y' :
Thomas White's avatar
Thomas White committed
344
			sym_str = strdup(optarg);
345
346
			break;

Thomas White's avatar
Thomas White committed
347
			case 2 :
Thomas White's avatar
Thomas White committed
348
349
350
			gpu_dev = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
351
			case 3 :
352
353
354
355
356
			min_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid minimum size.\n");
				return 1;
			}
357
358
359
			random_size++;
			break;

Thomas White's avatar
Thomas White committed
360
			case 4 :
361
362
363
364
365
			max_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid maximum size.\n");
				return 1;
			}
366
367
368
			random_size++;
			break;

Thomas White's avatar
Thomas White committed
369
			case 0 :
Thomas White's avatar
Thomas White committed
370
			break;
371

372
373
374
			case '?' :
			break;

Thomas White's avatar
Thomas White committed
375
			default :
376
377
			ERROR("Unhandled option '%c'\n", c);
			break;
Thomas White's avatar
Thomas White committed
378

Thomas White's avatar
Thomas White committed
379
		}
380
381
382

	}

383
	if ( config_random ) {
384
385
386
387
388
389
390
391
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

392
393
394
395
396
	if ( random_size == 1 ) {
		ERROR("You must specify both --min-size and --max-size.\n");
		return 1;
	}

397
398
399
400
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

401
402
403
404
405
406
407
408
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
409
410
411
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	/* sym_str is used below */
412

413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
	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
435
436
437
438
439
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

440
441
442
443
444
445
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

446
	if ( intfile == NULL ) {
447

448
449
450
451
452
453
		/* 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;
454
455
		flags = NULL;

456
	} else {
457

Thomas White's avatar
Thomas White committed
458
459
460
		RefList *reflections;

		reflections = read_reflections(intfile);
Thomas White's avatar
Thomas White committed
461
462
463
464
		if ( reflections == NULL ) {
			ERROR("Problem reading input file %s\n", intfile);
			return 1;
		}
Thomas White's avatar
Thomas White committed
465
		free(intfile);
466

467
		if ( grad == GRADIENT_PHASED ) {
Thomas White's avatar
Thomas White committed
468
			phases = phases_from_list(reflections);
469
470
471
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
472
473
474
		intensities = intensities_from_list(reflections);
		phases = phases_from_list(reflections);
		flags = flags_from_list(reflections);
475

Thomas White's avatar
Thomas White committed
476
		/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
477
		if ( check_list_symmetry(reflections, sym) ) {
Thomas White's avatar
Thomas White committed
478
			ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
479
			      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
480
481
482
			return 1;
		}

Thomas White's avatar
Thomas White committed
483
484
		reflist_free(reflections);

485
486
	}

487
488
489
490
491
492
493
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

494
495
496
497
498
499
	image.beam = get_beam_parameters(beamfile);
	if ( image.beam == NULL ) {
		ERROR("Failed to read beam parameters from '%s'\n", beamfile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
500
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
501
502
	image.width = image.det->max_fs + 1;
	image.height = image.det->max_ss + 1;
503
504
505
	if ( image.beam->photon_energy_from != NULL ) {
		ERROR("Photon energy must be specified, not taken from the"
		      " HDF5 file.  Please alter %s accordingly.\n", beamfile)
506
507
		return 1;
	} else {
508
509
		double wl = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
		image.lambda = wl;
510
	}
511
512
	image.bw = image.beam->bandwidth;
	image.div = image.beam->divergence;
513
	free(beamfile);
Thomas White's avatar
Thomas White committed
514
515
516
517

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
518
519
		exit(1);
	}
Thomas White's avatar
Thomas White committed
520
521

	/* Initialise stuff */
522
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
523
524
	image.features = NULL;
	image.flags = NULL;
Thomas White's avatar
Thomas White committed
525

526
527
	powder = calloc(image.width*image.height, sizeof(*powder));

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

531
532
	do {

Thomas White's avatar
Thomas White committed
533
		int na, nb, nc;
534
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
535
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
536

537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
		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
555

556
557
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
558
			orientation = random_quaternion();
559
		} else {
Thomas White's avatar
Thomas White committed
560
			orientation = read_quaternion();
561
562
		}

Thomas White's avatar
Thomas White committed
563
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
564
565
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
566

Thomas White's avatar
Thomas White committed
567
		if ( !quaternion_valid(orientation) ) {
568
			ERROR("Orientation modulus is not zero!\n");
569
570
571
			return 1;
		}

Thomas White's avatar
Thomas White committed
572
573
		cell = cell_rotate(input_cell, orientation);

574
575
576
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
577

578
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
Thomas White's avatar
Thomas White committed
579
580
		STATUS("Particle size = %i x %i x %i"
		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
581
582
	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);

583
		if ( config_gpu ) {
584
			if ( gctx == NULL ) {
585
				gctx = setup_gpu(config_nosfac,
Thomas White's avatar
Thomas White committed
586
				                 intensities, flags, sym_str,
587
				                 gpu_dev);
588
			}
589
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
590
		} else {
Thomas White's avatar
Thomas White committed
591
			get_diffraction(&image, na, nb, nc, intensities, phases,
Thomas White's avatar
Thomas White committed
592
			                flags, cell, grad, sym);
593
		}
594
		if ( image.data == NULL ) {
595
596
597
598
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

599
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
600

601
		if ( powder_fn != NULL ) {
602
603
604
605
606
607
608

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
609
				powder[x+w*y] += (double)image.data[x+w*y];
610
611
612
613
			}
			}

			if ( !(ndone % 10) ) {
614
				hdf5_write(powder_fn, powder,
615
				           image.width, image.height,
616
				           H5T_NATIVE_DOUBLE);
617
618
619
			}
		}

620
621
622
623
		if ( !config_noimages ) {

			char filename[1024];

624
625
626
627
628
629
630
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

631
			number++;
Thomas White's avatar
Thomas White committed
632

633
			/* Write the output file */
634
			hdf5_write_image(filename, &image);
635
636

		}
637

638
639
640
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
641
		cell_free(cell);
Thomas White's avatar
Thomas White committed
642

643
skip:
Thomas White's avatar
Thomas White committed
644
		ndone++;
645

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

648
	} while ( !done );
649

650
651
652
653
654
655
	if ( powder_fn != NULL ) {
		hdf5_write(powder_fn, powder,
		           image.width, image.height,
		           H5T_NATIVE_DOUBLE);
	}

656
657
658
659
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
660
661
	free(image.det->panels);
	free(image.det);
662
	free(image.beam);
663
	free(powder);
Thomas White's avatar
Thomas White committed
664
	cell_free(input_cell);
665
	free(intensities);
Thomas White's avatar
Thomas White committed
666
667
	free(outfile);
	free(filename);
Thomas White's avatar
Thomas White committed
668
669
	free(sym_str);
	free_symoplist(sym);
670

671
	return 0;
672
}