/*
 * geometry.c
 *
 * Geometrical (Lorentz) Corrections
 *
 * (c) 2006 Thomas White <taw27@cam.ac.uk>
 *  Synth2D - two-dimensional Fourier synthesis
 *
 */

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

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

#include "main.h"

#include "displaywindow.h"
#include "data.h"

typedef struct {

	GtkWidget *wn_entry;	/* Electon wavenumber */
	GtkWidget *cr_entry;	/* Laue Circle radius */
	GtkWidget *wn_label;
	GtkWidget *cr_label;
	GtkWidget *wn_nm_label;
	GtkWidget *cr_nm_label;
	GeometryCorrectionType correction_type;
	
} GeometryCorrectionWindow;

static int geometry_window_open = 0;

/* Apply "Lorentz" factor to reflections */
void geometry_correct(ReflectionList *reflections, GeometryCorrectionType correction_type, double wavenumber, double R) {

	unsigned int i;
	double a = data_a();
	double b = data_b();
	double c = data_c();
	
	for ( i=1; i<reflections->n_reflections; i++ ) {
	
		double corr = 1;
		double g, g_xy, g_z;
		signed int h = reflections->refs[i].h;
		signed int k = reflections->refs[i].k;
		signed int l = reflections->refs[i].l;
	
		/* Calculate g, g_xy and g_z */
		g = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)) + ((l*l)/(c*c)));
		g_xy = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)));
		g_z = l/c;
		
		switch ( correction_type ) {
		
			case GEOMETRY_CORRECTION_TYPE_VINCENT_MIDGLEY : {
		
				/* Vincent-Midgley */
				/* eqn (9) from Ultramicroscopy vol 53 1994 p271-282 */
				double corr1, corr2, corr3;

				corr1 = g/(2*R);
				corr2 = pow(corr1, 2);
				corr3 = 1 - corr2;
				corr = sqrt(corr3);
				
				break;
			
			}
		
			case GEOMETRY_CORRECTION_TYPE_GJONNES : {
			
				/* Gjonnes (without divergence) */
				/* eqn (9) from Ultramicroscopy 69 1997 p1-11 */
				double corr1, corr2, corr3, corr4;

				corr1 = 2*R*g_xy;	corr1 = pow(corr1, 2);
				corr2 = pow(g,2) - (2*wavenumber*g_z);	corr2 = pow(corr2, 2);
				corr3 = corr1 - corr2;
				corr4 = 2*wavenumber;
				corr = sqrt(corr3) / (2 * corr4);
				
				break;
			
			}

			case GEOMETRY_CORRECTION_TYPE_GJONNES_DIVERGENCE : {
		
				/* Gjønnes (with divergence) */
				/* eqn (37) from Ultramicroscopy 69 1997 p1-11 */
				double corr1, corr2, corr3, corr3b, corr4, corr4a, corr4b;

				corr1 = pow(R, 3) / pow(wavenumber, 2);
				corr2 = pow((g_xy / g), 2);
				corr3b = (pow(g, 2) - 2*wavenumber*g_z) / (2*R*g_xy);
				corr3 = 1 - pow(corr3b, 2);
				corr4a = pow(((pow(g, 2) - (2*wavenumber*g_z))/wavenumber), 2);
				corr4a = 1 - (corr4a / 2);
				corr4b = (wavenumber*g_z)/pow(g, 2);
				corr4b = 1 - (corr4b * 2);
				corr4 = corr4a / corr4b;
				corr = corr1 * corr2 * corr3 * corr4;
				
				//printf("corr1=%f, corr2=%f, corr3b=%f, corr3=%f, corr4a=%f, corr4b=%f, corr4=%f\n",
				//	corr1, corr2, corr3b, corr3, corr4a, corr4b, corr4);
				
				break;
			
			}
	
		}
		
		printf("%f %f\n", g, corr);
		
		if ( corr < 0 ) corr = 0;
		
		/* Square root here because we are working with amplitudes not intensities */
		reflections->refs[i].amplitude = sqrt(corr) * reflections->refs[i].amplitude;
		assert( reflections->refs[i].amplitude >= 0 );
		
	}
	
	displaywindow_switchview();

}

static gint geometry_window_response(GtkWidget *geometry_window, gint response, GeometryCorrectionWindow *gcw) {

	if ( response == GTK_RESPONSE_OK ) {

		const char *str;
		float wavenumber_f, circleradius_f;
		double wavenumber, circleradius;

		str = gtk_entry_get_text(GTK_ENTRY(gcw->wn_entry));
		sscanf(str, "%f", &wavenumber_f); wavenumber = wavenumber_f;

		str = gtk_entry_get_text(GTK_ENTRY(gcw->cr_entry));
		sscanf(str, "%f", &circleradius_f); circleradius = circleradius_f;
		
		main_geometry_correct(gcw->correction_type, wavenumber, circleradius);

	}
	
	geometry_window_open = 0;
	gtk_widget_destroy(geometry_window);
	free(gcw);
	
	return 0;

}

static gint geometry_select_vm(GtkWidget *widget, GeometryCorrectionWindow *gcw) {

	gcw->correction_type = GEOMETRY_CORRECTION_TYPE_VINCENT_MIDGLEY;
	
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_entry), FALSE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_label), FALSE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_nm_label), FALSE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_entry), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_nm_label), TRUE);
	
	return 0;
	
}

static gint geometry_select_gn(GtkWidget *widget, GeometryCorrectionWindow *gcw) {

	gcw->correction_type = GEOMETRY_CORRECTION_TYPE_GJONNES;
	
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_entry), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_nm_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_entry), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_nm_label), TRUE);
	
	return 0;

}

