/*
 * cflip.c
 *
 * Charge flipping iteration algorithms
 *
 * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
 * (c) 2008 Alex Eggeman <ase25@cam.ac.uk>
 *
 *  synth2d - two-dimensional Fourier synthesis
 *
 */

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

#include <gtk/gtk.h>
#include <fftw3.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>

#include "main.h"
#include "data.h"
#include "displaywindow.h"
#include "symmetry.h"
#include "gtk-symmetry.h"
#include "reflist.h"
#include "dpsynth.h"

static int cf_window_open = 0;

typedef struct struct_cflipdialog {

	/* Dialog box bits */
	GtkWidget *cf_window;
	GtkWidget *cf_go;
	GtkWidget *cf_stop;
	GtkWidget *cf_next;
	GtkWidget *cf_reset;
	GtkWidget *cf_iterations;

	/* Algorithm guts */
	ReflectionList *cf_reflections;		/* List of 'modulus constraints' */
	ReflectionList *cf_listout;		/* List of structure factors returned at the end of an iteration */
	Symmetry cf_symmetry;			/* Symmetry of the algorithm */
	fftw_complex *cf_drive_old;
	fftw_complex *cf_drive;
	fftw_complex *cf_in;			/* Input realspace array */
	fftw_complex *cf_out;			/* Array on the other end of the FFT of 'in' */
	unsigned int *cf_support;
	unsigned int cf_tangent;		/*Coehlo Constraints*/
	double *intensities;			/*Array of pixel intensities for auto-thresholding*/
	double *amplitudes;			/*Array of reflection intiensities for auto-thresholding*/
	fftw_plan cf_plan_i2d;
	fftw_plan cf_plan_d2i;
	unsigned int width;			/* Width of the array in pixels */
	unsigned int height;			/* Height of the array in pixels */
	double delta;
	//DPSynthWindow *dpsynth;

	/* Tracking */

	unsigned int running;
	unsigned int n_iterations;		/* Number of iterations performed so far */
	unsigned int period;			/* Number of iterations per symmetry application */

} CFlipDialog;

int compare_doubles (const void *xp,const void *yp) {

	const double *x, *y;
	int z;

	x=xp;
	y=yp;


	if (*x > *y) {
   		z = 1;
	} else if (*x < *y) {
      		z = -1;
      	} else {
   		z = 0;
  	}
  	return z;
}

