/*
 * main.c
 *
 * The Top Level Source File
 *
 * (c) 2006-2008 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2d - Two-Dimensional Crystallographic Fourier Synthesis
 *
 */

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

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

#include "displaywindow.h"
#include "data.h"
#include "main.h"
#include "reflist.h"
#include "normalise.h"
#include "dpsynth.h"
#include "argand.h"
#include "symmetry.h"
#include "gsf.h"
#include "cflip.h"
#include "cdm.h"
#include "clean.h"
#include "model.h"
#include "geometry.h"
#include "statistics.h"
#include "elements.h"
#include "refine.h"
#include "superlattice.h"
#include "options.h"


ReflectionList *main_reflections = NULL;

/* A rather shameless glue layer to avoid sharing information */
void main_show_patterson() { displaywindow_show_patterson(main_reflections); }
void main_show_pattersone() { displaywindow_show_pattersone(main_reflections); }
void main_show_knownphases() { displaywindow_show_knownphases(main_reflections); }
void main_show_calcphases() { displaywindow_show_calcphases(main_reflections); }
void main_show_difference() {
	ReflectionList *reflections = reflist_copy(main_reflections);
	model_calculate_difference_coefficients(reflections);
	displaywindow_show_difference(reflections);
	reflist_free(reflections);
}
void main_show_refsyn() {
	ReflectionList *reflections = reflist_copy(main_reflections);
	model_calculate_refinement_coefficients(reflections);
	displaywindow_show_refsyn(reflections);
	reflist_free(reflections);
}
void main_show_diffpatt() {
	ReflectionList *reflections = reflist_copy(main_reflections);
	model_calculate_difference_coefficients(reflections);
	displaywindow_show_diffpatt(reflections);
	reflist_free(reflections);
}
void main_wilsonplot() { normalise_wilsonplot(main_reflections); }
void main_falloffplot() { normalise_falloffplot(main_reflections); }
void main_dpsynth() { dpsynth_main_open(main_reflections); }
void main_dpsynth_update() { dpsynth_main_update(main_reflections); }
void main_argand() { argand_open(main_reflections); }
void main_argand_update() { argand_update(main_reflections); }
void main_dethermalise(double level) { normalise_dethermalise(main_reflections, level); }
void main_normalise_exponential(double a, double b, double c) { normalise_exponential(main_reflections, a, b, c); }
void main_normalise(double level) { normalise_execute(main_reflections, level); }
void main_gsf_initialise() { gsf_initialise(main_reflections); }
void main_gsf_reset() { gsf_reset(main_reflections); }
void main_aperture_open() { clean_aperture_open(main_reflections); }
unsigned int main_cdm_tangentexpansion(CDMContext *cdm) {
	cdm->reflections = main_reflections;
	return cdm_tangentexpansion(cdm);
}
void main_display_phasing_solution(CDMContext *cdm, PhasingSolution *sol) {
	cdm_display_phasing_solution(cdm, sol, main_reflections);
}
void main_symmetrise(Symmetry symmetry) {

	char string[256];
	double rsym;
	SymFlags flags = 0;

	if ( displaywindow_mode() == DWV_KNOWNPHASE ) {
		flags = flags | SYMFLAG_PHASES_KNOWN;
	} else if ( displaywindow_mode() == DWV_CALCPHASE ) {
		flags = flags | SYMFLAG_PHASES_CALC;
	} /* else amplitudes only */

	rsym = symmetry_symmetrise(main_reflections, symmetry, flags);

	main_displayr();
	main_dpsynth_update();

	displaywindow_switchview();
	snprintf(string, 255, "Rsym=%.2f%%", rsym*100);
	displaywindow_statusbar(string);

}
void main_geometry_correct(GeometryCorrectionType correction_type, double wavenumber, double circleradius) {
	geometry_correct(main_reflections, correction_type, wavenumber, circleradius);
}