static gint geometry_select_gd(GtkWidget *widget, GeometryCorrectionWindow *gcw) {

	gcw->correction_type = GEOMETRY_CORRECTION_TYPE_GJONNES_DIVERGENCE;
	
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_entry), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->wn_nm_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_entry), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_label), TRUE);
	gtk_widget_set_sensitive(GTK_WIDGET(gcw->cr_nm_label), TRUE);
	
	return 0;
	
}

/* Wavelength of an electron (in m) given accelerating potential (in V) */
static double lambda(double V) {

	double m = 9.110E-31;
	double h = 6.625E-34;
	double e = 1.60E-19;
	double c = 2.998E8;
	
	return h / sqrt(2*m*e*V*(1+((e*V) / (2*m*c*c))));

}

static void geometry_set_defaults(GeometryCorrectionWindow *gcw) {

	double wn;
	char string[32];
	int i;
	double r;
	ReflectionList *reflections;
	
	wn = 1/(lambda(300000)*1e9);	/* Answer in nm^-1 */
	snprintf(string, 31, "%f", wn);
	gtk_entry_set_text(GTK_ENTRY(gcw->wn_entry), string);
	
	reflections = main_reflist();
	r = 0;
	double a = data_a();
	double b = data_b();
	for ( i=1; i<reflections->n_reflections; i++ ) {
	
		double g;
		signed int h = reflections->refs[i].h;
		signed int k = reflections->refs[i].k;
		
		g = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)));
		
		if ( g > r ) r = g;
		
	}
	
	snprintf(string, 31, "%f", r);
	gtk_entry_set_text(GTK_ENTRY(gcw->cr_entry), string);

}

void geometry_dialog_open() {

	GtkWidget *geometry_window;
	GtkWidget *vbox;
	GtkWidget *hbox;
	GtkWidget *table;
	GtkWidget *geometry_type_vm;
	GtkWidget *geometry_type_gn;
	GtkWidget *geometry_type_gd;
	GeometryCorrectionWindow *gcw;
	
	if ( geometry_window_open ) {
		return;
	}
	geometry_window_open = 1;
	
	gcw = malloc(sizeof(GeometryCorrectionWindow));
	
	geometry_window = gtk_dialog_new_with_buttons("Geometrical Correction", GTK_WINDOW(displaywindow_gtkwindow()),
		GTK_DIALOG_DESTROY_WITH_PARENT,	GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE, GTK_STOCK_OK, GTK_RESPONSE_OK, NULL);
	
	vbox = gtk_vbox_new(FALSE, 0);
	hbox = gtk_hbox_new(TRUE, 0);
	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(geometry_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7);
	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5);
	
	table = gtk_table_new(5, 3, FALSE);
	gtk_table_set_row_spacings(GTK_TABLE(table), 5);
	gtk_table_set_col_spacings(GTK_TABLE(table), 5);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), FALSE, FALSE, 0);
	
	geometry_type_vm = gtk_radio_button_new_with_label(NULL, "Vincent-Midgley");
	geometry_type_gn = gtk_radio_button_new_with_label_from_widget(GTK_RADIO_BUTTON(geometry_type_vm),
							"Gjønnes (without divergence)");
	geometry_type_gd = gtk_radio_button_new_with_label_from_widget(GTK_RADIO_BUTTON(geometry_type_vm),
							"Gjønnes (with divergence)");
	gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(geometry_type_gd), TRUE);
	gcw->correction_type = GEOMETRY_CORRECTION_TYPE_GJONNES_DIVERGENCE;
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(geometry_type_vm), 1, 4, 1, 2);
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(geometry_type_gn), 1, 4, 2, 3);
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(geometry_type_gd), 1, 4, 3, 4);
	
	gcw->wn_label = gtk_label_new("Electron Wavenumber:");
	gtk_misc_set_alignment(GTK_MISC(gcw->wn_label), 1, 0.5);
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gcw->wn_label), 1, 2, 4, 5);
	gcw->wn_nm_label = gtk_label_new("nm^-1");
	gtk_label_set_markup(GTK_LABEL(gcw->wn_nm_label), "nm<sup>-1</sup>");
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gcw->wn_nm_label), 3, 4, 4, 5);
	
	gcw->wn_entry = gtk_entry_new();
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gcw->wn_entry), 2, 3, 4, 5);
	
	gcw->cr_label = gtk_label_new("Laue Circle Radius:");
	gtk_misc_set_alignment(GTK_MISC(gcw->cr_label), 1, 0.5);
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gcw->cr_label), 1, 2, 5, 6);
	gcw->cr_nm_label = gtk_label_new("nm^-1");
	gtk_label_set_markup(GTK_LABEL(gcw->cr_nm_label), "nm<sup>-1</sup>");
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gcw->cr_nm_label), 3, 4, 5, 6);
	
	gcw->cr_entry = gtk_entry_new();
	gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gcw->cr_entry), 2, 3, 5, 6);
	
	geometry_set_defaults(gcw);
		
	g_signal_connect(G_OBJECT(geometry_type_vm), "clicked", G_CALLBACK(geometry_select_vm), gcw);
	g_signal_connect(G_OBJECT(geometry_type_gn), "clicked", G_CALLBACK(geometry_select_gn), gcw);
	g_signal_connect(G_OBJECT(geometry_type_gd), "clicked", G_CALLBACK(geometry_select_gd), gcw);

	g_signal_connect(G_OBJECT(geometry_window), "response", G_CALLBACK(geometry_window_response), gcw);
	
	gtk_widget_show_all(geometry_window);
	
	gtk_widget_grab_focus(GTK_WIDGET(gcw->wn_entry));

}