pattern_sim.c 17.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
 * Copyright © 2012 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * 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/>.
22
23
24
25
26
27
28
29
30
31
32
33
34
 *
 */


#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
35
#include <getopt.h>
36

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


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


118
static void show_details()
Thomas White's avatar
Thomas White committed
119
{
Thomas White's avatar
Thomas White committed
120
121
122
123
124
	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"
125
126
"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
127
"\n"
128
129
130
"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
131
132
133
134
"\n"
"na = number of unit cells in 'a' direction (likewise nb, nc)\n"
" q = reciprocal vector (1/d convention, not 2pi/d)\n"
"\n"
135
136
137
"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
138
139
140
"\n"
"Expected intensities at the CCD are then calculated using:\n"
"\n"
141
"I(q) = I0 * r^2 * I_latt(q) * I_mol(q) * S\n"
Thomas White's avatar
Thomas White committed
142
143
144
145
"\n"
"I0 = number of photons per unit area in the incident beam\n"
" r = Thomson radius\n"
" S = solid angle of corresponding pixel\n"
146
147
148
"\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
149
150
"\n"
"Poisson counts are generated from the expected intensities using Knuth's\n"
Thomas White's avatar
Thomas White committed
151
152
153
"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"
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
215
216
217
218
219
220
221
222
223
224
225
226
227
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;

}


228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
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 {
255
			ERROR("Invalid rotation '%s'\n", line);
256
257
258
		}

	} while ( 1 );
259
260
261
}


262
263
264
265
266
267
268
269
270
271
272
273
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;
}


274
275
int main(int argc, char *argv[])
{
276
	int c;
277
	struct image image;
278
	struct gpu_context *gctx = NULL;
279
	double *powder;
280
281
	char *intfile = NULL;
	double *intensities;
282
	char *rval;
283
	double *phases;
284
	unsigned char *flags;
Thomas White's avatar
Thomas White committed
285
	int config_simdetails = 0;
286
	int config_randomquat = 0;
287
	int config_noimages = 0;
Thomas White's avatar
Thomas White committed
288
	int config_nonoise = 0;
289
	int config_nosfac = 0;
290
	int config_gpu = 0;
291
	int config_random = 0;
292
	char *powder_fn = NULL;
293
	char *filename = NULL;
294
	char *grad_str = NULL;
295
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
296
	char *geometry = NULL;
297
	char *beamfile = NULL;
298
	GradientMethod grad;
Thomas White's avatar
Thomas White committed
299
300
301
	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 */
302
	int done = 0;
Thomas White's avatar
Thomas White committed
303
304
	UnitCell *input_cell;
	struct quaternion orientation;
Thomas White's avatar
Thomas White committed
305
	int gpu_dev = -1;
306
307
308
	int random_size = 0;
	double min_size = 0.0;
	double max_size = 0.0;
Thomas White's avatar
Thomas White committed
309
310
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
311

312
	/* Long options */
Thomas White's avatar
Thomas White committed
313
	const struct option longopts[] = {
314
315
		{"help",               0, NULL,               'h'},
		{"simulation-details", 0, &config_simdetails,  1},
316
		{"gpu",                0, &config_gpu,         1},
317
318
		{"random-orientation", 0, NULL,               'r'},
		{"number",             1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
319
		{"no-images",          0, &config_noimages,    1},
Thomas White's avatar
Thomas White committed
320
		{"no-noise",           0, &config_nonoise,     1},
321
		{"intensities",        1, NULL,               'i'},
322
		{"symmetry",           1, NULL,               'y'},
323
		{"powder",             1, NULL,               'w'},
Thomas White's avatar
Thomas White committed
324
		{"gradients",          1, NULL,               't'},
325
		{"pdb",                1, NULL,               'p'},
326
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
327
		{"geometry",           1, NULL,               'g'},
328
		{"beam",               1, NULL,               'b'},
329
		{"really-random",      0, &config_random,      1},
Thomas White's avatar
Thomas White committed
330
		{"gpu-dev",            1, NULL,                2},
331
332
		{"min-size",           1, NULL,                3},
		{"max-size",           1, NULL,                4},
333
		{0, 0, NULL, 0}
Thomas White's avatar
Thomas White committed
334
	};
335

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

Thomas White's avatar
Thomas White committed
340
		switch (c) {
Thomas White's avatar
Thomas White committed
341
		case 'h' :
Thomas White's avatar
Thomas White committed
342
343
			show_help(argv[0]);
			return 0;
344

Thomas White's avatar
Thomas White committed
345
		case 'r' :
346
347
348
			config_randomquat = 1;
			break;

Thomas White's avatar
Thomas White committed
349
		case 'n' :
350
351
352
353
354
			n_images = strtol(optarg, &rval, 10);
			if ( *rval != '\0' ) {
				ERROR("Invalid number of images.\n");
				return 1;
			}
355
356
			break;

Thomas White's avatar
Thomas White committed
357
		case 'i' :
358
359
360
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
361
		case 't' :
362
363
364
			grad_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
365
		case 'p' :
366
367
368
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
369
		case 'o' :
370
371
372
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
373
		case 'w' :
374
375
376
			powder_fn = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
377
		case 'g' :
Thomas White's avatar
Thomas White committed
378
379
380
			geometry = strdup(optarg);
			break;

381
382
383
384
		case 'b' :
			beamfile = strdup(optarg);
			break;

385
		case 'y' :
Thomas White's avatar
Thomas White committed
386
			sym_str = strdup(optarg);
387
388
			break;

Thomas White's avatar
Thomas White committed
389
390
391
392
		case 2 :
			gpu_dev = atoi(optarg);
			break;

393
		case 3 :
394
395
396
397
398
			min_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid minimum size.\n");
				return 1;
			}
399
400
401
402
			random_size++;
			break;

		case 4 :
403
404
405
406
407
			max_size = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid maximum size.\n");
				return 1;
			}
408
409
410
			random_size++;
			break;

Thomas White's avatar
Thomas White committed
411
		case 0 :
Thomas White's avatar
Thomas White committed
412
			break;
413

Thomas White's avatar
Thomas White committed
414
		default :
Thomas White's avatar
Thomas White committed
415
416
			return 1;
		}
