/*
 * 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;
	
}

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;


}

/* 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;

}