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

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

#include "reflections.h"

static void reflection_addfirst(ReflectionContext *reflectionctx) {
	/* Create first items on lists - saves faffing later.  Corresponds to a central marker.
		Reflections are only stored if they have non-zero value. */
	reflectionctx->reflections = malloc(sizeof(Reflection));
	reflectionctx->reflections->next = NULL;
	reflectionctx->reflections->x = 0;
	reflectionctx->reflections->y = 0;
	reflectionctx->reflections->z = 0;
	reflectionctx->reflections->type = REFLECTION_CENTRAL;
	reflectionctx->last_reflection = reflectionctx->reflections;
	reflectionctx->n_reflections = 1;
	reflectionctx->list_capped = 0;
}

ReflectionContext *reflection_init() {

	ReflectionContext *reflectionctx = malloc(sizeof(ReflectionContext));
	reflection_addfirst(reflectionctx);
	reflectionctx->n_reflections = 0;
	reflectionctx->list_capped = 0;
	
	return reflectionctx;

}

void reflection_clear_markers(ReflectionContext *reflectionctx) {

	Reflection *reflection = reflectionctx->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 {
				reflectionctx->reflections = next;
			}
		} else {
			prev = reflection;
		}
		
		reflection = next;
	
	} while ( reflection );
	
	reflectionctx->n_reflections -= del;
	reflectionctx->last_reflection = prev;
	
}
void reflection_clear(ReflectionContext *reflectionctx) {

	Reflection *reflection = reflectionctx->reflections;
	do {
		Reflection *next = reflection->next;
		free(reflection);
		reflection = next;
	} while ( reflection );
	
	reflectionctx->n_reflections = 0;
	reflectionctx->list_capped = 0;
	reflection_addfirst(reflectionctx);

}

void reflection_free(ReflectionContext *reflectionctx) {
	reflection_clear(reflectionctx);
	free(reflectionctx->reflections);
	free(reflectionctx);
}

Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) {

	Reflection *new_reflection;
	
	if ( reflectionctx->list_capped ) return NULL;
	
	if ( reflectionctx->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");
		reflectionctx->list_capped = 1;
	}
	reflectionctx->n_reflections++;
	
	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;
	
	reflectionctx->last_reflection->next = new_reflection;
	reflectionctx->last_reflection = new_reflection;
	
	return new_reflection;
	
}

/* 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) {

	/* "Input" space */
	double d;
	
	/* 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;
	
	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 ( ctx->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 ( ctx->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;
	}
	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;
	
	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);

}

void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r) {
	r->next = NULL;
	rctx->last_reflection->next = r;
	rctx->last_reflection = r;
}