static int cflip_tangent (ReflectionList *listout, double emin, double gmin) {

	ReflectionList *strongest_reflections;
	TripletList *triplet_list;
	unsigned int i;
	unsigned int p;
	unsigned int q;
	unsigned int n_total_triplets;
	double E1;
	double E2;
	double E3;
	double G;
	double Th;
	double Bh;
	double *Mh;
	double maxMh;
	double *phase_tangent;
	double alpha_tan, phi_a, phi_tan;
	unsigned int r, t;
	signed int ph, pk, pl, pph, qh, qk, ql, qph, rh, rk, rl, rph;

	strongest_reflections = reflist_new_parent(listout);
	/* Initialise triplet list */
	triplet_list = malloc(sizeof(TripletList));
	triplet_list->n_triplets = 0;

	for ( i=0; i<listout->n_reflections; i++ ) {
		listout->refs[i].parent_index = i;
	}

	for ( i=1; i<listout->n_reflections; i++ ) {
		if ( (listout->refs[i].amplitude >= emin) && (listout->refs[i].h != 0) && (listout->refs[i].k != 0) ) {
			reflist_addref_parent(strongest_reflections, listout->refs[i].h, listout->refs[i].k, listout->refs[i].l,
					listout->refs[i].parent_index);
		}
	}

	Mh = malloc(strongest_reflections->n_reflections*sizeof(double));

	phase_tangent = malloc(strongest_reflections->n_reflections*sizeof(double));

	/* Iterate doubly over the reflections */
	n_total_triplets = 0;
	for ( p=1; p<strongest_reflections->n_reflections; p++ ) {
		for ( q=1; q<strongest_reflections->n_reflections; q++ ) {

			ph = strongest_reflections->refs[p].h;
			pk = strongest_reflections->refs[p].k;
			pl = strongest_reflections->refs[p].l;
			pph = strongest_reflections->refs[p].phase_known;
			qh = strongest_reflections->refs[q].h;
			qk = strongest_reflections->refs[q].k;
			ql = strongest_reflections->refs[q].l;
			qph = strongest_reflections->refs[q].phase_known;

			assert(ph || pk || pl);
			assert(qh || qk || ql);

			/* Determine the third reflection in the triplet */
			rh = -ph-qh;
			rk = -pk-qk;
			if ( pl != ql ) {
				printf("Coelho: pl != ql: %3i %3i\n", pl, ql);
				continue;
			}
			//assert(pl == ql);	/* This failure mode should have been weeded out ages ago */
			rl = pl;

			/* See if the third reflection is available */
			if ( (r = reflist_inlist(strongest_reflections, rh, rk, rl)) ) {

				/* Trust me on this ... */
				E1 = strongest_reflections->parent->refs[strongest_reflections->refs[p].parent_index].amplitude;
				E2 = strongest_reflections->parent->refs[strongest_reflections->refs[q].parent_index].amplitude;
				E3 = strongest_reflections->parent->refs[strongest_reflections->refs[r].parent_index].amplitude;

				G = E1 * E2 * E3;

				rph = strongest_reflections->parent->refs[strongest_reflections->refs[r].parent_index].phase_known;

				/* Convincing enough? */
				if ( G > gmin ) {
					triplet_list->triplets[triplet_list->n_triplets].p.h = ph;
					triplet_list->triplets[triplet_list->n_triplets].p.k = pk;
					triplet_list->triplets[triplet_list->n_triplets].p.l = pl;
					triplet_list->triplets[triplet_list->n_triplets].q.h = qh;
					triplet_list->triplets[triplet_list->n_triplets].q.k = qk;
					triplet_list->triplets[triplet_list->n_triplets].q.l = ql;
					triplet_list->triplets[triplet_list->n_triplets].G = G;
					triplet_list->n_triplets++;
					if ( triplet_list->n_triplets == MAX_TRIPLETS ) {
						printf("DM: Too many triplets!\n");
						return NULL;
					}

				}
				n_total_triplets++;

			}

		}
	}
	maxMh = 0;

	for ( i=0; i<strongest_reflections->n_reflections; i++ ) {

		signed int h = strongest_reflections->refs[i].h;
		signed int k = strongest_reflections->refs[i].k;

		for ( t=0; t<triplet_list->n_triplets; t++ ) {

			signed int ph = triplet_list->triplets[t].p.h;
			signed int pk = triplet_list->triplets[t].p.k;


				/* Does this triplet concern this reflection? */
				if ( (ph == h) && (pk == k) ) {

					signed int qh = triplet_list->triplets[t].q.h;
					signed int qk = triplet_list->triplets[t].q.k;
					signed int rh = -ph-qh;
					signed int rk = -pk-qk;
					double G = triplet_list->triplets[t].G;

					/* Are the other reflections in the triplet still in the list? */
					if ( reflist_inlist_2d(listout, qh, qk) && reflist_inlist_2d(listout, rh, rk) ) {
						Th += G * sin(( qph - rph));
						Bh += G * cos(( qph - rph));
					}


				}

		}
		phase_tangent[i] = atan( Th / Bh );
		Mh[i] = sqrt( (Th*Th) / (Bh*Bh) );
		if ( Mh[i] > maxMh ) {
			maxMh = Mh[i];
		}

	}

	for ( i=0; i<strongest_reflections->n_reflections; i++ ) {

		alpha_tan = Mh[i] / maxMh;

		phi_a = strongest_reflections->parent->refs[strongest_reflections->refs[i].parent_index].phase_known;

		phi_tan = phi_a + alpha_tan*(phase_tangent[i]-phi_a);

		listout->refs[strongest_reflections->refs[i].parent_index].phase_known = phi_tan;

	}

	free(strongest_reflections);
	free(triplet_list);

	return 0;

}

void cflip_mask(CFlipDialog *cflip_dialog, fftw_complex *in) {


	unsigned int width = cflip_dialog->width;
	unsigned int height = cflip_dialog->height;
	unsigned int i;
	fftw_complex *in_new = malloc(width*height*sizeof(fftw_complex));
	bzero(in_new, width*height*sizeof(fftw_complex));
	ReflectionList *reflections;

	reflections = cflip_dialog->cf_reflections;

	for ( i=1; i<reflections->n_reflections; i++ ) {

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

		if ( h < 0 ) { h = width+h; }
		if ( k < 0 ) { k = height+k; }

		in_new[k + height*h][0] = in[k + height*h][0];
		in_new[k + height*h][1] = in[k + height*h][1];

	}

	memcpy(in, in_new, width*height*sizeof(fftw_complex));
	free(in_new);

}

