aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2021-01-11 09:53:47 +0100
committerThomas White <taw@physics.org>2021-01-11 09:57:58 +0100
commit9d9e819bc535102ca52d6270b6a928e6e70539b6 (patch)
tree135655915dc6418356f72076b572fc1ed032159e
parent708590dd4bfac92bde90d9762f413c7496afbddf (diff)
UnitCell: Store all representations once they're calculated
Previously, the "getter" functions would re-calculate the requested representation every time they were called. This could mean doing a matrix inversion in the middle of a tight loop, wasting loads of time. Now, it stores the calculated values and returns them directly next time. Setting the parameters invalidates the values for all representations other than the one given. The cost of doing this is that the cell can no longer be "const" in the getter functions. This tracked through some other code, but nothing too severe.
-rw-r--r--libcrystfel/src/cell-utils.c4
-rw-r--r--libcrystfel/src/cell-utils.h2
-rw-r--r--libcrystfel/src/cell.c349
-rw-r--r--libcrystfel/src/cell.h10
-rw-r--r--libcrystfel/src/stream.c2
-rw-r--r--libcrystfel/src/stream.h2
-rw-r--r--src/post-refinement.c8
7 files changed, 178 insertions, 199 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,
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 09438a82..89b468c4 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -56,7 +56,7 @@ struct rf_alteration
struct rf_priv
{
- const Crystal *cr; /**< Original crystal (before any refinement) */
+ Crystal *cr; /**< Original crystal (before any refinement) */
const RefList *full;
int serial;
int scaleflags;
@@ -98,7 +98,7 @@ const char *str_prflag(enum prflag flag)
}
-static void rotate_cell_xy(const UnitCell *source, UnitCell *tgt,
+static void rotate_cell_xy(UnitCell *source, UnitCell *tgt,
double ang1, double ang2)
{
double asx, asy, asz;
@@ -150,7 +150,7 @@ static void rotate_cell_xy(const UnitCell *source, UnitCell *tgt,
}
-static void apply_parameters(const Crystal *cr_orig, Crystal *cr_tgt,
+static void apply_parameters(Crystal *cr_orig, Crystal *cr_tgt,
struct rf_alteration alter)
{
double R, lambda;
@@ -161,7 +161,7 @@ static void apply_parameters(const Crystal *cr_orig, Crystal *cr_tgt,
R = crystal_get_profile_radius(cr_orig);
lambda = crystal_get_image_const(cr_orig)->lambda;
- rotate_cell_xy(crystal_get_cell_const(cr_orig), crystal_get_cell(cr_tgt),
+ rotate_cell_xy(crystal_get_cell(cr_orig), crystal_get_cell(cr_tgt),
alter.rot_x, alter.rot_y);
crystal_set_profile_radius(cr_tgt, R+alter.delta_R);
image->lambda = lambda+alter.delta_wave;