diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/cell.c | 231 | ||||
| -rw-r--r-- | src/cell.h | 2 | ||||
| -rw-r--r-- | src/diffraction.c | 12 | ||||
| -rw-r--r-- | src/dirax.c | 14 | ||||
| -rw-r--r-- | src/image.h | 1 | ||||
| -rw-r--r-- | src/index.c | 9 | ||||
| -rw-r--r-- | src/indexamajig.c | 12 | ||||
| -rw-r--r-- | src/intensities.c | 6 | ||||
| -rw-r--r-- | src/intensities.h | 3 | ||||
| -rw-r--r-- | src/pattern_sim.c | 4 |
10 files changed, 259 insertions, 35 deletions
@@ -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 */ diff --git a/src/diffraction.c b/src/diffraction.c index eb67e4a8..84206908 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -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, diff --git a/src/dirax.c b/src/dirax.c index 2c0675da..4161c79c 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -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 */ diff --git a/src/image.h b/src/image.h index 0f64d248..59435b4e 100644 --- a/src/image.h +++ b/src/image.h @@ -78,6 +78,7 @@ struct image { struct rvec *qvecs; double *twotheta; struct molecule *molecule; + UnitCell *indexed_cell; struct quaternion orientation; diff --git a/src/index.c b/src/index.c index b3c0b7e7..a20ad3a5 100644 --- a/src/index.c +++ b/src/index.c @@ -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; } } diff --git a/src/indexamajig.c b/src/indexamajig.c index 1a94882b..245d1a81 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -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 ) { diff --git a/src/intensities.c b/src/intensities.c index 8b100d9f..d5b3604c 100644 --- a/src/intensities.c +++ b/src/intensities.c @@ -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++ ) { diff --git a/src/intensities.h b/src/intensities.h index 8f863287..10664ab2 100644 --- a/src/intensities.h +++ b/src/intensities.h @@ -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 */ diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 9e4abf90..7a27ceb6 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -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 ) { |