/* Perform a single CF iteration */
void cflip_iteration(CFlipDialog *cflip_dialog) {

	double re, im;
	double am, ph, max, i_max;
	signed int h, k;
	unsigned int j, indexj;
	unsigned int x, y;
	double fzero;
	double emin = 0;
	double gmin = 0;
	double mean, neg;
	double coehlo_thresh;
	unsigned int mean_n;
	double meansq;
	double sig;
	double delta;
	unsigned int i;
	unsigned int width = cflip_dialog->width;
	unsigned int height = cflip_dialog->height;
	ReflectionList *reflections = cflip_dialog->cf_reflections;
	ReflectionList *listout = cflip_dialog->cf_listout;

	cflip_dialog->intensities = malloc(width*height*sizeof(double));
	cflip_dialog->amplitudes = malloc(reflections->n_reflections*sizeof(double));

//	listout = malloc(sizeof(ReflectionList));
	listout = NULL;
	assert(cflip_dialog->cf_drive_old != NULL);
	assert(cflip_dialog->cf_drive != NULL);
	assert(cflip_dialog->cf_in != NULL);
	assert(cflip_dialog->cf_out != NULL);

	if ( !cflip_dialog->cf_tangent ) {

		/* Transform back to reciprocal space and scale.  This time, the driving function is
			being transformed. */
		fftw_execute(cflip_dialog->cf_plan_d2i);
		for ( x=0; x<width; x++ ) {
			for ( y=0; y<((height/2)+1); y++ ) {
				cflip_dialog->cf_in[y+((height/2)+1)*x][0] = cflip_dialog->cf_in[y+((height/2)+1)*x][0] / (width*height);	/* Re */
				cflip_dialog->cf_in[y+((height/2)+1)*x][1] = cflip_dialog->cf_in[y+((height/2)+1)*x][1] / (width*height);	/* Im */
			}
		}

		/* Impose observed intensities */
		for ( i=0; i<reflections->n_reflections; i++ ) {

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

			if ( h < 0 ) { h = width+h; }
			if ( k < 0 ) { k = height+k; }

			re = cflip_dialog->cf_in[k + height*h][0];
			im = cflip_dialog->cf_in[k + height*h][1];
			ph = atan2(im, re);
			am = reflections->refs[i].amplitude;

			cflip_dialog->cf_in[k + height*h][0] = am*cos(ph);
			cflip_dialog->cf_in[k + height*h][1] = am*sin(ph);

		}

		fzero = cflip_dialog->cf_in[0][0];
		cflip_mask(cflip_dialog, cflip_dialog->cf_in);

		/* Impose symmetry */
		if ( cflip_dialog->n_iterations % cflip_dialog->period == 0 ) {

			listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
			symmetry_symmetrise(listout, cflip_dialog->cf_symmetry, SYMFLAG_PHASES_KNOWN);
			reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
			//dpsynth_update(cflip_dialog->dpsynth, listout);
			reflist_free(listout);


		} else {

			listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
			symmetry_symmetrise(listout, cflip_dialog->cf_symmetry & SYMMETRY_FRIEDEL, SYMFLAG_PHASES_KNOWN);
			reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
			//dpsynth_update(cflip_dialog->dpsynth, listout);
			reflist_free(listout);

		}



		/* Transform to real space */
		fftw_execute(cflip_dialog->cf_plan_i2d);

		memcpy(cflip_dialog->cf_out, cflip_dialog->cf_drive, width*height*sizeof(fftw_complex));
		displaywindow_switchview();
		mean = 0; mean_n = 0; meansq = 0; sig = 0; delta = 0;
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				re = cflip_dialog->cf_drive[y+height*x][0];
				im = cflip_dialog->cf_drive[y+height*x][1];
				am = sqrt(re*re + im*im);
				mean += am;
				mean_n += 1;
				meansq += am * am;
			}
		}

		sig = sqrt((meansq/mean_n)-((mean/mean_n) * (mean/mean_n)));
		delta = sig * cflip_dialog->delta;

		/* Determine the support */
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				if ( cflip_dialog->cf_drive[y+height*x][0] > delta ) {
					cflip_dialog->cf_support[y+height*x] = 1;
				} else {
					cflip_dialog->cf_support[y+height*x] = 0;
				}
			}
		}


		/* Impose non-negativity, sharpen peaks and apply feedback */

		/*mean = 0; neg = 0;

		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				re = cflip_dialog->cf_drive[y+height*x][0];
				im = cflip_dialog->cf_drive[y+height*x][1];
				am = sqrt(re*re + im*im);

				mean += am;
				if ( re<0 ) {
					neg += am;
				}

			}
		}

		mean = mean /(width * height);
		
		double entropy;
		entropy = 0;
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				entropy += (am/mean) * log (am/mean);
			}
		}
		entropy = -entropy;
		printf("step %i; Entropy : %f, Negativity : %f \n", cflip_dialog->n_iterations, entropy, neg);*/
	

		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {

				if ( !cflip_dialog->cf_support[y+height*x] ) {

					cflip_dialog->cf_drive[y+height*x][0] = -cflip_dialog->cf_drive[y+height*x][0];
					cflip_dialog->cf_drive[y+height*x][1] = -cflip_dialog->cf_drive[y+height*x][1];

				}

			}
		}

	} else { 	/*Rework the Algorithm from here to incorporate the Coehlo constraints*/

		/* Transform back to reciprocal space and scale.  This time, the driving function is being transformed. */
		fftw_execute(cflip_dialog->cf_plan_d2i);
		for ( x=0; x<width; x++ ) {
			for ( y=0; y<((height/2)+1); y++ ) {
				cflip_dialog->cf_in[y+((height/2)+1)*x][0] = cflip_dialog->cf_in[y+((height/2)+1)*x][0] / (width*height);	/* Re */
				cflip_dialog->cf_in[y+((height/2)+1)*x][1] = cflip_dialog->cf_in[y+((height/2)+1)*x][1] / (width*height);	/* Im */
			}
		}

		/* Do Tangent Formula Wizardry */
		cflip_mask(cflip_dialog, cflip_dialog->cf_in);
		listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);

		i_max = 0;

		for ( i=0; i<listout->n_reflections; i++ ) {

			cflip_dialog->amplitudes[i]=listout->refs[i].amplitude;

			if ( listout->refs[i].amplitude > i_max ) {
				i_max = listout->refs[i].amplitude;
			}
		}

		emin = 0.5 * i_max;
		gmin = 0.2 * i_max;

		cflip_tangent (listout, emin, gmin);

		i_max = 0; j=0; max = 0; indexj=0;

		for ( i=0; i<reflections->n_reflections; i++ ) {

			cflip_dialog->amplitudes[i]=reflections->refs[i].amplitude;
			if ( reflections->refs[i].amplitude > i_max ) {
				i_max = reflections->refs[i].amplitude;
			}
		}

		qsort(cflip_dialog->amplitudes, listout->n_reflections, sizeof(double), compare_doubles);

		j = (int)(listout->n_reflections / 2);
		j = (int)(listout->n_reflections / 4);
		max = cflip_dialog->amplitudes[j];

		/* Impose observed intensities */
		//unsigned int n_sub = 0;
		for ( i=0; i<listout->n_reflections; i++ ) {

			h = listout->refs[i].h;
			k = listout->refs[i].k;

			indexj = reflist_inlist(reflections, h, k, 0);

			if ( h < 0 ) { h = width+h; }
			if ( k < 0 ) { k = height+k; }

			ph = listout->refs[i].phase_known;
			am = reflections->refs[indexj].amplitude;
			if (am < max ) {
				am = 0;
        			//n_sub += 1;
			}

			//printf("number of Zero reflections is %i.\n", n_sub);
			cflip_dialog->cf_in[k + height*h][0] = am*cos(ph);
			cflip_dialog->cf_in[k + height*h][1] = am*sin(ph);

		}

		fzero = cflip_dialog->cf_in[0][0];
		//free(listout);

		/* Impose symmetry */
		if ( cflip_dialog->n_iterations % cflip_dialog->period == 0 ) {

			listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
			symmetry_symmetrise(listout, cflip_dialog->cf_symmetry, SYMFLAG_PHASES_KNOWN);
			reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
			//dpsynth_update(cflip_dialog->dpsynth, listout);
			reflist_free(listout);


		} else {

			listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
			symmetry_symmetrise(listout, cflip_dialog->cf_symmetry & SYMMETRY_FRIEDEL, SYMFLAG_PHASES_KNOWN);
			reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
			//dpsynth_update(cflip_dialog->dpsynth, listout);
			reflist_free(listout);

		}

	        /* Transform to real space */
		fftw_execute(cflip_dialog->cf_plan_i2d);

	        max = 0;

		memcpy(cflip_dialog->cf_out, cflip_dialog->cf_drive, width*height*sizeof(fftw_complex));
		displaywindow_switchview();

		/* Find the Threshold Real-Space Intensity */
		mean = 0; mean_n = 0; sig = 0; delta = 0;
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				re = cflip_dialog->cf_drive[y+height*x][0];
				im = cflip_dialog->cf_drive[y+height*x][1];
				am = sqrt((re*re) + (im*im));
				if ( am > max ) {
				        max = am;
				}

			}
		}

		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				re = cflip_dialog->cf_drive[y+height*x][0];
				im = cflip_dialog->cf_drive[y+height*x][1];
				am = sqrt(re*re + im*im) / max;
				mean += am;
				cflip_dialog->intensities[y+height*x]=am;
			}
		}

		qsort(cflip_dialog->intensities, mean_n, sizeof(double), compare_doubles);

	        coehlo_thresh = 0;
	        i=0;
		while ( coehlo_thresh < (0.6 * mean) ) {

			coehlo_thresh += cflip_dialog->intensities[i];
			delta = cflip_dialog->intensities[i];
			i += 1;
		}

		/* Determine the support */
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				if ( cflip_dialog->cf_drive[y+height*x][0] > delta ) {
					cflip_dialog->cf_support[y+height*x] = 1;
				} else {
					cflip_dialog->cf_support[y+height*x] = 0;
				}
			}
		}

		/* Perform the Charge Flip */
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {

				if ( !cflip_dialog->cf_support[y+height*x] ) {

					cflip_dialog->cf_drive[y+height*x][0] = -cflip_dialog->cf_drive[y+height*x][0];
					cflip_dialog->cf_drive[y+height*x][1] = -cflip_dialog->cf_drive[y+height*x][1];

				} else {
					//cflip_dialog->cf_drive[y+height*x][0] = delta + sqrt((cflip_dialog->cf_drive[y+height*x][0])-delta);
					cflip_dialog->cf_drive[y+height*x][0] = cflip_dialog->cf_drive[y+height*x][0];
					cflip_dialog->cf_drive[y+height*x][1] = cflip_dialog->cf_drive[y+height*x][1];

				}
			}
		}
	}