void main_superlattice_split(unsigned int xc, unsigned int yc) {

	char string[256];

	main_substitutereflections(superlattice_split(main_reflections, xc, yc));
	displaywindow_switchview();
	snprintf(string, 255, "Split %ix%i superlattice", xc, yc);
	displaywindow_statusbar(string);

}

/* Cut out reflections with |g|<2.5 nm^-1 */
void main_hpfilter() {

	unsigned int i;
	double a, b, c;
	char string[256];
	unsigned int n_del, n_orig;

	a = data_a();  b = data_b();  c = data_c();
	n_del = 0;  n_orig = main_reflections->n_reflections;

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

		double g;
		signed int h = main_reflections->refs[i].h;
		signed int k = main_reflections->refs[i].k;
		signed int l = main_reflections->refs[i].l;

		g = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)) + ((l*l)/(c*c)));

		if ( g < 2.5 ) {
			reflist_delref(main_reflections, h, k, l);
			n_del++;
		}
	}

	displaywindow_switchview();
	snprintf(string, 255, "%0.2f%% of reflections filtered out (|g|<2.5nm^-1)", ((double)n_del/n_orig)*100);
	displaywindow_statusbar(string);

}

void main_savereflections(const char *filename) {

	FILE *fh;
	unsigned int i;
	int sim;

	fh = fopen(filename, "w");

	fprintf(fh, "a %f\n", data_a());
	fprintf(fh, "b %f\n", data_b());
	fprintf(fh, "c %f\n", data_c());
	fprintf(fh, "angle %f\n", 180*data_gamma()/M_PI);
	fprintf(fh, "scale %i\n", data_get_image_scale());

	switch ( displaywindow_mode() ) {
		case DWV_PATTERSON	: sim = 0; break;
		case DWV_PATTERSONE	: sim = 0; break;
		case DWV_KNOWNPHASE	: sim = 0; break;
		/* DWV_REALSPACE should never get here */
		case DWV_DIFFERENCE	: sim = 1; break;
		case DWV_DIFFPATT	: sim = 1; break;
		case DWV_REFSYN		: sim = 1; break;
		case DWV_MODEL		: sim = 1; break;
		case DWV_SIMPATT	: sim = 1; break;
		case DWV_SIMFOLZPATT	: sim = 1; break;
		case DWV_EXITWAVE	: sim = 1; break;
		default			: sim = 0; break;
	}

	if ( sim ) {

		ReflectionList *calc;

		switch ( displaywindow_mode() ) {
			case DWV_DIFFERENCE	: calc = reflist_copy(main_reflections);
						  model_calculate_difference_coefficients(calc); break;
			case DWV_DIFFPATT	: calc = reflist_copy(main_reflections);
						  model_calculate_difference_coefficients(calc); break;
			case DWV_REFSYN		: calc = reflist_copy(main_reflections);
						  model_calculate_refinement_coefficients(calc); break;
			case DWV_MODEL		: calc = model_calculate_f(NULL, NULL, 0); break;
			case DWV_SIMPATT	: calc = model_calculate_f(NULL, NULL, 0); break;
			case DWV_SIMFOLZPATT	: calc = model_calculate_f(NULL, NULL, 1); break;
			case DWV_EXITWAVE	: calc = model_calculate_f(NULL, NULL, 0); break;
			default			: calc = model_calculate_f(NULL, NULL, 0); break;
		}

		for ( i=1; i<calc->n_reflections; i++ ) {
			signed int h = calc->refs[i].h;
			signed int k = calc->refs[i].k;
			signed int l = calc->refs[i].l;
			double am = calc->refs[i].amplitude;
			if ( (displaywindow_mode() != DWV_DIFFPATT) && (displaywindow_mode() != DWV_SIMPATT) && (displaywindow_mode() != DWV_SIMFOLZPATT) ) {
				fprintf(fh, "%3i %3i %3i %8f %8f\n", h, k, l, am, calc->refs[i].phase_known);
			} else {
				fprintf(fh, "%3i %3i %3i %8f\n", h, k, l, am);
			}
		}
		reflist_free(calc);

	} else {

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

			signed int h = main_reflections->refs[i].h;
			signed int k = main_reflections->refs[i].k;
			signed int l = main_reflections->refs[i].l;
			double am = main_reflections->refs[i].amplitude;

			switch ( displaywindow_mode() ) {

				case DWV_PATTERSON : fprintf(fh, "%3i %3i %3i %8f\n", h, k, l, am); break;
				case DWV_PATTERSONE : fprintf(fh, "%3i %3i %3i %8f\n", h, k, l, am); break;
				case DWV_KNOWNPHASE : if ( main_reflections->refs[i].phase_known_set )
							fprintf(fh, "%3i %3i %3i %8f %8f\n", h, k, l, am, main_reflections->refs[i].phase_known);
							break;
				case DWV_CALCPHASE : if ( main_reflections->refs[i].phase_calc_set )
							fprintf(fh, "%3i %3i %3i %8f %8f\n", h, k, l, am, main_reflections->refs[i].phase_calc);
							break;
				default : break;

			}

		}

	}

	fclose(fh);

}

