get_hkl.c 2.95 KB
Newer Older
1
2
3
4
5
/*
 * get_hkl.c
 *
 * Small program to write out a list of h,k,l,I values given a structure
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2010 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>

#include "utils.h"
#include "sfac.h"
#include "reflections.h"


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Write idealised intensity lists.\n"
"\n"
Thomas White's avatar
Thomas White committed
35
36
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
37
38
39
"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
"      --poisson              Simulate Poisson samples.\n"
"  -o  --output=<filename>    Output filename (default: stdout).\n");
Thomas White's avatar
Thomas White committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
}


static double *template_reflections(double *ref, const char *filename)
{
	char *rval;
	double *out;
	FILE *fh;

	fh = fopen(filename, "r");
	if ( fh == NULL ) {
		return NULL;
	}

	out = new_list_intensity();

	do {

		char line[1024];
		double val;
		int r;
		signed int h, k, l;

		rval = fgets(line, 1023, fh);

		r = sscanf(line, "%i %i %i", &h, &k, &l);
		if ( r != 3 ) continue;

		val = lookup_intensity(ref, h, k, l);
		set_intensity(out, h, k, l, val);

	} while ( rval != NULL );

	fclose(fh);

	return out;
}


/* Apply Poisson noise to all reflections */
static void noisify_reflections(double *ref)
{
	signed int h, k, l;

	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {
	for ( l=-INDMAX; l<INDMAX; l++ ) {

		double val;
Thomas White's avatar
Thomas White committed
89
		int c;
Thomas White's avatar
Thomas White committed
90
91

		val = lookup_intensity(ref, h, k, l);
Thomas White's avatar
Thomas White committed
92
93
		c = poisson_noise(val);
		set_intensity(ref, h, k, l, c);
Thomas White's avatar
Thomas White committed
94
95
96

	}
	}
Thomas White's avatar
Thomas White committed
97
	progress_bar(h+INDMAX, 2*INDMAX, "Simulating noise");
Thomas White's avatar
Thomas White committed
98
	}
99
100
101
102
103
104
}


int main(int argc, char *argv[])
{
	int c;
Thomas White's avatar
Thomas White committed
105
	double *ref;
106
	struct molecule *mol;
Thomas White's avatar
Thomas White committed
107
108
	char *template = NULL;
	int config_noisify = 0;
Thomas White's avatar
Thomas White committed
109
	char *output = NULL;
110
111
112
113

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
114
115
		{"template",           1, NULL,               't'},
		{"poisson",            0, &config_noisify,     1},
Thomas White's avatar
Thomas White committed
116
		{"output",             1, NULL,               'o'},
117
118
119
120
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
121
	while ((c = getopt_long(argc, argv, "ht:o:", longopts, NULL)) != -1) {
122
123
124
125
126
127
128

		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}

Thomas White's avatar
Thomas White committed
129
130
131
132
133
		case 't' : {
			template = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
134
135
136
137
138
		case 'o' : {
			output = strdup(optarg);
			break;
		}

139
140
141
142
143
144
145
146
147
148
149
150
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

	mol = load_molecule();
151
	get_reflections_cached(mol, eV_to_J(1.8e3));
152
	ref = ideal_intensities(mol->reflections);
Thomas White's avatar
Thomas White committed
153
154

	if ( template != NULL ) {
Thomas White's avatar
Thomas White committed
155
156
157

		double *tref;

Thomas White's avatar
Thomas White committed
158
159
160
161
162
163
164
		tref = template_reflections(ref, template);
		if ( tref == NULL ) {
			ERROR("Couldn't read template file!\n");
			return 1;
		}
		free(ref);
		ref = tref;
Thomas White's avatar
Thomas White committed
165

Thomas White's avatar
Thomas White committed
166
167
168
169
	}

	if ( config_noisify ) noisify_reflections(ref);

Thomas White's avatar
Thomas White committed
170
	write_reflections(output, NULL, ref, 1, mol->cell);
171
172
173

	return 0;
}