free(cflip_dialog->intensities);
free(cflip_dialog->amplitudes);
}

static gint cflip_threshsup_deltachanged(GtkWidget *widget, CFlipDialog *cflip_dialog) {

	cflip_dialog->delta = gtk_range_get_value(GTK_RANGE(widget));
	return 0;
}

void cflip_finalise(CFlipDialog *cflip_dialog) {

	free(cflip_dialog->cf_drive);		cflip_dialog->cf_drive = NULL;
	free(cflip_dialog->cf_drive_old);	cflip_dialog->cf_drive_old = NULL;
	free(cflip_dialog->cf_support);		cflip_dialog->cf_support = NULL;
	free(cflip_dialog->cf_in);		cflip_dialog->cf_in = NULL;
	reflist_free(cflip_dialog->cf_reflections);	cflip_dialog->cf_reflections = NULL;	
	/* Don't free cf_out, since it is now backing up the display */

	fftw_destroy_plan(cflip_dialog->cf_plan_i2d);
	fftw_destroy_plan(cflip_dialog->cf_plan_d2i);

}


static gint cflip_window_periodchanged(GtkWidget *cfw_period, CFlipDialog *cflip_dialog) {

	const char *str;

	str = gtk_entry_get_text(GTK_ENTRY(cfw_period));
	cflip_dialog->period = atoi(str);
	return 0;
}