417
418
419

	}

420
	if ( config_random ) {
421
422
423
424
425
426
427
428
		FILE *fh;
		unsigned int seed;
		fh = fopen("/dev/urandom", "r");
		fread(&seed, sizeof(seed), 1, fh);
		fclose(fh);
		srand(seed);
	}

429
430
431
432
433
	if ( random_size == 1 ) {
		ERROR("You must specify both --min-size and --max-size.\n");
		return 1;
	}

434
435
436
437
	if ( filename == NULL ) {
		filename = strdup("molecule.pdb");
	}

438
439
440
441
442
443
444
445
	if ( outfile == NULL ) {
		if ( n_images == 1 ) {
			outfile = strdup("sim.h5");
		} else {
			outfile = strdup("sim");
		}
	}

Thomas White's avatar
Thomas White committed
446
447
448
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	/* sym_str is used below */
449

Thomas White's avatar
Thomas White committed
450
	if ( config_simdetails ) {
451
		show_details();
Thomas White's avatar
Thomas White committed
452
453
454
		return 0;
	}

455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
	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
477
478
479
480
481
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

482
483
484
485
486
487
	if ( beamfile == NULL ) {
		ERROR("You need to specify a beam parameter file"
		      " with --beam\n");
		return 1;
	}

488
	if ( intfile == NULL ) {
489

490
491
492
493
494
495
		/* 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;
496
497
		flags = NULL;

498
	} else {
499

Thomas White's avatar
Thomas White committed
500
501
502
		RefList *reflections;

		reflections = read_reflections(intfile);
Thomas White's avatar
Thomas White committed
503
504
505
506
		if ( reflections == NULL ) {
			ERROR("Problem reading input file %s\n", intfile);
			return 1;
		}
Thomas White's avatar
Thomas White committed
507
		free(intfile);
508

509
		if ( grad == GRADIENT_PHASED ) {
Thomas White's avatar
Thomas White committed
510
			phases = phases_from_list(reflections);
511
512
513
		} else {
			phases = NULL;
		}
Thomas White's avatar
Thomas White committed
514
515
516
		intensities = intensities_from_list(reflections);
		phases = phases_from_list(reflections);
		flags = flags_from_list(reflections);
517

Thomas White's avatar
Thomas White committed
518
		/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
519
		if ( check_list_symmetry(reflections, sym) ) {
Thomas White's avatar
Thomas White committed
520
			ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
521
			      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
522
523
524
			return 1;
		}

Thomas White's avatar
Thomas White committed
525
526
		reflist_free(reflections);

527
528
	}

529
530
531
532
533
534
535
	image.det = get_detector_geometry(geometry);
	if ( image.det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

536
537
538
539
540
541
542
	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
543
	/* Define image parameters */
Thomas White's avatar
Thomas White committed
544
545
	image.width = image.det->max_fs + 1;
	image.height = image.det->max_ss + 1;
546
	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
Thomas White's avatar
Thomas White committed
547
548
549
550

	/* Load unit cell */
	input_cell = load_cell_from_pdb(filename);
	if ( input_cell == NULL ) {
Thomas White's avatar
Thomas White committed
551
552
		exit(1);
	}
Thomas White's avatar
Thomas White committed
553
554

	/* Initialise stuff */
555
	image.filename = NULL;
Thomas White's avatar
Thomas White committed
556
557
	image.features = NULL;
	image.flags = NULL;
Thomas White's avatar
Thomas White committed
558
559
	image.i0 = 1.0;
	image.i0_available = 1;
Thomas White's avatar
Thomas White committed
560

561
562
	powder = calloc(image.width*image.height, sizeof(*powder));

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

566
567
	do {

Thomas White's avatar
Thomas White committed
568
		int na, nb, nc;
569
		double a, b, c, d;
Thomas White's avatar
Thomas White committed
570
		UnitCell *cell;
Thomas White's avatar
Thomas White committed
571

572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
		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
590

591
592
		/* Read quaternion from stdin */
		if ( config_randomquat ) {
Thomas White's avatar
Thomas White committed
593
			orientation = random_quaternion();
594
		} else {
Thomas White's avatar
Thomas White committed
595
			orientation = read_quaternion();
596
597
		}

Thomas White's avatar
Thomas White committed
598
		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
Thomas White's avatar
Thomas White committed
599
600
		       orientation.w, orientation.x,
		       orientation.y, orientation.z);
Thomas White's avatar
Thomas White committed
601

Thomas White's avatar
Thomas White committed
602
		if ( !quaternion_valid(orientation) ) {
603
			ERROR("Orientation modulus is not zero!\n");
604
605
606
			return 1;
		}

Thomas White's avatar
Thomas White committed
607
608
		cell = cell_rotate(input_cell, orientation);

609
610
611
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
612

613
		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
Thomas White's avatar
Thomas White committed
614
615
		STATUS("Particle size = %i x %i x %i"
		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
616
617
	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);

618
		if ( config_gpu ) {
619
			if ( gctx == NULL ) {
620
				gctx = setup_gpu(config_nosfac,
Thomas White's avatar
Thomas White committed
621
				                 intensities, flags, sym_str,
622
				                 gpu_dev);
623
			}
624
			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
625
		} else {
Thomas White's avatar
Thomas White committed
626
			get_diffraction(&image, na, nb, nc, intensities, phases,
Thomas White's avatar
Thomas White committed
627
			                flags, cell, grad, sym);
628
		}
629
		if ( image.data == NULL ) {
630
631
632
633
			ERROR("Diffraction calculation failed.\n");
			goto skip;
		}

634
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
635

636
		if ( powder_fn != NULL ) {
637
638
639
640
641
642
643

			int x, y, w;

			w = image.width;

			for ( x=0; x<image.width; x++ ) {
			for ( y=0; y<image.height; y++ ) {
644
				powder[x+w*y] += (double)image.data[x+w*y];
645
646
647
648
			}
			}

			if ( !(ndone % 10) ) {
649
				hdf5_write(powder_fn, powder,
650
				           image.width, image.height,
651
				           H5T_NATIVE_DOUBLE);
652
653
654
			}
		}

655
656
657
658
		if ( !config_noimages ) {

			char filename[1024];

659
660
661
662
663
664
665
			if ( n_images != 1 ) {
				snprintf(filename, 1023, "%s-%i.h5",
				         outfile, number);
			} else {
				strncpy(filename, outfile, 1023);
			}

666
			number++;
Thomas White's avatar
Thomas White committed
667

668
669
			/* Write the output file */
			hdf5_write(filename, image.data,
670
			           image.width, image.height, H5T_NATIVE_FLOAT);
671
672

		}
673

674
675
676
		/* Clean up */
		free(image.data);
		free(image.twotheta);
Thomas White's avatar
Thomas White committed
677
		cell_free(cell);
Thomas White's avatar
Thomas White committed
678

679
skip:
Thomas White's avatar
Thomas White committed
680
		ndone++;
681

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

684
	} while ( !done );
685

686
687
688
689
690
691
	if ( powder_fn != NULL ) {
		hdf5_write(powder_fn, powder,
		           image.width, image.height,
		           H5T_NATIVE_DOUBLE);
	}

692
693
694
695
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

Thomas White's avatar
Thomas White committed
696
697
	free(image.det->panels);
	free(image.det);
698
	free(image.beam);
699
	free(powder);
Thomas White's avatar
Thomas White committed
700
	cell_free(input_cell);
701
	free(intensities);
Thomas White's avatar
Thomas White committed
702
703
	free(outfile);
	free(filename);
Thomas White's avatar
Thomas White committed
704
705
	free(sym_str);
	free_symoplist(sym);
706

707
	return 0;
708
}