Check unit cell is sensible

4 files changed, 33 insertions, 13 deletions

diff --git a/src/cell.c b/src/cell.c
index caeab3c7..c31697bc 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1142,3 +1142,23 @@ UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot)
 
 	return out;
 }
+
+
+int cell_is_sensible(UnitCell *cell)
+{
+	double a, b, c, al, be, ga;
+
+	cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+	if (   al + be + ga >= 2.0*M_PI ) return 0;
+	if (   al + be - ga >= 2.0*M_PI ) return 0;
+	if (   al - be + ga >= 2.0*M_PI ) return 0;
+	if ( - al + be + ga >= 2.0*M_PI ) return 0;
+	if (   al + be + ga <= 0.0 ) return 0;
+	if (   al + be - ga <= 0.0 ) return 0;
+	if (   al - be + ga <= 0.0 ) return 0;
+	if ( - al + be + ga <= 0.0 ) return 0;
+	if ( isnan(al) ) return 0;
+	if ( isnan(be) ) return 0;
+	if ( isnan(ga) ) return 0;
+	return 1;
+}
diff --git a/src/cell.h b/src/cell.h
index 23b28da5..b5d31fc6 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -102,5 +102,6 @@ extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template);
 
 extern UnitCell *load_cell_from_pdb(const char *filename);
 
+extern int cell_is_sensible(UnitCell *cell);
 
 #endif	/* CELL_H */
diff --git a/src/geometry.c b/src/geometry.c
index 56bb2b0a..8d5ec6f8 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -236,23 +236,16 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
 	int hmax, kmax, lmax;
 	double mres;
 	signed int h, k, l;
-	double a, b, c, al, be, ga;
 
 	reflections = reflist_new();
 
 	/* Cell angle check from Foadi and Evans (2011) */
-	cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
-	if (   al + be + ga >= 2.0*M_PI ) return NULL;
-	if (   al + be - ga >= 2.0*M_PI ) return NULL;
-	if (   al - be + ga >= 2.0*M_PI ) return NULL;
-	if ( - al + be + ga >= 2.0*M_PI ) return NULL;
-	if (   al + be + ga <= 0.0 ) return NULL;
-	if (   al + be - ga <= 0.0 ) return NULL;
-	if (   al - be + ga <= 0.0 ) return NULL;
-	if ( - al + be + ga <= 0.0 ) return NULL;
-	if ( isnan(al) ) return NULL;
-	if ( isnan(be) ) return NULL;
-	if ( isnan(ga) ) return NULL;
+	if ( !cell_is_sensible(cell) ) {
+		ERROR("Invalid unit cell parameters given to"
+		      " find_intersections()\n");
+		cell_print(cell);
+		return NULL;
+	}
 
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
 	                          &bsx, &bsy, &bsz,
diff --git a/src/partial_sim.c b/src/partial_sim.c
index ea7be5e0..8fc2fb21 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -242,6 +242,12 @@ int main(int argc, char *argv[])
 	}
 	free(cellfile);
 
+	if ( !cell_is_sensible(cell) ) {
+		ERROR("Invalid unit cell parameters:\n");
+		cell_print(cell);
+		return 1;
+	}
+
 	/* Load geometry */
 	if ( geomfile == NULL ) {
 		ERROR("You need to give a geometry file.\n");