/*
 * 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(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 metres (in real space!) */
void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) {

	double nx, ny, nz;

	/* Real space */
	double L;
	double d;

	/* Shared */
	double theta;
	double psi;
	
	/* Reciprocal space */
	double r;
	double tilt;
	double omega;
	double x_temp, y_temp, z_temp;
	
	tilt = 2*M_PI*(tilt_degrees/360);	/* Convert to Radians */
	omega = 2*M_PI*(ctx->omega/360);	/* Likewise */
	
	d = sqrt((x*x) + (y*y));
	L = ctx->camera_length;
	theta = atan2(d, L);
	psi = atan2(y, x);
	
	r = 1/ctx->lambda;
	x_temp = r*sin(theta)*cos(psi);
	y_temp = -r*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
	z_temp = r- r*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;
	
	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);

}

/* Alternative formulation for when the input data is already in reciprocal space units 
	x and y in metres^-1 (in reciprocal space) */
void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) {

	double nx, ny, nz;
	
	/* Input (reciprocal space) */
	double dr;
	
	/* Shared */
	double theta;
	double psi;
	double r;
	
	/* Reciprocal space */
	double tilt;
	double omega;
	double x_temp, y_temp, z_temp;
	
	tilt = M_PI*(tilt_degrees/180);	/* Convert to Radians */
	omega = M_PI*(ctx->omega/180);	/* Likewise */
	
	dr = sqrt((x*x) + (y*y));
	r = 1/ctx->lambda;
	theta = atan2(dr, r);
	psi = atan2(y, x);
	
	x_temp = r*sin(theta)*cos(psi);
	y_temp = -r*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
	z_temp = r - r*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;
	
	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;
}