Commit 9c3d9caa authored by Thomas White's avatar Thomas White
Browse files

Match the unit cell to a model cell after indexing

parent d19a20b8
......@@ -22,6 +22,7 @@
#include "cell.h"
#include "utils.h"
#include "image.h"
/* Update the cartesian representation from the crystallographic one */
......@@ -190,6 +191,55 @@ UnitCell *cell_new_from_parameters(double a, double b, double c,
}
static UnitCell *cell_new_from_axes(struct rvec as, struct rvec bs,
struct rvec cs)
{
UnitCell *cell;
int s;
gsl_matrix *m;
gsl_matrix *inv;
gsl_permutation *perm;
cell = cell_new();
if ( !cell ) return NULL;
m = gsl_matrix_alloc(3, 3);
gsl_matrix_set(m, 0, 0, as.u);
gsl_matrix_set(m, 0, 1, as.v);
gsl_matrix_set(m, 0, 2, as.w);
gsl_matrix_set(m, 1, 0, bs.u);
gsl_matrix_set(m, 1, 1, bs.v);
gsl_matrix_set(m, 1, 2, bs.w);
gsl_matrix_set(m, 2, 0, cs.u);
gsl_matrix_set(m, 2, 1, cs.v);
gsl_matrix_set(m, 2, 2, cs.w);
/* Invert */
perm = gsl_permutation_alloc(m->size1);
inv = gsl_matrix_alloc(m->size1, m->size2);
gsl_linalg_LU_decomp(m, perm, &s);
gsl_linalg_LU_invert(m, perm, inv);
gsl_permutation_free(perm);
gsl_matrix_free(m);
/* Transpose */
gsl_matrix_transpose(inv);
cell->ax = gsl_matrix_get(inv, 0, 0);
cell->ay = gsl_matrix_get(inv, 0, 1);
cell->az = gsl_matrix_get(inv, 0, 2);
cell->bx = gsl_matrix_get(inv, 1, 0);
cell->by = gsl_matrix_get(inv, 1, 1);
cell->bz = gsl_matrix_get(inv, 1, 2);
cell->cx = gsl_matrix_get(inv, 2, 0);
cell->cy = gsl_matrix_get(inv, 2, 1);
cell->cz = gsl_matrix_get(inv, 2, 2);
cell_update_crystallographic(cell);
return cell;
}
void cell_get_cartesian(UnitCell *cell,
double *ax, double *ay, double *az,
double *bx, double *by, double *bz,
......@@ -247,6 +297,187 @@ void cell_get_reciprocal(UnitCell *cell,
}
static void cell_print(UnitCell *cell)
{
STATUS(" a b c alpha beta gamma\n");
STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n",
cell->a*1e9, cell->b*1e9, cell->c*1e9,
rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
}
#define MAX_CAND (1024)
static int within_tolerance(double a, double b, double percent)
{
double tol = a * (percent/100.0);
if ( fabs(b-a) < tol ) return 1;
return 0;
}
struct cvec {
struct rvec vec;
float na;
float nb;
float nc;
};
/* Attempt to make 'cell' fit into 'template' somehow */
UnitCell *match_cell(UnitCell *cell, UnitCell *template)
{
signed int n1l, n2l, n3l;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
int i, j;
double lengths[3];
double angles[3];
struct cvec *cand[3];
int ncand[3] = {0,0,0};
float ltl = 5.0;
float angtol = 5.0;
UnitCell *new_cell = NULL;
cell_get_reciprocal(template, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
lengths[0] = modulus(asx, asy, asz);
lengths[1] = modulus(bsx, bsy, bsz);
lengths[2] = modulus(csx, csy, csz);
angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
STATUS("Performing unit cell magic...\n");
/* Negative values mean 1/n, positive means n, zero means zero */
for ( n1l=-4; n1l<=4; n1l++ ) {
for ( n2l=-4; n2l<=4; n2l++ ) {
for ( n3l=-4; n3l<=4; n3l++ ) {
float n1, n2, n3;
signed int b1, b2, b3;
n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
/* 'bit' values can be +1 or -1 */
for ( b1=-1; b1<=1; b1+=2 ) {
for ( b2=-1; b2<=1; b2+=2 ) {
for ( b3=-1; b3<=1; b3+=2 ) {
double tx, ty, tz;
double tlen;
int i;
n1 *= b1; n2 *= b2; n3 *= b3;
tx = n1*asx + n2*asy + n3*asz;
ty = n1*bsx + n2*bsy + n3*bsz;
tz = n1*csx + n2*csy + n3*csz;
tlen = modulus(tx, ty, tz);
/* Test modulus for agreement with moduli of template */
for ( i=0; i<3; i++ ) {
if ( within_tolerance(lengths[i], tlen, ltl) ) {
STATUS("sought %e, found %e (%e %e) for %i\n",
lengths[i], tlen, 1.0/lengths[i],
1.0/tlen, i);
cand[i][ncand[i]].vec.u = tx;
cand[i][ncand[i]].vec.v = ty;
cand[i][ncand[i]].vec.w = tz;
cand[i][ncand[i]].na = n1;
cand[i][ncand[i]].nb = n2;
cand[i][ncand[i]].nc = n3;
ncand[i]++;
}
}
}
}
}
}
}
}
for ( i=0; i<ncand[0]; i++ ) {
for ( j=0; j<ncand[1]; j++ ) {
double ang;
int k;
/* Measure the angle between the ith candidate for axis 0
* and the jth candidate for axis 1 */
ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
cand[0][i].vec.w, cand[1][j].vec.u,
cand[1][j].vec.v, cand[1][j].vec.w);
/* Angle between axes 0 and 1 should be angle 2 */
if ( !within_tolerance(ang, angles[2], angtol) ) continue;
for ( k=0; k<ncand[2]; k++ ) {
/* Measure the angle between the current candidate for
* axis 0 and the kth candidate for axis 2 */
ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
cand[0][i].vec.w, cand[2][k].vec.u,
cand[2][k].vec.v, cand[2][k].vec.w);
/* ... it should be angle 1 ... */
if ( !within_tolerance(ang, angles[1], angtol) ) continue;
/* Finally, the angle between the current candidate for
* axis 1 and the kth candidate for axis 2 */
ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
cand[1][j].vec.w, cand[2][k].vec.u,
cand[2][k].vec.v, cand[2][k].vec.w);
/* ... it should be angle 0 ... */
if ( !within_tolerance(ang, angles[0], angtol) ) continue;
new_cell = cell_new_from_axes(cand[0][i].vec,
cand[1][j].vec,
cand[2][k].vec);
STATUS("Success! --------------- \n");
cell_print(new_cell);
STATUS("%f %f %f\n", cand[0][i].na, cand[0][i].nb,
cand[0][i].nc);
STATUS("%f %f %f\n", cand[1][j].na, cand[1][j].nb,
cand[1][j].nc);
STATUS("%f %f %f\n", cand[2][k].na, cand[2][k].nb,
cand[2][k].nc);
goto done;
}
}
}
done:
free(cand[0]);
free(cand[1]);
free(cand[2]);
return new_cell;
}
double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
{
const double a = cell->a;
......
......@@ -72,4 +72,6 @@ extern void cell_get_reciprocal(UnitCell *cell,
extern double resolution(UnitCell *cell,
signed int h, signed int k, signed int l);
extern UnitCell *match_cell(UnitCell *cell, UnitCell *template);
#endif /* CELL_H */
......@@ -140,24 +140,12 @@ void get_diffraction(struct image *image, int na, int nb, int nc)
double bx, by, bz;
double cx, cy, cz;
double a, b, c, d;
struct molecule *mtmp;
/* Generate the array of reciprocal space vectors in image->qvecs */
if ( image->qvecs == NULL ) {
get_ewald(image);
}
/* FIXME: Nasty */
mtmp = load_molecule();
if ( image->molecule == NULL ) {
image->molecule = mtmp;
} else {
int i;
for ( i=0; i<32; i++ ) {
image->molecule->species[i] = mtmp->species[i];
}
image->molecule->n_species = mtmp->n_species;
}
if ( image->molecule == NULL ) return;
cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
......
......@@ -69,13 +69,7 @@ static void dirax_parseline(const char *line, struct image *image)
if ( line[i] == 'R' ) rf = 1;
if ( (line[i] == 'D') && rf ) {
image->dirax_read_cell = 1;
if ( image->molecule == NULL ) {
image->molecule = malloc(sizeof(struct molecule));
} else if ( image->molecule->cell ) {
free(image->molecule->cell);
free(image->molecule);
}
image->molecule->cell = cell_new();
image->indexed_cell = cell_new();
return;
}
i++;
......@@ -86,7 +80,7 @@ static void dirax_parseline(const char *line, struct image *image)
/* First row of unit cell values */
float ax, ay, az, d;
sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &ax, &ay, &az);
cell_set_cartesian_a(image->molecule->cell,
cell_set_cartesian_a(image->indexed_cell,
ax*1e-10, ay*1e-10, az*1e-10);
image->dirax_read_cell++;
return;
......@@ -94,7 +88,7 @@ static void dirax_parseline(const char *line, struct image *image)
/* First row of unit cell values */
float bx, by, bz, d;
sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &bx, &by, &bz);
cell_set_cartesian_b(image->molecule->cell,
cell_set_cartesian_b(image->indexed_cell,
bx*1e-10, by*1e-10, bz*1e-10);
image->dirax_read_cell++;
return;
......@@ -102,7 +96,7 @@ static void dirax_parseline(const char *line, struct image *image)
/* First row of unit cell values */
float cx, cy, cz, d;
sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &cx, &cy, &cz);
cell_set_cartesian_c(image->molecule->cell,
cell_set_cartesian_c(image->indexed_cell,
cx*1e-10, cy*1e-10, cz*1e-10);
STATUS("Read a direct space unit cell from DirAx\n");
/* FIXME: Do something */
......
......@@ -78,6 +78,7 @@ struct image {
struct rvec *qvecs;
double *twotheta;
struct molecule *molecule;
UnitCell *indexed_cell;
struct quaternion orientation;
......
......@@ -96,7 +96,16 @@ void index_pattern(struct image *image, int no_index, int use_dirax)
}
if ( use_dirax ) {
UnitCell *new_cell;
run_dirax(image, no_index);
new_cell = match_cell(image->indexed_cell,
image->molecule->cell);
if ( new_cell != NULL ) image->indexed_cell = new_cell;
return;
}
}
......@@ -30,6 +30,7 @@
#include "peaks.h"
#include "diffraction.h"
#include "detector.h"
#include "sfac.h"
static void show_help(const char *s)
......@@ -132,8 +133,9 @@ int main(int argc, char *argv[])
chomp(line);
image.features = NULL;
image.molecule = NULL;
image.molecule = load_molecule();
image.data = NULL;
image.indexed_cell = NULL;
STATUS("Processing '%s'\n", line);
......@@ -162,10 +164,8 @@ int main(int argc, char *argv[])
if ( config_noindex ) goto done;
/* Calculate orientation matrix (by magic) */
index_pattern(&image, config_noindex,
config_dirax);
if ( image.molecule == NULL ) goto done;
index_pattern(&image, config_noindex, config_dirax);
if ( image.indexed_cell == NULL ) goto done;
if ( config_nearbragg || config_simulate ) {
......@@ -186,7 +186,7 @@ int main(int argc, char *argv[])
if ( config_nearbragg ) {
/* Read h,k,l,I */
output_intensities(&image);
output_intensities(&image, image.indexed_cell);
}
if ( config_simulate ) {
......
......@@ -49,7 +49,7 @@ static int sum_nearby_points(int16_t *data, int width, int x, int y)
}
void output_intensities(struct image *image)
void output_intensities(struct image *image, UnitCell *cell)
{
int x, y;
double ax, ay, az;
......@@ -59,9 +59,7 @@ void output_intensities(struct image *image)
int n_hits = 0;
int i;
cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
&bx, &by, &bz,
&cx, &cy, &cz);
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
......
......@@ -17,7 +17,8 @@
#define INTENSITIES_H
#include "image.h"
#include "cell.h"
extern void output_intensities(struct image *image);
extern void output_intensities(struct image *image, UnitCell *cell);
#endif /* INTENSITIES_H */
......@@ -218,7 +218,7 @@ int main(int argc, char *argv[])
image.camera_len = 0.05; /* 5 cm (front CCD can move from 5cm-20cm) */
image.resolution = 13333.3; /* 75 micron pixel size */
image.lambda = ph_en_to_lambda(eV_to_J(2.0e3)); /* Wavelength */
image.molecule = NULL;
image.molecule = load_molecule();
/* Splurge a few useful numbers */
STATUS("Wavelength is %f nm\n", image.lambda/1.0e-9);
......@@ -263,7 +263,7 @@ int main(int argc, char *argv[])
!config_nobloom);
if ( config_nearbragg ) {
output_intensities(&image);
output_intensities(&image, image.molecule->cell);
}
if ( !config_noimages ) {
......
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