6 files changed, 174 insertions, 195 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 93dea212..6b5d7717 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -333,8 +333,6 @@ void cell_print_full(UnitCell *cell)
 		       rad2deg(angle_between(asx, asy, asz, csx, csy, csz)),
 		       rad2deg(angle_between(asx, asy, asz, bsx, bsy, bsz)));
 
-		STATUS("Cell representation is %s.\n", cell_rep(cell));
-
 	}
 }
 
@@ -910,7 +908,7 @@ static int get_angle_rad(char **bits, int nbits, double *pl)
  * Writes \p cell to \p fh, in CrystFEL unit cell file format
  *
  */
-void write_cell(const UnitCell *cell, FILE *fh)
+void write_cell(UnitCell *cell, FILE *fh)
 {
 	double a, b, c, al, be, ga;
 	LatticeType lt;
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 0a110bf2..5789d606 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -58,7 +58,7 @@ extern void cell_print_full(UnitCell *cell);
 
 extern UnitCell *load_cell_from_pdb(const char *filename);
 extern UnitCell *load_cell_from_file(const char *filename);
-extern void write_cell(const UnitCell *cell, FILE *fh);
+extern void write_cell(UnitCell *cell, FILE *fh);
 
 extern int cell_is_sensible(UnitCell *cell);
 
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 23a5e46b..50a28470 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -35,6 +35,7 @@
 #endif
 
 #include <math.h>
+#include <assert.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
@@ -54,19 +55,14 @@
  */
 
 
-typedef enum {
-	CELL_REP_CRYST,
-	CELL_REP_CART,
-	CELL_REP_RECIP
-} CellRepresentation;
-
 struct _unitcell {
 
-	CellRepresentation rep;
-
-	int have_parameters;
+	LatticeType  lattice_type;
+	char         centering;
+	char         unique_axis;
 
 	/* Crystallographic representation */
+	int have_cryst;
 	double a;	/* m */
 	double b;	/* m */
 	double c;	/* m */
@@ -75,18 +71,16 @@ struct _unitcell {
 	double gamma;	/* Radians */
 
 	/* Cartesian representation */
+	int have_cart;
 	double ax;	double bx;	double cx;
 	double ay;	double by;	double cy;
 	double az;	double bz;	double cz;
 
 	/* Cartesian representation of reciprocal axes */
+	int have_recip;
 	double axs;	double bxs;	double cxs;
 	double ays;	double bys;	double cys;
 	double azs;	double bzs;	double czs;
-
-	LatticeType  lattice_type;
-	char         centering;
-	char         unique_axis;
 };
 
 typedef enum {
@@ -127,12 +121,13 @@ UnitCell *cell_new()
 	cell->beta = 0.0;
 	cell->gamma = 0.0;
 
-	cell->rep = CELL_REP_CRYST;
+	cell->have_cryst = 0;
+	cell->have_cart = 0;
+	cell->have_recip = 0;
 
 	cell->lattice_type = L_TRICLINIC;
 	cell->centering = 'P';
 	cell->unique_axis = '?';
-	cell->have_parameters = 0;
 
 	return cell;
 }
@@ -160,7 +155,9 @@ void cell_free(UnitCell *cell)
 int cell_has_parameters(const UnitCell *cell)
 {
 	if ( cell == NULL ) return 0;
-	return cell->have_parameters;
+	return (cell->have_cryst > 0)
+	    || (cell->have_cart > 0 )
+	    || (cell->have_recip > 0);
 }
 
 
@@ -176,8 +173,9 @@ void cell_set_parameters(UnitCell *cell, double a, double b, double c,
 	cell->beta = beta;
 	cell->gamma = gamma;
 
-	cell->rep = CELL_REP_CRYST;
-	cell->have_parameters = 1;
+	cell->have_cryst = 1;
+	cell->have_cart = 0;
+	cell->have_recip = 0;
 }
 
 
@@ -192,8 +190,9 @@ void cell_set_cartesian(UnitCell *cell,
 	cell->bx = bx;  cell->by = by;  cell->bz = bz;
 	cell->cx = cx;  cell->cy = cy;  cell->cz = cz;
 
-	cell->rep = CELL_REP_CART;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 1;
+	cell->have_recip = 0;
 }
 
 
@@ -223,8 +222,9 @@ UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
 	cell->bxs = bs.u;  cell->bys = bs.v;  cell->bzs = bs.w;
 	cell->cxs = cs.u;  cell->cys = cs.v;  cell->czs = cs.w;
 
-	cell->rep = CELL_REP_RECIP;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 0;
+	cell->have_recip = 1;
 
 	return cell;
 }
@@ -241,8 +241,9 @@ UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c)
 	cell->bx = b.u;  cell->by = b.v;  cell->bz = b.w;
 	cell->cx = c.u;  cell->cy = c.v;  cell->cz = c.w;
 
