aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2017-02-02 15:08:59 +0100
committerThomas White <taw@physics.org>2017-02-02 15:08:59 +0100
commit6ece93b3ab590da9da08fdc42e8177ca64056145 (patch)
tree3de5cd3bc4d193a40d28d5c19cc4ca85dd53c3b0 /libcrystfel
parent2cefd4a8807d87cbe459c83d4b944d32ee38d2af (diff)
parent80b54c66c36c46a247d6d015be60ea9e73aa025f (diff)
Merge branch 'master' into taketwo
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/asdf.c45
-rw-r--r--libcrystfel/src/geometry.c26
-rw-r--r--libcrystfel/src/index.c13
-rw-r--r--libcrystfel/src/index.h7
-rw-r--r--libcrystfel/src/stream.c58
-rw-r--r--libcrystfel/src/symmetry.c21
-rw-r--r--libcrystfel/src/symmetry.h8
7 files changed, 123 insertions, 55 deletions
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index df20f7ef..49917ad9 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -322,6 +322,15 @@ static int compare_doubles(const void *a, const void *b)
}
+static double max(double a, double b, double c)
+{
+ double m = a;
+ if ( m < b ) m = b;
+ if ( m < c ) m = c;
+ return m;
+}
+
+
/* Compares tvectors by length */
static int compare_tvectors(const void *a, const void *b)
{
@@ -889,12 +898,14 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit,
return 0;
}
+
void swap(int *a, int *b) {
int c = *a;
*a = *b;
*b = c;
}
+
/* Generate triplets of integers < N_reflections in random sequence */
static int **generate_triplets(int N_reflections, int N_triplets_max, int *N)
{
@@ -904,10 +915,13 @@ static int **generate_triplets(int N_reflections, int N_triplets_max, int *N)
int N_triplets = N_reflections * (N_reflections - 1) *
(N_reflections - 2) / 6;
- if ( N_triplets > N_triplets_max ) N_triplets = N_triplets_max;
+ if ( N_triplets > N_triplets_max || N_reflections > 1000 ) {
+ N_triplets = N_triplets_max;
+ }
*N = N_triplets;
int **triplets = malloc(N_triplets * sizeof(int *));
+
if (triplets == NULL) {
ERROR("Failed to allocate triplets in generate_triplets!\n");
return 0;
@@ -962,6 +976,7 @@ static int **generate_triplets(int N_reflections, int N_triplets_max, int *N)
return triplets;
}
+
static int index_refls(gsl_vector **reflections, int N_reflections,
double d_max, double volume_min, double volume_max,
double LevelFit, double IndexFit, int N_triplets_max,
@@ -1094,20 +1109,38 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
double LevelFit = 1./1000;
double IndexFit = 1./500;
- double d_max = 220.; // thrice the maximum expected axis length
+ double d_max = 1000.; // thrice the maximum expected axis length
double volume_min = 100.;
- double volume_max = 1000000.;
+ double volume_max = 100000000.;
int N_triplets_max = 10000; // maximum number of triplets
struct asdf_private *dp = (struct asdf_private *)ipriv;
- if ( dp->indm & INDEXING_CHECK_CELL_AXES ) {
+ if ( dp->indm & INDEXING_CHECK_CELL_AXES ||
+ dp->indm & INDEXING_CHECK_CELL_COMBINATIONS) {
+ double a, b, c, gamma, beta, alpha;
+ cell_get_parameters(dp->template, &a, &b, &c,
+ &alpha, &beta, &gamma);
+
+ d_max = max(a, b, c) * 3 * 1e10;
+
double volume = cell_get_volume(dp->template);
double vtol = (dp->ltl[0] + dp->ltl[1] + dp->ltl[2]) / 100 +
dp->ltl[3] / 180 * M_PI;
- volume_min = volume * (1 - vtol);
- volume_max = volume * (1 + vtol);
+
+ /* Divide volume constraints by number of lattice points per
+ * unit cell since asdf always finds primitive cell */
+ int latt_points_per_uc = 1;
+ char centering = cell_get_centering(dp->template);
+ if ( centering == 'A' ||
+ centering == 'B' ||
+ centering == 'C' ||
+ centering == 'I' ) latt_points_per_uc = 2;
+ else if ( centering == 'F' ) latt_points_per_uc = 4;
+
+ volume_min = volume * (1 - vtol)/latt_points_per_uc;
+ volume_max = volume * (1 + vtol)/latt_points_per_uc;
}
int n = image_feature_count(image->features);
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 266d6267..b968dc06 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -52,16 +52,15 @@ static int locate_peak_on_panel(double x, double y, double z, double k,
struct panel *p,
double *pfs, double *pss)
{
- double r2, ctt, tta, phi;
+ double ctt, tta, phi;
gsl_vector *v;
gsl_vector *t;
gsl_matrix *M;
double fs, ss, one_over_mu;
/* Calculate 2theta (scattering angle) and azimuth (phi) */
- r2 = x*x + y*y + z*z;
- ctt = 1.0 - r2 / (2.0*k*k); /* cos(2theta) */
- tta = acos(ctt);
+ tta = atan2(sqrt(x*x+y*y), k+z);
+ ctt = cos(tta);
phi = atan2(y, x);
/* Set up matrix equation */
@@ -637,25 +636,18 @@ void polarisation_correction(RefList *list, UnitCell *cell, struct image *image)
refl != NULL;
refl = next_refl(refl, iter) )
{
- double pol, pa, pb, phi, tt, ool;
+ double pol;
double intensity;
- double xl, yl, zl;
+ double xl, yl;
signed int h, k, l;
+ const double P = 1.0; /* degree of polarisation */
get_indices(refl, &h, &k, &l);
- /* Polarisation correction assuming 100% polarisation
- * along the x direction */
- xl = h*asx + k*bsx + l*csx;
- yl = h*asy + k*bsy + l*csy;
- zl = h*asz + k*bsz + l*csz;
+ xl = (h*asx + k*bsx + l*csx)*image->lambda;
+ yl = (h*asy + k*bsy + l*csy)*image->lambda;
- ool = 1.0 / image->lambda;
- tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool);
- phi = atan2(yl, xl);
- pa = pow(sin(phi)*sin(tt), 2.0);
- pb = pow(cos(tt), 2.0);
- pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb);
+ pol = P*(1.0 - xl*xl) + (1.0-P)*(1.0 - yl*yl);
intensity = get_intensity(refl);
set_intensity(refl, intensity / pol);
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index e0bebfa4..e7fd879d 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -3,12 +3,12 @@
*
* Perform indexing (somehow)
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2016 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
* 2010-2011 Richard Kirian <rkirian@asu.edu>
* 2012 Lorenzo Galli
* 2013 Cornelius Gati <cornelius.gati@cfel.de>
@@ -434,10 +434,16 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image)
}
}
-
void index_pattern(struct image *image,
IndexingMethod *indms, IndexingPrivate **iprivs)
{
+ index_pattern_2(image, indms, iprivs, NULL);
+}
+
+
+void index_pattern_2(struct image *image, IndexingMethod *indms,
+ IndexingPrivate **iprivs, int *ping)
+{
int n = 0;
ImageFeatureList *orig;
@@ -465,6 +471,7 @@ void index_pattern(struct image *image,
ntry++;
done = finished_retry(indms[n], r, image);
if ( ntry > 5 ) done = 1;
+ if ( ping != NULL ) (*ping)++;
} while ( !done );
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 23545506..33d568f6 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -3,13 +3,13 @@
*
* Perform indexing (somehow)
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2016 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
* 2010 Richard Kirian
* 2012 Lorenzo Galli
* 2015 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
@@ -161,6 +161,9 @@ extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
extern void index_pattern(struct image *image,
IndexingMethod *indms, IndexingPrivate **iprivs);
+extern void index_pattern_2(struct image *image, IndexingMethod *indms,
+ IndexingPrivate **iprivs, int *ping);
+
extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs);
#ifdef __cplusplus
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 789a9d7a..e858172e 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -65,9 +65,11 @@ struct _stream
int major_version;
int minor_version;
+
+ long long int ln;
};
-static int read_peaks(FILE *fh, struct image *image)
+static int read_peaks(Stream *st, struct image *image)
{
char *rval = NULL;
int first = 1;
@@ -82,7 +84,8 @@ static int read_peaks(FILE *fh, struct image *image)
struct panel *p = NULL;
float add_x, add_y;
- rval = fgets(line, 1023, fh);
+ rval = fgets(line, 1023, st->fh);
+ st->ln++;
if ( rval == NULL ) continue;
chomp(line);
@@ -126,7 +129,7 @@ static int read_peaks(FILE *fh, struct image *image)
}
-static int read_peaks_2_3(FILE *fh, struct image *image)
+static int read_peaks_2_3(Stream *st, struct image *image)
{
char *rval = NULL;
int first = 1;
@@ -142,7 +145,8 @@ static int read_peaks_2_3(FILE *fh, struct image *image)
struct panel *p = NULL;
float add_x, add_y;
- rval = fgets(line, 1023, fh);
+ rval = fgets(line, 1023, st->fh);
+ st->ln++;
if ( rval == NULL ) continue;
chomp(line);
@@ -269,7 +273,7 @@ static int write_peaks_2_3(struct image *image, FILE *ofh)
}
-static RefList *read_stream_reflections_2_3(FILE *fh, struct detector *det)
+static RefList *read_stream_reflections_2_3(Stream *st, struct detector *det)
{
char *rval = NULL;
int first = 1;
@@ -290,7 +294,8 @@ static RefList *read_stream_reflections_2_3(FILE *fh, struct detector *det)
int r;
Reflection *refl;
- rval = fgets(line, 1023, fh);
+ rval = fgets(line, 1023, st->fh);
+ st->ln++;
if ( rval == NULL ) continue;
chomp(line);
@@ -342,7 +347,7 @@ static RefList *read_stream_reflections_2_3(FILE *fh, struct detector *det)
}
-static RefList *read_stream_reflections_2_1(FILE *fh, struct detector *det)
+static RefList *read_stream_reflections_2_1(Stream *st, struct detector *det)
{
char *rval = NULL;
int first = 1;
@@ -364,7 +369,8 @@ static RefList *read_stream_reflections_2_1(FILE *fh, struct detector *det)
int r;
Reflection *refl;
- rval = fgets(line, 1023, fh);
+ rval = fgets(line, 1023, st->fh);
+ st->ln++;
if ( rval == NULL ) continue;
chomp(line);
@@ -426,7 +432,7 @@ static RefList *read_stream_reflections_2_1(FILE *fh, struct detector *det)
}
-static RefList *read_stream_reflections_2_2(FILE *fh, struct detector *det)
+static RefList *read_stream_reflections_2_2(Stream *st, struct detector *det)
{
char *rval = NULL;
int first = 1;
@@ -442,7 +448,8 @@ static RefList *read_stream_reflections_2_2(FILE *fh, struct detector *det)
int r;
Reflection *refl;
- rval = fgets(line, 1023, fh);
+ rval = fgets(line, 1023, st->fh);
+ st->ln++;
if ( rval == NULL ) continue;
chomp(line);
@@ -887,14 +894,15 @@ int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
}
-static int find_start_of_chunk(FILE *fh)
+static int find_start_of_chunk(Stream *st)
{
char *rval = NULL;
char line[1024];
do {
- rval = fgets(line, 1023, fh);
+ rval = fgets(line, 1023, st->fh);
+ st->ln++;
/* Trouble? */
if ( rval == NULL ) return 1;
@@ -942,6 +950,7 @@ static void read_crystal(Stream *st, struct image *image, StreamReadFlags srf)
char c;
rval = fgets(line, 1023, st->fh);
+ st->ln++;
/* Trouble? */
if ( rval == NULL ) break;
@@ -975,7 +984,8 @@ static void read_crystal(Stream *st, struct image *image, StreamReadFlags srf)
centering = c;
have_cen = 1;
} else {
- ERROR("Duplicate centering ignored.\n");
+ ERROR("Duplicate centering (line %lli) - "
+ "stream may be corrupted!\n", st->ln);
}
}
@@ -986,7 +996,8 @@ static void read_crystal(Stream *st, struct image *image, StreamReadFlags srf)
unique_axis = c;
have_ua = 1;
} else {
- ERROR("Duplicate unique axis ignored.\n");
+ ERROR("Duplicate unique axis (line %lli) - "
+ "stream may be corrupted!\n", st->ln);
}
}
@@ -997,7 +1008,8 @@ static void read_crystal(Stream *st, struct image *image, StreamReadFlags srf)
lattice_type = lattice_from_str(line+15);
have_latt = 1;
} else {
- ERROR("Duplicate lattice type ignored.\n");
+ ERROR("Duplicate lattice type (line %lli) - "
+ "stream may be corrupted!\n", st->ln);
}
}
@@ -1030,13 +1042,13 @@ static void read_crystal(Stream *st, struct image *image, StreamReadFlags srf)
/* The reflection list format in the stream diverges
* after 2.2 */
if ( AT_LEAST_VERSION(st, 2, 3) ) {
- reflist = read_stream_reflections_2_3(st->fh,
+ reflist = read_stream_reflections_2_3(st,
image->det);
} else if ( AT_LEAST_VERSION(st, 2, 2) ) {
- reflist = read_stream_reflections_2_2(st->fh,
+ reflist = read_stream_reflections_2_2(st,
image->det);
} else {
- reflist = read_stream_reflections_2_1(st->fh,
+ reflist = read_stream_reflections_2_1(st,
image->det);
}
if ( reflist == NULL ) {
@@ -1148,7 +1160,7 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf)
int have_filename = 0;
int have_ev = 0;
- if ( find_start_of_chunk(st->fh) ) return 1;
+ if ( find_start_of_chunk(st) ) return 1;
image->lambda = -1.0;
image->features = NULL;
@@ -1167,6 +1179,7 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf)
float div, bw;
rval = fgets(line, 1023, st->fh);
+ st->ln++;
/* Trouble? */
if ( rval == NULL ) break;
@@ -1257,9 +1270,9 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf)
int fail;
if ( AT_LEAST_VERSION(st, 2, 3) ) {
- fail = read_peaks_2_3(st->fh, image);
+ fail = read_peaks_2_3(st, image);
} else {
- fail = read_peaks(st->fh, image);
+ fail = read_peaks(st, image);
}
if ( fail ) {
ERROR("Failed while reading peaks\n");
@@ -1358,6 +1371,8 @@ Stream *open_stream_for_read(const char *filename)
return NULL;
}
+ st->ln = 1;
+
return st;
}
@@ -1656,6 +1671,7 @@ void write_geometry_file(Stream *st, const char *geom_filename) {
*/
int rewind_stream(Stream *st)
{
+ st->ln = 0;
return fseek(st->fh, 0, SEEK_SET);
}
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 503078a7..d675fb30 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -3,12 +3,12 @@
*
* Symmetry
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
- * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ * 2010-2014,2016 Thomas White <taw@physics.org>
+ * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
*
* This file is part of CrystFEL.
*
@@ -1088,6 +1088,21 @@ SymOpList *get_pointgroup(const char *sym)
}
+void pointgroup_warning(const char *sym)
+{
+ if ( (strcmp(sym, "m") == 0)
+ || (strcmp(sym, "2/m") == 0)
+ || (strcmp(sym, "2") == 0) )
+ {
+ ERROR("WARNING: You have specified a monoclinic point group "
+ "without a unique axis. The default unique axis is 'c'. "
+ "If you want unique axis b, append '_uab' to your point "
+ "group.\n");
+ }
+
+}
+
+
static void do_op(const IntegerMatrix *op,
signed int h, signed int k, signed int l,
signed int *he, signed int *ke, signed int *le)
diff --git a/libcrystfel/src/symmetry.h b/libcrystfel/src/symmetry.h
index c29164ef..a61ca96f 100644
--- a/libcrystfel/src/symmetry.h
+++ b/libcrystfel/src/symmetry.h
@@ -3,12 +3,12 @@
*
* Symmetry
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
- * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ * 2010-2014,2016 Thomas White <taw@physics.org>
+ * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
*
* This file is part of CrystFEL.
*
@@ -90,6 +90,8 @@ extern void describe_symmetry(const SymOpList *s);
extern int is_centric(signed int h, signed int k, signed int l,
const SymOpList *ops);
+extern void pointgroup_warning(const char *sym);
+
extern void add_symop(SymOpList *ops, IntegerMatrix *m);
extern SymOpList *parse_symmetry_operations(const char *s);
extern char *get_matrix_name(const IntegerMatrix *m, int row);