pattern_sim.c 17.7 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
560
	powder = calloc(image.width*image.height, sizeof(*powder));

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

564
565
	do {

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

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

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

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

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

Thomas White's avatar
Thomas White committed
605
606
		cell = cell_rotate(input_cell, orientation);

607
608
609
		/* Ensure no residual information */
		image.data = NULL;
		image.twotheta = NULL;
610

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

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

632
		record_image(&image, !config_nonoise);
Thomas White's avatar
Thomas White committed
633

634
		if ( powder_fn != NULL ) {
635
636
637
638
639
640
641

			int x, y, w;

			w = image.width;

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

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

653
654
655
656
		if ( !config_noimages ) {

			char filename[1024];

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

664
			number++;
Thomas White's avatar
Thomas White committed
665

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

		}
671

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

677
skip:
Thomas White's avatar
Thomas White committed
678
		ndone++;
679

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

682
	} while ( !done );
683

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

690
691
692
693
	if ( gctx != NULL ) {
		cleanup_gpu(gctx);
	}

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

705
	return 0;
706
}