-	cell->rep = CELL_REP_CART;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 1;
+	cell->have_recip = 0;
 
 	return cell;
 }
@@ -268,8 +269,9 @@ void cell_set_reciprocal(UnitCell *cell,
 	cell->bxs = bsx;  cell->bys = bsy;  cell->bzs = bsz;
 	cell->cxs = csx;  cell->cys = csy;  cell->czs = csz;
 
-	cell->rep = CELL_REP_RECIP;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 0;
+	cell->have_recip = 1;
 }
 
 
@@ -291,31 +293,25 @@ void cell_set_unique_axis(UnitCell *cell, char unique_axis)
 }
 
 
-/************************* Getter helper functions ****************************/
+/************************* Conversion functions ****************************/
 
-static int cell_crystallographic_to_cartesian(const UnitCell *cell,
-                                              double *ax, double *ay, double *az,
-                                              double *bx, double *by, double *bz,
-                                              double *cx, double *cy, double *cz)
+static void crystallographic_to_cartesian(UnitCell *cell)
 {
 	double tmp, V, cosalphastar, cstar;
 
-	if ( !cell->have_parameters ) {
-		ERROR("Unit cell has unspecified parameters.\n");
-		return 1;
-	}
+	assert(cell->have_cryst == 1);
 
 	/* Firstly: Get a in terms of x, y and z
 	 * +a (cryst) is defined to lie along +x (cart) */
-	*ax = cell->a;
-	*ay = 0.0;
-	*az = 0.0;
+	cell->ax = cell->a;
+	cell->ay = 0.0;
+	cell->az = 0.0;
 
 	/* b in terms of x, y and z
 	 * b (cryst) is defined to lie in the xy (cart) plane */
-	*bx = cell->b*cos(cell->gamma);
-	*by = cell->b*sin(cell->gamma);
-	*bz = 0.0;
+	cell->bx = cell->b*cos(cell->gamma);
+	cell->by = cell->b*sin(cell->gamma);
+	cell->bz = 0.0;
 
 	tmp = cos(cell->alpha)*cos(cell->alpha)
 		+ cos(cell->beta)*cos(cell->beta)
@@ -329,21 +325,41 @@ static int cell_crystallographic_to_cartesian(const UnitCell *cell,
 	cstar = (cell->a * cell->b * sin(cell->gamma))/V;
 
 	/* c in terms of x, y and z */
-	*cx = cell->c*cos(cell->beta);
-	*cy = -cell->c*sin(cell->beta)*cosalphastar;
-	*cz = 1.0/cstar;
+	cell->cx = cell->c*cos(cell->beta);
+	cell->cy = -cell->c*sin(cell->beta)*cosalphastar;
+	cell->cz = 1.0/cstar;
 
-	return 0;
+	cell->have_cart = 1;
+}
+
+
+static void cartesian_to_crystallographic(UnitCell *cell)
+{
+	assert(cell->have_cart);
+
+	/* Convert cartesian -> crystallographic */
+	cell->a = modulus(cell->ax, cell->ay, cell->az);
+	cell->b = modulus(cell->bx, cell->by, cell->bz);
+	cell->c = modulus(cell->cx, cell->cy, cell->cz);
+
+	cell->alpha = angle_between(cell->bx, cell->by, cell->bz,
+	                            cell->cx, cell->cy, cell->cz);
+	cell->beta = angle_between(cell->ax, cell->ay, cell->az,
+	                           cell->cx, cell->cy, cell->cz);
+	cell->gamma = angle_between(cell->ax, cell->ay, cell->az,
+	                            cell->bx, cell->by, cell->bz);
+
+	cell->have_cryst = 1;
 }
 
 
 /* Why yes, I do enjoy long argument lists...! */
-static int cell_invert(double ax, double ay, double az,
-                       double bx, double by, double bz,
-                       double cx, double cy, double cz,
-                       double *asx, double *asy, double *asz,
-                       double *bsx, double *bsy, double *bsz,
-                       double *csx, double *csy, double *csz)
+static int invert(double ax, double ay, double az,
+                  double bx, double by, double bz,
+                  double cx, double cy, double cz,
+                  double *asx, double *asy, double *asz,
+                  double *bsx, double *bsy, double *bsz,
+                  double *csx, double *csy, double *csz)
 {
 	int s;
 	gsl_matrix *m;
@@ -413,157 +429,143 @@ static int cell_invert(double ax, double ay, double az,
 }
 
 
+static int reciprocal_to_cartesian(UnitCell *cell)
+{
+	assert(cell->have_recip);
+
+	if ( invert(cell->axs, cell->ays, cell->azs,
+	            cell->bxs, cell->bys, cell->bzs,
+	            cell->cxs, cell->cys, cell->czs,
+	            &cell->ax, &cell->ay, &cell->az,
+	            &cell->bx, &cell->by, &cell->bz,
+	            &cell->cx, &cell->cy, &cell->cz) ) return 1;
+
+	cell->have_cart = 1;
+	return 0;
+}
+
+
+static int cartesian_to_reciprocal(UnitCell *cell)
+{
+	assert(cell->have_cart);
+
+	if ( invert(cell->ax, cell->ay, cell->az,
+	            cell->bx, cell->by, cell->bz,
+	            cell->cx, cell->cy, cell->cz,
+	            &cell->axs, &cell->ays, &cell->azs,
+	            &cell->bxs, &cell->bys, &cell->bzs,
+	            &cell->cxs, &cell->cys, &cell->czs) ) return 1;
+
+	cell->have_recip = 1;
+	return 0;
+}
+
+
 /********************************** Getters ***********************************/
 
-int cell_get_parameters(const UnitCell *cell, double *a, double *b, double *c,
+int cell_get_parameters(UnitCell *cell,
+                        double *a, double *b, double *c,
                         double *alpha, double *beta, double *gamma)
 {
-	double ax, ay, az, bx, by, bz, cx, cy, cz;
-
 	if ( cell == NULL ) return 1;
 
-	if ( !cell->have_parameters ) {
+	if ( cell->have_cryst ) {
+
+		/* Nothing to do */
+
+	} else if ( cell->have_cart ) {
+
+		cartesian_to_crystallographic(cell);
+
+	} else if ( cell->have_recip ) {
+
+		if ( reciprocal_to_cartesian(cell) ) return 1;
+		cartesian_to_crystallographic(cell);
+
+	} else {
+
 		ERROR("Unit cell has unspecified parameters.\n");
 		return 1;
-	}
 
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		/* Direct response */
-		*a = cell->a;
-		*b = cell->b;
-		*c = cell->c;
-		*alpha = cell->alpha;
-		*beta = cell->beta;
-		*gamma = cell->gamma;
-		return 0;
-
-		case CELL_REP_CART:
-		/* Convert cartesian -> crystallographic */
-		*a = modulus(cell->ax, cell->ay, cell->az);
-		*b = modulus(cell->bx, cell->by, cell->bz);
-		*c = modulus(cell->cx, cell->cy, cell->cz);
-
-		*alpha = angle_between(cell->bx, cell->by, cell->bz,
-		                       cell->cx, cell->cy, cell->cz);
-		*beta = angle_between(cell->ax, cell->ay, cell->az,
-		                      cell->cx, cell->cy, cell->cz);
-		*gamma = angle_between(cell->ax, cell->ay, cell->az,
-		                       cell->bx, cell->by, cell->bz);
-		return 0;
-
-		case CELL_REP_RECIP:
-		/* Convert reciprocal -> crystallographic.
-                 * Start by converting reciprocal -> cartesian */
-		if ( cell_invert(cell->axs, cell->ays, cell->azs,
-		                 cell->bxs, cell->bys, cell->bzs,
-		                         cell->cxs, cell->cys, cell->czs,
-		                         &ax, &ay, &az,
-		                         &bx, &by, &bz,
-		                         &cx, &cy, &cz) ) return 1;
-
-		/* Now convert cartesian -> crystallographic */
-		*a = modulus(ax, ay, az);
-		*b = modulus(bx, by, bz);
-		*c = modulus(cx, cy, cz);
-
-		*alpha = angle_between(bx, by, bz, cx, cy, cz);
-		*beta = angle_between(ax, ay, az, cx, cy, cz);
-		*gamma = angle_between(ax, ay, az, bx, by, bz);
-		return 0;
 	}
 
-	return 1;
+	*a = cell->a;
+	*b = cell->b;
+	*c = cell->c;
+	*alpha = cell->alpha;
+	*beta = cell->beta;
+	*gamma = cell->gamma;
+	return 0;
 }
 
 
-int cell_get_cartesian(const UnitCell *cell,
+int cell_get_cartesian(UnitCell *cell,
                        double *ax, double *ay, double *az,
                        double *bx, double *by, double *bz,
                        double *cx, double *cy, double *cz)
 {
 	if ( cell == NULL ) return 1;
 
-	if ( !cell->have_parameters ) {
+	if ( cell->have_cart ) {
+
+		/* Nothing to do */
+
+	} else if ( cell->have_recip ) {
+
+		/* NB recip->cart has priority over
+		 * cryst->cart, to preserve orientation */
+		reciprocal_to_cartesian(cell);
+
+	} else if ( cell->have_cryst ) {
+
+		crystallographic_to_cartesian(cell);
+
+	} else {
+
 		ERROR("Unit cell has unspecified parameters.\n");
 		return 1;
-	}
-
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		/* Convert crystallographic -> cartesian. */
-		return cell_crystallographic_to_cartesian(cell,
-		                                          ax, ay, az,
-		                                          bx, by, bz,
-		                                          cx, cy, cz);
-
-		case CELL_REP_CART:
-		/* Direct response */
-		*ax = cell->ax;  *ay = cell->ay;  *az = cell->az;
-		*bx = cell->bx;  *by = cell->by;  *bz = cell->bz;
-		*cx = cell->cx;  *cy = cell->cy;  *cz = cell->cz;
-		return 0;
-
-		case CELL_REP_RECIP:
-		/* Convert reciprocal -> cartesian */
-		return cell_invert(cell->axs, cell->ays, cell->azs,
-		                   cell->bxs, cell->bys, cell->bzs,
-		                   cell->cxs, cell->cys, cell->czs,
-		                   ax, ay, az, bx, by, bz, cx, cy, cz);
 
 	}
 
-	return 1;
+	*ax = cell->ax;  *ay = cell->ay;  *az = cell->az;
+	*bx = cell->bx;  *by = cell->by;  *bz = cell->bz;
+	*cx = cell->cx;  *cy = cell->cy;  *cz = cell->cz;
+	return 0;
 }
 
 
-int cell_get_reciprocal(const UnitCell *cell,
+int cell_get_reciprocal(UnitCell *cell,
                         double *asx, double *asy, double *asz,
                         double *bsx, double *bsy, double *bsz,
                         double *csx, double *csy, double *csz)
 {
-	int r;
-	double ax, ay, az, bx, by, bz, cx, cy, cz;
-
 	if ( cell == NULL ) return 1;
 
-	if ( !cell->have_parameters ) {
+	if ( cell->have_recip ) {
+
+		/* Nothing to do */
+
+	} else if ( cell->have_cart ) {
+
+		/* NB cart->recip has priority over cryst->recip */
+		cartesian_to_reciprocal(cell);
+
+	} else if ( cell->have_cryst ) {
+
+		crystallographic_to_cartesian(cell);
+		cartesian_to_reciprocal(cell);
+
+	} else {
+
 		ERROR("Unit cell has unspecified parameters.\n");
 		return 1;
-	}
-
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		/* Convert crystallographic -> reciprocal */
-		r = cell_crystallographic_to_cartesian(cell,
-		                                       &ax, &ay, &az,
-		                                       &bx, &by, &bz,
-		                                       &cx, &cy, &cz);
-		if ( r ) return r;
-		return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz,
-		                   asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
-		case CELL_REP_CART:
-		/* Convert cartesian -> reciprocal */
-		cell_invert(cell->ax, cell->ay, cell->az,
-		            cell->bx, cell->by, cell->bz,
-		            cell->cx, cell->cy, cell->cz,
-		            asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-		return 0;
-
-		case CELL_REP_RECIP:
-		/* Direct response */
-		*asx = cell->axs;  *asy = cell->ays;  *asz = cell->azs;
-		*bsx = cell->bxs;  *bsy = cell->bys;  *bsz = cell->bzs;
-		*csx = cell->cxs;  *csy = cell->cys;  *csz = cell->czs;
-		return 0;
 
 	}
 
-	return 1;
+	*asx = cell->axs;  *asy = cell->ays;  *asz = cell->azs;
+	*bsx = cell->bxs;  *bsy = cell->bys;  *bsz = cell->bzs;
+	*csx = cell->cxs;  *csy = cell->cys;  *csz = cell->czs;
+	return 0;
 }
 
 
@@ -579,7 +581,7 @@ LatticeType cell_get_lattice_type(const UnitCell *cell)
 }
 
 
-struct g6 cell_get_G6(const UnitCell *cell)
+struct g6 cell_get_G6(UnitCell *cell)
 {
 	double a, b, c, al, be, ga;
 	struct g6 g;
@@ -600,25 +602,6 @@ char cell_get_unique_axis(const UnitCell *cell)
 }
 
 
-const char *cell_rep(UnitCell *cell)
-{
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		return "crystallographic, direct space";
-
-		case CELL_REP_CART:
-		return "cartesian, direct space";
-
-		case CELL_REP_RECIP:
-		return "cartesian, reciprocal space";
-
-	}
-
-	return "unknown";
-}
-
-
 UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m)
 {
 	gsl_matrix *c;
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index bf5d689c..064dacd0 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -107,15 +107,15 @@ extern void cell_set_cartesian(UnitCell *cell,
 extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
 				double alpha, double beta, double gamma);
 
-extern int cell_get_parameters(const UnitCell *cell, double *a, double *b, double *c,
+extern int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
                                double *alpha, double *beta, double *gamma);
 
-extern int cell_get_cartesian(const UnitCell *cell,
+extern int cell_get_cartesian(UnitCell *cell,
                               double *ax, double *ay, double *az,
                               double *bx, double *by, double *bz,
                               double *cx, double *cy, double *cz);
 
-extern int cell_get_reciprocal(const UnitCell *cell,
+extern int cell_get_reciprocal(UnitCell *cell,
                                double *asx, double *asy, double *asz,
                                double *bsx, double *bsy, double *bsz,
                                double *csx, double *csy, double *csz);
@@ -138,7 +138,7 @@ struct g6
 	double F;
 };
 
-extern struct g6 cell_get_G6(const UnitCell *cell);
+extern struct g6 cell_get_G6(UnitCell *cell);
 
 extern char cell_get_centering(const UnitCell *cell);
 extern void cell_set_centering(UnitCell *cell, char centering);
@@ -146,8 +146,6 @@ extern void cell_set_centering(UnitCell *cell, char centering);
 extern char cell_get_unique_axis(const UnitCell *cell);
 extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
 
-extern const char *cell_rep(UnitCell *cell);
-
 extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m);
 
 extern UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 08f0e332..afaf40a5 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -1250,7 +1250,7 @@ Stream *stream_open_fd_for_write(int fd, const DataTemplate *dtempl)
 }
 
 
-void stream_write_target_cell(Stream *st, const UnitCell *cell)
+void stream_write_target_cell(Stream *st, UnitCell *cell)
 {
 	if ( cell == NULL ) return;
 	fprintf(st->fh, STREAM_CELL_START_MARKER"\n");
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index 3d57c7b6..f503f391 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -98,7 +98,7 @@ extern void stream_close(Stream *st);
 extern void stream_write_geometry_file(Stream *st,
                                        const char *geom_filename);
 extern void stream_write_target_cell(Stream *st,
-                                     const UnitCell *cell);
+                                     UnitCell *cell);
 extern void stream_write_commandline_args(Stream *st,
                                           int argc, char *argv[]);
 extern void stream_write_indexing_methods(Stream *st,