aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/cell.c231
-rw-r--r--src/cell.h2
-rw-r--r--src/diffraction.c12
-rw-r--r--src/dirax.c14
-rw-r--r--src/image.h1
-rw-r--r--src/index.c9
-rw-r--r--src/indexamajig.c12
-rw-r--r--src/intensities.c6
-rw-r--r--src/intensities.h3
-rw-r--r--src/pattern_sim.c4
10 files changed, 259 insertions, 35 deletions
diff --git a/src/cell.c b/src/cell.c
index 86be7cc4..a3fdc72a 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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;
diff --git a/src/cell.h b/src/cell.h
index d1b49d22..687c7ef6 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -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 ) {
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 18:50:33 +0000