diff options
-rw-r--r-- | libcrystfel/src/cell-utils.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/cell.c | 349 | ||||
-rw-r--r-- | libcrystfel/src/cell.h | 10 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 2 | ||||
-rw-r--r-- | libcrystfel/src/stream.h | 2 | ||||
-rw-r--r-- | src/post-refinement.c | 8 |
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; |