Merge branch 'master' into taketwo

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);