path: root/src/reflist.c

/*
 * 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);
		}

	}

}