static void main_reset() {

	unsigned int i;
	unsigned int phases_known;

	if ( main_reflections ) {
		free(main_reflections);
	}
	main_reflections = reflist_new();
	memcpy(main_reflections, data_getreflections(), sizeof(ReflectionList));
	data_free();

	/* Decide whether to start with known phases or Patterson */
	phases_known = 0;
	for ( i=0; i<main_reflections->n_reflections; i++ ) {
		if ( main_reflections->refs[i].phase_known_set ) {
			phases_known++;;
		}
	}
	printf("MA: Have phase values for %i of %i reflections\n", phases_known, main_reflections->n_reflections);
	if ( phases_known == 0 ) {
		displaywindow_show_patterson(main_reflections);
		printf("MA: No phase values known - displaying Patterson map instead.\n");
		displaywindow_statusbar("No phase values known - displaying Patterson map instead");
	} else {
		displaywindow_show_knownphases(main_reflections);
		displaywindow_statusbar("Displaying known phases from file");
	}

	displaywindow_brightness_auto(NULL, NULL);

}

int main(int argc, char *argv[]) {

	/* Sort out configuration */
	options_load();

	g_thread_init(NULL);
	gtk_init(&argc, &argv);

	if ( argc != 2 ) {
		fprintf(stderr, "Syntax: %s <data file>\n", argv[0]);
		return 1;
	}

	gsl_set_error_handler_off();

	/* Read data */
	if ( data_read(argv[1]) ) {
		return 1;
	}

	elements_initialise();
	model_default();

	/* Open main window */
	displaywindow_open(argv[1]);

	/* Execute first transform */
	main_reset();

	gtk_main();

	options_save();

	return 0;

}

void main_substitutereflections(ReflectionList *new) {
	memcpy(main_reflections, new, sizeof(ReflectionList));
	displaywindow_switchview();
}

void main_displayr() {

	ReflectionList *model_reflections;
	char r_string[256];

	model_reflections = model_calculate_f(main_reflections, NULL, 69);
	dpsynth_simdp_update(model_reflections);
	if ( !model_current_is_blank() ) {
		snprintf(r_string, 255, "R1=%.2f%%, R2=%.2f%%", 100*stat_r(main_reflections, model_reflections),
								100*stat_r2(main_reflections, model_reflections));
		displaywindow_statusbar(r_string);
	}
	reflist_free(model_reflections);

}

/* Caution! */
ReflectionList *main_reflist() { return main_reflections; }

unsigned int main_max_h() {
	unsigned int max, i;
	max = 0;
	for ( i=1; i<main_reflections->n_reflections; i++ ) {
		if ( abs(main_reflections->refs[i].h) > max ) max = abs(main_reflections->refs[i].h);
	}
	return max;
}

