Generate vectors for ReAx search

1 files changed, 36 insertions, 0 deletions

diff --git a/src/reax.c b/src/reax.c
index 44dc52a6..60784da7 100644
--- a/src/reax.c
+++ b/src/reax.c
@@ -38,18 +38,54 @@ struct dvec
 struct reax_private
 {
 	IndexingPrivate base;
+	struct dvec *directions;
 };
 
 
 IndexingPrivate *reax_prepare()
 {
 	struct reax_private *priv;
+	int ui, vi;
+	int samp;
+	double angular_inc;
 
 	priv = calloc(1, sizeof(*priv));
 	if ( priv == NULL ) return NULL;
 
 	priv->base.indm = INDEXING_REAX;
 
+	/* Decide on sampling interval */
+	angular_inc = 0.03;  /* From Steller (1997) */
+	samp = (2.0 * M_PI) / angular_inc;
+
+	priv->directions = malloc(samp*samp*sizeof(struct dvec));
+	if ( priv == NULL) {
+		free(priv);
+		return NULL;
+	}
+
+	for ( ui=0; ui<samp; ui++ ) {
+	for ( vi=0; vi<samp; vi++ ) {
+
+		double u, v;
+		double th, ph;
+		struct dvec *dir;
+
+		u = (double)ui/samp;
+		v = (double)vi/samp;
+
+		th = 2.0 * M_PI * u;
+		ph = acos(2.0*v - 1.0);
+
+		dir = &priv->directions[ui + vi*samp];
+
+		dir->x = cos(th) * sin(ph);
+		dir->y = sin(th) * sin(th);
+		dir->z = cos(ph);
+
+	}
+	}
+
 	return (IndexingPrivate *)priv;
 }