static void cflip_window_close(GtkWidget *widget, gint arg1, CFlipDialog *cflip_dialog) {

	if  ( cflip_dialog->running ) {
		cflip_dialog->running = 0;
	}
	cf_window_open = 0;

	cflip_finalise(cflip_dialog);

	gtk_widget_destroy(widget);

}

static void cflip_update_iterations(CFlipDialog *cflip_dialog) {

	char *text = malloc(40);
	snprintf(text, 39, "Iterations performed so far: %i", cflip_dialog->n_iterations);
	gtk_label_set_text(GTK_LABEL(cflip_dialog->cf_iterations), text);
	free(text);

}

void cflip_reset(CFlipDialog *cflip_dialog) {

	reflist_free(cflip_dialog->cf_reflections);
	cflip_dialog->cf_reflections = reflist_copy(main_reflist());

	unsigned int width = cflip_dialog->width;
	unsigned int height = cflip_dialog->height;
	cflip_dialog->n_iterations = 0;
	cflip_update_iterations(cflip_dialog);
	ReflectionList *reflections = cflip_dialog->cf_reflections;

	unsigned int i;

	for ( i=0; i<reflections->n_reflections; i++ ) {

		double am, ph;
		signed int h, k;

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

		if ( h < 0 ) { h = width+h; }
		if ( k < 0 ) { k = height+k; }

		am = reflections->refs[i].amplitude;
		ph = (((double)random())/RAND_MAX) * (2*M_PI);

		cflip_dialog->cf_in[k + height*h][0] = am*cos(ph);
		cflip_dialog->cf_in[k + height*h][1] = am*sin(ph);

	}
	symmetry_symmetrise_array(cflip_dialog->cf_in, width, height, cflip_dialog->cf_symmetry & SYMMETRY_FRIEDEL);
	fftw_execute(cflip_dialog->cf_plan_i2d);

	memcpy(cflip_dialog->cf_out, cflip_dialog->cf_drive, width*height*sizeof(fftw_complex));
	displaywindow_switchview();

	//dpsynth_update(cflip_dialog->dpsynth, reflections);

}

