1 files changed, 287 insertions, 0 deletions
diff --git a/src/reflist.c b/src/reflist.c
new file mode 100644
index 0000000..e0224e9
--- /dev/null
+++ b/src/reflist.c
@@ -0,0 +1,287 @@
+/*
+ * reflist.c
+ *
+ * Reflection-handling code
+ *
+ * (c) 2006-2008 Thomas White <taw27@cam.ac.uk>
+ *
+ *  synth2d - two-dimensional Fourier synthesis
+ *
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+
+#include "reflist.h"
+
+ReflectionList *reflist_new() {
+
+	ReflectionList *new = malloc(sizeof(ReflectionList));
+	assert(new != NULL);
+
+	new->n_reflections = 0;
+
+	/* DC component makes reflist_inlist() work */
+	reflist_addref(new, 0, 0, 0);
+	new->n_reflections = 1;
+
+	return new;
+
+}
+
+ReflectionList *reflist_new_parent(ReflectionList *parent) {
+
+	ReflectionList *new = reflist_new();
+
+	new->parent = parent;
+
+	return new;
+
+}
+
+ReflectionList *reflist_copy(ReflectionList *list) {
+	ReflectionList *new = malloc(sizeof(ReflectionList));
+	memcpy(new, list, sizeof(ReflectionList));
+	return new;
+}
+
+void reflist_free(ReflectionList *list) {
+	free(list);
+}
+
+/* Add a reflection to a list, including all associated data */
+static int reflist_addref_all(ReflectionList *list, signed int h, signed int k, signed int l,
+				double am, double phase, double phase_set, double alpha, double delta_theta, signed int multiplier,
+				unsigned int parent_index) {
+
+	unsigned int i;
+
+	if ( list->n_reflections >= MAX_REFLECTIONS ) {
+		fprintf(stderr, "Too many reflections (%i): increase MAX_REFLECTIONS in reflist.h (addref)\n",
+					list->n_reflections);
+		return -1;
+	}
+
+	//if ( (i = reflist_inlist(list, h, k, l)) ) {
+		//printf("RL: Duplicated reflection %i %i %i\n", h, k, l);
+	//} else {
+		i = list->n_reflections;
+	//}
+	if ( (h==0) && (k==0) && (l==0) ) i = 0;
+
+	list->refs[i].h = h;
+	list->refs[i].k = k;
+	list->refs[i].l = l;
+	list->refs[i].amplitude = am;
+	list->refs[i].am_error = 1;
+	list->refs[i].weight = 1;
+	list->refs[i].phase_known = phase;
+	list->refs[i].phase_known_set = phase_set;
+	list->refs[i].alpha = alpha;
+	list->refs[i].delta_theta = delta_theta;
+	list->refs[i].multiplier = multiplier;
+	list->refs[i].parent_index = parent_index;
+
+	if ( !((h==0) && (k==0) && (l==0)) ) list->n_reflections++;
+
+	return i;
+
+}
+
+/* Add a reflection to a list */
+int reflist_addref(ReflectionList *list, signed int h, signed int k, signed int l) {
+	return reflist_addref_all(list, h, k, l, 0, 0, 0, 0, 0, 1, 0);
+}
+
+/* Add a reflection to a list */
+int reflist_addref_am(ReflectionList *list, signed int h, signed int k, signed int l, double am) {
+	return reflist_addref_all(list, h, k, l, am, 0, 0, 0, 0, 1, 0);
+}
+
+int reflist_addref_am_ph(ReflectionList *list, signed int h, signed int k, signed int l, double am, double ph) {
+	return reflist_addref_all(list, h, k, l, am, ph, 1, 0, 0, 1, 0);
+}
+
+/* Add a reflection to a list, including alpha value */
+int reflist_addref_alpha_parent(ReflectionList *list, signed int h, signed int k, signed int l, double alpha, unsigned int parent_index) {
+	return reflist_addref_all(list, h, k, l, 0, 0, 0, alpha, 0, 1, parent_index);
+}
+
+/* Add a reflection to a list, including phase value */
+int reflist_addref_phase(ReflectionList *list, signed int h, signed int k, signed int l, double phase) {
+	return reflist_addref_all(list, h, k, l, 0, phase, 1, 0, 0, 1, 0);
+}
+
+/* Add a reflection to a list, including "delta-phase" value */
+int reflist_addref_deltatheta(ReflectionList *list, signed int h, signed int k, signed int l, double delta_theta, signed int multiplier) {
+	if ( !reflist_inlist(list, h, k, l) ) {
+		return reflist_addref_all(list, h, k, l, 0, 0, 0, 0, delta_theta, multiplier, 0);
+	} else {
+		//printf("RL: Not re-adding reflection %3i %3i %3i\n", h, k, l);
+		return -1;
+	}
+}
+
+int reflist_addref_parent(ReflectionList *list, signed int h, signed int k, signed int l, unsigned int parent_index) {
+	return reflist_addref_all(list, h, k, l, 0, 0, 0, 0, 0, 1, parent_index);
+}
+
+int reflist_addref_am_parent(ReflectionList *list, signed int h, signed int k, signed int l, double am, unsigned int parent_index) {
+	return reflist_addref_all(list, h, k, l, am, 0, 0, 0, 0, 1, parent_index);
+}
+
+/* Delete a reflection from a list */
+void reflist_delref(ReflectionList *list, signed int h, signed int k, signed int l) {
+
+	unsigned int i;
+
+	if ( list->n_reflections >= MAX_REFLECTIONS ) {
+		fprintf(stderr, "Too many reflections (%i): increase MAX_REFLECTIONS in reflist.h (delref)\n",
+					list->n_reflections);
+		return;
+	}
+
+	if ( list->n_reflections == 0 ) { return; }
+
+	for ( i=1; i<list->n_reflections; i++ ) {
+		if ( (list->refs[i].h == h) && (list->refs[i].k == k) ) {
+
+			if ( i < list->n_reflections ) {
+				unsigned int j;
+				/* Shove everything up one place */
+				for ( j=i+1; j<list->n_reflections; j++ ) {
+					list->refs[j-1].h = list->refs[j].h;
+					list->refs[j-1].k = list->refs[j].k;
+					list->refs[j-1].l = list->refs[j].l;
+					list->refs[j-1].amplitude = list->refs[j].amplitude;
+					list->refs[j-1].phase_known = list->refs[j].phase_known;
+					list->refs[j-1].phase_known_set = list->refs[j].phase_known_set;
+					list->refs[j-1].phase_calc = list->refs[j].phase_calc;
+					list->refs[j-1].phase_calc_set = list->refs[j].phase_calc_set;
+					list->refs[j-1].alpha = list->refs[j].alpha;
+					list->refs[j-1].delta_theta = list->refs[j].delta_theta;
+					list->refs[j-1].multiplier = list->refs[j].multiplier;
+					list->refs[j-1].parent_index = list->refs[j].parent_index;
+				}
+			} /* else nothing to do */
+
+			list->n_reflections--;
+
+		}
+	}
+
+}
+
+/* See if a reflection is in a list */
+int reflist_inlist(const ReflectionList *list, signed int h, signed int k, signed int l) {
+
+	unsigned int i;
+
+	if ( list->n_reflections >= MAX_REFLECTIONS ) {
+		fprintf(stderr, "Too many reflections (%i): increase MAX_REFLECTIONS in reflist.h (refinlist)\n",
+					list->n_reflections);
+		return 0;
+	}
+
+	if ( list->n_reflections == 0 ) { return 0; }
+
+	for ( i=0; i<list->n_reflections; i++ ) {
+		if ( (list->refs[i].h == h) && (list->refs[i].k == k) && (list->refs[i].l == l) ) return i;
+	}
+
+	return 0;
+
+}
+
+/* As above, but ignoring "l" */
+int reflist_inlist_2d(const ReflectionList *list, signed int h, signed int k) {
+
+	unsigned int i;
+
+	if ( list->n_reflections >= MAX_REFLECTIONS ) {
+		fprintf(stderr, "Too many reflections (%i): increase MAX_REFLECTIONS in reflist.h (refinlist)\n",
+					list->n_reflections);
+		return 0;
+	}
+
+	if ( list->n_reflections == 0 ) { return 0; }
+
+	for ( i=0; i<list->n_reflections; i++ ) {
+		if ( (list->refs[i].h == h) && (list->refs[i].k == k) ) return i;
+	}
+
+	return 0;
+
+}
+
+void reflist_set_components(ReflectionList *list, signed int h, signed int k, signed int l, double re, double im) {
+
+	unsigned int idx = reflist_inlist(list, h, k, l);
+	if ( !idx ) return;
+	list->refs[idx].re = re;
+	list->refs[idx].im = im;
+
+}
+
+ReflectionList *reflist_new_from_array(fftw_complex *array, int width, int height) {
+
+	ReflectionList *reflections;
+	signed int h, k;
+
+	reflections = reflist_new();
+
+	for ( h=-(width-1)/2; h<=(width-1)/2; h++ ) {
+		for ( k=-(height-1)/2; k<=(height-1)/2; k++ ) {
+
+			double re, im, am, ph;
+			signed int hc, kc;
+
+			hc = h; kc = k;
+			if ( h < 0 ) { hc = width+hc; }
+			if ( k < 0 ) { kc = height+kc; }
+			re = array[kc+height*hc][0];	im = array[kc+height*hc][1];
+			am = sqrt(re*re + im*im);	ph = atan2(im, re);
+
+			if ( am > 0 ) {
+				reflist_addref_am_ph(reflections, h, k, 0, am, ph);
+				//printf("Got %3i %3i am=%f ph=%f\n", h, k, am, ph);
+			}
+
+		}
+	}
+
+	return reflections;
+
+}
+
+void reflist_fill_array(fftw_complex *array, const ReflectionList *reflections, int width, int height) {
+
+	signed int h, k;
+	unsigned int i;
+
+	bzero(array, width*height*sizeof(fftw_complex));
+
+	for ( i=0; i<reflections->n_reflections; i++ ) {
+
+		h = reflections->refs[i].h;
+		k = reflections->refs[i].k;
+
+		if ( abs(h) > width/2 ) {
+			printf("Index %i %i (%i) is above the Nyquist frequency!\n", h, k, reflections->refs[i].l);
+			continue;
+		}
+
+		if ( h < 0 ) { h = width+h; }
+		if ( k < 0 ) { k = height+k; }
+
+		if ( reflections->refs[i].phase_known_set ) {
+			array[k + height*h][0] = reflections->refs[i].amplitude * cos(reflections->refs[i].phase_known);
+			array[k + height*h][1] = reflections->refs[i].amplitude * sin(reflections->refs[i].phase_known);
+		}
+
+	}
+
+}