/*
 * reflections.c
 *
 * Data structures in 3D space
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

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

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

#include "reflections.h"
#include "utils.h"

static void reflectionlist_init(ReflectionList *reflectionlist) {
	reflectionlist->n_reflections = 0;
	reflectionlist->list_capped = 0;
	reflectionlist->reflections = NULL;
	reflectionlist->last_reflection = NULL;
}

ReflectionList *reflectionlist_new() {

	ReflectionList *reflectionlist = malloc(sizeof(ReflectionList));
	
	reflectionlist_init(reflectionlist);
	
	return reflectionlist;

}

void reflectionlist_clear_markers(ReflectionList *reflectionlist) {

	Reflection *reflection = reflectionlist->reflections;
	Reflection *prev = NULL;
	int del = 0;
	
	do {
		Reflection *next = reflection->next;
		
		if ( (reflection->type == REFLECTION_MARKER) || (reflection->type == REFLECTION_GENERATED)
		 || (reflection->type == REFLECTION_VECTOR_MARKER_1) || (reflection->type == REFLECTION_VECTOR_MARKER_2)
		 || (reflection->type == REFLECTION_VECTOR_MARKER_3) ) {
			free(reflection);
			del++;
			if ( prev ) {
				prev->next = next;
			} else {
				reflectionlist->reflections = next;
			}
		} else {
			prev = reflection;
		}
		
		reflection = next;
	
	} while ( reflection );
	
	reflectionlist->n_reflections -= del;
	reflectionlist->last_reflection = prev;
	
}
void reflectionlist_clear(ReflectionList *reflectionlist) {

	Reflection *reflection = reflectionlist->reflections;
	do {
		Reflection *next = reflection->next;
		free(reflection);
		reflection = next;
	} while ( reflection );
	
	reflectionlist_init(reflectionlist);

}

void reflectionlist_free(ReflectionList *reflectionlist) {
	reflectionlist_clear(reflectionlist);
	free(reflectionlist);
}

Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, double z, double intensity, ReflectionType type) {

	Reflection *new_reflection;
	Reflection *nearest;
	
	if ( reflectionlist->list_capped ) return NULL;
	
	if ( reflectionlist->n_reflections > 1e7 ) {
		fprintf(stderr, "More than 10 million reflections on list.  I think this is silly.\n");
		fprintf(stderr, "No further reflections will be stored.  Go and fix the peak detection.\n");
		reflectionlist->list_capped = 1;
	}
	
	nearest = reflectionlist_find_nearest_type(reflectionlist, x, y, z, type);
	if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL;
	
	new_reflection = malloc(sizeof(Reflection));
	new_reflection->next = NULL;
	new_reflection->x = x;
	new_reflection->y = y;
	new_reflection->z = z;
	new_reflection->intensity = intensity;
	new_reflection->type = type;
	new_reflection->found = 0;
	
	if ( reflectionlist->last_reflection ) {
		reflectionlist->last_reflection->next = new_reflection;
		reflectionlist->last_reflection = new_reflection;
	} else {
		reflectionlist->reflections = new_reflection;
		reflectionlist->last_reflection = new_reflection;
	}
	reflectionlist->n_reflections++;
	
	return new_reflection;
	
}

int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {

	/* "Input" space */
	double d, x, y;
	ImageRecord *imagerecord;
	
	/* Angular description of reflection */
	double theta, psi, k;
	
	/* Reciprocal space */
	double tilt;
	double omega;
	
	double x_temp, y_temp, z_temp;
	double nx, ny, nz;
	
	imagerecord = refl->parent;
	x = refl->x;	y = refl->y;
	tilt = 2*M_PI*(imagerecord->tilt/360);	/* Convert to Radians */
	omega = 2*M_PI*(imagerecord->omega/360);	/* Likewise */
	k = 1/imagerecord->lambda;
	
	/* Calculate an angular description of the reflection */
	if ( imagerecord->fmode == FORMULATION_CLEN ) {
		x /= imagerecord->resolution;
		y /= imagerecord->resolution;	/* Convert pixels to metres */
		d = sqrt((x*x) + (y*y));
		theta = atan2(d, imagerecord->camera_len);
	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
		x *= imagerecord->pixel_size;
		y *= imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
		d = sqrt((x*x) + (y*y));
		theta = atan2(d, k);
	} else {
		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
		return -1;
	}
	psi = atan2(y, x);
	
	x_temp = k*sin(theta)*cos(psi);
	y_temp = k*sin(theta)*sin(psi);
	z_temp = k - k*cos(theta);
	
	/* Apply the rotations...
		First: rotate image clockwise until tilt axis is aligned horizontally. */
	nx = x_temp*cos(omega) + y_temp*sin(omega);
	ny = -x_temp*sin(omega) + y_temp*cos(omega);
	nz = z_temp;
	
	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
	x_temp = nx; y_temp = ny; z_temp = nz;	
	nx = x_temp;
	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
	
	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
	x_temp = nx; y_temp = ny; z_temp = nz;	
	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
	nz = z_temp;
	
	*ddx = nx;
	*ddy = ny;
	*ddz = nz;
	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
	
	return 0;

}