void cflip_initialise(CFlipDialog *cflip_dialog) {

	ReflectionList *reflections = reflist_copy(main_reflist());
	Symmetry symmetry;
	cflip_dialog->width = data_width();
	cflip_dialog->height = data_height();
	cflip_dialog->n_iterations = 0;
	cflip_dialog->period = 1;
	unsigned int width = cflip_dialog->width;
	unsigned int height = cflip_dialog->height;
	symmetry = cflip_dialog->cf_symmetry;

	cflip_dialog->cf_drive_old = fftw_malloc(width * height * sizeof(fftw_complex));
	cflip_dialog->cf_drive = fftw_malloc(width * height * sizeof(fftw_complex));
	cflip_dialog->cf_support = fftw_malloc(width * height * sizeof(unsigned int));
	cflip_dialog->cf_in = fftw_malloc(width * height * sizeof(fftw_complex));
	memset(cflip_dialog->cf_in, 0, width * height * sizeof(fftw_complex));
	cflip_dialog->cf_out = fftw_malloc(width * height * sizeof(fftw_complex));

	cflip_dialog->cf_plan_i2d = fftw_plan_dft_2d(width, height, cflip_dialog->cf_in, cflip_dialog->cf_drive, FFTW_BACKWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT);
	cflip_dialog->cf_plan_d2i = fftw_plan_dft_2d(width, height, cflip_dialog->cf_drive, cflip_dialog->cf_in, FFTW_FORWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT);

	cflip_dialog->cf_reflections = reflections;

	displaywindow_set_realspace(cflip_dialog->cf_out, DWR_GSF);
	//cflip_dialog->dpsynth = dpsynth_open(reflections, "Current Diffraction Pattern", 1);

	cflip_reset(cflip_dialog);

}

static gint cflip_window_go(GtkWidget *widget, CFlipDialog *cflip_dialog) {


	if ( cflip_dialog->running ) {
		return 0;
	}
	cflip_dialog->running = 1;
	while ( cflip_dialog->running ) {
		cflip_dialog->n_iterations++;
		cflip_update_iterations(cflip_dialog);
		cflip_iteration(cflip_dialog);
		while ( gtk_events_pending() ) gtk_main_iteration();
	}
	return 0;

}

