From 189da15810deabd739d7c11c6e95fea55739fe60 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 1 Aug 2020 15:13:49 +0200 Subject: Initial import from archive --- src/reflist.c | 287 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 287 insertions(+) create mode 100644 src/reflist.c (limited to 'src/reflist.c') 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 + * + * synth2d - two-dimensional Fourier synthesis + * + */ + +#include +#include +#include +#include +#include + +#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; in_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; jn_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; in_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; in_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; in_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); + } + + } + +} -- cgit v1.2.3