aboutsummaryrefslogtreecommitdiff
path: root/src/index.c
blob: 470d9e4c3bbf49f4a23915ebf7f9ded387350762 (plain)
1
2
3
4
5
6
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
35
36
37
38
39
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
/*
 * index.c
 *
 * Perform indexing (somehow)
 *
 * (c) 2006-2010 Thomas White <taw@physics.org>
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


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

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <assert.h>

#include "image.h"
#include "utils.h"
#include "peaks.h"
#include "dirax.h"
#include "sfac.h"
#include "detector.h"
#include "index.h"



static void write_drx(struct image *image)
{
	FILE *fh;
	int i;

	STATUS("Writing xfel.drx file.  Remember that it uses units of "
	       "reciprocal Angstroms!\n");

	fh = fopen("xfel.drx", "w");
	if ( !fh ) {
		ERROR("Couldn't open temporary file xfel.drx\n");
		return;
	}
	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  */

	for ( i=0; i<image_feature_count(image->features); i++ ) {

		struct imagefeature *f;

		f = image_get_feature(image->features, i);
		if ( f == NULL ) continue;

		fprintf(fh, "%10f %10f %10f %8f\n",
		        f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0);

	}
	fclose(fh);
}


void index_pattern(struct image *image, IndexingMethod indm, int no_match,
                   int verbose)
{
	int i;
	UnitCell *new_cell = NULL;
	int nc = 0;

	/* Map positions to 3D */
	for ( i=0; i<image_feature_count(image->features); i++ ) {

		struct imagefeature *f;
		int c;

		f = image_get_feature(image->features, i);
		if ( f == NULL ) continue;

		c = map_position(image, f->x, f->y, &f->rx, &f->ry, &f->rz);
		if ( c != 0 ) nc++;

	}
	if ( nc ) {
		ERROR("Failed to map %i reflections\n", nc);
	}

	write_drx(image);

	/* Index (or not) as appropriate */
	if ( indm == INDEXING_NONE ) return;
	if ( indm == INDEXING_DIRAX ) run_dirax(image);

	if ( image->indexed_cell == NULL ) {
		STATUS("No cell found.\n");
		return;
	} else if ( verbose ) {
		STATUS("--------------------\n");
		STATUS("The indexed cell (before matching):\n");
		cell_print(image->indexed_cell);
		STATUS("--------------------\n");
	}

	if ( !no_match ) {
		new_cell = match_cell(image->indexed_cell,
			              image->molecule->cell, verbose);
		free(image->indexed_cell);
		image->indexed_cell = new_cell;
	}
}