Commit 1c2e32b7 authored by Thomas White's avatar Thomas White
Browse files

Check unit cell is sensible

parent 0dc96534
......@@ -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;
}
......@@ -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 */
......@@ -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,
......
......@@ -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");
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment