diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/geometry.c | 155 | ||||
| -rw-r--r-- | src/geometry.h | 27 |
2 files changed, 182 insertions, 0 deletions
diff --git a/src/geometry.c b/src/geometry.c new file mode 100644 index 00000000..16914994 --- /dev/null +++ b/src/geometry.c @@ -0,0 +1,155 @@ +/* + * geometry.h + * + * Geometry of diffraction + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdlib.h> + +#include "utils.h" +#include "cell.h" +#include "image.h" +#include "peaks.h" + + +#define MAX_HITS (1024) + + +struct reflhit *find_intersections(struct image *image, UnitCell *cell, + double divergence, double bandwidth, + int *n, int output) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + struct reflhit *hits; + int np = 0; + int hmax, kmax, lmax; + double mres; + signed int h, k, l; + + hits = malloc(sizeof(struct reflhit)*MAX_HITS); + if ( hits == NULL ) { + *n = 0; + return NULL; + } + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + mres = 1.0 / 8.0e-10; /* 8 Angstroms */ + hmax = mres / modulus(asx, asy, asz); + kmax = mres / modulus(bsx, bsy, bsz); + lmax = mres / modulus(csx, csy, csz); + + for ( h=-hmax; h<hmax; h++ ) { + for ( k=-kmax; k<kmax; k++ ) { + for ( l=-lmax; l<lmax; l++ ) { + + double xl, yl, zl; + double ds_sq, dps_sq; + double delta, divfact; + double llow, lhigh; + double xd, yd, cl; + double xda, yda; + int p; + int found = 0; + + if ( (h==0) && (k==0) && (l==0) ) continue; + + llow = image->lambda - image->lambda*bandwidth/2.0; + lhigh = image->lambda + image->lambda*bandwidth/2.0; + + /* Get the coordinates of the reciprocal lattice point */ + zl = h*asz + k*bsz + l*csz; + if ( zl < 0.0 ) continue; /* Do this check very early */ + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + + ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */ + delta = divergence/image->lambda; + dps_sq = ds_sq + pow(delta, 2.0); /* d'*^2 */ + + /* In range? */ + divfact = 2.0 * delta * sqrt(xl*xl + yl*yl); + if ( ds_sq - 2.0*zl/llow > 0.0 ) continue; + if ( ds_sq - 2.0*zl/lhigh < 0.0 ) continue; + + /* Work out which panel this peak would fall on */ + for ( p=0; p<image->det->n_panels; p++ ) { + + /* Camera length for this panel */ + cl = image->det->panels[p].clen; + + /* Coordinates of peak relative to central beam, in m */ + xd = cl*xl / (ds_sq/(2.0*zl) - zl); + yd = cl*yl / (ds_sq/(2.0*zl) - zl); + + /* Convert to pixels */ + xd *= image->det->panels[p].res; + yd *= image->det->panels[p].res; + + /* Add the coordinates of the central beam */ + xda = xd + image->det->panels[p].cx; + yda = yd + image->det->panels[p].cy; + + /* Now, is this on this panel? */ + if ( xda < image->det->panels[p].min_x ) continue; + if ( xda > image->det->panels[p].max_x ) continue; + if ( yda < image->det->panels[p].min_y ) continue; + if ( yda > image->det->panels[p].max_y ) continue; + + /* Woohoo! */ + found = 1; + break; + + } + + if ( !found ) continue; + + hits[np].h = h; + hits[np].k = k; + hits[np].l = l; + np++; + + if ( output ) { + printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda); + } + + } + } + } + + *n = np; + return hits; +} + + +double integrate_all(struct image *image, struct reflhit *hits, int n) +{ + double itot = 0.0; + int i; + + for ( i=0; i<n; i++ ) { + + float x, y, intensity; + + if ( integrate_peak(image, hits[i].x, hits[i].y, &x, &y, + &intensity, 0, 0) ) continue; + + itot += intensity; + } + + return itot; +} diff --git a/src/geometry.h b/src/geometry.h new file mode 100644 index 00000000..a1e73127 --- /dev/null +++ b/src/geometry.h @@ -0,0 +1,27 @@ +/* + * geometry.h + * + * Geometry of diffraction + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef GEOMETRY_H +#define GEOMETRY_H + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +extern struct reflhit *find_intersections(struct image *image, UnitCell *cell, + double divergence, double bandwidth, + int *n, int output); + +extern double integrate_all(struct image *image, struct reflhit *hits, int n); + + +#endif /* GEOMETRY_H */ |