aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-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;