/*
 * reproject.c
 *
 * Synthesize diffraction patterns
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

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

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

#define MAX_IMAGE_REFLECTIONS 8*1024

ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx, ControlContext *ctx) {

	ImageReflection *refl;
	Reflection *reflection;
	int i;
	double smax = 0.5e9;
	double tilt, omega;
	
	tilt = deg2rad(image.tilt);
	omega = deg2rad(image.omega);
	
	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
	
	i = 0;
	
	reflection = rctx->reflections;
	do {
	
		double xl, yl, zl;
		double nx, ny, nz;
		double xt, yt, zt;
		double a, b, c;
		double s1, s2, s;
		
		/* Get the coordinates of the reciprocal lattice point */
		xl = reflection->x;
		yl = reflection->y;
		zl = reflection->z;
		
		/* Now calculate the (normalised) incident electron wavevector */
		xt = 0;
		yt = - sin(tilt);
		zt = - cos(tilt);
		nx = xt*cos(omega) + yt*-sin(omega);
		ny = xt*sin(omega) + yt*cos(omega);
		nz = zt;
		
		/* Next, solve the relrod equation to calculate the excitation error */
		a = 1.0;
		b = 2.0*(xl*nx + yl*ny + zl*nz - nz/image.lambda);
		c = xl*xl + yl*yl + zl*zl - 2.0*zl/image.lambda;
		s1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
		s2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
		
		/* Skip this reflection if s is large */
		if ( fabs(s) <= smax ) {
		
			double xi, yi, zi;
			double gx, gy, gz;
			double kx, ky, kz;
			double cx, cy, cz;
			double ux, uy, uz;
			double uxt, uyt, uzt;
			double theta, psi;
			double x, y;
			
			/* Determine the intersection point */
			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
			reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_CENTRAL);
			reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER);
			
			/* Calculate Bragg angle and azimuth of point in image */
			kx = nx / image.lambda;
			ky = ny / image.lambda;
			kz = nz / image.lambda;	/* This is the centre of the Ewald sphere */
			gx = xi - kx;
			gy = yi - ky;
			gz = zi - kz;	/* This is the vector from the centre of the sphere to the intersection */
			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
			cx = nz*gy - ny*gz;
			cy = nz*gx - nx*gz;
			cz = ny*gx - nx*gy;
			uxt = 0;
			uyt = cos(tilt);
			uzt = -sin(tilt);
			ux = uyt*-sin(omega);
			uy = uyt*cos(omega);
			uz = uzt;
			psi = angle_between(cx, cy, cz, ux, uy, uz) - M_PI_2;
			
			if ( image.fmode == FORMULATION_CLEN ) {
				x = image.camera_len*sin(theta)*cos(psi);
				y = image.camera_len*sin(theta)*sin(psi);
				x *= image.resolution;
				y *= image.resolution;
			} else if ( image.fmode == FORMULATION_PIXELSIZE ) {
				x = sin(theta)*cos(psi) / image.lambda;
				y = sin(theta)*sin(psi) / image.lambda;
				x /= image.pixel_size;
				y /= image.pixel_size;
			} else {
				fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n");
				return NULL;
			}
			
			/* Adjust centre */
			x += image.x_centre;
			y += image.y_centre;
			
			/* Sanity check */
			if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) {
			
				/* Record the reflection */
				refl[i].x = x;
				refl[i].y = y;
				i++;
				
				if ( i > MAX_IMAGE_REFLECTIONS ) {
					fprintf(stderr, "Too many reflections\n");
					break;
				}
				
				//printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
			
			} else {
				//fprintf(stderr, "Reflection failed sanity test (x=%f, y=%f)\n", x, y);
			}
		
		}
		
		reflection = reflection->next;

	} while ( reflection );
	
	//printf("Found %i reflections in image\n", i);
	*n = i;
	return refl;

}