/*
 * luzzatti.c
 *
 * Plot of R-factor against resolution
 *
 * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2d - Two-Dimensional Crystallographic Fourier Synthesis
 *
 */

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

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

#include "data.h"
#include "main.h"
#include "statistics.h"
#include "displaywindow.h"
#include "model.h"

#define LUZZATTI_BINS 50

static void luzzatti_gpwrite(FILE *gnuplot, const char *string) {
	fwrite(string, strlen(string), 1, gnuplot);
}

void luzzatti_show() {

	FILE *fh;
	float a = data_a();
	float b = data_b();
	float c = data_c();
	unsigned int i;
	FILE *gnuplot;
	double scale;
	ReflectionList *reflections;
	ReflectionList *model_reflections;
	
	reflections = main_reflist();
	model_reflections = model_calculate_f(reflections, NULL, 69);
	
	scale = stat_scale(reflections, model_reflections);
	fh = fopen("synth2d-luzzatti.dat", "w");
		
	for ( i=1; i<reflections->n_reflections; i++ ) {

		double s, residual;
		signed int h, k, l;
		
		h = reflections->refs[i].h;
		k = reflections->refs[i].k;
		l = reflections->refs[i].l;
		
		s = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)) + ((l*l)/(c*c)));
		residual = 100*fabs((reflections->refs[i].amplitude - scale*model_reflections->refs[i].amplitude)/reflections->refs[i].amplitude);
		printf("%3i %3i %3i: obs=%f calc=%f => residual=%f\n", h, k, l, reflections->refs[i].amplitude,
										scale*model_reflections->refs[i].amplitude, residual);
		fprintf(fh, "%f %f\n", s, residual);
	}
	
	fclose(fh);
	
	gnuplot = popen("gnuplot -persist -", "w");
	if ( !gnuplot ) {
		error_report("Couldn't invoke gnuplot.  Please check your PATH.");
		return;
	}
	
	luzzatti_gpwrite(gnuplot, "set autoscale\n");
	luzzatti_gpwrite(gnuplot, "unset log\n");
	luzzatti_gpwrite(gnuplot, "unset label\n");
	luzzatti_gpwrite(gnuplot, "set xtic auto\n");
	luzzatti_gpwrite(gnuplot, "set ytic auto\n");
	luzzatti_gpwrite(gnuplot, "set grid\n");
	luzzatti_gpwrite(gnuplot, "set ylabel 'R (%)' font \"Helvetica,10\"\n");
	luzzatti_gpwrite(gnuplot, "set xlabel '1/d / nm^(-1)' font \"Helvetica,10\"\n");
	luzzatti_gpwrite(gnuplot, "set title 'Luzzatti Plot' font \"Helvetica,10\"\n");
	luzzatti_gpwrite(gnuplot, "plot 'synth2d-luzzatti.dat'\n");
	
	if ( pclose(gnuplot) == -1 ) {
		error_report("gnuplot returned an error code.");
		return;
	}
	
}