aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-07-29 16:07:16 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:35 +0100
commitbca7a90bc0bb8e5aecbbdd946e950bcaa24f1c77 (patch)
tree23e914f4e512c3b8fd6aeb66192c409eab4a6855 /src
parentd7ee0d1aff5be156ba5f71d03d7591f364a2aa04 (diff)
More ReAx stuff
Diffstat (limited to 'src')
-rw-r--r--src/reax.c122
1 files changed, 114 insertions, 8 deletions
diff --git a/src/reax.c b/src/reax.c
index 77ca124f..20e8c94c 100644
--- a/src/reax.c
+++ b/src/reax.c
@@ -18,6 +18,8 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
+#include <assert.h>
+#include <fftw3.h>
#include "image.h"
#include "utils.h"
@@ -39,9 +41,116 @@ struct reax_private
{
IndexingPrivate base;
struct dvec *directions;
+ int n_dir;
};
+static double check_dir(struct dvec *dir, ImageFeatureList *flist, double modv,
+ double pmax, double *fft_in, fftw_complex *fft_out,
+ fftw_plan plan)
+{
+ int n, i, v;
+ double *vals;
+
+ n = image_feature_count(flist);
+
+ vals = malloc(n*sizeof(double));
+ v = 0;
+ for ( i=0; i<n; i++ ) {
+
+ struct imagefeature *f;
+
+ f = image_get_feature(flist, i);
+ if ( f == NULL ) continue;
+
+ vals[v] = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z;
+
+ v++; /* Yuk, yuk, yuk */
+
+ }
+
+
+
+
+ free(vals);
+
+ return 0.0;
+}
+
+
+void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
+{
+ int i;
+ struct reax_private *p;
+ double fom;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double mod_as;
+ double dx, dy, dz;
+ int nel, n;
+ double pmax;
+ double *fft_in;
+ fftw_complex *fft_out;
+ fftw_plan plan;
+ const double cellmax = 50.0e9; /* 50 nm max cell size */
+
+ assert(pp->indm == INDEXING_REAX);
+ p = (struct reax_private *)pp;
+
+ cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ mod_as = modulus(asx, asy, asz);
+
+ n = image_feature_count(image->features);
+ pmax = 0.0;
+ for ( i=0; i<n; i++ ) {
+
+ struct imagefeature *f;
+ double val;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ val = modulus(f->rx, f->ry, f->rz);
+
+ if ( val > pmax ) pmax = val;
+
+ }
+ nel = 2.0*pmax*5.0*cellmax;
+
+ fft_in = fftw_malloc(nel*sizeof(double));
+ fft_out = fftw_malloc((nel/2 + 1)*sizeof(fftw_complex));
+
+ plan = fftw_plan_dft_r2c_1d(nel, fft_in, fft_out, FFTW_ESTIMATE);
+
+ /* Search for a* */
+ fom = 0.0; dx = 0.0; dy = 0.0; dz = 0.0;
+ for ( i=0; i<p->n_dir; i++ ) {
+
+ double new_fom;
+
+ new_fom = check_dir(&p->directions[i], image->features, mod_as,
+ pmax, fft_in, fft_out, plan);
+ if ( new_fom > fom ) {
+ fom = new_fom;
+ dx = p->directions[i].x;
+ dy = p->directions[i].x;
+ dz = p->directions[i].x;
+ }
+
+ }
+
+ fftw_destroy_plan(plan);
+ fftw_free(fft_in);
+ fftw_free(fft_out);
+
+ /* No improvement from zero? */
+ if ( fom == 0.0 ) return;
+}
+
+
IndexingPrivate *reax_prepare()
{
struct reax_private *priv;
@@ -58,12 +167,14 @@ IndexingPrivate *reax_prepare()
angular_inc = 0.03; /* From Steller (1997) */
samp = (2.0 * M_PI) / angular_inc;
- priv->directions = malloc(samp*samp*sizeof(struct dvec));
+ priv->n_dir = samp*samp;
+ priv->directions = malloc(priv->n_dir*sizeof(struct dvec));
if ( priv == NULL) {
free(priv);
return NULL;
}
+ /* Generate vectors for 1D Fourier transforms */
for ( ui=0; ui<samp; ui++ ) {
for ( vi=0; vi<samp; vi++ ) {
@@ -74,8 +185,9 @@ IndexingPrivate *reax_prepare()
u = (double)ui/samp;
v = (double)vi/samp;
+ /* Uniform sampling of a hemisphere */
th = 2.0 * M_PI * u;
- ph = acos(2.0*v - 1.0);
+ ph = acos(v);
dir = &priv->directions[ui + vi*samp];
@@ -88,9 +200,3 @@ IndexingPrivate *reax_prepare()
return (IndexingPrivate *)priv;
}
-
-
-void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell)
-{
-
-}