static gint cflip_window_next(GtkWidget *widget, CFlipDialog *cflip_dialog) {

		cflip_dialog->n_iterations++;
		cflip_update_iterations(cflip_dialog);
		cflip_iteration(cflip_dialog);

	return 0;
}

static gint cflip_window_stop(GtkWidget *widget, CFlipDialog *cflip_dialog) {

	cflip_dialog->running = 0;
	return 0;

}

static gint cflip_window_reset(GtkWidget *widget, CFlipDialog *cflip_dialog) {

	cflip_dialog->running = 0;
	cflip_reset(cflip_dialog);
	return 0;
}

static gint cflip_symmetry_changed(GtkWidget *widget, CFlipDialog *cflip_dialog) {

	cflip_dialog->cf_symmetry = gtk_symmetry_get_symmetry(GTK_SYMMETRY(widget));
	return 0;
}

GtkWidget *cfw_tangent;
void cflip_tangent_toggled(GtkWidget *widget, CFlipDialog *cflip_dialog) {

	cflip_dialog->cf_tangent = 1-cflip_dialog->cf_tangent;

}

void cflip_dialog_open() {

	GtkWidget *cfw_window;
	GtkWidget *cfw_period;
	GtkWidget *cfw_period_hbox;
	GtkWidget *cfw_threshsup_delta;
	GtkWidget *cfw_threshsup_delta_label;
	GtkWidget *cfw_threshsup_hbox;
	GtkWidget *vbox;
	GtkWidget *hbox;
	GtkWidget *sym_define;
	GtkWidget *cfw_support_label;
	GtkWidget *cfw_process_label;
	GtkWidget *cfw_period_label;
	GtkWidget *cfw_gostop_hbox;
	GtkWidget *cfw_tangent;
	char *text;
	CFlipDialog *cflip_dialog;

	cflip_dialog = malloc(sizeof(CFlipDialog));

	if ( cf_window_open ) {
		return;
	}
	cf_window_open = 1;

	cflip_dialog->running = 0;
	cflip_dialog->cf_tangent = 0;
	cflip_dialog->n_iterations = 0;
	cflip_dialog->delta = 0;

	cfw_window = gtk_dialog_new_with_buttons("Charge Flipping", GTK_WINDOW(displaywindow_gtkwindow()),
		GTK_DIALOG_DESTROY_WITH_PARENT, GTK_STOCK_CLOSE, GTK_RESPONSE_CLOSE, NULL);

	vbox = gtk_vbox_new(FALSE, 0);
	hbox = gtk_hbox_new(TRUE, 0);
	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(cfw_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7);
	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5);

	cflip_dialog->cf_iterations = gtk_label_new("");
	text = malloc(40);
	snprintf(text, 39, "Iterations performed so far: %i", cflip_dialog->n_iterations);
	gtk_label_set_text(GTK_LABEL(cflip_dialog->cf_iterations), text);
	free(text);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cflip_dialog->cf_iterations), FALSE, FALSE, 5);

	sym_define = gtk_symmetry_new(2, 2, TRUE);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(sym_define), FALSE, FALSE, 5);
	g_signal_connect(G_OBJECT(sym_define), "changed", G_CALLBACK(cflip_symmetry_changed), cflip_dialog);
	cflip_dialog->cf_symmetry = PLANEGROUP_P1;

	cfw_period_label = gtk_label_new("Iterations to symmetrize:");
	gtk_misc_set_alignment(GTK_MISC(cfw_period_label), 1, 0.5);
	cfw_period = gtk_entry_new();
	gtk_entry_set_width_chars(GTK_ENTRY(cfw_period), 5);
	gtk_entry_set_text(GTK_ENTRY(cfw_period), "1");
	cfw_period_hbox = gtk_hbox_new(FALSE, 0);
	gtk_box_pack_start(GTK_BOX(cfw_period_hbox), GTK_WIDGET(cfw_period_label), FALSE, FALSE, 5);
	gtk_box_pack_start(GTK_BOX(cfw_period_hbox), GTK_WIDGET(cfw_period), FALSE, FALSE, 5);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_period_hbox), FALSE, FALSE, 5);
	g_signal_connect(G_OBJECT(cfw_period), "activate", G_CALLBACK(cflip_window_periodchanged), cflip_dialog);

	cfw_support_label = gtk_label_new("");
	gtk_label_set_markup(GTK_LABEL(cfw_support_label), "<span weight=\"bold\">Threshold</span>");
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_support_label), FALSE, FALSE, 10);

	cfw_tangent = gtk_check_button_new_with_label("Coehlo Constraints");
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_tangent), FALSE, FALSE, 5);
	if ( cflip_dialog->cf_tangent ) {
		gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(cfw_tangent), TRUE);
	} else {
		gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(cfw_tangent), FALSE);
	}

	cfw_process_label = gtk_label_new("Maximum unflipped intensity");
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_process_label), FALSE, FALSE, 5);

	cfw_threshsup_hbox = gtk_hbox_new(FALSE, 0);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_threshsup_hbox), FALSE, FALSE, 5);
	cfw_threshsup_delta_label = gtk_label_new("δ:");
	gtk_box_pack_start(GTK_BOX(cfw_threshsup_hbox), GTK_WIDGET(cfw_threshsup_delta_label), FALSE, FALSE, 5);
	cfw_threshsup_delta = gtk_hscale_new_with_range(-5, 5, 0.1);
	gtk_scale_set_value_pos(GTK_SCALE(cfw_threshsup_delta), GTK_POS_RIGHT);
	gtk_box_pack_start(GTK_BOX(cfw_threshsup_hbox), GTK_WIDGET(cfw_threshsup_delta), TRUE, TRUE, 5);
	gtk_range_set_value(GTK_RANGE(cfw_threshsup_delta), cflip_dialog->delta);

	gtk_widget_set_sensitive(cfw_threshsup_delta, TRUE);
	gtk_widget_set_sensitive(cfw_threshsup_delta_label, TRUE);

	cfw_gostop_hbox = gtk_hbox_new(FALSE, 0);
	cflip_dialog->cf_go = gtk_button_new_with_label("Run");
	gtk_button_set_image(GTK_BUTTON(cflip_dialog->cf_go), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PLAY, GTK_ICON_SIZE_BUTTON));
	cflip_dialog->cf_next = gtk_button_new_with_label("Next");
	gtk_button_set_image(GTK_BUTTON(cflip_dialog->cf_next), gtk_image_new_from_stock(GTK_STOCK_MEDIA_NEXT, GTK_ICON_SIZE_BUTTON));
	cflip_dialog->cf_stop = gtk_button_new_from_stock(GTK_STOCK_MEDIA_STOP);
	cflip_dialog->cf_reset = gtk_button_new_with_label("Reset");
	gtk_button_set_image(GTK_BUTTON(cflip_dialog->cf_reset), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PREVIOUS, GTK_ICON_SIZE_BUTTON));
	gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_reset), FALSE, FALSE, 5);
	gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_stop), FALSE, FALSE, 5);
	gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_next), FALSE, FALSE, 5);
	gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_go), FALSE, FALSE, 5);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_gostop_hbox), FALSE, FALSE, 5);
	g_signal_connect(G_OBJECT(cflip_dialog->cf_go), "clicked", G_CALLBACK(cflip_window_go), cflip_dialog);
	g_signal_connect(G_OBJECT(cflip_dialog->cf_stop), "clicked", G_CALLBACK(cflip_window_stop), cflip_dialog);
	g_signal_connect(G_OBJECT(cflip_dialog->cf_next), "clicked", G_CALLBACK(cflip_window_next), cflip_dialog);
	g_signal_connect(G_OBJECT(cflip_dialog->cf_reset), "clicked", G_CALLBACK(cflip_window_reset), cflip_dialog);
	g_signal_connect(G_OBJECT(cfw_tangent), "toggled", G_CALLBACK(cflip_tangent_toggled), cflip_dialog);
	g_signal_connect(G_OBJECT(cfw_window), "response", G_CALLBACK(cflip_window_close), cflip_dialog);
	g_signal_connect(G_OBJECT(cfw_threshsup_delta), "value-changed", G_CALLBACK(cflip_threshsup_deltachanged), cflip_dialog);
	g_signal_connect(G_OBJECT(displaywindow_gtkwindow()), "delete-event", G_CALLBACK(cflip_window_stop), cflip_dialog);

	cflip_initialise(cflip_dialog);
	displaywindow_forceview(DWV_REALSPACE);

	gtk_widget_show_all(cfw_window);

}