Move a load more stuff into libcrystfel

1 files changed, 0 insertions, 341 deletions

diff --git a/src/geometry.c b/src/geometry.c
deleted file mode 100644
index 485abba3..00000000
--- a/src/geometry.c
+++ /dev/null
@@ -1,341 +0,0 @@
-/*
- * 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);
-}