/*
 * symmetry.c
 *
 * Symmetry stuff
 *
 * (c) 2006-2008 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2d - Two-Dimensional Crystallographic Fourier Synthesis
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>

#include "symmetry.h"
#include "reflist.h"
#include "model.h"

#define is_odd(a) ((a)%2==1)

/* Return a list of atoms containing the given atom and all its symmetry equivalents */
AtomicModel *symmetry_generate_equivalent_atoms(AtomicModel *model, size_t j, ModelTarget target) {

	AtomicModel *list;

	list = model_new();
	list->atoms[list->n_atoms] = model->atoms[j];  list->n_atoms++;

	if ( target == MODEL_TARGET_DISPLAY ) {

		if ( model->atoms[j].x == 0 ) {
			list->atoms[list->n_atoms] = model->atoms[j];
			list->atoms[list->n_atoms].x = 1.0;
			list->n_atoms++;
		}

		if ( model->atoms[j].y == 0 ) {
			list->atoms[list->n_atoms] = model->atoms[j];
			list->atoms[list->n_atoms].y = 1.0;
			list->n_atoms++;
		}

		if ( (model->atoms[j].x == 0) && (model->atoms[j].y == 0) ) {
			list->atoms[list->n_atoms] = model->atoms[j];
			list->atoms[list->n_atoms].x = 1.0;
			list->atoms[list->n_atoms].y = 1.0;
			list->n_atoms++;
		}

	}

	if ( model->sym & SYMMETRY_CENTRE ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = 1-list->atoms[list->n_atoms].x;
		list->atoms[list->n_atoms].y = 1-list->atoms[list->n_atoms].y;
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_MIRROR_HORIZONTAL ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].y = 1-list->atoms[list->n_atoms].y;
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_MIRROR_VERTICAL ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = 1-list->atoms[list->n_atoms].x;
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_GLIDE_HORIZONTAL ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = fmod(0.5 + list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = 1-list->atoms[list->n_atoms].y;
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_GLIDE_VERTICAL ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = 1-list->atoms[list->n_atoms].x;
		list->atoms[list->n_atoms].y = fmod(0.5 + list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = fmod(0.5 + list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = fmod(1.5 - list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_GLIDE_VERTICAL_QUARTER ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = fmod(1.5 - list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = fmod(0.5 + list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	}

	/* Fix cm */
	if ( (model->sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER) && (model->sym & SYMMETRY_MIRROR_HORIZONTAL) ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = fmod(0.5 + list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = 1 - fmod(1.5 - list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	} else if ( (model->sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (model->sym & SYMMETRY_MIRROR_VERTICAL) ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = 1 - fmod(1.5 - list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = fmod(0.5 + list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	}

	/* Fix c2mm */
	if ( (model->sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (model->sym & SYMMETRY_MIRROR_HORIZONTAL) ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = fmod(1.5 - list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = 1 - fmod(0.5 + list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_MIRROR_HORIZONTAL_QUARTER ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = list->atoms[list->n_atoms].x;
		list->atoms[list->n_atoms].y = fmod(1.5 - list->atoms[list->n_atoms].y, 1);
		list->n_atoms++;
	}

	if ( model->sym & SYMMETRY_MIRROR_VERTICAL_QUARTER ) {
		list->atoms[list->n_atoms] = model->atoms[j];
		list->atoms[list->n_atoms].x = fmod(1.5 - list->atoms[list->n_atoms].x, 1);
		list->atoms[list->n_atoms].y = list->atoms[list->n_atoms].y;
		list->n_atoms++;
	}

	return list;

}

ReflectionList *symmetry_generate_equivalent_reflections(Symmetry sym, signed int h, signed int k, signed int l) {

	ReflectionList *result;

	result = reflist_new();

	/* Remove the minus in front of 'l' to make even HOLZ conditional transforms real.
	 *	This is generally the wrong thing to do. */
	if ( sym & SYMMETRY_FRIEDEL ) reflist_addref_deltatheta(result, -h, -k, -l, 0, -1);

	if ( sym & SYMMETRY_MIRROR_HORIZONTAL ) reflist_addref_deltatheta(result, h, -k, l, 0, 1);
	if ( sym & SYMMETRY_MIRROR_VERTICAL ) reflist_addref_deltatheta(result, -h, k, l, 0, 1);
	if ( sym & SYMMETRY_MIRROR_DIAGONAL ) reflist_addref_deltatheta(result, -h, -k, l, 0, 1);

	if ( sym & SYMMETRY_GLIDE_HORIZONTAL ) reflist_addref_deltatheta(result, h, -k, l, M_PI, 1);
	if ( sym & SYMMETRY_GLIDE_VERTICAL ) reflist_addref_deltatheta(result, -h, k, l, M_PI, 1);

	if ( sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER ) {
		if ( is_odd(abs(h)+abs(k)) ) {
			reflist_addref_deltatheta(result, h, -k, l, M_PI, 1);
		} else {
			reflist_addref_deltatheta(result, h, -k, l, 0, 1);
		}
	}

	if ( sym & SYMMETRY_GLIDE_VERTICAL_QUARTER ) {
		if ( is_odd(abs(h)+abs(k)) ) {
			reflist_addref_deltatheta(result, -h, k, l, M_PI, 1);
		} else {
			reflist_addref_deltatheta(result, -h, k, l, 0, 1);
		}
	}

	if ( sym & SYMMETRY_ROTATION_4 ) {
		if ( sym & SYMMETRY_GLIDE_DIAGONAL ) {
			if ( is_odd(abs(h)+abs(k)) ) {
				reflist_addref_deltatheta(result, h, -k, l, M_PI, 1);
				reflist_addref_deltatheta(result, -h, k, l, M_PI, 1);
				reflist_addref_deltatheta(result, -h, -k, l, 0, 1);
				reflist_addref_deltatheta(result, k, h, l, 0, 1);
				reflist_addref_deltatheta(result, -k, h, l, M_PI, 1);
				reflist_addref_deltatheta(result, k, -h, l, M_PI, 1);
				reflist_addref_deltatheta(result, -k, -h, l, 0, 1);
			} else {
				reflist_addref_deltatheta(result, h, -k, l, 0, 1);
				reflist_addref_deltatheta(result, -h, k, l, 0, 1);
				reflist_addref_deltatheta(result, -h, -k, l, 0, 1);
				reflist_addref_deltatheta(result, k, h, l, 0, 1);
				reflist_addref_deltatheta(result, -k, h, l, 0, 1);
				reflist_addref_deltatheta(result, k, -h, l, 0, 1);
				reflist_addref_deltatheta(result, -k, -h, l, 0, 1);
			}
		} else {
			reflist_addref_deltatheta(result, -h, -k, l, 0, 1);
			reflist_addref_deltatheta(result, -h, k, l, 0, 1);
			reflist_addref_deltatheta(result, h, -k, l, 0, 1);
		}
	}

	return result;

}

static double symmetry_realphase(double phk) {

	if ( (phk >= 0) && (phk < M_PI_2) ) phk = 0;
	else if ( (phk > M_PI_2) && (phk < M_PI) ) phk = M_PI;
	else if ( (phk < 0) && (phk > -M_PI/2) ) phk = 0;
	else if ( (phk <= -M_PI_2) && (phk > -M_PI) ) phk = M_PI;

	return phk;

}

void symmetry_centricity(ReflectionList *reflections, unsigned int i, Symmetry sym, SymFlags flags) {

	double am, ph;
	signed int h, k;

	ph = 0;
	h = reflections->refs[i].h;
	k = reflections->refs[i].k;

	am = reflections->refs[i].amplitude;
	if ( flags & SYMFLAG_PHASES_KNOWN ) ph = fmod(reflections->refs[i].phase_known, 2*M_PI);
	if ( flags & SYMFLAG_PHASES_CALC ) ph = fmod(reflections->refs[i].phase_calc, 2*M_PI);

	if ( (sym & SYMMETRY_MIRROR_HORIZONTAL) && (sym & SYMMETRY_FRIEDEL) && (h == 0) ) {
		ph = symmetry_realphase(ph);
	}

	if ( (sym & SYMMETRY_MIRROR_VERTICAL) && (sym & SYMMETRY_FRIEDEL) && (k == 0) ) {
		ph = symmetry_realphase(ph);
	}

	if ( (sym & SYMMETRY_CENTRE) && (sym & SYMMETRY_FRIEDEL) ) {
		ph = symmetry_realphase(ph);
	}

	if ((	( (sym & SYMMETRY_GLIDE_HORIZONTAL) && (sym & SYMMETRY_MIRROR_VERTICAL) )
	     || ( (sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER) && (sym & SYMMETRY_MIRROR_VERTICAL) )
	     || ( (sym & SYMMETRY_GLIDE_VERTICAL) && (sym & SYMMETRY_MIRROR_HORIZONTAL) )
	     || ( (sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (sym & SYMMETRY_MIRROR_HORIZONTAL) ) )
	     && (is_odd(abs(h)+abs(k))) ) {
		if ( (ph >= 0) && (ph < M_PI) ) ph = M_PI_2;
		else if ( (ph >= M_PI_2) && (ph < 2*M_PI) ) ph = -M_PI_2;
		else if ( (ph < 0) && (ph > -M_PI) ) ph = -M_PI_2;
		else if ( (ph <= -M_PI) ) ph = M_PI_2;
	}

	if ( flags & SYMFLAG_PHASES_KNOWN ) reflections->refs[i].phase_known = ph;
	if ( flags & SYMFLAG_PHASES_CALC ) reflections->refs[i].phase_calc = ph;

}

unsigned int symmetry_reflection_allowed(signed int h, signed int k, signed int l, Symmetry sym) {

	if ( (sym & SYMMETRY_GLIDE_HORIZONTAL) && (h == 0) && (k % 2) ) return 0;
	if ( (sym & SYMMETRY_GLIDE_HORIZONTAL_QUARTER) && (h == 0) && (k % 2) ) return 0;
	if ( (sym & SYMMETRY_GLIDE_VERTICAL) && (k == 0) && (h % 2) ) return 0;
	if ( (sym & SYMMETRY_GLIDE_VERTICAL_QUARTER) && (k == 0) && (h % 2) ) return 0;
	if ( (sym & SYMMETRY_ROTATION_4 ) && ( sym & SYMMETRY_GLIDE_DIAGONAL ) && (k == 0) && (h % 2) ) return 0;
	if ( (sym & SYMMETRY_ROTATION_4 ) && ( sym & SYMMETRY_GLIDE_DIAGONAL ) && (h == 0) && (k % 2) ) return 0;

	return 1;

}

/* Symmetrise a list of reflections (expanding all equivalents to P1) */
double symmetry_symmetrise(ReflectionList *reflections, Symmetry sym, SymFlags flags) {

	unsigned int i = 0;
	double r_sym = 0;
	double r_sym_norm = 0;
	unsigned int new_ref = 0;
	unsigned int elim_ref = 0;
	double am_elim = 0;
	unsigned int n = reflections->n_reflections;

	for ( i=1; i<n; i++ ) {	/* Don't symmetrise 000 */

		signed int h, k, l;

		h = reflections->refs[i].h;
		k = reflections->refs[i].k;
		l = reflections->refs[i].l;

		if ( (h==0) && (k==0) && (l==0) ) continue;
		if ( symmetry_reflection_allowed(h, k, l, sym) ) {

			double ph = 69;
			unsigned int j;
			double av;
			unsigned int av_n;
			ReflectionList *equivalents;

			equivalents = symmetry_generate_equivalent_reflections(sym, h, k, l);
			symmetry_centricity(reflections, i, sym, flags);

			/* First pass: calculate average amplitude, determine phase */
			av = reflections->refs[i].amplitude; av_n = 1;
			for ( j=1; j<equivalents->n_reflections; j++ ) {
				unsigned int p;
				p = reflist_inlist(reflections, equivalents->refs[j].h, equivalents->refs[j].k, equivalents->refs[j].l);
				if ( p ) {
					av += reflections->refs[p].amplitude;
					av_n++;
				}
			}
			av = av / av_n;
			r_sym += fabs(av - reflections->refs[i].amplitude);
			r_sym_norm += reflections->refs[i].amplitude;
			if ( flags & SYMFLAG_PHASES_KNOWN ) ph = fmod(reflections->refs[i].phase_known, 2*M_PI);
			if ( flags & SYMFLAG_PHASES_CALC ) ph = fmod(reflections->refs[i].phase_calc, 2*M_PI);
			//printf("Phase of %3i %3i %3i is %f\n", h, k, l, ph);

			/* Second pass: set all equivalent reflections to the average amplitude */
			reflections->refs[i].amplitude = av;
			for ( j=1; j<equivalents->n_reflections; j++ ) {

				unsigned int p;
				double phn;
				p = reflist_inlist(reflections, equivalents->refs[j].h, equivalents->refs[j].k, equivalents->refs[j].l);
				if ( p ) {
					r_sym += fabs(av - reflections->refs[p].amplitude);
					r_sym_norm += reflections->refs[p].amplitude;
					reflections->refs[p].amplitude = av;
				} else {
					p = reflist_addref_am(reflections, equivalents->refs[j].h, equivalents->refs[j].k, equivalents->refs[j].l, av);
					//printf("SY: Generating %3i %3i %3i am=%f\n", equivalents->refs[j].h,
					//			equivalents->refs[j].k, equivalents->refs[j].l, av);
					new_ref++;
				}
				phn = (ph + equivalents->refs[j].delta_theta) * equivalents->refs[j].multiplier;
				if ( flags & SYMFLAG_PHASES_KNOWN ) {
					reflections->refs[p].phase_known = phn;
					reflections->refs[p].phase_known_set = 1;
				}
				if ( flags & SYMFLAG_PHASES_CALC ) {
					reflections->refs[p].phase_calc = phn;
					reflections->refs[p].phase_calc_set = 1;
				}
				//printf("Set phase of %3i %3i %3i to %f (dt=%f, mul=%i)\n", equivalents->refs[j].h,
							//equivalents->refs[j].k, equivalents->refs[j].l,
							//phn, equivalents->refs[j].delta_theta, equivalents->refs[j].multiplier);


			}

			reflist_free(equivalents);

		} else {

			am_elim += reflections->refs[i].amplitude;
			reflist_delref(reflections, h, k, l);
			//printf("SY: Eliminating systematically absent reflection %i %i %i\n", h, k, l);
			elim_ref++;
			i--;

		}

	}

	if ( r_sym > 0 ) {
		r_sym = r_sym / r_sym_norm;
	}
	//printf("SY: R_sym = %.2f%%, %i reflections generated, %i reflections eliminated (total amplitude %f)\n", r_sym*100, new_ref, elim_ref, am_elim);

	return r_sym;

}

Symmetry symmetry_encode(const char *symmetry) {

	if ( strcmp(symmetry, "p1") == 0 ) return PLANEGROUP_P1;
	if ( strcmp(symmetry, "p2") == 0 ) return PLANEGROUP_P2;
	if ( strcmp(symmetry, "pm (m // x)") == 0 ) return PLANEGROUP_PM_X;
	if ( strcmp(symmetry, "pm (m // y)") == 0 ) return PLANEGROUP_PM_Y;
	if ( strcmp(symmetry, "pg (g // x)") == 0 ) return PLANEGROUP_PG_X;
	if ( strcmp(symmetry, "pg (g // y)") == 0 ) return PLANEGROUP_PG_Y;
	if ( strcmp(symmetry, "cm (m // x)") == 0 ) return PLANEGROUP_CM_X;
	if ( strcmp(symmetry, "cm (m // y)") == 0 ) return PLANEGROUP_CM_Y;
	if ( strcmp(symmetry, "p2mm") == 0 ) return PLANEGROUP_P2MM;
	if ( strcmp(symmetry, "p2mg (m // x)") == 0 ) return PLANEGROUP_P2MG_X;
	if ( strcmp(symmetry, "p2mg (m // y)") == 0 ) return PLANEGROUP_P2MG_Y;
	if ( strcmp(symmetry, "p2gg") == 0 ) return PLANEGROUP_P2GG;
	if ( strcmp(symmetry, "c2mm") == 0 ) return PLANEGROUP_C2MM;
	if ( strcmp(symmetry, "p4") == 0 ) return PLANEGROUP_P4;
	if ( strcmp(symmetry, "p4mm") == 0 ) return PLANEGROUP_P4MM;
	if ( strcmp(symmetry, "p4gm") == 0 ) return PLANEGROUP_P4GM;
	if ( strcmp(symmetry, "p3") == 0 ) return PLANEGROUP_P3;
	if ( strcmp(symmetry, "p3m1") == 0 ) return PLANEGROUP_P3M1;
	if ( strcmp(symmetry, "p31m") == 0 ) return PLANEGROUP_P31M;
	if ( strcmp(symmetry, "p6") == 0 ) return PLANEGROUP_P6;
	if ( strcmp(symmetry, "p6mm") == 0 ) return PLANEGROUP_P6MM;

	fprintf(stderr, "Unrecognised symmetry identifier '%s'\n", symmetry);
	return SYMMETRY_IDENTITY;

}

const char *symmetry_decode(Symmetry sym) {

	if ( sym == PLANEGROUP_P1 ) return "p1";
	if ( sym == PLANEGROUP_P2 ) return "p2";
	if ( sym == PLANEGROUP_PM_X ) return "pm (m // x)";
	if ( sym == PLANEGROUP_PM_Y ) return "pm (m // y)";
	if ( sym == PLANEGROUP_PG_X ) return "pg (g // x)";
	if ( sym == PLANEGROUP_PG_Y ) return "pg (g // y)";
	if ( sym == PLANEGROUP_CM_X ) return "cm (m // x)";
	if ( sym == PLANEGROUP_CM_Y ) return "cm (m // y)";
	if ( sym == PLANEGROUP_P2MM ) return "p2mm";
	if ( sym == PLANEGROUP_P2MG_X ) return "p2mg (m // x)";
	if ( sym == PLANEGROUP_P2MG_Y ) return "p2mg (m // y)";
	if ( sym == PLANEGROUP_P2GG ) return "p2gg";
	if ( sym == PLANEGROUP_C2MM ) return "c2mm";
	if ( sym == PLANEGROUP_P4 ) return "p4";
	if ( sym == PLANEGROUP_P4MM ) return "p4mm";
	if ( sym == PLANEGROUP_P4GM ) return "p4gm";
	if ( sym == PLANEGROUP_P3 ) return "p3";
	if ( sym == PLANEGROUP_P3M1 ) return "p3m1";
	if ( sym == PLANEGROUP_P31M ) return "p31m";
	if ( sym == PLANEGROUP_P6 ) return "p6";
	if ( sym == PLANEGROUP_P6MM ) return "p6mm";
	if ( sym == (PLANEGROUP_P1 | SYMMETRY_FRIEDEL) ) return "p1 & Friedel";
	if ( sym == (PLANEGROUP_P2 | SYMMETRY_FRIEDEL) ) return "p2 & Friedel";
	if ( sym == (PLANEGROUP_PM_X | SYMMETRY_FRIEDEL) ) return "pm (m // x) & Friedel";
	if ( sym == (PLANEGROUP_PM_Y | SYMMETRY_FRIEDEL) ) return "pm (m // y) & Friedel";
	if ( sym == (PLANEGROUP_PG_X | SYMMETRY_FRIEDEL) ) return "pg (g // x) & Friedel";
	if ( sym == (PLANEGROUP_PG_Y | SYMMETRY_FRIEDEL) ) return "pg (g // y) & Friedel";
	if ( sym == (PLANEGROUP_CM_X | SYMMETRY_FRIEDEL) ) return "cm (m // x) & Friedel";
	if ( sym == (PLANEGROUP_CM_Y | SYMMETRY_FRIEDEL) ) return "cm (m // y) & Friedel";
	if ( sym == (PLANEGROUP_P2MM | SYMMETRY_FRIEDEL) ) return "p2mm & Friedel";
	if ( sym == (PLANEGROUP_P2MG_X | SYMMETRY_FRIEDEL) ) return "p2mg (m // x) & Friedel";
	if ( sym == (PLANEGROUP_P2MG_Y | SYMMETRY_FRIEDEL) ) return "p2mg (m // y) & Friedel";
	if ( sym == (PLANEGROUP_P2GG | SYMMETRY_FRIEDEL) ) return "p2gg & Friedel";
	if ( sym == (PLANEGROUP_C2MM | SYMMETRY_FRIEDEL) ) return "c2mm & Friedel";
	if ( sym == (PLANEGROUP_P4 | SYMMETRY_FRIEDEL) ) return "p4 & Friedel";
	if ( sym == (PLANEGROUP_P4MM | SYMMETRY_FRIEDEL) ) return "p4mm & Friedel";
	if ( sym == (PLANEGROUP_P4GM | SYMMETRY_FRIEDEL) ) return "p4gm & Friedel";
	if ( sym == (PLANEGROUP_P3 | SYMMETRY_FRIEDEL) ) return "p3 & Friedel";
	if ( sym == (PLANEGROUP_P3M1 | SYMMETRY_FRIEDEL) ) return "p3m1 & Friedel";
	if ( sym == (PLANEGROUP_P31M | SYMMETRY_FRIEDEL) ) return "p31m & Friedel";
	if ( sym == (PLANEGROUP_P6 | SYMMETRY_FRIEDEL) ) return "p6 & Friedel";
	if ( sym == (PLANEGROUP_P6MM | SYMMETRY_FRIEDEL) ) return "p6mm & Friedel";

	return "(unknown symmetry)";

}

void symmetry_symmetrise_array(fftw_complex *in, signed int width, signed int height, Symmetry sym) {

	ReflectionList *reflections;

	reflections = reflist_new_from_array(in, width, height);
	symmetry_symmetrise(reflections, sym, SYMFLAG_PHASES_KNOWN);
	reflist_fill_array(in, reflections, width, height);
	reflist_free(reflections);

}