unsigned int main_max_k() {
	unsigned int max, i;
	max = 0;
	for ( i=1; i<main_reflections->n_reflections; i++ ) {
		if ( abs(main_reflections->refs[i].k) > max ) max = abs(main_reflections->refs[i].k);
	}
	return max;
}

void main_stripzero() {

	int i;

	i = reflist_inlist(main_reflections, 0, 0, 0);
	assert(i == 0);
	main_reflections->refs[i].amplitude = 0.0;

}

/* Anti-alias the pattern by restricting resolution to a circle */
void main_antialias() {

	unsigned int i;
	double max;
	double ph, pk, mh, mk;

	/* Find the resolution in each direction */
	max = 0.0;
	for ( i=1; i<main_reflections->n_reflections; i++ ) {
		double res;
		/* Find reflections on centre line */
		if ( main_reflections->refs[i].k != 0 ) continue;
		if ( main_reflections->refs[i].h < 0 ) continue;
		res = resolution(main_reflections->refs[i].h,
				 main_reflections->refs[i].k,
				 main_reflections->refs[i].l,
				 data_a(), data_b(), data_c(), data_gamma());
		/* Find highest resolution in this direction */
		if ( res > max ) max = res;
	}
	ph = max;

	max = 0.0;
	for ( i=1; i<main_reflections->n_reflections; i++ ) {
		double res;
		/* Find reflections on centre line */
		if ( main_reflections->refs[i].k != 0 ) continue;
		if ( main_reflections->refs[i].h > 0 ) continue;
		res = resolution(main_reflections->refs[i].h,
				 main_reflections->refs[i].k,
				 main_reflections->refs[i].l,
				 data_a(), data_b(), data_c(), data_gamma());
		/* Find highest resolution in this direction */
		if ( res > max ) max = res;
	}
	mh = max;

	max = 0.0;
	for ( i=1; i<main_reflections->n_reflections; i++ ) {
		double res;
		/* Find reflections on centre line */
		if ( main_reflections->refs[i].h != 0 ) continue;
		if ( main_reflections->refs[i].k < 0 ) continue;
		res = resolution(main_reflections->refs[i].h,
				 main_reflections->refs[i].k,
				 main_reflections->refs[i].l,
				 data_a(), data_b(), data_c(), data_gamma());
		/* Find highest resolution in this direction */
		if ( res > max ) max = res;
	}
	pk = max;

	max = 0.0;
	for ( i=1; i<main_reflections->n_reflections; i++ ) {
		double res;
		/* Find reflections on centre line */
		if ( main_reflections->refs[i].h != 0 ) continue;
		if ( main_reflections->refs[i].k > 0 ) continue;
		res = resolution(main_reflections->refs[i].h,
				 main_reflections->refs[i].k,
				 main_reflections->refs[i].l,
				 data_a(), data_b(), data_c(), data_gamma());
		/* Find highest resolution in this direction */
		if ( res > max ) max = res;
	}
	mk = max;

	/* Find the smallest */
	double hm, km, m;
	if ( mh < ph ) hm = mh; else hm = ph;
	if ( mk < pk ) km = mk; else km = pk;
	if ( hm < km ) m = hm; else m=km;

	printf("Data resolution is %f nm^-1\n", ph);

	i = 1;
	while ( i < main_reflections->n_reflections ) {

		double res;
		signed int h, k, l;

		h = main_reflections->refs[i].h;
		k = main_reflections->refs[i].k;
		l = main_reflections->refs[i].l;
		res = resolution(h, k, l, data_a(), data_b(), data_c(),
					  data_gamma());
		//printf("Resolution of %3i %3i %3i is %f - ", h, k, l, res);
		if ( res >= m ) {
			reflist_delref(main_reflections, h, k, l);
		//	printf("eliminated\n");
		} else {
		//	main_reflections->refs[i].amplitude = 1000000000;
		//	printf("maximised\n");
			i++;
		}

	}

}