Move a load more stuff into libcrystfel

1 files changed, 341 insertions, 0 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
new file mode 100644
index 00000000..485abba3
--- /dev/null
+++ b/libcrystfel/src/geometry.c
@@ -0,0 +1,341 @@
+/*
+ * geometry.c
+ *
+ * Geometry of diffraction
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <assert.h>
+
+#include "utils.h"
+#include "cell.h"
+#include "image.h"
+#include "peaks.h"
+#include "beam-parameters.h"
+#include "reflist.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+
+
+static signed int locate_peak(double x, double y, double z, double k,
+                              struct detector *det, double *xdap, double *ydap)
+{
+	int i;
+	signed int found = -1;
+	const double den = k + z;
+
+	*xdap = -1;  *ydap = -1;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		double xd, yd;
+		double fs, ss, plx, ply;
+		struct panel *p;
+
+		p = &det->panels[i];
+
+		/* Coordinates of peak relative to central beam, in m */
+		xd = p->clen * x / den;
+		yd = p->clen * y / den;
+
+		/* Convert to pixels */
+		xd *= p->res;
+		yd *= p->res;
+
+		/* Convert to relative to the panel corner */
+		plx = xd - p->cnx;
+		ply = yd - p->cny;
+
+		fs = p->xfs*plx + p->yfs*ply;
+		ss = p->xss*plx + p->yss*ply;
+
+		fs += p->min_fs;
+		ss += p->min_ss;
+
+		/* Now, is this on this panel? */
+		if ( fs < p->min_fs ) continue;
+		if ( fs > p->max_fs ) continue;
+		if ( ss < p->min_ss ) continue;
+		if ( ss > p->max_ss ) continue;
+
+		/* If peak appears on multiple panels, reject it */
+		if ( found != -1 ) return -1;
+
+		/* Woohoo! */
+		found = i;
+		*xdap = fs;
+		*ydap = ss;
+
+	}
+
+	return found;
+}
+
+
+static double excitation_error(double xl, double yl, double zl,
+                               double ds, double k, double divergence,
+                               double tt)
+{
+	double al;
+	double r;
+	double delta;
+
+	al = M_PI_2 - asin(-zl/ds);
+
+	r = ( ds * sin(al) / sin(tt) ) - k;
+
+	delta = sqrt(2.0 * pow(ds, 2.0) * (1.0-cos(divergence)));
+	if ( divergence > 0.0 ) {
+		r += delta;
+	} else {
+		r -= delta;
+	}
+
+	return r;
+}
+
+
+static double partiality(double r1, double r2, double r)
+{
+	double q1, q2;
+	double p1, p2;
+
+	/* Calculate degrees of penetration */
+	q1 = (r1 + r)/(2.0*r);
+	q2 = (r2 + r)/(2.0*r);
+
+	/* Convert to partiality */
+	p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
+	p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
+
+	return p2 - p1;
+}
+
+
+static Reflection *check_reflection(struct image *image,
+                                    signed int h, signed int k, signed int l,
+                                    double asx, double asy, double asz,
+                                    double bsx, double bsy, double bsz,
+                                    double csx, double csy, double csz)
+{
+	const int output = 0;
+	double xl, yl, zl;
+	double ds, ds_sq;
+	double rlow, rhigh;     /* "Excitation error" */
+	signed int p;           /* Panel number */
+	double xda, yda;        /* Position on detector */
+	int close, inside;
+	double part;            /* Partiality */
+	int clamp_low = 0;
+	int clamp_high = 0;
+	double bandwidth = image->bw;
+	double divergence = image->div;
+	double lambda = image->lambda;
+	double klow, kcen, khigh;    /* Wavenumber */
+	Reflection *refl;
+	double tt;
+
+	/* "low" gives the largest Ewald sphere,
+	 * "high" gives the smallest Ewald sphere. */
+	klow = 1.0/(lambda - lambda*bandwidth/2.0);
+	kcen = 1.0/lambda;
+	khigh = 1.0/(lambda + lambda*bandwidth/2.0);
+
+	/* Get the coordinates of the reciprocal lattice point */
+	zl = h*asz + k*bsz + l*csz;
+	/* Throw out if it's "in front".  A tiny bit "in front" is OK. */
+	if ( zl > image->profile_radius ) return NULL;
+	xl = h*asx + k*bsx + l*csx;
+	yl = h*asy + k*bsy + l*csy;
+
+	tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+kcen);
+	if ( tt > deg2rad(90.0) ) return NULL;
+
+	ds_sq = modulus_squared(xl, yl, zl);  /* d*^2 */
+	ds = sqrt(ds_sq);
+
+	/* Calculate excitation errors */
+	rlow = excitation_error(xl, yl, zl, ds, klow, -divergence/2.0, tt);
+	rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence/2.0, tt);
+
+	/* Is the reciprocal lattice point close to either extreme of
+	 * the sphere, maybe just outside the "Ewald volume"? */
+	close = (fabs(rlow) < image->profile_radius)
+	     || (fabs(rhigh) < image->profile_radius);
+
+	/* Is the reciprocal lattice point somewhere between the
+	 * extremes of the sphere, i.e. inside the "Ewald volume"? */
+	inside = signbit(rlow) ^ signbit(rhigh);
+
+	/* Can't be both inside and close */
+	if ( inside ) close = 0;
+
+	/* Neither?  Skip it. */
+	if ( !(close || inside) ) return NULL;
+
+	/* If the "lower" Ewald sphere is a long way away, use the
+	 * position at which the Ewald sphere would just touch the
+	 * reflection. */
+	if ( rlow < -image->profile_radius ) {
+		rlow = -image->profile_radius;
+		clamp_low = -1;
+	}
+	if ( rlow > +image->profile_radius ) {
+		rlow = +image->profile_radius;
+		clamp_low = +1;
+	}
+	/* Likewise the "higher" Ewald sphere */
+	if ( rhigh < -image->profile_radius ) {
+		rhigh = -image->profile_radius;
+		clamp_high = -1;
+	}
+	if ( rhigh > +image->profile_radius ) {
+		rhigh = +image->profile_radius;
+		clamp_high = +1;
+	}
+	assert(clamp_low <= clamp_high);
+	/* The six possible combinations of clamp_{low,high} (including
+	 * zero) correspond to the six situations in Table 3 of Rossmann
+	 * et al. (1979). */
+
+	/* Calculate partiality */
+	part = partiality(rlow, rhigh, image->profile_radius);
+
+	/* Locate peak on detector. */
+	p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
+	if ( p == -1 ) return NULL;
+
+	/* Add peak to list */
+	refl = reflection_new(h, k, l);
+	set_detector_pos(refl, 0.0, xda, yda);
+	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	set_symmetric_indices(refl, h, k, l);
+	set_redundancy(refl, 1);
+
+	if ( output ) {
+		printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
+		       h, k, l, 0.0, xda, yda, part);
+	}
+
+	return refl;
+}
+
+
+RefList *find_intersections(struct image *image, UnitCell *cell)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	RefList *reflections;
+	int hmax, kmax, lmax;
+	double mres;
+	signed int h, k, l;
+
+	reflections = reflist_new();
+
+	/* Cell angle check from Foadi and Evans (2011) */
+	if ( !cell_is_sensible(cell) ) {
+		ERROR("Invalid unit cell parameters given to"
+		      " find_intersections()\n");
+		cell_print(cell);
+		return NULL;
+	}
+
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+
+	/* We add a horrific 20% fudge factor because bandwidth, divergence
+	 * and so on mean reflections appear beyond the largest q */
+	mres = 1.2 * largest_q(image);
+
+	hmax = mres / modulus(asx, asy, asz);
+	kmax = mres / modulus(bsx, bsy, bsz);
+	lmax = mres / modulus(csx, csy, csz);
+
+	if ( (hmax >= 256) || (kmax >= 256) || (lmax >= 256) ) {
+		ERROR("Unit cell is stupidly large.\n");
+		cell_print(cell);
+		if ( hmax >= 256 ) hmax = 255;
+		if ( kmax >= 256 ) kmax = 255;
+		if ( lmax >= 256 ) lmax = 255;
+	}
+
+	for ( h=-hmax; h<=hmax; h++ ) {
+	for ( k=-kmax; k<=kmax; k++ ) {
+	for ( l=-lmax; l<=lmax; l++ ) {
+
+		Reflection *refl;
+
+		refl = check_reflection(image, h, k, l,
+		                        asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+
+		if ( refl != NULL ) {
+			refl = add_refl_to_list(refl, reflections);
+		}
+
+	}
+	}
+	}
+
+	return reflections;
+}
+
+
+/* Calculate partialities and apply them to the image's reflections */
+void update_partialities(struct image *image)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *predicted;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+
+	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
+
+	/* Scratch list to give check_reflection() something to add to */
+	predicted = reflist_new();
+
+	for ( refl = first_refl(image->reflections, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		Reflection *vals;
+		double r1, r2, p, x, y;
+		signed int h, k, l;
+		int clamp1, clamp2;
+
+		get_symmetric_indices(refl, &h, &k, &l);
+
+		vals = check_reflection(image, h, k, l,
+		                        asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+
+		if ( vals == NULL ) {
+			set_redundancy(refl, 0);
+			continue;
+		}
+		set_redundancy(refl, 1);
+
+		/* Transfer partiality stuff */
+		get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2);
+		set_partial(refl, r1, r2, p, clamp1, clamp2);
+
+		/* Transfer detector location */
+		get_detector_pos(vals, &x, &y);
+		set_detector_pos(refl, 0.0, x, y);
+	}
+
+	reflist_free(predicted);
+}