/* x and y in pixels, measured from centre of image */
void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord *imagerecord, double intensity) {

	double nx, ny, nz, twotheta;
	ImageReflection refl;
	
	refl.parent = imagerecord;
	refl.x = x;
	refl.y = y;
	refl.intensity = intensity;
	
	if ( !reflection_map_to_space(&refl, &nx, &ny, &nz, &twotheta) ) {
		reflection_add(ctx->reflectionlist, nx, ny, nz, intensity, REFLECTION_NORMAL);
	}

}

void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r) {
	if ( reflectionlist->last_reflection ) {
		reflectionlist->last_reflection->next = r;
		reflectionlist->last_reflection = r;
	} else {
		reflectionlist->reflections = r;
		reflectionlist->last_reflection = r;
	}
	r->next = NULL;
	reflectionlist->n_reflections++;
}

double reflectionlist_largest_g(ReflectionList *reflectionlist) {

	double max = 0.0;
	Reflection *reflection;
	
	reflection = reflectionlist->reflections;
	while ( reflection ) {
		if ( reflection->type == REFLECTION_NORMAL ) {
			double mod;
			mod = modulus(reflection->x, reflection->y, reflection->z);
			if ( mod > max ) max = mod;
		}
		reflection = reflection->next;
	};
	
	return max;

}

Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x, double y, double z) {

	double max = +INFINITY;
	Reflection *reflection;
	Reflection *best = NULL;
	
	reflection = reflectionlist->reflections;
	while ( reflection ) {
		if ( reflection->type == REFLECTION_NORMAL ) {
			double mod;
			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
			if ( mod < max ) {
				max = mod;
				best = reflection;
			}
		}
		reflection = reflection->next;
	};
	
	return best;

}

Reflection *reflectionlist_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) {

	double max = +INFINITY;
	Reflection *reflection;
	Reflection *best = NULL;
	
	reflection = reflectionlist->reflections;
	while ( reflection ) {
		if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) {
			double mod;
			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
			if ( (mod < max) && (mod >= min_distance) ) {
				max = mod;
				best = reflection;
			}
		}
		reflection = reflection->next;
	};
	
	return best;

}

Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type) {

	double max = +INFINITY;
	Reflection *reflection;
	Reflection *best = NULL;
	
	reflection = reflectionlist->reflections;
	while ( reflection ) {
		if ( reflection->type == type ) {
			double mod;
			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
			if ( mod < max ) {
				max = mod;
				best = reflection;
			}
		}
		reflection = reflection->next;
	};
	
	return best;


}

/* This destroys the lfom values in the input list */
ReflectionList *reflectionlist_sort_lfom(ReflectionList *in) {

	ReflectionList *out;
	Reflection *best;
	
	out = reflectionlist_new();
	
	do {
	
		Reflection *reflection;
		int lfom = 0;
		
		reflection = in->reflections;
		best = NULL;
		while ( reflection ) {
			if ( reflection->lfom > lfom ) {
				best = reflection;
				lfom = reflection->lfom;
			}
			reflection = reflection->next;
		};
		
		if ( best ) {
			Reflection *new;
			new = reflection_add(out, best->x, best->y, best->z, best->intensity, best->type);
			new->lfom = best->lfom;
			best->lfom = 0;
		}
		
	} while ( best );
	
	return out;

}

/* Generate a list of reflections from a unit cell */
ReflectionList *reflection_list_from_cell(Basis *basis) {

	ReflectionList *ordered;
	double max_res;
	signed int h, k, l;
	int max_order_a, max_order_b, max_order_c;
	
	
	max_res = 20e9;
	
	ordered = reflectionlist_new();
	
	max_order_a = max_res/modulus(basis->a.x, basis->a.y, basis->a.z);
	max_order_b = max_res/modulus(basis->b.x, basis->b.y, basis->b.z);
	max_order_c = max_res/modulus(basis->c.x, basis->c.y, basis->c.z);
	for ( h=-max_order_a; h<=max_order_a; h++ ) {
		for ( k=-max_order_b; k<=max_order_b; k++ ) {
			for ( l=-max_order_c; l<=max_order_c; l++ ) {

				double x, y, z;
				
				x = h*basis->a.x + k*basis->b.x + l*basis->c.x;
				y = h*basis->a.y + k*basis->b.y + l*basis->c.y;
				z = h*basis->a.z + k*basis->b.z + l*basis->c.z;
				
				if ( ( x*x + y*y + z*z ) <= max_res*max_res ) {
					Reflection *ref;
					ref = reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
					if ( ref ) {	/* Check it's not NULL */
						ref->h = h;	ref->k = k;	ref->l = l;
					}
					if ( (h!=0) || (k!=0) || (l!=0) ) {
					//	reflection_add(ctx->reflectionlist, x, y, z, 1, REFLECTION_GENERATED);
						reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
					}
				}
				
			}
		}
	}
	
	return ordered;

}