From 0dcc9e2e1ea4fcbd37db5ab7ac74146098c2f4d7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 17 Aug 2010 18:14:44 +0200 Subject: Generation of templates (needs debugging) --- src/facetron.c | 135 +-------------------------------------------------------- 1 file changed, 1 insertion(+), 134 deletions(-) (limited to 'src/facetron.c') diff --git a/src/facetron.c b/src/facetron.c index 0fa6b643..0c2b4d56 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -26,10 +26,7 @@ #include "cell.h" #include "hdf5-file.h" #include "detector.h" -#include "peaks.h" - - -#define MAX_HITS (1024) +#include "geometry.h" static void show_help(const char *s) @@ -116,136 +113,6 @@ static int find_chunk(FILE *fh, UnitCell **cell, char **filename) } -static 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; hlambda - 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; pdet->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; -} - - -static double integrate_all(struct image *image, struct reflhit *hits, int n) -{ - double itot = 0.0; - int i; - - for ( i=0; i