From 469efb904b59f137ac9e85e5ff23edd0c113de5c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 12:17:59 +0100 Subject: Move a load more stuff into libcrystfel --- src/geometry.c | 341 --------------------------------------------------------- 1 file changed, 341 deletions(-) delete mode 100644 src/geometry.c (limited to 'src/geometry.c') 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 - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include - -#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; in_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); -} -- cgit v1.2.3