aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/beam-parameters.c133
-rw-r--r--src/beam-parameters.h40
-rw-r--r--src/cell.c1164
-rw-r--r--src/cell.h107
-rw-r--r--src/detector.c1234
-rw-r--r--src/detector.h131
-rw-r--r--src/hdf5-file.c817
-rw-r--r--src/hdf5-file.h55
-rw-r--r--src/image.h8
-rw-r--r--src/list_tmp.h106
-rw-r--r--src/partial_sim.c2
-rw-r--r--src/reflist.c1048
-rw-r--r--src/reflist.h112
-rw-r--r--src/scaling-report.c2
-rw-r--r--src/thread-pool.c283
-rw-r--r--src/thread-pool.h71
-rw-r--r--src/utils.c488
-rw-r--r--src/utils.h236
18 files changed, 6 insertions, 6031 deletions
diff --git a/src/beam-parameters.c b/src/beam-parameters.c
deleted file mode 100644
index ac30ad04..00000000
--- a/src/beam-parameters.c
+++ /dev/null
@@ -1,133 +0,0 @@
-/*
- * beam-parameters.c
- *
- * Beam parameters
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#include "beam-parameters.h"
-#include "utils.h"
-
-
-struct beam_params *get_beam_parameters(const char *filename)
-{
- FILE *fh;
- struct beam_params *b;
- char *rval;
- int reject;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) return NULL;
-
- b = calloc(1, sizeof(struct beam_params));
- if ( b == NULL ) return NULL;
-
- b->fluence = -1.0;
- b->beam_radius = -1.0;
- b->photon_energy = -1.0;
- b->bandwidth = -1.0;
- b->divergence = -1.0;
- b->dqe = -1.0;
- b->adu_per_photon = -1.0;
-
- do {
-
- int n1;
- char line[1024];
- int i;
- char **bits;
-
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) break;
- chomp(line);
-
- n1 = assplode(line, " \t", &bits, ASSPLODE_NONE);
- if ( n1 < 3 ) {
- for ( i=0; i<n1; i++ ) free(bits[i]);
- free(bits);
- continue;
- }
-
- if ( bits[1][0] != '=' ) {
- for ( i=0; i<n1; i++ ) free(bits[i]);
- free(bits);
- continue;
- }
-
- if ( strcmp(bits[0], "beam/fluence") == 0 ) {
- b->fluence = atof(bits[2]);
- } else if ( strcmp(bits[0], "beam/radius") == 0 ) {
- b->beam_radius = atof(bits[2]);
- } else if ( strcmp(bits[0], "beam/photon_energy") == 0 ) {
- b->photon_energy = atof(bits[2]);
- } else if ( strcmp(bits[0], "beam/bandwidth") == 0 ) {
- b->bandwidth = atof(bits[2]);
- } else if ( strcmp(bits[0], "beam/divergence") == 0 ) {
- b->divergence = atof(bits[2]);
- } else if ( strcmp(bits[0], "detector/dqe") == 0 ) {
- b->dqe = atof(bits[2]);
- } else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) {
- b->adu_per_photon = atof(bits[2]);
- } else {
- ERROR("Unrecognised field '%s'\n", bits[0]);
- }
-
- for ( i=0; i<n1; i++ ) free(bits[i]);
- free(bits);
-
- } while ( rval != NULL );
- fclose(fh);
-
- reject = 0;
-
- if ( b->fluence < 0.0 ) {
- ERROR("Invalid or unspecified value for 'beam/fluence'.\n");
- reject = 1;
- }
- if ( b->beam_radius < 0.0 ) {
- ERROR("Invalid or unspecified value for 'beam/radius'.\n");
- reject = 1;
- }
- if ( b->photon_energy < 0.0 ) {
- ERROR("Invalid or unspecified value for"
- " 'beam/photon_energy'.\n");
- reject = 1;
- }
- if ( b->bandwidth < 0.0 ) {
- ERROR("Invalid or unspecified value for 'beam/bandwidth'.\n");
- reject = 1;
- }
- if ( b->divergence < 0.0 ) {
- ERROR("Invalid or unspecified value for 'beam/divergence'.\n");
- reject = 1;
- }
- if ( b->dqe < 0.0 ) {
- ERROR("Invalid or unspecified value for 'detector/dqe'.\n");
- reject = 1;
- }
- if ( b->adu_per_photon < 0.0 ) {
- ERROR("Invalid or unspecified value for"
- " 'detector/adu_per_photon'.\n");
- reject = 1;
- }
-
- if ( reject ) {
- ERROR("Please fix the above problems with the beam"
- " parameters file and try again.\n");
- return NULL;
- }
-
- return b;
-}
diff --git a/src/beam-parameters.h b/src/beam-parameters.h
deleted file mode 100644
index 411ab449..00000000
--- a/src/beam-parameters.h
+++ /dev/null
@@ -1,40 +0,0 @@
-/*
- * beam-parameters.h
- *
- * Beam parameters
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifndef BEAM_PARAMETERS_H
-#define BEAM_PARAMETERS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-struct beam_params
-{
- double fluence; /* photons per pulse */
- double beam_radius; /* metres */
- double photon_energy; /* eV per photon */
- double bandwidth; /* FWHM(wavelength) over wavelength.
- * Note: current simulation code just uses
- * a rectangular distribution with this as
- * its (full) width. */
- double divergence; /* divergence (radians) */
-
- double dqe; /* Detector DQE (fraction) */
- double adu_per_photon; /* Detector "gain" */
-};
-
-
-extern struct beam_params *get_beam_parameters(const char *filename);
-
-
-#endif /* BEAM_PARAMETERS_H */
diff --git a/src/cell.c b/src/cell.c
deleted file mode 100644
index c31697bc..00000000
--- a/src/cell.c
+++ /dev/null
@@ -1,1164 +0,0 @@
-/*
- * cell.c
- *
- * Unit Cell Calculations
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_linalg.h>
-
-#include "cell.h"
-#include "utils.h"
-#include "image.h"
-
-
-/**
- * SECTION:unitcell
- * @short_description: Unit cell
- * @title: UnitCell
- * @section_id:
- * @see_also:
- * @include: "cell.h"
- * @Image:
- *
- * This structure represents a unit cell.
- */
-
-
-
-/* Weighting factor of lengths relative to angles */
-#define LWEIGHT (10.0e-9)
-
-typedef enum {
- CELL_REP_CRYST,
- CELL_REP_CART,
- CELL_REP_RECIP
-} CellRepresentation;
-
-struct _unitcell {
-
- CellRepresentation rep;
-
- /* Crystallographic representation */
- double a; /* m */
- double b; /* m */
- double c; /* m */
- double alpha; /* Radians */
- double beta; /* Radians */
- double gamma; /* Radians */
-
- /* Cartesian representation */
- double ax; double bx; double cx;
- double ay; double by; double cy;
- double az; double bz; double cz;
-
- /* Cartesian representation of reciprocal axes */
- double axs; double bxs; double cxs;
- double ays; double bys; double cys;
- double azs; double bzs; double czs;
-
- char *pointgroup;
- char *spacegroup;
-};
-
-
-/************************** Setters and Constructors **************************/
-
-
-/**
- * cell_new:
- *
- * Create a new %UnitCell.
- *
- * Returns: the new unit cell, or NULL on failure.
- *
- */
-UnitCell *cell_new()
-{
- UnitCell *cell;
-
- cell = malloc(sizeof(UnitCell));
- if ( cell == NULL ) return NULL;
-
- cell->a = 1.0;
- cell->b = 1.0;
- cell->c = 1.0;
- cell->alpha = M_PI_2;
- cell->beta = M_PI_2;
- cell->gamma = M_PI_2;
-
- cell->rep = CELL_REP_CRYST;
-
- cell->pointgroup = strdup("1");
- cell->spacegroup = strdup("P 1");
-
- return cell;
-}
-
-
-/**
- * cell_free:
- * @cell: A %UnitCell to free.
- *
- * Frees a %UnitCell, and all internal resources concerning that cell.
- *
- */
-void cell_free(UnitCell *cell)
-{
- if ( cell == NULL ) return;
- free(cell->pointgroup);
- free(cell->spacegroup);
- free(cell);
-}
-
-
-void cell_set_parameters(UnitCell *cell, double a, double b, double c,
- double alpha, double beta, double gamma)
-{
- if ( cell == NULL ) return;
-
- cell->a = a;
- cell->b = b;
- cell->c = c;
- cell->alpha = alpha;
- cell->beta = beta;
- cell->gamma = gamma;
-
- cell->rep = CELL_REP_CRYST;
-}
-
-
-void cell_set_cartesian(UnitCell *cell,
- double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz)
-{
- if ( cell == NULL ) return;
-
- cell->ax = ax; cell->ay = ay; cell->az = az;
- cell->bx = bx; cell->by = by; cell->bz = bz;
- cell->cx = cx; cell->cy = cy; cell->cz = cz;
-
- cell->rep = CELL_REP_CART;
-}
-
-
-void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az)
-{
- if ( cell == NULL ) return;
- cell->ax = ax; cell->ay = ay; cell->az = az;
- cell->rep = CELL_REP_CART;
-}
-
-
-void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz)
-{
- if ( cell == NULL ) return;
- cell->bx = bx; cell->by = by; cell->bz = bz;
- cell->rep = CELL_REP_CART;
-}
-
-
-void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz)
-{
- if ( cell == NULL ) return;
- cell->cx = cx; cell->cy = cy; cell->cz = cz;
- cell->rep = CELL_REP_CART;
-}
-
-
-UnitCell *cell_new_from_parameters(double a, double b, double c,
- double alpha, double beta, double gamma)
-{
- UnitCell *cell;
-
- cell = cell_new();
- if ( cell == NULL ) return NULL;
-
- cell_set_parameters(cell, a, b, c, alpha, beta, gamma);
-
- return cell;
-}
-
-
-UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
- struct rvec cs)
-{
- UnitCell *cell;
-
- cell = cell_new();
- if ( cell == NULL ) return NULL;
-
- cell->axs = as.u; cell->ays = as.v; cell->azs = as.w;
- cell->bxs = bs.u; cell->bys = bs.v; cell->bzs = bs.w;
- cell->cxs = cs.u; cell->cys = cs.v; cell->czs = cs.w;
-
- cell->rep = CELL_REP_RECIP;
-
- return cell;
-}
-
-
-UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c)
-{
- UnitCell *cell;
-
- cell = cell_new();
- if ( cell == NULL ) return NULL;
-
- cell->ax = a.u; cell->ay = a.v; cell->az = a.w;
- cell->bx = b.u; cell->by = b.v; cell->bz = b.w;
- cell->cx = c.u; cell->cy = c.v; cell->cz = c.w;
-
- cell->rep = CELL_REP_CART;
-
- return cell;
-}
-
-
-UnitCell *cell_new_from_cell(UnitCell *orig)
-{
- UnitCell *new;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
-
- new = cell_new();
-
- cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz);
- cell_set_pointgroup(new, orig->pointgroup);
- cell_set_spacegroup(new, orig->spacegroup);
-
- return new;
-}
-
-
-void cell_set_reciprocal(UnitCell *cell,
- double asx, double asy, double asz,
- double bsx, double bsy, double bsz,
- double csx, double csy, double csz)
-{
- if ( cell == NULL ) return;
-
- cell->axs = asx; cell->ays = asy; cell->azs = asz;
- cell->bxs = bsx; cell->bys = bsy; cell->bzs = bsz;
- cell->cxs = csx; cell->cys = csy; cell->czs = csz;
-
- cell->rep = CELL_REP_RECIP;
-}
-
-
-void cell_set_spacegroup(UnitCell *cell, const char *sym)
-{
- free(cell->spacegroup);
- cell->spacegroup = strdup(sym);
-}
-
-
-void cell_set_pointgroup(UnitCell *cell, const char *sym)
-{
- free(cell->pointgroup);
- cell->pointgroup = strdup(sym);
-}
-
-
-/************************* Getter helper functions ****************************/
-
-static int cell_crystallographic_to_cartesian(UnitCell *cell,
- double *ax, double *ay, double *az,
- double *bx, double *by, double *bz,
- double *cx, double *cy, double *cz)
-{
- double tmp, V, cosalphastar, cstar;
-
- /* Firstly: Get a in terms of x, y and z
- * +a (cryst) is defined to lie along +x (cart) */
- *ax = cell->a;
- *ay = 0.0;
- *az = 0.0;
-
- /* b in terms of x, y and z
- * b (cryst) is defined to lie in the xy (cart) plane */
- *bx = cell->b*cos(cell->gamma);
- *by = cell->b*sin(cell->gamma);
- *bz = 0.0;
-
- tmp = cos(cell->alpha)*cos(cell->alpha)
- + cos(cell->beta)*cos(cell->beta)
- + cos(cell->gamma)*cos(cell->gamma)
- - 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma);
- V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp);
-
- cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha);
- cosalphastar /= sin(cell->beta)*sin(cell->gamma);
-
- cstar = (cell->a * cell->b * sin(cell->gamma))/V;
-
- /* c in terms of x, y and z */
- *cx = cell->c*cos(cell->beta);
- *cy = -cell->c*sin(cell->beta)*cosalphastar;
- *cz = 1.0/cstar;
-
- return 0;
-}
-
-
-/* Why yes, I do enjoy long argument lists...! */
-static int cell_invert(double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz,
- double *asx, double *asy, double *asz,
- double *bsx, double *bsy, double *bsz,
- double *csx, double *csy, double *csz)
-{
- int s;
- gsl_matrix *m;
- gsl_matrix *inv;
- gsl_permutation *perm;
-
- m = gsl_matrix_alloc(3, 3);
- if ( m == NULL ) {
- ERROR("Couldn't allocate memory for matrix\n");
- return 1;
- }
- gsl_matrix_set(m, 0, 0, ax);
- gsl_matrix_set(m, 0, 1, bx);
- gsl_matrix_set(m, 0, 2, cx);
- gsl_matrix_set(m, 1, 0, ay);
- gsl_matrix_set(m, 1, 1, by);
- gsl_matrix_set(m, 1, 2, cy);
- gsl_matrix_set(m, 2, 0, az);
- gsl_matrix_set(m, 2, 1, bz);
- gsl_matrix_set(m, 2, 2, cz);
-
- /* Invert */
- perm = gsl_permutation_alloc(m->size1);
- if ( perm == NULL ) {
- ERROR("Couldn't allocate permutation\n");
- gsl_matrix_free(m);
- return 1;
- }
- inv = gsl_matrix_alloc(m->size1, m->size2);
- if ( inv == NULL ) {
- ERROR("Couldn't allocate inverse\n");
- gsl_matrix_free(m);
- gsl_permutation_free(perm);
- return 1;
- }
- if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
- ERROR("Couldn't decompose matrix\n");
- gsl_matrix_free(m);
- gsl_permutation_free(perm);
- return 1;
- }
- if ( gsl_linalg_LU_invert(m, perm, inv) ) {
- ERROR("Couldn't invert matrix\n");
- gsl_matrix_free(m);
- gsl_permutation_free(perm);
- return 1;
- }
- gsl_permutation_free(perm);
- gsl_matrix_free(m);
-
- /* Transpose */
- gsl_matrix_transpose(inv);
-
- *asx = gsl_matrix_get(inv, 0, 0);
- *bsx = gsl_matrix_get(inv, 0, 1);
- *csx = gsl_matrix_get(inv, 0, 2);
- *asy = gsl_matrix_get(inv, 1, 0);
- *bsy = gsl_matrix_get(inv, 1, 1);
- *csy = gsl_matrix_get(inv, 1, 2);
- *asz = gsl_matrix_get(inv, 2, 0);
- *bsz = gsl_matrix_get(inv, 2, 1);
- *csz = gsl_matrix_get(inv, 2, 2);
-
- gsl_matrix_free(inv);
-
- return 0;
-}
-
-
-/********************************** Getters ***********************************/
-
-int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
- double *alpha, double *beta, double *gamma)
-{
- double ax, ay, az, bx, by, bz, cx, cy, cz;
-
- if ( cell == NULL ) return 1;
-
- switch ( cell->rep ) {
-
- case CELL_REP_CRYST:
- /* Direct response */
- *a = cell->a;
- *b = cell->b;
- *c = cell->c;
- *alpha = cell->alpha;
- *beta = cell->beta;
- *gamma = cell->gamma;
- return 0;
-
- case CELL_REP_CART:
- /* Convert cartesian -> crystallographic */
- *a = modulus(cell->ax, cell->ay, cell->az);
- *b = modulus(cell->bx, cell->by, cell->bz);
- *c = modulus(cell->cx, cell->cy, cell->cz);
-
- *alpha = angle_between(cell->bx, cell->by, cell->bz,
- cell->cx, cell->cy, cell->cz);
- *beta = angle_between(cell->ax, cell->ay, cell->az,
- cell->cx, cell->cy, cell->cz);
- *gamma = angle_between(cell->ax, cell->ay, cell->az,
- cell->bx, cell->by, cell->bz);
- return 0;
-
- case CELL_REP_RECIP:
- /* Convert reciprocal -> crystallographic.
- * Start by converting reciprocal -> cartesian */
- cell_invert(cell->axs, cell->ays, cell->azs,
- cell->bxs, cell->bys, cell->bzs,
- cell->cxs, cell->cys, cell->czs,
- &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
- /* Now convert cartesian -> crystallographic */
- *a = modulus(ax, ay, az);
- *b = modulus(bx, by, bz);
- *c = modulus(cx, cy, cz);
-
- *alpha = angle_between(bx, by, bz, cx, cy, cz);
- *beta = angle_between(ax, ay, az, cx, cy, cz);
- *gamma = angle_between(ax, ay, az, bx, by, bz);
- return 0;
- }
-
- return 1;
-}
-
-
-int cell_get_cartesian(UnitCell *cell,
- double *ax, double *ay, double *az,
- double *bx, double *by, double *bz,
- double *cx, double *cy, double *cz)
-{
- if ( cell == NULL ) return 1;
-
- switch ( cell->rep ) {
-
- case CELL_REP_CRYST:
- /* Convert crystallographic -> cartesian. */
- return cell_crystallographic_to_cartesian(cell,
- ax, ay, az,
- bx, by, bz,
- cx, cy, cz);
-
- case CELL_REP_CART:
- /* Direct response */
- *ax = cell->ax; *ay = cell->ay; *az = cell->az;
- *bx = cell->bx; *by = cell->by; *bz = cell->bz;
- *cx = cell->cx; *cy = cell->cy; *cz = cell->cz;
- return 0;
-
- case CELL_REP_RECIP:
- /* Convert reciprocal -> cartesian */
- return cell_invert(cell->axs, cell->ays, cell->azs,
- cell->bxs, cell->bys, cell->bzs,
- cell->cxs, cell->cys, cell->czs,
- ax, ay, az, bx, by, bz, cx, cy, cz);
-
- }
-
- return 1;
-}
-
-
-int cell_get_reciprocal(UnitCell *cell,
- double *asx, double *asy, double *asz,
- double *bsx, double *bsy, double *bsz,
- double *csx, double *csy, double *csz)
-{
- int r;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
- if ( cell == NULL ) return 1;
-
- switch ( cell->rep ) {
-
- case CELL_REP_CRYST:
- /* Convert crystallographic -> reciprocal */
- r = cell_crystallographic_to_cartesian(cell,
- &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
- if ( r ) return r;
- return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz,
- asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
- case CELL_REP_CART:
- /* Convert cartesian -> reciprocal */
- cell_invert(cell->ax, cell->ay, cell->az,
- cell->bx, cell->by, cell->bz,
- cell->cx, cell->cy, cell->cz,
- asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
- return 0;
-
- case CELL_REP_RECIP:
- /* Direct response */
- *asx = cell->axs; *asy = cell->ays; *asz = cell->azs;
- *bsx = cell->bxs; *bsy = cell->bys; *bsz = cell->bzs;
- *csx = cell->cxs; *csy = cell->cys; *csz = cell->czs;
- return 0;
-
- }
-
- return 1;
-}
-
-
-const char *cell_get_pointgroup(UnitCell *cell)
-{
- return cell->pointgroup;
-}
-
-
-const char *cell_get_spacegroup(UnitCell *cell)
-{
- return cell->spacegroup;
-}
-
-
-
-
-
-/********************************* Utilities **********************************/
-
-static const char *cell_rep(UnitCell *cell)
-{
- switch ( cell->rep ) {
- case CELL_REP_CRYST:
- return "crystallographic, direct space";
- case CELL_REP_CART:
- return "cartesian, direct space";
- case CELL_REP_RECIP:
- return "cartesian, reciprocal space";
- }
-
- return "unknown";
-}
-
-
-/**
- * cell_rotate:
- * @in: A %UnitCell to rotate
- * @quat: A %quaternion
- *
- * Rotate a %UnitCell using a %quaternion.
- *
- * Returns: a newly allocated rotated copy of @in.
- *
- */
-UnitCell *cell_rotate(UnitCell *in, struct quaternion quat)
-{
- struct rvec a, b, c;
- struct rvec an, bn, cn;
- UnitCell *out = cell_new_from_cell(in);
-
- cell_get_cartesian(in, &a.u, &a.v, &a.w,
- &b.u, &b.v, &b.w,
- &c.u, &c.v, &c.w);
-
- an = quat_rot(a, quat);
- bn = quat_rot(b, quat);
- cn = quat_rot(c, quat);
-
- cell_set_cartesian(out, an.u, an.v, an.w,
- bn.u, bn.v, bn.w,
- cn.u, cn.v, cn.w);
-
- return out;
-}
-
-
-void cell_print(UnitCell *cell)
-{
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double a, b, c, alpha, beta, gamma;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
-
- cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
-
- STATUS(" a b c alpha beta gamma\n");
- STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n",
- a*1e9, b*1e9, c*1e9,
- rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
-
- cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
- STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az);
- STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz);
- STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz);
-
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
-
- STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
- asx, asy, asz, modulus(asx, asy, asz));
- STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
- bsx, bsy, bsz, modulus(bsx, bsy, bsz));
- STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
- csx, csy, csz, modulus(csx, csy, csz));
-
- STATUS("Point group: %s\n", cell_get_pointgroup(cell));
- STATUS("Cell representation is %s.\n", cell_rep(cell));
-}
-
-
-#define MAX_CAND (1024)
-
-static int right_handed(struct rvec a, struct rvec b, struct rvec c)
-{
- struct rvec aCb;
- double aCb_dot_c;
-
- /* "a" cross "b" */
- aCb.u = a.v*b.w - a.w*b.v;
- aCb.v = - (a.u*b.w - a.w*b.u);
- aCb.w = a.u*b.v - a.v*b.u;
-
- /* "a cross b" dot "c" */
- aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w;
-
- if ( aCb_dot_c > 0.0 ) return 1;
- return 0;
-}
-
-
-struct cvec {
- struct rvec vec;
- float na;
- float nb;
- float nc;
- float fom;
-};
-
-
-static int same_vector(struct cvec a, struct cvec b)
-{
- if ( a.na != b.na ) return 0;
- if ( a.nb != b.nb ) return 0;
- if ( a.nc != b.nc ) return 0;
- return 1;
-}
-
-
-/* Attempt to make 'cell' fit into 'template' somehow */
-UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
- int reduce)
-{
- signed int n1l, n2l, n3l;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- int i, j;
- double lengths[3];
- double angles[3];
- struct cvec *cand[3];
- UnitCell *new_cell = NULL;
- float best_fom = +999999999.9; /* Large number.. */
- int ncand[3] = {0,0,0};
- signed int ilow, ihigh;
- float ltl = 5.0; /* percent */
- float angtol = deg2rad(1.5);
-
- if ( cell_get_reciprocal(template, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz) ) {
- ERROR("Couldn't get reciprocal cell for template.\n");
- return NULL;
- }
-
- lengths[0] = modulus(asx, asy, asz);
- lengths[1] = modulus(bsx, bsy, bsz);
- lengths[2] = modulus(csx, csy, csz);
-
- angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
- angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
- angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
-
- cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
- cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
- cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
-
- if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz) ) {
- ERROR("Couldn't get reciprocal cell.\n");
- return NULL;
- }
-
- if ( reduce ) {
- ilow = -2; ihigh = 4;
- } else {
- ilow = 0; ihigh = 1;
- }
-
- /* Negative values mean 1/n, positive means n, zero means zero */
- for ( n1l=ilow; n1l<=ihigh; n1l++ ) {
- for ( n2l=ilow; n2l<=ihigh; n2l++ ) {
- for ( n3l=ilow; n3l<=ihigh; n3l++ ) {
-
- float n1, n2, n3;
- signed int b1, b2, b3;
-
- n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
- n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
- n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
-
- if ( !reduce ) {
- if ( n1l + n2l + n3l > 1 ) continue;
- }
-
- /* 'bit' values can be +1 or -1 */
- for ( b1=-1; b1<=1; b1+=2 ) {
- for ( b2=-1; b2<=1; b2+=2 ) {
- for ( b3=-1; b3<=1; b3+=2 ) {
-
- double tx, ty, tz;
- double tlen;
- int i;
-
- n1 *= b1; n2 *= b2; n3 *= b3;
-
- tx = n1*asx + n2*bsx + n3*csx;
- ty = n1*asy + n2*bsy + n3*csy;
- tz = n1*asz + n2*bsz + n3*csz;
- tlen = modulus(tx, ty, tz);
-
- /* Test modulus for agreement with moduli of template */
- for ( i=0; i<3; i++ ) {
-
- if ( !within_tolerance(lengths[i], tlen, ltl) )
- continue;
-
- cand[i][ncand[i]].vec.u = tx;
- cand[i][ncand[i]].vec.v = ty;
- cand[i][ncand[i]].vec.w = tz;
- cand[i][ncand[i]].na = n1;
- cand[i][ncand[i]].nb = n2;
- cand[i][ncand[i]].nc = n3;
- cand[i][ncand[i]].fom = fabs(lengths[i] - tlen);
- if ( ncand[i] == MAX_CAND ) {
- ERROR("Too many candidates\n");
- } else {
- ncand[i]++;
- }
- }
-
- }
- }
- }
-
- }
- }
- }
-
- if ( verbose ) {
- STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]);
- }
-
- for ( i=0; i<ncand[0]; i++ ) {
- for ( j=0; j<ncand[1]; j++ ) {
-
- double ang;
- int k;
- float fom1;
-
- if ( same_vector(cand[0][i], cand[1][j]) ) continue;
-
- /* Measure the angle between the ith candidate for axis 0
- * and the jth candidate for axis 1 */
- ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
- cand[0][i].vec.w, cand[1][j].vec.u,
- cand[1][j].vec.v, cand[1][j].vec.w);
-
- /* Angle between axes 0 and 1 should be angle 2 */
- if ( fabs(ang - angles[2]) > angtol ) continue;
-
- fom1 = fabs(ang - angles[2]);
-
- for ( k=0; k<ncand[2]; k++ ) {
-
- float fom2, fom3;
-
- if ( same_vector(cand[1][j], cand[2][k]) ) continue;
-
- /* Measure the angle between the current candidate for
- * axis 0 and the kth candidate for axis 2 */
- ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
- cand[0][i].vec.w, cand[2][k].vec.u,
- cand[2][k].vec.v, cand[2][k].vec.w);
-
- /* ... it should be angle 1 ... */
- if ( fabs(ang - angles[1]) > angtol ) continue;
-
- fom2 = fom1 + fabs(ang - angles[1]);
-
- /* Finally, the angle between the current candidate for
- * axis 1 and the kth candidate for axis 2 */
- ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
- cand[1][j].vec.w, cand[2][k].vec.u,
- cand[2][k].vec.v, cand[2][k].vec.w);
-
- /* ... it should be angle 0 ... */
- if ( fabs(ang - angles[0]) > angtol ) continue;
-
- /* Unit cell must be right-handed */
- if ( !right_handed(cand[0][i].vec, cand[1][j].vec,
- cand[2][k].vec) ) continue;
-
- fom3 = fom2 + fabs(ang - angles[0]);
- fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom
- + cand[2][k].fom);
-
- if ( fom3 < best_fom ) {
- if ( new_cell != NULL ) free(new_cell);
- new_cell = cell_new_from_reciprocal_axes(
- cand[0][i].vec, cand[1][j].vec,
- cand[2][k].vec);
- best_fom = fom3;
- }
-
- }
-
- }
- }
-
- free(cand[0]);
- free(cand[1]);
- free(cand[2]);
-
- return new_cell;
-}
-
-
-
-UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- int i;
- double lengths[3];
- int used[3];
- struct rvec real_a, real_b, real_c;
- struct rvec params[3];
- double alen, blen, clen;
- float ltl = 5.0; /* percent */
- int have_real_a;
- int have_real_b;
- int have_real_c;
-
- /* Get the lengths to match */
- if ( cell_get_cartesian(template, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz) )
- {
- ERROR("Couldn't get cell for template.\n");
- return NULL;
- }
- alen = modulus(ax, ay, az);
- blen = modulus(bx, by, bz);
- clen = modulus(cx, cy, cz);
-
- /* Get the lengths from the cell and turn them into anonymous vectors */
- if ( cell_get_cartesian(cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz) )
- {
- ERROR("Couldn't get cell.\n");
- return NULL;
- }
- lengths[0] = modulus(ax, ay, az);
- lengths[1] = modulus(bx, by, bz);
- lengths[2] = modulus(cx, cy, cz);
- used[0] = 0; used[1] = 0; used[2] = 0;
- params[0].u = ax; params[0].v = ay; params[0].w = az;
- params[1].u = bx; params[1].v = by; params[1].w = bz;
- params[2].u = cx; params[2].v = cy; params[2].w = cz;
-
- real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0;
- real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0;
- real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0;
-
- /* Check each vector against a and b */
- have_real_a = 0;
- have_real_b = 0;
- for ( i=0; i<3; i++ ) {
- if ( within_tolerance(lengths[i], alen, ltl)
- && !used[i] && !have_real_a )
- {
- used[i] = 1;
- memcpy(&real_a, &params[i], sizeof(struct rvec));
- have_real_a = 1;
- }
- if ( within_tolerance(lengths[i], blen, ltl)
- && !used[i] && !have_real_b )
- {
- used[i] = 1;
- memcpy(&real_b, &params[i], sizeof(struct rvec));
- have_real_b = 1;
- }
- }
-
- /* Have we matched both a and b? */
- if ( !(have_real_a && have_real_b) ) return NULL;
-
- /* "c" is "the other one" */
- have_real_c = 0;
- for ( i=0; i<3; i++ ) {
- if ( !used[i] ) {
- memcpy(&real_c, &params[i], sizeof(struct rvec));
- have_real_c = 1;
- }
- }
-
- if ( !have_real_c ) {
- ERROR("Huh? Couldn't find the third vector.\n");
- ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]);
- return NULL;
- }
-
- /* Flip c if not right-handed */
- if ( !right_handed(real_a, real_b, real_c) ) {
- real_c.u = -real_c.u;
- real_c.v = -real_c.v;
- real_c.w = -real_c.w;
- }
-
- return cell_new_from_direct_axes(real_a, real_b, real_c);
-}
-
-
-/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */
-double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
-{
- double a, b, c, alpha, beta, gamma;
-
- cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
-
- const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha)
- - cos(beta)*cos(beta)
- - cos(gamma)*cos(gamma)
- + 2*cos(alpha)*cos(beta)*cos(gamma) );
-
- const double S11 = b*b*c*c*sin(alpha)*sin(alpha);
- const double S22 = a*a*c*c*sin(beta)*sin(beta);
- const double S33 = a*a*b*b*sin(gamma)*sin(gamma);
- const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma));
- const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha));
- const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta));
-
- const double brackets = S11*h*h + S22*k*k + S33*l*l
- + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
- const double oneoverdsq = brackets / Vsq;
- const double oneoverd = sqrt(oneoverdsq);
-
- return oneoverd / 2;
-}
-
-
-static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym)
-{
- char *new = NULL;
-
- if ( strcmp(sym, "P 1") == 0 ) new = "1";
- if ( strcmp(sym, "P 63") == 0 ) new = "6";
- if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222";
- if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222";
- if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422";
-
- if ( new != NULL ) {
- cell_set_pointgroup(cell, new);
- } else {
- ERROR("Can't determine point group for '%s'\n", sym);
- }
-}
-
-
-UnitCell *load_cell_from_pdb(const char *filename)
-{
- FILE *fh;
- char *rval;
- UnitCell *cell = NULL;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) {
- ERROR("Couldn't open '%s'\n", filename);
- return NULL;
- }
-
- do {
-
- char line[1024];
-
- rval = fgets(line, 1023, fh);
-
- if ( strncmp(line, "CRYST1", 6) == 0 ) {
-
- float a, b, c, al, be, ga;
- char as[10], bs[10], cs[10];
- char als[8], bes[8], gas[8];
- char *sym;
- int r;
-
- memcpy(as, line+6, 9); as[9] = '\0';
- memcpy(bs, line+15, 9); bs[9] = '\0';
- memcpy(cs, line+24, 9); cs[9] = '\0';
- memcpy(als, line+33, 7); als[7] = '\0';
- memcpy(bes, line+40, 7); bes[7] = '\0';
- memcpy(gas, line+47, 7); gas[7] = '\0';
-
- r = sscanf(as, "%f", &a);
- r += sscanf(bs, "%f", &b);
- r += sscanf(cs, "%f", &c);
- r += sscanf(als, "%f", &al);
- r += sscanf(bes, "%f", &be);
- r += sscanf(gas, "%f", &ga);
-
- if ( r != 6 ) {
- STATUS("Couldn't understand CRYST1 line.\n");
- continue;
- }
-
- cell = cell_new_from_parameters(a*1e-10,
- b*1e-10, c*1e-10,
- deg2rad(al),
- deg2rad(be),
- deg2rad(ga));
-
- if ( strlen(line) > 65 ) {
- sym = strndup(line+55, 10);
- notrail(sym);
- cell_set_pointgroup_from_pdb(cell, sym);
- cell_set_spacegroup(cell, sym);
- free(sym);
- } else {
- cell_set_pointgroup_from_pdb(cell, "P 1");
- cell_set_spacegroup(cell, "P 1");
- ERROR("CRYST1 line without space group.\n");
- }
-
- break; /* Done */
- }
-
- } while ( rval != NULL );
-
- fclose(fh);
-
- return cell;
-}
-
-
-#ifdef GSL_FUDGE
-/* Force the linker to bring in CBLAS to make GSL happy */
-void cell_fudge_gslcblas()
-{
- STATUS("%p\n", cblas_sgemm);
-}
-#endif
-
-
-UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot)
-{
- UnitCell *out;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double xnew, ynew, znew;
-
- cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
- &bsz, &csx, &csy, &csz);
-
- /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */
- xnew = asx*cos(omega) + asy*sin(omega);
- ynew = -asx*sin(omega) + asy*cos(omega);
- znew = asz;
- asx = xnew; asy = ynew; asz = znew;
- xnew = bsx*cos(omega) + bsy*sin(omega);
- ynew = -bsx*sin(omega) + bsy*cos(omega);
- znew = bsz;
- bsx = xnew; bsy = ynew; bsz = znew;
- xnew = csx*cos(omega) + csy*sin(omega);
- ynew = -csx*sin(omega) + csy*cos(omega);
- znew = csz;
- csx = xnew; csy = ynew; csz = znew;
-
- /* Rotate by "phi" about +x (not parallel to anything specific) */
- xnew = asx;
- ynew = asy*cos(phi) + asz*sin(phi);
- znew = -asy*sin(phi) + asz*cos(phi);
- asx = xnew; asy = ynew; asz = znew;
- xnew = bsx;
- ynew = bsy*cos(phi) + bsz*sin(phi);
- znew = -bsy*sin(phi) + bsz*cos(phi);
- bsx = xnew; bsy = ynew; bsz = znew;
- xnew = csx;
- ynew = csy*cos(phi) + csz*sin(phi);
- znew = -csy*sin(phi) + csz*cos(phi);
- csx = xnew; csy = ynew; csz = znew;
-
- /* Rotate by "rot" about the new +z (in-plane rotation) */
- xnew = asx*cos(rot) + asy*sin(rot);
- ynew = -asx*sin(rot) + asy*cos(rot);
- znew = asz;
- asx = xnew; asy = ynew; asz = znew;
- xnew = bsx*cos(rot) + bsy*sin(rot);
- ynew = -bsx*sin(rot) + bsy*cos(rot);
- znew = bsz;
- bsx = xnew; bsy = ynew; bsz = znew;
- xnew = csx*cos(rot) + csy*sin(rot);
- ynew = -csx*sin(rot) + csy*cos(rot);
- znew = csz;
- csx = xnew; csy = ynew; csz = znew;
-
- out = cell_new_from_cell(in);
- cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
- return out;
-}
-
-
-int cell_is_sensible(UnitCell *cell)
-{
- double a, b, c, al, be, ga;
-
- cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
- if ( al + be + ga >= 2.0*M_PI ) return 0;
- if ( al + be - ga >= 2.0*M_PI ) return 0;
- if ( al - be + ga >= 2.0*M_PI ) return 0;
- if ( - al + be + ga >= 2.0*M_PI ) return 0;
- if ( al + be + ga <= 0.0 ) return 0;
- if ( al + be - ga <= 0.0 ) return 0;
- if ( al - be + ga <= 0.0 ) return 0;
- if ( - al + be + ga <= 0.0 ) return 0;
- if ( isnan(al) ) return 0;
- if ( isnan(be) ) return 0;
- if ( isnan(ga) ) return 0;
- return 1;
-}
diff --git a/src/cell.h b/src/cell.h
deleted file mode 100644
index b5d31fc6..00000000
--- a/src/cell.h
+++ /dev/null
@@ -1,107 +0,0 @@
-/*
- * cell.h
- *
- * Unit Cell Calculations
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifndef CELL_H
-#define CELL_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "utils.h"
-
-/* A 3D vector in reciprocal space */
-struct rvec
-{
- double u;
- double v;
- double w;
-};
-
-
-/**
- * UnitCell:
- *
- * This data structure is opaque. You must use the available accessor functions
- * to read and write its contents.
- **/
-typedef struct _unitcell UnitCell;
-
-extern UnitCell *cell_new(void);
-extern UnitCell *cell_new_from_cell(UnitCell *orig);
-extern void cell_free(UnitCell *cell);
-
-/* Lengths in m, angles in radians */
-extern UnitCell *cell_new_from_parameters(double a, double b, double c,
- double alpha, double beta, double gamma);
-
-extern UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
- struct rvec cs);
-
-extern UnitCell *cell_new_from_direct_axes(struct rvec as, struct rvec bs,
- struct rvec cs);
-
-extern void cell_set_cartesian(UnitCell *cell,
- double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz);
-
-extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
- double alpha, double beta, double gamma);
-
-extern void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az);
-extern void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz);
-extern void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz);
-extern void cell_set_spacegroup(UnitCell *cell, const char *sym);
-extern void cell_set_pointgroup(UnitCell *cell, const char *sym);
-
-
-extern int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
- double *alpha, double *beta, double *gamma);
-
-extern int cell_get_cartesian(UnitCell *cell,
- double *ax, double *ay, double *az,
- double *bx, double *by, double *bz,
- double *cx, double *cy, double *cz);
-
-extern int cell_get_reciprocal(UnitCell *cell,
- double *asx, double *asy, double *asz,
- double *bsx, double *bsy, double *bsz,
- double *csx, double *csy, double *csz);
-
-extern void cell_set_reciprocal(UnitCell *cell,
- double asx, double asy, double asz,
- double bsx, double bsy, double bsz,
- double csx, double csy, double csz);
-
-extern const char *cell_get_pointgroup(UnitCell *cell);
-
-extern const char *cell_get_spacegroup(UnitCell *cell);
-
-extern double resolution(UnitCell *cell,
- signed int h, signed int k, signed int l);
-
-extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
-extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
- double rot);
-
-extern void cell_print(UnitCell *cell);
-
-extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
- int reduce);
-
-extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template);
-
-extern UnitCell *load_cell_from_pdb(const char *filename);
-
-extern int cell_is_sensible(UnitCell *cell);
-
-#endif /* CELL_H */
diff --git a/src/detector.c b/src/detector.c
deleted file mode 100644
index 7f29ec82..00000000
--- a/src/detector.c
+++ /dev/null
@@ -1,1234 +0,0 @@
-/*
- * detector.c
- *
- * Detector properties
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- * (c) 2011 Rick Kirian <rkirian@asu.edu>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <assert.h>
-#include <ctype.h>
-
-#include "image.h"
-#include "utils.h"
-#include "diffraction.h"
-#include "detector.h"
-#include "beam-parameters.h"
-#include "hdf5-file.h"
-
-
-/**
- * SECTION:detector
- * @short_description: Detector geometry
- * @title: Detector
- * @section_id:
- * @see_also:
- * @include: "detector.h"
- * @Image:
- *
- * This structure represents the detector geometry
- */
-
-
-static int atob(const char *a)
-{
- if ( strcasecmp(a, "true") == 0 ) return 1;
- if ( strcasecmp(a, "false") == 0 ) return 0;
- return atoi(a);
-}
-
-
-static int assplode_algebraic(const char *a_orig, char ***pbits)
-{
- int len, i;
- int nexp;
- char **bits;
- char *a;
- int idx, istr;
-
- len = strlen(a_orig);
-
- /* Add plus at start if no sign there already */
- if ( (a_orig[0] != '+') && (a_orig[0] != '-') ) {
- len += 1;
- a = malloc(len+1);
- snprintf(a, len+1, "+%s", a_orig);
- a[len] = '\0';
-
- } else {
- a = strdup(a_orig);
- }
-
- /* Count the expressions */
- nexp = 0;
- for ( i=0; i<len; i++ ) {
- if ( (a[i] == '+') || (a[i] == '-') ) nexp++;
- }
-
- bits = calloc(nexp, sizeof(char *));
-
- /* Break the string up */
- idx = -1;
- istr = 0;
- assert((a[0] == '+') || (a[0] == '-'));
- for ( i=0; i<len; i++ ) {
-
- char ch;
-
- ch = a[i];
-
- if ( (ch == '+') || (ch == '-') ) {
- if ( idx >= 0 ) bits[idx][istr] = '\0';
- idx++;
- bits[idx] = malloc(len+1);
- istr = 0;
- }
-
- if ( !isdigit(ch) && (ch != '.') && (ch != 'x') && (ch != 'y')
- && (ch != '+') && (ch != '-') )
- {
- ERROR("Invalid character '%C' found.\n", ch);
- return 0;
- }
-
- assert(idx >= 0);
- bits[idx][istr++] = ch;
-
- }
- if ( idx >= 0 ) bits[idx][istr] = '\0';
-
- *pbits = bits;
- free(a);
-
- return nexp;
-}
-
-
-/* Parses the scan directions (accounting for possible rotation)
- * Assumes all white spaces have been already removed */
-static int dir_conv(const char *a, double *sx, double *sy)
-{
- int n;
- char **bits;
- int i;
-
- *sx = 0.0; *sy = 0.0;
-
- n = assplode_algebraic(a, &bits);
-
- if ( n == 0 ) {
- ERROR("Invalid direction '%s'\n", a);
- return 1;
- }
-
- for ( i=0; i<n; i++ ) {
-
- int len;
- double val;
- char axis;
- int j;
-
- len = strlen(bits[i]);
- assert(len != 0);
- axis = bits[i][len-1];
- if ( (axis != 'x') && (axis != 'y') ) {
- ERROR("Invalid symbol '%C' - must be x or y.\n", axis);
- return 1;
- }
-
- /* Chop off the symbol now it's dealt with */
- bits[i][len-1] = '\0';
-
- /* Check for anything that isn't part of a number */
- for ( j=0; j<strlen(bits[i]); j++ ) {
- if ( isdigit(bits[i][j]) ) continue;
- if ( bits[i][j] == '+' ) continue;
- if ( bits[i][j] == '-' ) continue;
- if ( bits[i][j] == '.' ) continue;
- ERROR("Invalid coefficient '%s'\n", bits[i]);
- }
-
- if ( strlen(bits[i]) == 0 ) {
- val = 1.0;
- } else {
- val = atof(bits[i]);
- }
- if ( strlen(bits[i]) == 1 ) {
- if ( bits[i][0] == '+' ) val = 1.0;
- if ( bits[i][0] == '-' ) val = -1.0;
- }
- if ( axis == 'x' ) {
- *sx += val;
- } else if ( axis == 'y' ) {
- *sy += val;
- }
-
- free(bits[i]);
-
- }
- free(bits);
-
- //STATUS("'%s' -> %5.2fx + %5.2fy\n", a, *sx, *sy);
-
- return 0;
-}
-
-
-struct rvec get_q_for_panel(struct panel *p, double fs, double ss,
- double *ttp, double k)
-{
- struct rvec q;
- double twotheta, r, az;
- double rx, ry;
- double xs, ys;
-
- /* Convert xs and ys, which are in fast scan/slow scan coordinates,
- * to x and y */
- xs = fs*p->fsx + ss*p->ssx;
- ys = fs*p->fsy + ss*p->ssy;
-
- rx = (xs + p->cnx) / p->res;
- ry = (ys + p->cny) / p->res;
-
- /* Calculate q-vector for this sub-pixel */
- r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
-
- twotheta = atan2(r, p->clen);
- az = atan2(ry, rx);
- if ( ttp != NULL ) *ttp = twotheta;
-
- q.u = k * sin(twotheta)*cos(az);
- q.v = k * sin(twotheta)*sin(az);
- q.w = k * (cos(twotheta) - 1.0);
-
- return q;
-}
-
-
-struct rvec get_q(struct image *image, double fs, double ss,
- double *ttp, double k)
-{
- struct panel *p;
- const unsigned int fsi = fs;
- const unsigned int ssi = ss; /* Explicit rounding */
-
- /* Determine which panel to use */
- p = find_panel(image->det, fsi, ssi);
- assert(p != NULL);
-
- return get_q_for_panel(p, fs-(double)p->min_fs, ss-(double)p->min_ss,
- ttp, k);
-}
-
-
-int in_bad_region(struct detector *det, double fs, double ss)
-{
- double rx, ry;
- struct panel *p;
- double xs, ys;
- int i;
-
- /* Determine which panel to use */
- p = find_panel(det, fs, ss);
-
- /* No panel found -> definitely bad! */
- if ( p == NULL ) return 1;
-
- /* Convert xs and ys, which are in fast scan/slow scan coordinates,
- * to x and y */
- xs = (fs-(double)p->min_fs)*p->fsx + (ss-(double)p->min_ss)*p->ssx;
- ys = (fs-(double)p->min_fs)*p->fsy + (ss-(double)p->min_ss)*p->ssy;
-
- rx = xs + p->cnx;
- ry = ys + p->cny;
-
- for ( i=0; i<det->n_bad; i++ ) {
- struct badregion *b = &det->bad[i];
- if ( rx < b->min_x ) continue;
- if ( rx > b->max_x ) continue;
- if ( ry < b->min_y ) continue;
- if ( ry > b->max_y ) continue;
- return 1;
- }
-
- return 0;
-}
-
-
-double get_tt(struct image *image, double fs, double ss)
-{
- double r, rx, ry;
- struct panel *p;
- double xs, ys;
-
- p = find_panel(image->det, fs, ss);
-
- /* Convert xs and ys, which are in fast scan/slow scan coordinates,
- * to x and y */
- xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx;
- ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy;
-
- rx = (xs + p->cnx) / p->res;
- ry = (ys + p->cny) / p->res;
-
- r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
-
- return atan2(r, p->clen);
-}
-
-
-void record_image(struct image *image, int do_poisson)
-{
- int x, y;
- double total_energy, energy_density;
- double ph_per_e;
- double area;
- double max_tt = 0.0;
- int n_inf1 = 0;
- int n_neg1 = 0;
- int n_nan1 = 0;
- int n_inf2 = 0;
- int n_neg2 = 0;
- int n_nan2 = 0;
-
- /* How many photons are scattered per electron? */
- area = M_PI*pow(image->beam->beam_radius, 2.0);
- total_energy = image->beam->fluence * ph_lambda_to_en(image->lambda);
- energy_density = total_energy / area;
- ph_per_e = (image->beam->fluence /area) * pow(THOMSON_LENGTH, 2.0);
- STATUS("Fluence = %8.2e photons, "
- "Energy density = %5.3f kJ/cm^2, "
- "Total energy = %5.3f microJ\n",
- image->beam->fluence, energy_density/1e7, total_energy*1e6);
-
- for ( x=0; x<image->width; x++ ) {
- for ( y=0; y<image->height; y++ ) {
-
- double counts;
- double cf;
- double intensity, sa;
- double pix_area, Lsq;
- double xs, ys, rx, ry;
- double dsq, proj_area;
- struct panel *p;
-
- intensity = (double)image->data[x + image->width*y];
- if ( isinf(intensity) ) n_inf1++;
- if ( intensity < 0.0 ) n_neg1++;
- if ( isnan(intensity) ) n_nan1++;
-
- p = find_panel(image->det, x, y);
-
- /* Area of one pixel */
- pix_area = pow(1.0/p->res, 2.0);
- Lsq = pow(p->clen, 2.0);
-
- /* Area of pixel as seen from crystal (approximate) */
- proj_area = pix_area * cos(image->twotheta[x + image->width*y]);
-
- /* Calculate distance from crystal to pixel */
- xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx;
- ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy;
- rx = (xs + p->cnx) / p->res;
- ry = (ys + p->cny) / p->res;
- dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
-
- /* Projected area of pixel divided by distance squared */
- sa = proj_area / (dsq + Lsq);
-
- if ( do_poisson ) {
- counts = poisson_noise(intensity * ph_per_e
- * sa * image->beam->dqe );
- } else {
- cf = intensity * ph_per_e * sa * image->beam->dqe;
- counts = cf;
- }
-
- image->data[x + image->width*y] = counts
- * image->beam->adu_per_photon;
-
- /* Sanity checks */
- if ( isinf(image->data[x+image->width*y]) ) n_inf2++;
- if ( isnan(image->data[x+image->width*y]) ) n_nan2++;
- if ( image->data[x+image->width*y] < 0.0 ) n_neg2++;
-
- if ( image->twotheta[x + image->width*y] > max_tt ) {
- max_tt = image->twotheta[x + image->width*y];
- }
-
- }
- progress_bar(x, image->width-1, "Post-processing");
- }
-
- STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
- rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
-
- double tt_side = image->twotheta[(image->width/2)+image->width*0];
- STATUS("At middle of bottom edge: %.2f deg, min d = %.2f nm\n",
- rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9);
-
- tt_side = image->twotheta[0+image->width*(image->height/2)];
- STATUS("At middle of left edge: %.2f deg, min d = %.2f nm\n",
- rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9);
-
- STATUS("Halve the d values to get the voxel size for a synthesis.\n");
-
- if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) {
- ERROR("WARNING: The raw calculation produced %i negative"
- " values, %i infinities and %i NaNs.\n",
- n_neg1, n_inf1, n_nan1);
- ERROR("WARNING: After processing, there were %i negative"
- " values, %i infinities and %i NaNs.\n",
- n_neg2, n_inf2, n_nan2);
- }
-}
-
-
-struct panel *find_panel(struct detector *det, double fs, double ss)
-{
- int p;
-
- for ( p=0; p<det->n_panels; p++ ) {
- if ( (fs >= det->panels[p].min_fs)
- && (fs <= det->panels[p].max_fs)
- && (ss >= det->panels[p].min_ss)
- && (ss <= det->panels[p].max_ss) ) {
- return &det->panels[p];
- }
- }
-
- return NULL;
-}
-
-
-int find_panel_number(struct detector *det, int fs, int ss)
-{
- int p;
-
- for ( p=0; p<det->n_panels; p++ ) {
- if ( (fs >= det->panels[p].min_fs)
- && (fs <= det->panels[p].max_fs)
- && (ss >= det->panels[p].min_ss)
- && (ss <= det->panels[p].max_ss) ) return p;
- }
-
- return -1;
-}
-
-
-void fill_in_values(struct detector *det, struct hdfile *f)
-{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
-
- struct panel *p = &det->panels[i];
-
- if ( p->clen_from != NULL ) {
- p->clen = get_value(f, p->clen_from) * 1.0e-3;
- free(p->clen_from);
- p->clen_from = NULL;
- }
-
- p->clen += p->coffset;
-
- }
-}
-
-
-static struct panel *new_panel(struct detector *det, const char *name)
-{
- struct panel *new;
-
- det->n_panels++;
- det->panels = realloc(det->panels, det->n_panels*sizeof(struct panel));
-
- new = &det->panels[det->n_panels-1];
- memcpy(new, &det->defaults, sizeof(struct panel));
-
- /* Create a new copy of the camera length location if needed */
- if ( new->clen_from != NULL ) {
- new->clen_from = strdup(new->clen_from);
- }
- strcpy(new->name, name);
-
- return new;
-}
-
-
-static struct badregion *new_bad_region(struct detector *det, const char *name)
-{
- struct badregion *new;
-
- det->n_bad++;
- det->bad = realloc(det->bad, det->n_bad*sizeof(struct badregion));
-
- new = &det->bad[det->n_bad-1];
- new->min_x = NAN;
- new->max_x = NAN;
- new->min_y = NAN;
- new->max_y = NAN;
- strcpy(new->name, name);
-
- return new;
-}
-
-
-struct panel *find_panel_by_name(struct detector *det, const char *name)
-{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
- if ( strcmp(det->panels[i].name, name) == 0 ) {
- return &det->panels[i];
- }
- }
-
- return NULL;
-}
-
-
-static struct badregion *find_bad_region_by_name(struct detector *det,
- const char *name)
-{
- int i;
-
- for ( i=0; i<det->n_bad; i++ ) {
- if ( strcmp(det->bad[i].name, name) == 0 ) {
- return &det->bad[i];
- }
- }
-
- return NULL;
-}
-
-
-static char *find_or_add_rg(struct detector *det, const char *name)
-{
- int i;
- char **new;
- char *tmp;
-
- for ( i=0; i<det->num_rigid_groups; i++ ) {
-
- if ( strcmp(det->rigid_groups[i], name) == 0 ) {
- return det->rigid_groups[i];
- }
-
- }
-
- new = realloc(det->rigid_groups,
- (1+det->num_rigid_groups)*sizeof(char *));
- if ( new == NULL ) return NULL;
-
- det->rigid_groups = new;
-
- tmp = strdup(name);
- det->rigid_groups[det->num_rigid_groups] = tmp;
-
- det->num_rigid_groups++;
-
- return tmp;
-}
-
-
-static int parse_field_for_panel(struct panel *panel, const char *key,
- const char *val, struct detector *det)
-{
- int reject = 0;
-
- if ( strcmp(key, "min_fs") == 0 ) {
- panel->min_fs = atof(val);
- } else if ( strcmp(key, "max_fs") == 0 ) {
- panel->max_fs = atof(val);
- } else if ( strcmp(key, "min_ss") == 0 ) {
- panel->min_ss = atof(val);
- } else if ( strcmp(key, "max_ss") == 0 ) {
- panel->max_ss = atof(val);
- } else if ( strcmp(key, "corner_x") == 0 ) {
- panel->cnx = atof(val);
- } else if ( strcmp(key, "corner_y") == 0 ) {
- panel->cny = atof(val);
- } else if ( strcmp(key, "rigid_group") == 0 ) {
- panel->rigid_group = find_or_add_rg(det, val);
- } else if ( strcmp(key, "clen") == 0 ) {
-
- char *end;
- double v = strtod(val, &end);
- if ( end == val ) {
- /* This means "fill in later" */
- panel->clen = -1.0;
- panel->clen_from = strdup(val);
- } else {
- panel->clen = v;
- panel->clen_from = NULL;
- }
-
- } else if ( strcmp(key, "coffset") == 0) {
- panel->coffset = atof(val);
- } else if ( strcmp(key, "res") == 0 ) {
- panel->res = atof(val);
- } else if ( strcmp(key, "peak_sep") == 0 ) {
- panel->peak_sep = atof(val);
- } else if ( strcmp(key, "integr_radius") == 0 ) {
- panel->integr_radius = atof(val);
- } else if ( strcmp(key, "badrow_direction") == 0 ) {
- panel->badrow = val[0]; /* First character only */
- if ( (panel->badrow != 'x') && (panel->badrow != 'y')
- && (panel->badrow != 'f') && (panel->badrow != 's')
- && (panel->badrow != '-') ) {
- ERROR("badrow_direction must be x, y, f, s or '-'\n");
- ERROR("Assuming '-'\n.");
- panel->badrow = '-';
- }
- if ( panel->badrow == 'x' ) panel->badrow = 'f';
- if ( panel->badrow == 'y' ) panel->badrow = 's';
- } else if ( strcmp(key, "no_index") == 0 ) {
- panel->no_index = atob(val);
- } else if ( strcmp(key, "fs") == 0 ) {
- if ( dir_conv(val, &panel->fsx, &panel->fsy) != 0 ) {
- ERROR("Invalid fast scan direction '%s'\n", val);
- reject = 1;
- }
- } else if ( strcmp(key, "ss") == 0 ) {
- if ( dir_conv(val, &panel->ssx, &panel->ssy) != 0 ) {
- ERROR("Invalid slow scan direction '%s'\n", val);
- reject = 1;
- }
- } else {
- ERROR("Unrecognised field '%s'\n", key);
- }
-
- return reject;
-}
-
-
-static int parse_field_bad(struct badregion *panel, const char *key,
- const char *val)
-{
- int reject = 0;
-
- if ( strcmp(key, "min_x") == 0 ) {
- panel->min_x = atof(val);
- } else if ( strcmp(key, "max_x") == 0 ) {
- panel->max_x = atof(val);
- } else if ( strcmp(key, "min_y") == 0 ) {
- panel->min_y = atof(val);
- } else if ( strcmp(key, "max_y") == 0 ) {
- panel->max_y = atof(val);
- } else {
- ERROR("Unrecognised field '%s'\n", key);
- }
-
- return reject;
-}
-
-
-static void parse_toplevel(struct detector *det, const char *key,
- const char *val)
-{
- if ( strcmp(key, "mask") == 0 ) {
-
- det->mask = strdup(val);
-
- } else if ( strcmp(key, "mask_bad") == 0 ) {
-
- char *end;
- double v = strtod(val, &end);
-
- if ( end != val ) {
- det->mask_bad = v;
- }
-
- } else if ( strcmp(key, "mask_good") == 0 ) {
-
- char *end;
- double v = strtod(val, &end);
-
- if ( end != val ) {
- det->mask_good = v;
- }
-
- } else if ( strcmp(key, "peak_sep") == 0 ) {
- det->defaults.peak_sep = atof(val);
- } else if ( strcmp(key, "integr_radius") == 0 ) {
- det->defaults.integr_radius = atof(val);
- } else if ( parse_field_for_panel(&det->defaults, key, val, det) ) {
- ERROR("Unrecognised top level field '%s'\n", key);
- }
-}
-
-
-struct detector *get_detector_geometry(const char *filename)
-{
- FILE *fh;
- struct detector *det;
- char *rval;
- char **bits;
- int i;
- int reject = 0;
- int x, y, max_fs, max_ss;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) return NULL;
-
- det = calloc(1, sizeof(struct detector));
- if ( det == NULL ) {
- fclose(fh);
- return NULL;
- }
-
- det->n_panels = 0;
- det->panels = NULL;
- det->n_bad = 0;
- det->bad = NULL;
- det->mask_good = 0;
- det->mask_bad = 0;
- det->mask = NULL;
- det->num_rigid_groups = 0;
- det->rigid_groups = NULL;
-
- /* The default defaults... */
- det->defaults.min_fs = -1;
- det->defaults.min_ss = -1;
- det->defaults.max_fs = -1;
- det->defaults.max_ss = -1;
- det->defaults.cnx = NAN;
- det->defaults.cny = NAN;
- det->defaults.clen = -1.0;
- det->defaults.coffset = 0.0;
- det->defaults.res = -1.0;
- det->defaults.badrow = '-';
- det->defaults.no_index = 0;
- det->defaults.peak_sep = 50.0;
- det->defaults.integr_radius = 3.0;
- det->defaults.fsx = 1.0;
- det->defaults.fsy = 0.0;
- det->defaults.ssx = 0.0;
- det->defaults.ssy = 1.0;
- det->defaults.rigid_group = NULL;
- strncpy(det->defaults.name, "", 1023);
-
- do {
-
- int n1, n2;
- char **path;
- char line[1024];
- struct badregion *badregion = NULL;
- struct panel *panel = NULL;
- char *key;
- char wholeval[1024];
-
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) break;
- chomp(line);
-
- if ( line[0] == ';' ) continue;
-
- n1 = assplode(line, " \t", &bits, ASSPLODE_NONE);
- if ( n1 < 3 ) {
- for ( i=0; i<n1; i++ ) free(bits[i]);
- free(bits);
- continue;
- }
-
- /* Stitch the pieces of the "value" back together */
- wholeval[0] = '\0'; /* Empty string */
- for ( i=2; i<n1; i++ ) {
- if ( bits[i][0] == ';' ) break; /* Stop on comment */
- strncat(wholeval, bits[i], 1023);
- }
-
- if ( bits[1][0] != '=' ) {
- for ( i=0; i<n1; i++ ) free(bits[i]);
- free(bits);
- continue;
- }
-
- n2 = assplode(bits[0], "/\\.", &path, ASSPLODE_NONE);
- if ( n2 < 2 ) {
- /* This was a top-level option, not handled above. */
- parse_toplevel(det, bits[0], bits[2]);
- for ( i=0; i<n1; i++ ) free(bits[i]);
- free(bits);
- for ( i=0; i<n2; i++ ) free(path[i]);
- free(path);
- continue;
- }
-
- if ( strncmp(path[0], "bad", 3) == 0 ) {
- badregion = find_bad_region_by_name(det, path[0]);
- if ( badregion == NULL ) {
- badregion = new_bad_region(det, path[0]);
- }
- } else {
- panel = find_panel_by_name(det, path[0]);
- if ( panel == NULL ) {
- panel = new_panel(det, path[0]);
- }
- }
-
- key = path[1];
-
- if ( panel != NULL ) {
- if ( parse_field_for_panel(panel, path[1],
- wholeval, det) )
- {
- reject = 1;
- }
- } else {
- if ( parse_field_bad(badregion, path[1], wholeval) ) {
- reject = 1;
- }
- }
-
- for ( i=0; i<n1; i++ ) free(bits[i]);
- for ( i=0; i<n2; i++ ) free(path[i]);
- free(bits);
- free(path);
-
- } while ( rval != NULL );
-
- if ( det->n_panels == -1 ) {
- ERROR("No panel descriptions in geometry file.\n");
- fclose(fh);
- free(det);
- return NULL;
- }
-
- max_fs = 0;
- max_ss = 0;
- for ( i=0; i<det->n_panels; i++ ) {
-
- if ( det->panels[i].min_fs < 0 ) {
- ERROR("Please specify the minimum FS coordinate for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( det->panels[i].max_fs < 0 ) {
- ERROR("Please specify the maximum FS coordinate for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( det->panels[i].min_ss < 0 ) {
- ERROR("Please specify the minimum SS coordinate for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( det->panels[i].max_ss < 0 ) {
- ERROR("Please specify the maximum SS coordinate for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( isnan(det->panels[i].cnx) ) {
- ERROR("Please specify the corner X coordinate for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( isnan(det->panels[i].cny) ) {
- ERROR("Please specify the corner Y coordinate for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( (det->panels[i].clen < 0.0)
- && (det->panels[i].clen_from == NULL) ) {
- ERROR("Please specify the camera length for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- if ( det->panels[i].res < 0 ) {
- ERROR("Please specify the resolution for"
- " panel %s\n", det->panels[i].name);
- reject = 1;
- }
- /* It's OK if the badrow direction is '0' */
- /* It's not a problem if "no_index" is still zero */
- /* The default peak_sep is OK (maybe) */
- /* The default transformation matrix is at least valid */
-
- if ( det->panels[i].max_fs > max_fs ) {
- max_fs = det->panels[i].max_fs;
- }
- if ( det->panels[i].max_ss > max_ss ) {
- max_ss = det->panels[i].max_ss;
- }
-
- }
-
- for ( i=0; i<det->n_bad; i++ ) {
-
- if ( isnan(det->bad[i].min_x) ) {
- ERROR("Please specify the minimum x coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- if ( isnan(det->bad[i].min_y) ) {
- ERROR("Please specify the minimum y coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- if ( isnan(det->bad[i].max_x) ) {
- ERROR("Please specify the maximum x coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- if ( isnan(det->bad[i].max_y) ) {
- ERROR("Please specify the maximum y coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- }
-
- for ( x=0; x<=max_fs; x++ ) {
- for ( y=0; y<=max_ss; y++ ) {
- if ( find_panel(det, x, y) == NULL ) {
- ERROR("Detector geometry invalid: contains gaps.\n");
- reject = 1;
- goto out;
- }
- }
- }
-
-out:
- det->max_fs = max_fs;
- det->max_ss = max_ss;
-
- /* Calculate matrix inverse */
- for ( i=0; i<det->n_panels; i++ ) {
-
- struct panel *p;
- double d;
-
- p = &det->panels[i];
-
- if ( p->fsx*p->ssy == p->ssx*p->fsy ) {
- ERROR("Panel %i transformation singular.\n", i);
- reject = 1;
- }
-
- d = (double)p->fsx*p->ssy - p->ssx*p->fsy;
- p->xfs = p->ssy / d;
- p->yfs = -p->ssx / d;
- p->xss = -p->fsy / d;
- p->yss = p->fsx / d;
-
- }
-
- if ( reject ) return NULL;
-
- fclose(fh);
-
- return det;
-}
-
-
-void free_detector_geometry(struct detector *det)
-{
- int i;
-
- for ( i=0; i<det->num_rigid_groups; i++ ) {
- free(det->rigid_groups[i]);
- }
- free(det->rigid_groups);
-
- free(det->panels);
- free(det->bad);
- free(det->mask);
- free(det);
-}
-
-
-struct detector *copy_geom(const struct detector *in)
-{
- struct detector *out;
- int i;
-
- out = malloc(sizeof(struct detector));
- memcpy(out, in, sizeof(struct detector));
-
- if ( in->mask != NULL ) {
- out->mask = strdup(in->mask);
- } else {
- out->mask = NULL; /* = in->mask */
- }
-
- out->panels = malloc(out->n_panels * sizeof(struct panel));
- memcpy(out->panels, in->panels, out->n_panels * sizeof(struct panel));
-
- out->bad = malloc(out->n_bad * sizeof(struct badregion));
- memcpy(out->bad, in->bad, out->n_bad * sizeof(struct badregion));
-
- if ( in->rigid_groups != NULL ) {
-
- out->rigid_groups = malloc(out->num_rigid_groups*sizeof(char *));
- memcpy(out->rigid_groups, in->rigid_groups,
- out->num_rigid_groups*sizeof(char *));
-
- for ( i=0; i<in->num_rigid_groups; i++ ) {
- out->rigid_groups[i] = strdup(in->rigid_groups[i]);
- }
-
- }
-
- for ( i=0; i<out->n_panels; i++ ) {
-
- struct panel *p;
-
- p = &out->panels[i];
-
- if ( p->clen_from != NULL ) {
- /* Make a copy of the clen_from fields unique to this
- * copy of the structure. */
- p->clen_from = strdup(p->clen_from);
- }
-
- }
-
- for ( i=0; i<in->num_rigid_groups; i++ ) {
-
- int j;
- char *rg = in->rigid_groups[i];
- char *rgn = out->rigid_groups[i];
-
- for ( j=0; j<in->n_panels; j++ ) {
-
- if ( in->panels[j].rigid_group == rg ) {
- out->panels[j].rigid_group = rgn;
- }
-
- }
-
- }
-
- return out;
-}
-
-
-struct detector *simple_geometry(const struct image *image)
-{
- struct detector *geom;
-
- geom = calloc(1, sizeof(struct detector));
-
- geom->n_panels = 1;
- geom->panels = calloc(1, sizeof(struct panel));
-
- geom->panels[0].min_fs = 0;
- geom->panels[0].max_fs = image->width-1;
- geom->panels[0].min_ss = 0;
- geom->panels[0].max_ss = image->height-1;
- geom->panels[0].cnx = -image->width / 2.0;
- geom->panels[0].cny = -image->height / 2.0;
- geom->panels[0].rigid_group = NULL;
-
- geom->panels[0].fsx = 1;
- geom->panels[0].fsy = 0;
- geom->panels[0].ssx = 0;
- geom->panels[0].ssy = 1;
-
- geom->panels[0].xfs = 1;
- geom->panels[0].xss = 0;
- geom->panels[0].yfs = 0;
- geom->panels[0].yss = 1;
-
- return geom;
-}
-
-
-int reverse_2d_mapping(double x, double y, double *pfs, double *pss,
- struct detector *det)
-{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
-
- struct panel *p = &det->panels[i];
- double cx, cy, fs, ss;
-
- /* Get position relative to corner */
- cx = x - p->cnx;
- cy = y - p->cny;
-
- /* Reverse the transformation matrix */
- fs = cx*p->xfs + cy*p->yfs;
- ss = cx*p->xss + cy*p->yss;
-
- /* In range? */
- if ( fs < 0 ) continue;
- if ( ss < 0 ) continue;
- if ( fs > (p->max_fs-p->min_fs+1) ) continue;
- if ( ss > (p->max_ss-p->min_ss+1) ) continue;
-
- *pfs = fs + p->min_fs;
- *pss = ss + p->min_ss;
- return 0;
-
- }
-
- return 1;
-}
-
-
-static void check_extents(struct panel p, double *min_x, double *min_y,
- double *max_x, double *max_y, double fs, double ss)
-{
- double xs, ys, rx, ry;
-
- xs = fs*p.fsx + ss*p.ssx;
- ys = fs*p.fsy + ss*p.ssy;
-
- rx = xs + p.cnx;
- ry = ys + p.cny;
-
- if ( rx > *max_x ) *max_x = rx;
- if ( ry > *max_y ) *max_y = ry;
- if ( rx < *min_x ) *min_x = rx;
- if ( ry < *min_y ) *min_y = ry;
-}
-
-
-double largest_q(struct image *image)
-{
- int fs, ss;
- double ttm = 0.0;
- double qmax = 0.0;
-
- for ( fs=0; fs<image->width; fs++ ) {
- for ( ss=0; ss<image->height; ss++ ) {
-
- struct rvec q;
- double tt;
-
- q = get_q(image, fs, ss, &tt, 1.0/image->lambda);
-
- if ( tt > ttm ) {
- qmax = modulus(q.u, q.v, q.w);
- ttm = tt;
- }
-
- }
- }
-
- return qmax;
-}
-
-
-double smallest_q(struct image *image)
-{
- int fs, ss;
- double ttm = +INFINITY;
- double qmin = +INFINITY;
- for ( fs=0; fs<image->width; fs++ ) {
- for ( ss=0; ss<image->height; ss++ ) {
-
- struct rvec q;
- double tt;
-
- q = get_q(image, fs, ss, &tt, 1.0/image->lambda);
-
- if ( tt < ttm ) {
- qmin = modulus(q.u, q.v, q.w);
- ttm = tt;
- }
-
- }
- }
-
- return qmin;
-}
-
-
-void get_pixel_extents(struct detector *det,
- double *min_x, double *min_y,
- double *max_x, double *max_y)
-{
- int i;
-
- *min_x = 0.0;
- *max_x = 0.0;
- *min_y = 0.0;
- *max_y = 0.0;
-
- /* To determine the maximum extents of the detector, put all four
- * corners of each panel through the transformations and watch for the
- * biggest */
-
- for ( i=0; i<det->n_panels; i++ ) {
-
- check_extents(det->panels[i], min_x, min_y, max_x, max_y,
- 0.0,
- 0.0);
-
- check_extents(det->panels[i], min_x, min_y, max_x, max_y,
- 0.0,
- det->panels[i].max_ss-det->panels[i].min_ss+1);
-
- check_extents(det->panels[i], min_x, min_y, max_x, max_y,
- det->panels[i].max_fs-det->panels[i].min_fs+1,
- 0.0);
-
- check_extents(det->panels[i], min_x, min_y, max_x, max_y,
- det->panels[i].max_fs-det->panels[i].min_fs+1,
- det->panels[i].max_ss-det->panels[i].min_ss+1);
-
-
- }
-}
-
-
-int write_detector_geometry(const char *filename, struct detector *det)
-{
- struct panel *p;
- int pi;
- FILE *fh;
-
- if ( filename == NULL ) return 2;
- if ( det->n_panels < 1 ) return 3;
-
- fh = fopen(filename, "w");
- if ( fh == NULL ) return 1;
-
- for ( pi=0; pi<det->n_panels; pi++) {
-
- p = &(det->panels[pi]);
-
- if ( p == NULL ) return 4;
-
- if ( pi > 0 ) fprintf(fh, "\n");
-
- fprintf(fh, "%s/min_fs = %d\n", p->name, p->min_fs);
- fprintf(fh, "%s/min_ss = %d\n", p->name, p->min_ss);
- fprintf(fh, "%s/max_fs = %d\n", p->name, p->max_fs);
- fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss);
- fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow);
- fprintf(fh, "%s/res = %g\n", p->name, p->res);
- fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep);
- fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from);
- fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy);
- fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy);
- fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx);
- fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny);
-
- if ( p->no_index ) {
- fprintf(fh, "%s/no_index = 1\n", p->name);
- } /* else don't clutter up the file */
-
- if ( p->rigid_group != NULL ) {
- fprintf(fh, "%s/rigid_group = %s\n",
- p->name, p->rigid_group);
- }
-
- }
- fclose(fh);
-
- return 0;
-}
diff --git a/src/detector.h b/src/detector.h
deleted file mode 100644
index dd5dccec..00000000
--- a/src/detector.h
+++ /dev/null
@@ -1,131 +0,0 @@
-/*
- * detector.h
- *
- * Detector properties
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- * (c) 2011 Rick Kirian <rkirian@asu.edu>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DETECTOR_H
-#define DETECTOR_H
-
-struct image;
-struct hdfile;
-
-#include "hdf5-file.h"
-#include "image.h"
-
-
-struct panel
-{
- char name[1024]; /* Name for this panel */
-
- int min_fs; /* Smallest FS value considered to be in the panel */
- int max_fs; /* Largest FS value considered to be in this panel */
- int min_ss; /* ... and so on */
- int max_ss;
- double cnx; /* Location of corner (min_fs,min_ss) in pixels */
- double cny;
- double coffset;
- double clen; /* Camera length in metres */
- char *clen_from;
- double res; /* Resolution in pixels per metre */
- char badrow; /* 'x' or 'y' */
- int no_index; /* Don't index peaks in this panel if non-zero */
- double peak_sep; /* Characteristic peak separation */
- double integr_radius; /* Peak integration radius */
- char *rigid_group; /* Rigid group, or -1 for none */
-
- double fsx;
- double fsy;
- double ssx;
- double ssy;
-
- double xfs;
- double yfs;
- double xss;
- double yss;
-};
-
-
-struct badregion
-{
- char name[1024];
- double min_x;
- double max_x;
- double min_y;
- double max_y;
-};
-
-
-struct detector
-{
- struct panel *panels;
- int n_panels;
-
- int max_fs;
- int max_ss; /* Size of overall array needed, minus 1 */
-
- struct badregion *bad;
- int n_bad;
-
- char *mask;
- unsigned int mask_bad;
- unsigned int mask_good;
-
- char **rigid_groups;
- int num_rigid_groups;
-
- struct panel defaults;
-};
-
-
-extern struct rvec get_q(struct image *image, double fs, double ss,
- double *ttp, double k);
-
-extern struct rvec get_q_for_panel(struct panel *p, double fs, double ss,
- double *ttp, double k);
-
-extern double get_tt(struct image *image, double xs, double ys);
-
-extern int in_bad_region(struct detector *det, double fs, double ss);
-
-extern void record_image(struct image *image, int do_poisson);
-
-extern struct panel *find_panel(struct detector *det, double fs, double ss);
-
-extern int find_panel_number(struct detector *det, int fs, int ss);
-
-extern struct detector *get_detector_geometry(const char *filename);
-
-extern void free_detector_geometry(struct detector *det);
-
-extern struct detector *simple_geometry(const struct image *image);
-
-extern void get_pixel_extents(struct detector *det,
- double *min_x, double *min_y,
- double *max_x, double *max_y);
-
-extern void fill_in_values(struct detector *det, struct hdfile *f);
-
-extern struct detector *copy_geom(const struct detector *in);
-
-extern int reverse_2d_mapping(double x, double y, double *pfs, double *pss,
- struct detector *det);
-
-extern double largest_q(struct image *image);
-
-extern double smallest_q(struct image *image);
-
-extern struct panel *find_panel_by_name(struct detector *det, const char *name);
-
-
-#endif /* DETECTOR_H */
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
deleted file mode 100644
index 355e97f1..00000000
--- a/src/hdf5-file.c
+++ /dev/null
@@ -1,817 +0,0 @@
-/*
- * hdf5-file.c
- *
- * Read/write HDF5 data files
- *
- * (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 <stdio.h>
-#include <stdint.h>
-#include <hdf5.h>
-#include <assert.h>
-
-#include "image.h"
-#include "hdf5-file.h"
-#include "utils.h"
-
-
-struct hdfile {
-
- const char *path; /* Current data path */
-
- size_t nx; /* Image width */
- size_t ny; /* Image height */
-
- hid_t fh; /* HDF file handle */
- hid_t dh; /* Dataset handle */
-
- int data_open; /* True if dh is initialised */
-};
-
-
-struct hdfile *hdfile_open(const char *filename)
-{
- struct hdfile *f;
-
- f = malloc(sizeof(struct hdfile));
- if ( f == NULL ) return NULL;
-
- /* Please stop spamming my terminal */
- H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
-
- f->fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
- if ( f->fh < 0 ) {
- ERROR("Couldn't open file: %s\n", filename);
- free(f);
- return NULL;
- }
-
- f->data_open = 0;
-
- return f;
-}
-
-
-int hdfile_set_image(struct hdfile *f, const char *path)
-{
- hsize_t size[2];
- hsize_t max_size[2];
- hid_t sh;
-
- f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT);
- if ( f->dh < 0 ) {
- ERROR("Couldn't open dataset\n");
- return -1;
- }
- f->data_open = 1;
-
- sh = H5Dget_space(f->dh);
- if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
- ERROR("Dataset is not two-dimensional\n");
- return -1;
- }
- H5Sget_simple_extent_dims(sh, size, max_size);
- H5Sclose(sh);
-
- f->nx = size[0];
- f->ny = size[1];
-
- return 0;
-}
-
-
-int get_peaks(struct image *image, struct hdfile *f, const char *p)
-{
- hid_t dh, sh;
- hsize_t size[2];
- hsize_t max_size[2];
- int i;
- float *buf;
- herr_t r;
- int tw;
-
- dh = H5Dopen2(f->fh, p, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Peak list (%s) not found.\n", p);
- return 1;
- }
-
- sh = H5Dget_space(dh);
- if ( sh < 0 ) {
- H5Dclose(dh);
- ERROR("Couldn't get dataspace for peak list.\n");
- return 1;
- }
-
- if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
- ERROR("Peak list has the wrong dimensionality (%i).\n",
- H5Sget_simple_extent_ndims(sh));
- H5Sclose(sh);
- H5Dclose(dh);
- return 1;
- }
-
- H5Sget_simple_extent_dims(sh, size, max_size);
-
- tw = size[1];
- if ( (tw != 3) && (tw != 4) ) {
- H5Sclose(sh);
- H5Dclose(dh);
- ERROR("Peak list has the wrong dimensions.\n");
- return 1;
- }
-
- buf = malloc(sizeof(float)*size[0]*size[1]);
- if ( buf == NULL ) {
- H5Sclose(sh);
- H5Dclose(dh);
- ERROR("Couldn't reserve memory for the peak list.\n");
- return 1;
- }
- r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
- if ( r < 0 ) {
- ERROR("Couldn't read peak list.\n");
- free(buf);
- return 1;
- }
-
- if ( image->features != NULL ) {
- image_feature_list_free(image->features);
- }
- image->features = image_feature_list_new();
-
- for ( i=0; i<size[0]; i++ ) {
-
- float fs, ss, val;
- struct panel *p;
-
- fs = buf[tw*i+0];
- ss = buf[tw*i+1];
- val = buf[tw*i+2];
-
- p = find_panel(image->det, fs, ss);
- if ( p == NULL ) continue;
- if ( p->no_index ) continue;
-
- image_add_feature(image->features, fs, ss, image, val, NULL);
-
- }
-
- free(buf);
- H5Sclose(sh);
- H5Dclose(dh);
-
- return 0;
-}
-
-
-static void cleanup(hid_t fh)
-{
- int n_ids, i;
- hid_t ids[256];
-
- n_ids = H5Fget_obj_ids(fh, H5F_OBJ_ALL, 256, ids);
- for ( i=0; i<n_ids; i++ ) {
-
- hid_t id;
- H5I_type_t type;
-
- id = ids[i];
- type = H5Iget_type(id);
-
- if ( type == H5I_GROUP ) H5Gclose(id);
- if ( type == H5I_DATASET ) H5Dclose(id);
- if ( type == H5I_DATATYPE ) H5Tclose(id);
- if ( type == H5I_DATASPACE ) H5Sclose(id);
- if ( type == H5I_ATTR ) H5Aclose(id);
-
- }
-}
-
-
-void hdfile_close(struct hdfile *f)
-{
- if ( f->data_open ) {
- H5Dclose(f->dh);
- }
- cleanup(f->fh);
-
- H5Fclose(f->fh);
- free(f);
-}
-
-
-int hdf5_write(const char *filename, const void *data,
- int width, int height, int type)
-{
- hid_t fh, gh, sh, dh; /* File, group, dataspace and data handles */
- hid_t ph; /* Property list */
- herr_t r;
- hsize_t size[2];
- hsize_t max_size[2];
-
- fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
- if ( fh < 0 ) {
- ERROR("Couldn't create file: %s\n", filename);
- return 1;
- }
-
- gh = H5Gcreate2(fh, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if ( gh < 0 ) {
- ERROR("Couldn't create group\n");
- H5Fclose(fh);
- return 1;
- }
-
- /* Note the "swap" here, according to section 3.2.5,
- * "C versus Fortran Dataspaces", of the HDF5 user's guide. */
- size[0] = height;
- size[1] = width;
- max_size[0] = height;
- max_size[1] = width;
- sh = H5Screate_simple(2, size, max_size);
-
- /* Set compression */
- ph = H5Pcreate(H5P_DATASET_CREATE);
- H5Pset_chunk(ph, 2, size);
- H5Pset_deflate(ph, 3);
-
- dh = H5Dcreate2(gh, "data", type, sh,
- H5P_DEFAULT, ph, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't create dataset\n");
- H5Fclose(fh);
- return 1;
- }
-
- /* Muppet check */
- H5Sget_simple_extent_dims(sh, size, max_size);
-
- r = H5Dwrite(dh, type, H5S_ALL,
- H5S_ALL, H5P_DEFAULT, data);
- if ( r < 0 ) {
- ERROR("Couldn't write data\n");
- H5Dclose(dh);
- H5Fclose(fh);
- return 1;
- }
-
- H5Pclose(ph);
- H5Gclose(gh);
- H5Dclose(dh);
- H5Fclose(fh);
-
- return 0;
-}
-
-
-static double get_wavelength(struct hdfile *f)
-{
- herr_t r;
- hid_t dh;
- double lambda;
- int nm = 1;
-
- dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_nm", H5P_DEFAULT);
- if ( dh < 0 ) {
- dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_A", H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't get photon wavelength from HDF5 file.\n");
- return -1.0;
- }
- nm = 0;
- }
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &lambda);
- H5Dclose(dh);
-
- if ( r < 0 ) return -1.0;
- if ( isnan(lambda) ) return -1.0;
-
- /* Convert nm -> m */
- if ( nm ) return lambda / 1.0e9;
- return lambda / 1.0e10;
-}
-
-
-static double get_i0(struct hdfile *f)
-{
- herr_t r;
- hid_t dh;
- double i0;
-
- dh = H5Dopen2(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT);
- if ( dh < 0 ) return -1.0;
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &i0);
- H5Dclose(dh);
- if ( r < 0 ) return -1.0;
-
- return i0;
-}
-
-
-static void debodge_saturation(struct hdfile *f, struct image *image)
-{
- hid_t dh, sh;
- hsize_t size[2];
- hsize_t max_size[2];
- int i;
- float *buf;
- herr_t r;
-
- dh = H5Dopen2(f->fh, "/processing/hitfinder/peakinfo_saturated",
- H5P_DEFAULT);
-
- if ( dh < 0 ) {
- /* This isn't an error */
- return;
- }
-
- sh = H5Dget_space(dh);
- if ( sh < 0 ) {
- H5Dclose(dh);
- ERROR("Couldn't get dataspace for saturation table.\n");
- return;
- }
-
- if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
- H5Sclose(sh);
- H5Dclose(dh);
- return;
- }
-
- H5Sget_simple_extent_dims(sh, size, max_size);
-
- if ( size[1] != 3 ) {
- H5Sclose(sh);
- H5Dclose(dh);
- ERROR("Saturation table has the wrong dimensions.\n");
- return;
- }
-
- buf = malloc(sizeof(float)*size[0]*size[1]);
- if ( buf == NULL ) {
- H5Sclose(sh);
- H5Dclose(dh);
- ERROR("Couldn't reserve memory for saturation table.\n");
- return;
- }
- r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
- if ( r < 0 ) {
- ERROR("Couldn't read saturation table.\n");
- free(buf);
- return;
- }
-
- for ( i=0; i<size[0]; i++ ) {
-
- unsigned int x, y;
- float val;
-
- x = buf[3*i+0];
- y = buf[3*i+1];
- val = buf[3*i+2];
-
- image->data[x+image->width*y] = val / 5.0;
- image->data[x+1+image->width*y] = val / 5.0;
- image->data[x-1+image->width*y] = val / 5.0;
- image->data[x+image->width*(y+1)] = val / 5.0;
- image->data[x+image->width*(y-1)] = val / 5.0;
-
- }
-
- free(buf);
- H5Sclose(sh);
- H5Dclose(dh);
-}
-
-
-int hdf5_read(struct hdfile *f, struct image *image, int satcorr)
-{
- herr_t r;
- float *buf;
- uint16_t *flags;
- hid_t mask_dh;
-
- /* Note the "swap" here, according to section 3.2.5,
- * "C versus Fortran Dataspaces", of the HDF5 user's guide. */
- image->width = f->ny;
- image->height = f->nx;
-
- buf = malloc(sizeof(float)*f->nx*f->ny);
-
- r = H5Dread(f->dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, buf);
- if ( r < 0 ) {
- ERROR("Couldn't read data\n");
- free(buf);
- return 1;
- }
- image->data = buf;
-
- if ( (image->det != NULL) && (image->det->mask != NULL) ) {
-
- mask_dh = H5Dopen2(f->fh, image->det->mask, H5P_DEFAULT);
- if ( mask_dh <= 0 ) {
- ERROR("Couldn't open flags\n");
- image->flags = NULL;
- } else {
- flags = malloc(sizeof(uint16_t)*f->nx*f->ny);
- r = H5Dread(mask_dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, flags);
- if ( r < 0 ) {
- ERROR("Couldn't read flags\n");
- free(flags);
- image->flags = NULL;
- } else {
- image->flags = flags;
- }
- H5Dclose(mask_dh);
- }
-
- }
-
- /* Read wavelength from file */
- image->lambda = get_wavelength(f);
-
- image->i0 = get_i0(f);
- if ( image->i0 < 0.0 ) {
- ERROR("Couldn't read incident intensity - using 1.0.\n");
- image->i0 = 1.0;
- image->i0_available = 0;
- } else {
- image->i0_available = 1;
- }
-
- if ( satcorr ) debodge_saturation(f, image);
-
- return 0;
-}
-
-
-static int looks_like_image(hid_t h)
-{
- hid_t sh;
- hsize_t size[2];
- hsize_t max_size[2];
-
- sh = H5Dget_space(h);
- if ( sh < 0 ) return 0;
-
- if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
- return 0;
- }
-
- H5Sget_simple_extent_dims(sh, size, max_size);
-
- if ( ( size[0] > 64 ) && ( size[1] > 64 ) ) return 1;
-
- return 0;
-}
-
-
-double get_value(struct hdfile *f, const char *name)
-{
- hid_t dh;
- hid_t sh;
- hsize_t size;
- hsize_t max_size;
- hid_t type;
- hid_t class;
- herr_t r;
- double buf;
-
- dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Couldn't open data\n");
- return 0.0;
- }
-
- type = H5Dget_type(dh);
- class = H5Tget_class(type);
-
- if ( class != H5T_FLOAT ) {
- ERROR("Not a floating point value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 0.0;
- }
-
- sh = H5Dget_space(dh);
- if ( H5Sget_simple_extent_ndims(sh) != 1 ) {
- ERROR("Not a scalar value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 0.0;
- }
-
- H5Sget_simple_extent_dims(sh, &size, &max_size);
- if ( size != 1 ) {
- ERROR("Not a scalar value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 0.0;
- }
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &buf);
- if ( r < 0 ) {
- ERROR("Couldn't read value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 0.0;
- }
-
- return buf;
-}
-
-
-struct copy_hdf5_field
-{
- char **fields;
- int n_fields;
- int max_fields;
-};
-
-
-struct copy_hdf5_field *new_copy_hdf5_field_list()
-{
- struct copy_hdf5_field *n;
-
- n = calloc(1, sizeof(struct copy_hdf5_field));
- if ( n == NULL ) return NULL;
-
- n->max_fields = 32;
- n->fields = malloc(n->max_fields*sizeof(char *));
- if ( n->fields == NULL ) {
- free(n);
- return NULL;
- }
-
- return n;
-}
-
-
-void free_copy_hdf5_field_list(struct copy_hdf5_field *n)
-{
- int i;
- for ( i=0; i<n->n_fields; i++ ) {
- free(n->fields[i]);
- }
- free(n->fields);
- free(n);
-}
-
-
-void add_copy_hdf5_field(struct copy_hdf5_field *copyme,
- const char *name)
-{
- assert(copyme->n_fields < copyme->max_fields); /* FIXME */
-
- copyme->fields[copyme->n_fields] = strdup(name);
- if ( copyme->fields[copyme->n_fields] == NULL ) {
- ERROR("Failed to add field for copying '%s'\n", name);
- return;
- }
-
- copyme->n_fields++;
-}
-
-
-void copy_hdf5_fields(struct hdfile *f, const struct copy_hdf5_field *copyme,
- FILE *fh)
-{
- int i;
-
- if ( copyme == NULL ) return;
-
- for ( i=0; i<copyme->n_fields; i++ ) {
-
- char *val;
- char *field;
-
- field = copyme->fields[i];
- val = hdfile_get_string_value(f, field);
-
- if ( field[0] == '/' ) {
- fprintf(fh, "hdf5%s = %s\n", field, val);
- } else {
- fprintf(fh, "hdf5/%s = %s\n", field, val);
- }
-
- free(val);
-
- }
-}
-
-
-char *hdfile_get_string_value(struct hdfile *f, const char *name)
-{
- hid_t dh;
- hid_t sh;
- hsize_t size;
- hsize_t max_size;
- hid_t type;
- hid_t class;
-
- dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
- if ( dh < 0 ) return NULL;
-
- type = H5Dget_type(dh);
- class = H5Tget_class(type);
-
- if ( class == H5T_STRING ) {
-
- herr_t r;
- char *tmp;
- hid_t sh;
-
- size = H5Tget_size(type);
- tmp = malloc(size+1);
-
- sh = H5Screate(H5S_SCALAR);
-
- r = H5Dread(dh, type, sh, sh, H5P_DEFAULT, tmp);
- if ( r < 0 ) goto fail;
-
- /* Two possibilities:
- * String is already zero-terminated
- * String is not terminated.
- * Make sure things are done properly... */
- tmp[size] = '\0';
- chomp(tmp);
-
- return tmp;
-
- }
-
- sh = H5Dget_space(dh);
- if ( H5Sget_simple_extent_ndims(sh) != 1 ) goto fail;
-
- H5Sget_simple_extent_dims(sh, &size, &max_size);
- if ( size != 1 ) {
- H5Dclose(dh);
- goto fail;
- }
-
- switch ( class ) {
- case H5T_FLOAT : {
-
- herr_t r;
- double buf;
- char *tmp;
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &buf);
- if ( r < 0 ) goto fail;
-
- tmp = malloc(256);
- snprintf(tmp, 255, "%f", buf);
-
- return tmp;
-
- }
- case H5T_INTEGER : {
-
- herr_t r;
- int buf;
- char *tmp;
-
- r = H5Dread(dh, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &buf);
- if ( r < 0 ) goto fail;
-
- tmp = malloc(256);
- snprintf(tmp, 255, "%d", buf);
-
- return tmp;
- }
- default : {
- goto fail;
- }
- }
-
-fail:
- H5Tclose(type);
- H5Dclose(dh);
- return NULL;
-}
-
-
-char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
- int **p_is_group, int **p_is_image)
-{
- hid_t gh;
- hsize_t num;
- char **res;
- int i;
- int *is_group;
- int *is_image;
- H5G_info_t ginfo;
-
- gh = H5Gopen2(f->fh, parent, H5P_DEFAULT);
- if ( gh < 0 ) {
- *n = 0;
- return NULL;
- }
-
- if ( H5Gget_info(gh, &ginfo) < 0 ) {
- /* Whoopsie */
- *n = 0;
- return NULL;
- }
- num = ginfo.nlinks;
- *n = num;
- if ( num == 0 ) return NULL;
-
- res = malloc(num*sizeof(char *));
- is_image = malloc(num*sizeof(int));
- is_group = malloc(num*sizeof(int));
- *p_is_image = is_image;
- *p_is_group = is_group;
-
- for ( i=0; i<num; i++ ) {
-
- char buf[256];
- hid_t dh;
- H5I_type_t type;
-
- H5Lget_name_by_idx(gh, ".", H5_INDEX_NAME, H5_ITER_NATIVE,
- i, buf, 255, H5P_DEFAULT);
- res[i] = malloc(256);
- if ( strlen(parent) > 1 ) {
- snprintf(res[i], 255, "%s/%s", parent, buf);
- } else {
- snprintf(res[i], 255, "%s%s", parent, buf);
- } /* ick */
-
- is_image[i] = 0;
- is_group[i] = 0;
- dh = H5Oopen(gh, buf, H5P_DEFAULT);
- if ( dh < 0 ) continue;
- type = H5Iget_type(dh);
-
- if ( type == H5I_GROUP ) {
- is_group[i] = 1;
- } else if ( type == H5I_DATASET ) {
- is_image[i] = looks_like_image(dh);
- }
- H5Oclose(dh);
-
- }
-
- return res;
-}
-
-
-int hdfile_set_first_image(struct hdfile *f, const char *group)
-{
- char **names;
- int *is_group;
- int *is_image;
- int n, i, j;
-
- names = hdfile_read_group(f, &n, group, &is_group, &is_image);
- if ( n == 0 ) return 1;
-
- for ( i=0; i<n; i++ ) {
-
- if ( is_image[i] ) {
- hdfile_set_image(f, names[i]);
- for ( j=0; j<n; j++ ) free(names[j]);
- free(is_image);
- free(is_group);
- free(names);
- return 0;
- } else if ( is_group[i] ) {
- if ( !hdfile_set_first_image(f, names[i]) ) {
- for ( j=0; j<n; j++ ) free(names[j]);
- free(is_image);
- free(is_group);
- free(names);
- return 0;
- }
- }
-
- }
-
- for ( j=0; j<n; j++ ) free(names[j]);
- free(is_image);
- free(is_group);
- free(names);
-
- return 1;
-}
diff --git a/src/hdf5-file.h b/src/hdf5-file.h
deleted file mode 100644
index 385e919c..00000000
--- a/src/hdf5-file.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/*
- * hdf5-file.h
- *
- * Read/write HDF5 data files
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef HDF5_H
-#define HDF5_H
-
-#include <stdint.h>
-#include <hdf5.h>
-
-#include "image.h"
-
-struct hdfile;
-
-struct copy_hdf5_field;
-
-
-extern int hdf5_write(const char *filename, const void *data,
- int width, int height, int type);
-
-extern int hdf5_read(struct hdfile *f, struct image *image, int satcorr);
-
-extern struct hdfile *hdfile_open(const char *filename);
-extern int hdfile_set_image(struct hdfile *f, const char *path);
-extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile,
- int binning, int16_t *maxp);
-extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
- int **p_is_group, int **p_is_image);
-extern int hdfile_set_first_image(struct hdfile *f, const char *group);
-extern void hdfile_close(struct hdfile *f);
-
-extern char *hdfile_get_string_value(struct hdfile *f, const char *name);
-extern int get_peaks(struct image *image, struct hdfile *f, const char *p);
-extern double get_value(struct hdfile *f, const char *name);
-
-extern struct copy_hdf5_field *new_copy_hdf5_field_list(void);
-extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f);
-extern void copy_hdf5_fields(struct hdfile *f,
- const struct copy_hdf5_field *copyme, FILE *fh);
-extern void add_copy_hdf5_field(struct copy_hdf5_field *copyme,
- const char *name);
-
-
-#endif /* HDF5_H */
diff --git a/src/image.h b/src/image.h
index 4d4ac2d7..e21ebe85 100644
--- a/src/image.h
+++ b/src/image.h
@@ -21,10 +21,10 @@
#include <complex.h>
#include <sys/types.h>
-#include "utils.h"
-#include "cell.h"
-#include "detector.h"
-#include "reflist.h"
+#include <utils.h>
+#include <cell.h>
+#include <detector.h>
+#include <reflist.h>
#define MAX_CELL_CANDIDATES (32)
diff --git a/src/list_tmp.h b/src/list_tmp.h
deleted file mode 100644
index a524b2f9..00000000
--- a/src/list_tmp.h
+++ /dev/null
@@ -1,106 +0,0 @@
-/*
- * Template for creating indexed 3D lists of a given type, usually indexed
- * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX.
- *
- * These are used, for example, for:
- * - a list of 'double complex' values for storing structure factors,
- * - a list of 'double' values for storing reflection intensities,
- * - a list of 'unsigned int' values for counts of some sort.
- *
- * When LABEL and TYPE are #defined appropriately, including this header
- * defines functions such as:
- * - new_list_<LABEL>(), for creating a new list of the given type,
- * - set_<LABEL>(), for setting a value in a list,
- * - lookup_<LABEL>(), for retrieving values from a list,
- * ..and so on.
- *
- * See src/utils.h for more information.
- *
- */
-
-#include <stdio.h>
-
-#define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
-
-static inline void LABEL(integrate)(TYPE *ref, signed int h,
- signed int k, signed int l,
- TYPE i)
-{
- int idx;
-
- if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
- ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
- ERROR_T("You need to re-configure INDMAX and re-run.\n");
- exit(1);
- }
-
- if ( h < 0 ) h += IDIM;
- if ( k < 0 ) k += IDIM;
- if ( l < 0 ) l += IDIM;
-
- idx = h + (IDIM*k) + (IDIM*IDIM*l);
- ref[idx] += i;
-}
-
-
-static inline void LABEL(set)(TYPE *ref, signed int h,
- signed int k, signed int l,
- TYPE i)
-{
- int idx;
-
- if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
- ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
- ERROR_T("You need to re-configure INDMAX and re-run.\n");
- exit(1);
- }
-
- if ( h < 0 ) h += IDIM;
- if ( k < 0 ) k += IDIM;
- if ( l < 0 ) l += IDIM;
-
- idx = h + (IDIM*k) + (IDIM*IDIM*l);
- ref[idx] = i;
-}
-
-
-static inline TYPE LABEL(lookup)(const TYPE *ref, signed int h,
- signed int k, signed int l)
-{
- int idx;
-
- if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
- ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
- ERROR_T("You need to re-configure INDMAX and re-run.\n");
- exit(1);
- }
-
- if ( h < 0 ) h += IDIM;
- if ( k < 0 ) k += IDIM;
- if ( l < 0 ) l += IDIM;
-
- idx = h + (IDIM*k) + (IDIM*IDIM*l);
- return ref[idx];
-}
-
-
-static inline TYPE *LABEL(new_list)(void)
-{
- TYPE *r;
- size_t r_size;
- r_size = IDIM*IDIM*IDIM*sizeof(TYPE);
- r = malloc(r_size);
- memset(r, 0, r_size);
- return r;
-}
-
-
-static inline void LABEL(zero_list)(TYPE *ref)
-{
- memset(ref, 0, IDIM*IDIM*IDIM*sizeof(TYPE));
-}
-
-
-#undef LABEL
-#undef TYPE
-#undef ERROR_T
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 9f6a4015..c339dba9 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -425,7 +425,7 @@ int main(int argc, char *argv[])
image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
image.div = beam->divergence;
image.bw = beam->bandwidth;
- image.profile_radius = 0.003e9;
+ image.profile_radius = 0.01e9;
image.i0_available = 0;
image.filename = malloc(256);
image.copyme = NULL;
diff --git a/src/reflist.c b/src/reflist.c
deleted file mode 100644
index 18856c67..00000000
--- a/src/reflist.c
+++ /dev/null
@@ -1,1048 +0,0 @@
-/*
- * reflist.c
- *
- * Fast reflection/peak list
- *
- * (c) 2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <assert.h>
-#include <stdio.h>
-#include <pthread.h>
-
-#include "reflist.h"
-#include "utils.h"
-
-/**
- * SECTION:reflist
- * @short_description: The fast reflection list
- * @title: RefList
- * @section_id:
- * @see_also:
- * @include: "reflist.h"
- * @Image:
- *
- * The fast reflection list stores reflections in an RB-tree indexed
- * by the Miller indices h, k and l. Any reflection can be found in a
- * length of time which scales logarithmically with the number of reflections in
- * the list.
- *
- * A RefList can contain any number of reflections, and can store more than
- * one reflection with a given set of indices, for example when two distinct
- * reflections are to be stored according to their asymmetric indices.
- *
- * There are getters and setters which can be used to get and set values for an
- * individual reflection. The reflection list does not calculate any values,
- * only stores what it was given earlier. As such, you will need to carefully
- * examine which fields your prior processing steps have filled in.
- */
-
-
-struct _refldata {
-
- /* Symmetric indices (i.e. the "real" indices) */
- signed int hs;
- signed int ks;
- signed int ls;
-
- /* Partiality and related geometrical stuff */
- double r1; /* First excitation error */
- double r2; /* Second excitation error */
- double p; /* Partiality */
- int clamp1; /* Clamp status for r1 */
- int clamp2; /* Clamp status for r2 */
-
- /* Location in image */
- double fs;
- double ss;
-
- /* The distance from the exact Bragg position to the coordinates
- * given above. */
- double excitation_error;
-
- /* Non-zero if this reflection can be used for scaling */
- int scalable;
-
- /* Non-zero if this reflection should be used as a "guide star" for
- * post refinement */
- int refinable;
-
- /* Intensity */
- double intensity;
- double esd_i;
-
- /* Phase */
- double phase;
- int have_phase;
-
- /* Redundancy */
- int redundancy;
-
- /* User-specified temporary values */
- double temp1;
- double temp2;
-};
-
-
-enum _nodecol {
- RED,
- BLACK
-};
-
-
-struct _reflection {
-
- /* Listy stuff */
- unsigned int serial; /* Unique serial number, key */
- struct _reflection *child[2]; /* Child nodes */
- struct _reflection *next; /* Next and previous in doubly linked */
- struct _reflection *prev; /* list of duplicate reflections */
- enum _nodecol col; /* Colour (red or black) */
-
- /* Payload */
- pthread_mutex_t lock; /* Protects the contents of "data" */
- struct _refldata data;
-};
-
-
-struct _reflist {
-
- struct _reflection *head;
- struct _reflection *tail;
-
-};
-
-
-#define SERIAL(h, k, l) ((((h)+256)<<18) + (((k)+256)<<9) + ((l)+256))
-#define GET_H(serial) ((((serial) & 0xfffc0000)>>18)-256)
-#define GET_K(serial) ((((serial) & 0x0003fe00)>>9)-256)
-#define GET_L(serial) (((serial) & 0x000001ff)-256)
-
-/**************************** Creation / deletion *****************************/
-
-static Reflection *new_node(unsigned int serial)
-{
- Reflection *new;
-
- new = calloc(1, sizeof(struct _reflection));
- new->serial = serial;
- new->next = NULL;
- new->prev = NULL;
- new->child[0] = NULL;
- new->child[1] = NULL;
- new->col = RED;
- pthread_mutex_init(&new->lock, NULL);
-
- return new;
-}
-
-
-/**
- * reflist_new:
- *
- * Creates a new reflection list.
- *
- * Returns: the new reflection list, or NULL on error.
- */
-RefList *reflist_new()
-{
- RefList *new;
-
- new = malloc(sizeof(struct _reflist));
- if ( new == NULL ) return NULL;
-
- new->head = NULL;
-
- return new;
-}
-
-
-/**
- * reflection_new:
- * @h: The h index of the new reflection
- * @k: The k index of the new reflection
- * @l: The l index of the new reflection
- *
- * Creates a new individual reflection. You'll probably want to use
- * add_refl_to_list() at some later point.
- */
-Reflection *reflection_new(signed int h, signed int k, signed int l)
-{
- return new_node(SERIAL(h, k, l));
-}
-
-
-/**
- * reflection_free:
- * @refl: The reflection to free.
- *
- * Destroys an individual reflection.
- */
-void reflection_free(Reflection *refl)
-{
- pthread_mutex_destroy(&refl->lock);
- free(refl);
-}
-
-
-static void recursive_free(Reflection *refl)
-{
- if ( refl->child[0] != NULL ) recursive_free(refl->child[0]);
- if ( refl->child[1] != NULL ) recursive_free(refl->child[1]);
-
- while ( refl != NULL ) {
- Reflection *next = refl->next;
- reflection_free(refl);
- refl = next;
- }
-}
-
-
-/**
- * reflist_free:
- * @list: The reflection list to free.
- *
- * Destroys a reflection list.
- */
-void reflist_free(RefList *list)
-{
- if ( list == NULL ) return;
- if ( list->head != NULL ) {
- recursive_free(list->head);
- } /* else empty list */
- free(list);
-}
-
-
-/********************************** Search ************************************/
-
-/**
- * find_refl:
- * @list: The reflection list to search in
- * @h: The 'h' index to search for
- * @k: The 'k' index to search for
- * @l: The 'l' index to search for
- *
- * This function finds the first reflection in 'list' with the given indices.
- *
- * Since a %RefList can contain multiple reflections with the same indices, you
- * may need to use next_found_refl() to get the other reflections.
- *
- * Returns: The found reflection, or NULL if no reflection with the given
- * indices could be found.
- **/
-Reflection *find_refl(const RefList *list,
- signed int h, signed int k, signed int l)
-{
- unsigned int search = SERIAL(h, k, l);
- Reflection *refl;
-
- if ( list->head == NULL ) return NULL;
-
- /* Indices greater than or equal to 256 are filtered out when
- * reflections are added, so don't even bother looking.
- * (also, looking for such reflections causes trouble because the search
- * serial number would be invalid) */
- if ( abs(h) >= 256 ) return NULL;
- if ( abs(k) >= 256 ) return NULL;
- if ( abs(l) >= 256 ) return NULL;
-
- refl = list->head;
-
- while ( refl != NULL ) {
-
- if ( refl->serial == search ) {
-
- assert(search == refl->serial);
- assert(h == GET_H(refl->serial));
- assert(k == GET_K(refl->serial));
- assert(l == GET_L(refl->serial));
- return refl;
-
- } else {
-
- int dir = search > refl->serial;
- if ( refl->child[dir] != NULL ) {
- refl = refl->child[dir];
- } else {
- /* Hit the bottom of the tree */
- return NULL;
- }
-
- }
-
- }
-
- return NULL;
-}
-
-
-/**
- * next_found_refl:
- * @refl: A reflection returned by find_refl() or next_found_refl()
- *
- * This function returns the next reflection in @refl's list with the same
- * indices.
- *
- * Returns: The found reflection, or NULL if there are no more reflections with
- * the same indices.
- **/
-Reflection *next_found_refl(Reflection *refl)
-{
- if ( refl->next != NULL ) assert(refl->serial == refl->next->serial);
-
- return refl->next; /* Well, that was easy... */
-}
-
-
-/********************************** Getters ***********************************/
-
-/**
- * get_excitation_error:
- * @refl: A %Reflection
- *
- * Returns: The excitation error for the reflection.
- **/
-double get_excitation_error(const Reflection *refl)
-{
- return refl->data.excitation_error;
-}
-
-
-/**
- * get_detector_pos:
- * @refl: A %Reflection
- * @fs: Location at which to store the fast scan offset of the reflection
- * @ss: Location at which to store the slow scan offset of the reflection
- *
- **/
-void get_detector_pos(const Reflection *refl, double *fs, double *ss)
-{
- *fs = refl->data.fs;
- *ss = refl->data.ss;
-}
-
-
-/**
- * get_indices:
- * @refl: A %Reflection
- * @h: Location at which to store the 'h' index of the reflection
- * @k: Location at which to store the 'k' index of the reflection
- * @l: Location at which to store the 'l' index of the reflection
- *
- **/
-void get_indices(const Reflection *refl,
- signed int *h, signed int *k, signed int *l)
-{
- *h = GET_H(refl->serial);
- *k = GET_K(refl->serial);
- *l = GET_L(refl->serial);
-}
-
-
-/**
- * get_symmetric_indices:
- * @refl: A %Reflection
- * @hs: Location at which to store the 'h' index of the reflection
- * @ks: Location at which to store the 'k' index of the reflection
- * @ls: Location at which to store the 'l' index of the reflection
- *
- * This function gives the symmetric indices, that is, the "real" indices before
- * squashing down to the asymmetric reciprocal unit. This may be useful if the
- * list is indexed according to the asymmetric indices, but you still need
- * access to the symmetric version. This happens during post-refinement.
- *
- **/
-void get_symmetric_indices(const Reflection *refl,
- signed int *hs, signed int *ks,
- signed int *ls)
-{
- *hs = refl->data.hs;
- *ks = refl->data.ks;
- *ls = refl->data.ls;
-}
-
-
-/**
- * get_partiality:
- * @refl: A %Reflection
- *
- * Returns: The partiality of the reflection.
- **/
-double get_partiality(const Reflection *refl)
-{
- return refl->data.p;
-}
-
-
-/**
- * get_intensity:
- * @refl: A %Reflection
- *
- * Returns: The intensity of the reflection.
- **/
-double get_intensity(const Reflection *refl)
-{
- return refl->data.intensity;
-}
-
-
-/**
- * get_partial:
- * @refl: A %Reflection
- * @r1: Location at which to store the first excitation error
- * @r2: Location at which to store the second excitation error
- * @p: Location at which to store the partiality
- * @clamp_low: Location at which to store the first clamp status
- * @clamp_high: Location at which to store the second clamp status
- *
- * This function is used during post refinement (in conjunction with
- * set_partial()) to get access to the details of the partiality calculation.
- *
- **/
-void get_partial(const Reflection *refl, double *r1, double *r2, double *p,
- int *clamp_low, int *clamp_high)
-{
- *r1 = refl->data.r1;
- *r2 = refl->data.r2;
- *p = get_partiality(refl);
- *clamp_low = refl->data.clamp1;
- *clamp_high = refl->data.clamp2;
-}
-
-
-/**
- * get_scalable:
- * @refl: A %Reflection
- *
- * Returns: non-zero if this reflection can be scaled.
- *
- **/
-int get_scalable(const Reflection *refl)
-{
- return refl->data.scalable;
-}
-
-
-/**
- * get_refinable:
- * @refl: A %Reflection
- *
- * Returns: non-zero if this reflection can be used for post refinement.
- *
- **/
-int get_refinable(const Reflection *refl)
-{
- return refl->data.refinable;
-}
-
-
-/**
- * get_redundancy:
- * @refl: A %Reflection
- *
- * The redundancy of the reflection is the number of measurements that have been
- * made of it. Note that a redundancy of zero may have a special meaning, such
- * as that the reflection was impossible to integrate. Note further that each
- * reflection in the list has its own redundancy, even if there are multiple
- * copies of the reflection in the list. The total number of reflection
- * measurements should always be the sum of the redundancies in the entire list.
- *
- * Returns: the number of measurements of this reflection.
- *
- **/
-int get_redundancy(const Reflection *refl)
-{
- return refl->data.redundancy;
-}
-
-
-/**
- * get_esd_intensity:
- * @refl: A %Reflection
- *
- * Returns: the standard error in the intensity measurement (as returned by
- * get_intensity()) for this reflection.
- *
- **/
-double get_esd_intensity(const Reflection *refl)
-{
- return refl->data.esd_i;
-}
-
-
-/**
- * get_phase:
- * @refl: A %Reflection
- * @have_phase: Place to store a non-zero value if the phase is set, or NULL.
- *
- * Returns: the phase for this reflection.
- *
- **/
-double get_phase(const Reflection *refl, int *have_phase)
-{
- if ( have_phase != NULL ) *have_phase = refl->data.have_phase;
- return refl->data.phase;
-}
-
-
-/**
- * get_temp1:
- * @refl: A %Reflection
- *
- * The temporary values can be used according to the needs of the calling
- * program.
- *
- * Returns: the first temporary value for this reflection.
- *
- **/
-double get_temp1(const Reflection *refl)
-{
- return refl->data.temp1;
-}
-
-
-/**
- * get_temp2:
- * @refl: A %Reflection
- *
- * The temporary values can be used according to the needs of the calling
- * program.
- *
- * Returns: the second temporary value for this reflection.
- *
- **/
-double get_temp2(const Reflection *refl)
-{
- return refl->data.temp2;
-}
-
-
-/********************************** Setters ***********************************/
-
-/**
- * copy_data:
- * @to: %Reflection to copy data into
- * @from: %Reflection to copy data from
- *
- * This function is used to copy the data (which is everything listed above in
- * the list of getters and setters, apart from the indices themselves) from one
- * reflection to another. This might be used when creating a new list from an
- * old one, perhaps using the asymmetric indices instead of the raw indices for
- * the new list.
- *
- **/
-void copy_data(Reflection *to, const Reflection *from)
-{
- memcpy(&to->data, &from->data, sizeof(struct _refldata));
-}
-
-
-/**
- * set_detector_pos:
- * @refl: A %Reflection
- * @exerr: The excitation error for this reflection
- * @fs: The fast scan offset of the reflection
- * @ss: The slow scan offset of the reflection
- *
- **/
-void set_detector_pos(Reflection *refl, double exerr, double fs, double ss)
-{
- refl->data.excitation_error = exerr;
- refl->data.fs = fs;
- refl->data.ss = ss;
-}
-
-
-/**
- * set_partial:
- * @refl: A %Reflection
- * @r1: The first excitation error
- * @r2: The second excitation error
- * @p: The partiality
- * @clamp_low: The first clamp status
- * @clamp_high: The second clamp status
- *
- * This function is used during post refinement (in conjunction with
- * get_partial()) to get access to the details of the partiality calculation.
- *
- **/
-void set_partial(Reflection *refl, double r1, double r2, double p,
- double clamp_low, double clamp_high)
-{
- refl->data.r1 = r1;
- refl->data.r2 = r2;
- refl->data.p = p;
- refl->data.clamp1 = clamp_low;
- refl->data.clamp2 = clamp_high;
-}
-
-
-/**
- * set_int:
- * @refl: A %Reflection
- * @intensity: The intensity for the reflection.
- *
- * Set the intensity for the reflection. Note that retrieval is done with
- * get_intensity().
- **/
-void set_int(Reflection *refl, double intensity)
-{
- refl->data.intensity = intensity;
-}
-
-
-/**
- * set_scalable:
- * @refl: A %Reflection
- * @scalable: Non-zero if this reflection should be scaled.
- *
- **/
-void set_scalable(Reflection *refl, int scalable)
-{
- refl->data.scalable = scalable;
-}
-
-
-/**
- * set_refinable:
- * @refl: A %Reflection
- * @refinable: Non-zero if this reflection can be used for post refinement.
- *
- **/
-void set_refinable(Reflection *refl, int refinable)
-{
- refl->data.refinable = refinable;
-}
-
-
-/**
- * set_redundancy:
- * @refl: A %Reflection
- * @red: New redundancy for the reflection
- *
- * The redundancy of the reflection is the number of measurements that have been
- * made of it. Note that a redundancy of zero may have a special meaning, such
- * as that the reflection was impossible to integrate. Note further that each
- * reflection in the list has its own redundancy, even if there are multiple
- * copies of the reflection in the list. The total number of reflection
- * measurements should always be the sum of the redundancies in the entire list.
- *
- **/
-void set_redundancy(Reflection *refl, int red)
-{
- refl->data.redundancy = red;
-}
-
-
-/**
- * set_esd_intensity:
- * @refl: A %Reflection
- * @esd: New standard error for this reflection's intensity measurement
- *
- **/
-void set_esd_intensity(Reflection *refl, double esd)
-{
- refl->data.esd_i = esd;
-}
-
-
-/**
- * set_ph:
- * @refl: A %Reflection
- * @phase: New phase for the reflection
- *
- **/
-void set_ph(Reflection *refl, double phase)
-{
- refl->data.phase = phase;
- refl->data.have_phase = 1;
-}
-
-
-/**
- * set_symmetric_indices:
- * @refl: A %Reflection
- * @hs: The 'h' index of the reflection
- * @ks: The 'k' index of the reflection
- * @ls: The 'l' index of the reflection
- *
- * This function gives the symmetric indices, that is, the "real" indices before
- * squashing down to the asymmetric reciprocal unit. This may be useful if the
- * list is indexed according to the asymmetric indices, but you still need
- * access to the symmetric version. This happens during post-refinement.
- *
- **/
-void set_symmetric_indices(Reflection *refl,
- signed int hs, signed int ks, signed int ls)
-{
- refl->data.hs = hs;
- refl->data.ks = ks;
- refl->data.ls = ls;
-}
-
-
-/**
- * set_temp1
- * @refl: A %Reflection
- * @temp: New temporary value for the reflection
- *
- * The temporary values can be used according to the needs of the calling
- * program.
- *
- **/
-void set_temp1(Reflection *refl, double temp)
-{
- refl->data.temp1 = temp;
-}
-
-
-/**
- * set_temp2
- * @refl: A %Reflection
- * @temp: New temporary value for the reflection
- *
- * The temporary values can be used according to the needs of the calling
- * program.
- *
- **/
-void set_temp2(Reflection *refl, double temp)
-{
- refl->data.temp2 = temp;
-}
-
-
-/********************************* Insertion **********************************/
-
-static Reflection *rotate_once(Reflection *refl, int dir)
-{
- Reflection *s = refl->child[!dir];
-
- refl->child[!dir] = s->child[dir];
- s->child[dir] = refl;
-
- refl->col = RED;
- s->col = BLACK;
-
- return s;
-}
-
-
-static Reflection *rotate_twice(Reflection *refl, int dir)
-{
- refl->child[!dir] = rotate_once(refl->child[!dir], !dir);
- return rotate_once(refl, dir);
-}
-
-
-static int is_red(Reflection *refl)
-{
- return (refl != NULL) && (refl->col == RED);
-}
-
-
-static Reflection *insert_node(Reflection *refl, Reflection *new)
-{
- if ( refl == NULL ) {
-
- refl = new;
-
- } else {
-
- int dir;
- Reflection *rcd;
-
- assert(new->serial != refl->serial);
- dir = new->serial > refl->serial;
- refl->child[dir] = insert_node(refl->child[dir], new);
-
- rcd = refl->child[dir];
- if ( is_red(rcd) ) {
-
- if ( is_red(refl->child[!dir]) ) {
-
- refl->col = RED;
- refl->child[0]->col = BLACK;
- refl->child[1]->col = BLACK;
-
- } else {
-
- if ( is_red(rcd->child[dir] ) ) {
- refl = rotate_once(refl, !dir);
- } else if ( is_red(rcd->child[!dir] ) ) {
- refl = rotate_twice(refl, !dir);
- }
-
- }
- }
-
- }
-
- return refl;
-}
-
-
-/**
- * add_refl
- * @list: A %RefList
- * @h: The 'h' index of the reflection
- * @k: The 'k' index of the reflection
- * @l: The 'l' index of the reflection
- *
- * Adds a new reflection to @list. Note that the implementation allows there to
- * be multiple reflections with the same indices in the list, so this function
- * should succeed even if the given indices already feature in the list.
- *
- * Returns: The newly created reflection, or NULL on failure.
- *
- **/
-Reflection *add_refl(RefList *list, signed int h, signed int k, signed int l)
-{
- Reflection *new;
- Reflection *f;
-
- assert(abs(h)<256);
- assert(abs(k)<256);
- assert(abs(l)<256);
-
- new = new_node(SERIAL(h, k, l));
- if ( new == NULL ) return NULL;
-
- f = find_refl(list, h, k, l);
- if ( f == NULL ) {
-
- list->head = insert_node(list->head, new);
- list->head->col = BLACK;
-
- } else {
-
- /* New reflection is identical to a previous one */
- while ( f->next != NULL ) {
- f = f->next;
- }
- f->next = new;
- new->prev = f;
- }
-
- return new;
-}
-
-
-/**
- * add_refl_to_list
- * @refl: A %Reflection
- * @list: A %RefList
- *
- * Adds a reflection to @list. The reflection that actually gets added will be
- * a newly created one, and all the data will be copied across. The original
- * reflection will be destroyed and the new reflection returned.
- *
- * Returns: The newly created reflection, or NULL on failure.
- *
- **/
-Reflection *add_refl_to_list(Reflection *refl, RefList *list)
-{
- signed int h, k, l;
- Reflection *r_added;
-
- get_indices(refl, &h, &k, &l);
- r_added = add_refl(list, h, k, l);
- if ( r_added == NULL ) return NULL;
-
- copy_data(r_added, refl);
- reflection_free(refl);
-
- return r_added;
-}
-
-
-/********************************* Iteration **********************************/
-
-struct _reflistiterator {
-
- int stack_size;
- int stack_ptr;
- Reflection **stack;
-
-};
-
-
-/**
- * first_refl:
- * @list: A %RefList to iterate over
- * @piter: Address at which to store a %RefListIterator
- *
- * This function sets up the state required for iteration over the entire list,
- * and then returns the first reflection in the list. An iterator object will
- * be created and its address stored at the location given in piter.
- *
- * Returns: the first reflection in the list.
- *
- **/
-Reflection *first_refl(RefList *list, RefListIterator **piter)
-{
- RefListIterator *iter;
-
- iter = malloc(sizeof(struct _reflistiterator));
- iter->stack_size = 32;
- iter->stack = malloc(iter->stack_size*sizeof(Reflection *));
- iter->stack_ptr = 0;
- *piter = iter;
-
- Reflection *refl = list->head;
-
- do {
-
- if ( refl != NULL ) {
- iter->stack[iter->stack_ptr++] = refl;
- if ( iter->stack_ptr == iter->stack_size ) {
- iter->stack_size += 32;
- iter->stack = realloc(iter->stack,
- iter->stack_size*sizeof(Reflection *));
- }
- refl = refl->child[0];
- continue;
- }
-
- if ( iter->stack_ptr == 0 ) {
- free(iter->stack);
- free(iter);
- return NULL;
- }
-
- refl = iter->stack[--iter->stack_ptr];
-
- return refl;
-
- } while ( 1 );
-}
-
-
-/**
- * next_refl:
- * @refl: A reflection
- * @iter: A %RefListIterator
- *
- * This function looks up the next reflection in the list that was given earlier
- * to first_refl().
- *
- * Returns: the next reflection in the list, or NULL if no more.
- *
- **/
-Reflection *next_refl(Reflection *refl, RefListIterator *iter)
-{
- int returned = 1;
-
- do {
-
- if ( returned ) refl = refl->child[1];
- returned = 0;
-
- if ( refl != NULL ) {
-
- iter->stack[iter->stack_ptr++] = refl;
- if ( iter->stack_ptr == iter->stack_size ) {
- iter->stack_size += 32;
- iter->stack = realloc(iter->stack,
- iter->stack_size*sizeof(Reflection *));
- }
- refl = refl->child[0];
- continue;
-
- }
- if ( iter->stack_ptr == 0 ) {
- free(iter->stack);
- free(iter);
- return NULL;
- }
-
- return iter->stack[--iter->stack_ptr];
-
- } while ( 1 );
-}
-
-
-/*********************************** Voodoo ***********************************/
-
-static int recursive_depth(Reflection *refl)
-{
- int depth_left, depth_right;
-
- if ( refl == NULL ) return 0;
-
- depth_left = recursive_depth(refl->child[0]);
- depth_right = recursive_depth(refl->child[1]);
-
- return 1 + biggest(depth_left, depth_right);
-}
-
-
-static int recursive_count(Reflection *refl)
-{
- int count_left, count_right;
-
- if ( refl == NULL ) return 0;
-
- count_left = recursive_count(refl->child[0]);
- count_right = recursive_count(refl->child[1]);
-
- return 1 + count_left + count_right;
-}
-
-
-/**
- * num_reflections:
- * @list: A %RefList
- *
- * Returns: the number of reflections in @list.
- *
- **/
-int num_reflections(RefList *list)
-{
- return recursive_count(list->head);
-}
-
-
-/**
- * tree_depth:
- * @list: A %RefList
- *
- * If the depth of the tree is more than about 20, access to the list will be
- * slow. This should never happen.
- *
- * Returns: the depth of the RB-tree used internally to represent @list.
- *
- **/
-int tree_depth(RefList *list)
-{
- return recursive_depth(list->head);
-}
-
-
-/**
- * lock_reflection:
- * @refl: A %Reflection
- *
- * Acquires a lock on the reflection.
- */
-void lock_reflection(Reflection *refl)
-{
- pthread_mutex_lock(&refl->lock);
-}
-
-
-/**
- * unlock_reflection:
- * @refl: A %Reflection
- *
- * Releases a lock on the reflection.
- */
-void unlock_reflection(Reflection *refl)
-{
- pthread_mutex_unlock(&refl->lock);
-}
diff --git a/src/reflist.h b/src/reflist.h
deleted file mode 100644
index 955db69b..00000000
--- a/src/reflist.h
+++ /dev/null
@@ -1,112 +0,0 @@
-/*
- * reflist.h
- *
- * Fast reflection/peak list
- *
- * (c) 2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifndef REFLIST_H
-#define REFLIST_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-/**
- * RefList:
- *
- * A %RefList represents a list of Bragg reflections.
- *
- * This data structure is opaque. You must use the available accessor functions
- * to read and write its contents.
- *
- **/
-typedef struct _reflist RefList;
-
-/**
- * Reflection:
- *
- * A %Reflection represents a single Bragg reflection.
- *
- * This data structure is opaque. You must use the available accessor functions
- * to read and write its contents.
- *
- **/
-typedef struct _reflection Reflection;
-
-/**
- * RefListIterator:
- *
- * A %RefListIterator is an opaque data type used when iterating over a
- * %RefList.
- *
- **/
-typedef struct _reflistiterator RefListIterator;
-
-/* Creation/deletion */
-extern RefList *reflist_new(void);
-extern void reflist_free(RefList *list);
-extern Reflection *reflection_new(signed int h, signed int k, signed int l);
-extern void reflection_free(Reflection *refl);
-
-/* Search */
-extern Reflection *find_refl(const RefList *list, signed int h, signed int k, signed int l);
-extern Reflection *next_found_refl(Reflection *refl);
-
-/* Get */
-extern double get_excitation_error(const Reflection *refl);
-extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
-extern double get_partiality(const Reflection *refl);
-extern void get_indices(const Reflection *refl,
- signed int *h, signed int *k, signed int *l);
-extern void get_symmetric_indices(const Reflection *refl,
- signed int *hs, signed int *ks,
- signed int *ls);
-extern double get_intensity(const Reflection *refl);
-extern void get_partial(const Reflection *refl, double *r1, double *r2,
- double *p, int *clamp_low, int *clamp_high);
-extern int get_scalable(const Reflection *refl);
-extern int get_refinable(const Reflection *refl);
-extern int get_redundancy(const Reflection *refl);
-extern double get_temp1(const Reflection *refl);
-extern double get_temp2(const Reflection *refl);
-extern double get_esd_intensity(const Reflection *refl);
-extern double get_phase(const Reflection *refl, int *have_phase);
-
-/* Set */
-extern void copy_data(Reflection *to, const Reflection *from);
-extern void set_detector_pos(Reflection *refl, double exerr,
- double fs, double ss);
-extern void set_partial(Reflection *refl, double r1, double r2, double p,
- double clamp_low, double clamp_high);
-extern void set_int(Reflection *refl, double intensity);
-extern void set_scalable(Reflection *refl, int scalable);
-extern void set_refinable(Reflection *refl, int refinable);
-extern void set_redundancy(Reflection *refl, int red);
-extern void set_temp1(Reflection *refl, double temp);
-extern void set_temp2(Reflection *refl, double temp);
-extern void set_esd_intensity(Reflection *refl, double esd);
-extern void set_ph(Reflection *refl, double phase);
-extern void set_symmetric_indices(Reflection *refl,
- signed int hs, signed int ks, signed int ls);
-
-/* Insertion */
-extern Reflection *add_refl(RefList *list,
- signed int h, signed int k, signed int l);
-extern Reflection *add_refl_to_list(Reflection *refl, RefList *list);
-
-/* Iteration */
-extern Reflection *first_refl(RefList *list, RefListIterator **piter);
-extern Reflection *next_refl(Reflection *refl, RefListIterator *iter);
-
-/* Misc */
-extern int num_reflections(RefList *list);
-extern int tree_depth(RefList *list);
-extern void lock_reflection(Reflection *refl);
-extern void unlock_reflection(Reflection *refl);
-
-#endif /* REFLIST_H */
diff --git a/src/scaling-report.c b/src/scaling-report.c
index b499e59d..32791bbc 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -234,7 +234,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
Ipart = get_intensity(refl);
Ifull = get_intensity(f);
- if ( Ifull < 10 ) continue; /* FIXME: Ugh */
+ //if ( Ifull < 10 ) continue; /* FIXME: Ugh */
pobs = Ipart/(images[i].osf*Ifull);
pcalc = get_partiality(refl);
diff --git a/src/thread-pool.c b/src/thread-pool.c
deleted file mode 100644
index 0ae26e17..00000000
--- a/src/thread-pool.c
+++ /dev/null
@@ -1,283 +0,0 @@
-/*
- * thread-pool.c
- *
- * A thread pool implementation
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <pthread.h>
-#include <assert.h>
-
-#ifdef HAVE_CPU_AFFINITY
-#include <sched.h>
-#endif
-
-
-#include "utils.h"
-
-
-/**
- * SECTION:thread-pool
- * @short_description: The thread pool
- * @title: The thread pool
- * @section_id:
- * @see_also:
- * @include: "thread-pool.h"
- * @Image:
- *
- * The thread pool helps when running many tasks in parallel. It takes care of
- * starting and stopping threads, and presents a relatively simple interface to
- * the individual programs.
- */
-
-/* ------------------------------ CPU affinity ------------------------------ */
-
-#ifdef HAVE_CPU_AFFINITY
-
-static void set_affinity(int n, int cpu_num, int cpu_groupsize, int cpu_offset)
-{
- cpu_set_t c;
- int group;
- int n_cpu_groups;
- int i;
-
- if ( cpu_num == 0 ) return;
-
- CPU_ZERO(&c);
-
- /* Work out which group this thread belongs to */
- group = (n / cpu_groupsize) + cpu_offset;
-
- /* Work out which CPUs should be used for this group */
- n_cpu_groups = cpu_num / cpu_groupsize;
- group = group % n_cpu_groups;
-
- /* Set flags */
- for ( i=0; i<cpu_groupsize; i++ ) {
-
- int cpu = cpu_groupsize*group + i;
-
- CPU_SET(cpu, &c);
-
- }
-
- if ( sched_setaffinity(0, sizeof(cpu_set_t), &c) ) {
-
- /* Cannot use ERROR() just yet */
- fprintf(stderr, "%i: Failed to set CPU affinity.\n", n);
-
- }
-}
-
-#else /* HAVE_CPU_AFFINITY */
-
-static void set_affinity(int n, int cpu_num, int cpu_groupsize, int cpu_offset)
-{
- /* Do absolutely nothing */
-}
-
-#endif /* HAVE_CPU_AFFINITY */
-
-
-/* --------------------------- Status label stuff --------------------------- */
-
-static int use_status_labels = 0;
-static pthread_key_t status_label_key;
-pthread_mutex_t stderr_lock = PTHREAD_MUTEX_INITIALIZER;
-
-struct worker_args
-{
- struct task_queue_range *tqr;
- struct task_queue *tq;
- int id;
- int cpu_num;
- int cpu_groupsize;
- int cpu_offset;
-};
-
-
-signed int get_status_label()
-{
- int *cookie;
-
- if ( !use_status_labels ) {
- return -1;
- }
-
- cookie = pthread_getspecific(status_label_key);
- return *cookie;
-}
-
-
-struct task_queue
-{
- pthread_mutex_t lock;
-
- int n_started;
- int n_completed;
- int max;
-
- void *(*get_task)(void *);
- void (*finalise)(void *, void *);
- void *queue_args;
- void (*work)(void *, int);
-};
-
-
-static void *task_worker(void *pargsv)
-{
- struct worker_args *w = pargsv;
- struct task_queue *q = w->tq;
- int *cookie;
-
- set_affinity(w->id, w->cpu_num, w->cpu_groupsize, w->cpu_offset);
-
- cookie = malloc(sizeof(int));
- *cookie = w->id;
- pthread_setspecific(status_label_key, cookie);
-
- free(w);
-
- do {
-
- void *task;
- int cookie;
-
- /* Get a task */
- pthread_mutex_lock(&q->lock);
- if ( (q->max) && (q->n_started >= q->max) ) {
- pthread_mutex_unlock(&q->lock);
- break;
- }
- task = q->get_task(q->queue_args);
-
- /* No more tasks? */
- if ( task == NULL ) {
- pthread_mutex_unlock(&q->lock);
- break;
- }
-
- q->n_started++;
- pthread_mutex_unlock(&q->lock);
-
- cookie = *(int *)pthread_getspecific(status_label_key);
- q->work(task, cookie);
-
- /* Update totals etc */
- pthread_mutex_lock(&q->lock);
- q->n_completed++;
- if ( q->finalise ) {
- q->finalise(q->queue_args, task);
- }
- pthread_mutex_unlock(&q->lock);
-
- } while ( 1 );
-
- free(cookie);
-
- return NULL;
-}
-
-
-/**
- * run_threads:
- * @n_threads: The number of threads to run in parallel
- * @work: The function to be called to do the work
- * @get_task: The function which will determine the next unassigned task
- * @final: The function which will be called to clean up after a task
- * @queue_args: A pointer to any data required to determine the next task
- * @max: Stop calling get_task after starting this number of jobs
- * @cpu_num: The number of CPUs in the system
- * @cpu_groupsize: The group size into which the CPUs are grouped
- * @cpu_offset: The CPU group number at which to start pinning threads
- *
- * 'get_task' will be called every time a worker is idle. It returns either
- * NULL, indicating that no further work is available, or a pointer which will
- * be passed to 'work'.
- *
- * 'final' will be called once per image, and will be given both queue_args
- * and the last task pointer.
- *
- * 'get_task' and 'final' will be called only under lock, and so do NOT need to
- * be re-entrant or otherwise thread safe. 'work', of course, needs to be
- * thread safe.
- *
- * Work will stop after 'max' tasks have been processed whether get_task
- * returned NULL or not. If "max" is zero, all tasks will be processed.
- *
- * Returns: The number of tasks completed.
- **/
-int run_threads(int n_threads, TPWorkFunc work,
- TPGetTaskFunc get_task, TPFinalFunc final,
- void *queue_args, int max,
- int cpu_num, int cpu_groupsize, int cpu_offset)
-{
- pthread_t *workers;
- int i;
- struct task_queue q;
-
- pthread_key_create(&status_label_key, NULL);
-
- workers = malloc(n_threads * sizeof(pthread_t));
-
- pthread_mutex_init(&q.lock, NULL);
- q.work = work;
- q.get_task = get_task;
- q.finalise = final;
- q.queue_args = queue_args;
- q.n_started = 0;
- q.n_completed = 0;
- q.max = max;
-
- /* Now it's safe to start using the status labels */
- if ( n_threads > 1 ) use_status_labels = 1;
-
- /* Start threads */
- for ( i=0; i<n_threads; i++ ) {
-
- struct worker_args *w;
-
- w = malloc(sizeof(struct worker_args));
-
- w->tq = &q;
- w->tqr = NULL;
- w->id = i;
- w->cpu_num = cpu_num;
- w->cpu_groupsize = cpu_groupsize;
- w->cpu_offset = cpu_offset;
-
- if ( pthread_create(&workers[i], NULL, task_worker, w) ) {
- /* Not ERROR() here */
- fprintf(stderr, "Couldn't start thread %i\n", i);
- n_threads = i;
- break;
- }
-
- }
-
- /* Join threads */
- for ( i=0; i<n_threads; i++ ) {
- pthread_join(workers[i], NULL);
- }
-
- use_status_labels = 0;
-
- free(workers);
-
- return q.n_completed;
-}
diff --git a/src/thread-pool.h b/src/thread-pool.h
deleted file mode 100644
index a99a7ade..00000000
--- a/src/thread-pool.h
+++ /dev/null
@@ -1,71 +0,0 @@
-/*
- * thread-pool.h
- *
- * A thread pool implementation
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifndef THREAD_POOL_H
-#define THREAD_POOL_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <pthread.h>
-
-extern pthread_mutex_t stderr_lock;
-extern signed int get_status_label(void);
-
-
-/**
- * TPGetTaskFunc:
- * @qargs: The queue_args pointer which was given to run_threads().
- * Returns: A pointer which will be passed to the worker function.
- *
- * This function is called, non-reentrantly, to get a new work item to give to
- * your work function. The stuff you need to generate the new work item should
- * have been stored in @qargs which was passed to run_threads().
- *
- **/
-typedef void *(*TPGetTaskFunc)(void *qargs);
-
-
-/**
- * TPWorkFunc:
- * @work: The queue_args pointer which was given to run_threads().
- * @cookie: A small integral number which is guaranteed to be unique among all
- * currently running threads.
- *
- * This function is called, reentrantly, for each work item.
- *
- **/
-typedef void (*TPWorkFunc)(void *work, int cookie);
-
-
-/**
- * TPFinalFunc:
- * @qargs: The queue_args pointer which was given to run_threads().
- * @work: The pointer which was returned by your get_task function.
- *
- * This function is called, non-reentrantly, after each work item has been
- * completed. A typical use might be to update some counters inside @qargs
- * according to fields withing @work which were filled by your 'work' function.
- *
- **/
-typedef void (*TPFinalFunc)(void *qargs, void *work);
-
-
-extern int run_threads(int n_threads, TPWorkFunc work,
- TPGetTaskFunc get_task, TPFinalFunc final,
- void *queue_args, int max,
- int cpu_num, int cpu_groupsize, int cpu_offset);
-
-
-#endif /* THREAD_POOL_H */
diff --git a/src/utils.c b/src/utils.c
deleted file mode 100644
index 6b29627f..00000000
--- a/src/utils.c
+++ /dev/null
@@ -1,488 +0,0 @@
-/*
- * utils.c
- *
- * Utility stuff
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#include <libgen.h>
-#include <math.h>
-#include <string.h>
-#include <stdio.h>
-#include <unistd.h>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
-
-#include "utils.h"
-#include "image.h"
-
-
-/**
- * SECTION:utils
- * @short_description: Miscellaneous utilities
- * @title: Utilities
- * @section_id:
- * @see_also:
- * @include: "utils.h"
- * @Image:
- *
- * Wibble
- */
-
-/**
- * show_matrix_eqn:
- * @M: A matrix
- * @v: A vector
- * @r: The number of elements in @v and the side length of @M.
- *
- * Displays a matrix equation of the form @M.a = @v.
- *
- * @M must be square.
- **/
-void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
-{
- int i, j;
-
- for ( i=0; i<r; i++ ) {
- STATUS("[ ");
- for ( j=0; j<r; j++ ) {
- STATUS("%+9.3e ", gsl_matrix_get(M, i, j));
- }
- STATUS("][ a%2i ] = [ %+9.3e ]\n", i, gsl_vector_get(v, i));
- }
-}
-
-
-size_t notrail(char *s)
-{
- size_t i;
- size_t munched = 0;
-
- for ( i=strlen(s)-1; i>=0; i-- ) {
- if ( (s[i] == ' ') || (s[i] == '\t') ) {
- s[i] = '\0';
- munched++;
- } else {
- return munched;
- }
- }
-
- return munched;
-}
-
-
-void chomp(char *s)
-{
- size_t i;
-
- if ( !s ) return;
-
- for ( i=0; i<strlen(s); i++ ) {
- if ( (s[i] == '\n') || (s[i] == '\r') ) {
- s[i] = '\0';
- return;
- }
- }
-}
-
-
-void progress_bar(int val, int total, const char *text)
-{
- double frac;
- int n, i;
- char s[1024];
- const int width = 50;
-
- if ( total == 0 ) return;
-
- if ( !isatty(STDERR_FILENO) ) return;
- if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
-
- frac = (double)val/total;
- n = (int)(frac*width);
-
- for ( i=0; i<n; i++ ) s[i] = '=';
- for ( i=n; i<width; i++ ) s[i] = '.';
- s[width] = '\0';
-
- pthread_mutex_lock(&stderr_lock);
- fprintf(stderr, "\r%s: |%s|", text, s);
- if ( val == total ) fprintf(stderr, "\n");
- pthread_mutex_unlock(&stderr_lock);
-
- fflush(stdout);
-}
-
-
-double random_flat(double max)
-{
- return max * (double)random()/RAND_MAX;
-}
-
-
-double flat_noise(double expected, double width)
-{
- double noise = random_flat(2.0*width);
- return expected+noise-width;
-}
-
-
-double gaussian_noise(double expected, double stddev)
-{
- double x1, x2, noise;
-
- /* A uniformly distributed random number between 0 and 1 */
- x1 = ((double)random()/RAND_MAX);
- x2 = ((double)random()/RAND_MAX);
-
- noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
-
- return expected + noise*stddev;
-}
-
-
-static int fake_poisson_noise(double expected)
-{
- double rf = gaussian_noise(expected, sqrt(expected));
- return (int)rf;
-}
-
-
-int poisson_noise(double expected)
-{
- double L;
- int k = 0;
- double p = 1.0;
-
- /* For large values of the mean, we get big problems with arithmetic.
- * In such cases, fall back on a Gaussian with the right variance. */
- if ( expected > 100.0 ) return fake_poisson_noise(expected);
-
- L = exp(-expected);
-
- do {
-
- double r;
-
- k++;
- r = (double)random()/(double)RAND_MAX;
- p *= r;
-
- } while ( p > L );
-
- return k - 1;
-}
-
-
-/**
- * SECTION:quaternion
- * @short_description: Simple quaternion handling
- * @title: Quaternion
- * @section_id:
- * @see_also:
- * @include: "utils.h"
- * @Image:
- *
- * There is a simple quaternion structure in CrystFEL. At the moment, it is
- * only used when simulating patterns, as an argument to cell_rotate() to
- * orient the unit cell.
- */
-
-/**
- * quaternion_modulus:
- * @q: A %quaternion
- *
- * If a quaternion represents a pure rotation, its modulus should be unity.
- *
- * Returns: the modulus of the given quaternion.
- **/
-double quaternion_modulus(struct quaternion q)
-{
- return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
-}
-
-
-/**
- * normalise_quaternion:
- * @q: A %quaternion
- *
- * Rescales the quaternion such that its modulus is unity.
- *
- * Returns: the normalised version of @q
- **/
-struct quaternion normalise_quaternion(struct quaternion q)
-{
- double mod;
- struct quaternion r;
-
- mod = quaternion_modulus(q);
-
- r.w = q.w / mod;
- r.x = q.x / mod;
- r.y = q.y / mod;
- r.z = q.z / mod;
-
- return r;
-}
-
-
-/**
- * random_quaternion:
- *
- * Returns: a randomly generated, normalised, quaternion.
- **/
-struct quaternion random_quaternion()
-{
- struct quaternion q;
-
- q.w = 2.0*(double)random()/RAND_MAX - 1.0;
- q.x = 2.0*(double)random()/RAND_MAX - 1.0;
- q.y = 2.0*(double)random()/RAND_MAX - 1.0;
- q.z = 2.0*(double)random()/RAND_MAX - 1.0;
- q = normalise_quaternion(q);
-
- return q;
-}
-
-
-/**
- * quaternion_valid:
- * @q: A %quaternion
- *
- * Checks if the given quaternion is normalised.
- *
- * This function performs a nasty floating point comparison of the form
- * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be
- * relied upon to spot anything other than the most obvious input error.
- *
- * Returns: 1 if the quaternion is normalised, 0 if not.
- **/
-int quaternion_valid(struct quaternion q)
-{
- double qmod;
-
- qmod = quaternion_modulus(q);
-
- /* Modulus = 1 to within some tolerance?
- * Nasty allowance for floating-point accuracy follows... */
- if ( (qmod > 0.999) && (qmod < 1.001) ) return 1;
-
- return 0;
-}
-
-
-/**
- * quat_rot
- * @q: A vector (in the form of an %rvec)
- * @z: A %quaternion
- *
- * Rotates a vector according to a quaternion.
- *
- * Returns: A rotated version of @p.
- **/
-struct rvec quat_rot(struct rvec q, struct quaternion z)
-{
- struct rvec res;
- double t01, t02, t03, t11, t12, t13, t22, t23, t33;
-
- t01 = z.w*z.x;
- t02 = z.w*z.y;
- t03 = z.w*z.z;
- t11 = z.x*z.x;
- t12 = z.x*z.y;
- t13 = z.x*z.z;
- t22 = z.y*z.y;
- t23 = z.y*z.z;
- t33 = z.z*z.z;
-
- res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
- + (2.0 * (t12 + t03)) * q.v
- + (2.0 * (t13 - t02)) * q.w;
-
- res.v = (2.0 * (t12 - t03)) * q.u
- + (1.0 - 2.0 * (t11 + t33)) * q.v
- + (2.0 * (t01 + t23)) * q.w;
-
- res.w = (2.0 * (t02 + t13)) * q.u
- + (2.0 * (t23 - t01)) * q.v
- + (1.0 - 2.0 * (t11 + t22)) * q.w;
-
- return res;
-}
-
-
-/* Return non-zero if c is in delims */
-static int assplode_isdelim(const char c, const char *delims)
-{
- size_t i;
- for ( i=0; i<strlen(delims); i++ ) {
- if ( c == delims[i] ) return 1;
- }
- return 0;
-}
-
-
-static int assplode_extract(char ***pbits, int n, size_t n_captured,
- size_t start, const char *a)
-{
- char **bits = *pbits;
- bits = realloc(bits, sizeof(char *)*(n+1));
- bits[n] = malloc(n_captured+1);
- memcpy(bits[n], a+start, n_captured);
- bits[n][n_captured] = '\0';
- n++;
- *pbits = bits;
- return n;
-}
-
-
-/* Split the string 'a' using 'delims' as a zero-terminated list of
- * deliminators.
- * Store each segment in bits[0...n] where n is the number of segments and is
- * the return value. pbits = &bits
- * Each segment needs to be freed with free() when finished with.
- * The array of bits also needs to be freed with free() when finished with,
- * unless n=0 in which case bits==NULL
- */
-int assplode(const char *a, const char *delims, char ***pbits,
- AssplodeFlag flags)
-{
- size_t i, start, n_captured;
- int n, last_was_delim;
- char **bits;
-
- n = 0;
- i = 0;
- n_captured = 0;
- start = 0;
- last_was_delim = 0;
- bits = NULL;
- while ( i < strlen(a) ) {
-
- if ( assplode_isdelim(a[i], delims) ) {
-
- if ( n_captured > 0 ) {
- /* This is a deliminator after a sequence of
- * non-deliminator chars */
- n = assplode_extract(&bits, n, n_captured,
- start, a);
- }
-
- n_captured = 0;
- if ( (flags & ASSPLODE_DUPS) && last_was_delim ) {
- n = assplode_extract(&bits, n, 0, start, a);
- }
- last_was_delim = 1;
-
- } else {
-
- if ( n_captured == 0 ) {
- /* No characters currently found, so this is the
- * start */
- start = i;
- }
- n_captured++;
- last_was_delim = 0;
-
- }
-
- i++;
-
- }
- /* Left over characters at the end? */
- if ( n_captured > 0 ) {
- n = assplode_extract(&bits, n, n_captured, start, a);
- }
-
- *pbits = bits;
- return n;
-}
-
-
-char *check_prefix(char *prefix)
-{
- int r;
- struct stat statbuf;
- char *new;
- size_t len;
-
- /* Is "prefix" a directory? */
- r = stat(prefix, &statbuf);
- if ( r != 0 ) {
- /* "prefix" probably doesn't exist. This is fine - assume
- * the user knows what they're doing, and that "prefix"
- * suffixed with the actual filename will produce something
- * sensible. */
- return prefix;
- }
-
- if ( !S_ISDIR(statbuf.st_mode) ) {
- /* Also fine, as above. */
- return prefix;
- }
-
- /* Does the prefix end in a slash? */
- if ( prefix[strlen(prefix)-1] == '/' ) {
- /* This looks sensible. */
- return prefix;
- }
-
- STATUS("Your prefix ('%s') is a directory, but doesn't end"
- " with a slash. I'm going to add it for you.\n", prefix);
- STATUS("If this isn't what you want, run with --no-check-prefix.\n");
- len = strlen(prefix)+2;
- new = malloc(len);
- snprintf(new, len, "%s/", prefix);
- free(prefix);
- return new;
-}
-
-
-char *safe_basename(const char *in)
-{
- int i;
- char *cpy;
- char *res;
-
- cpy = strdup(in);
-
- /* Get rid of any trailing slashes */
- for ( i=strlen(cpy)-1; i>0; i-- ) {
- if ( cpy[i] == '/' ) {
- cpy[i] = '\0';
- } else {
- break;
- }
- }
-
- /* Find the base name */
- for ( i=strlen(cpy)-1; i>0; i-- ) {
- if ( cpy[i] == '/' ) {
- i++;
- break;
- }
- }
-
- res = strdup(cpy+i);
- /* If we didn't find a previous slash, i==0 so res==cpy */
-
- free(cpy);
-
- return res;
-}
-
-
-#ifdef GSL_FUDGE
-/* Force the linker to bring in CBLAS to make GSL happy */
-void utils_fudge_gslcblas()
-{
- STATUS("%p\n", cblas_sgemm);
-}
-#endif
diff --git a/src/utils.h b/src/utils.h
deleted file mode 100644
index 1179a57f..00000000
--- a/src/utils.h
+++ /dev/null
@@ -1,236 +0,0 @@
-/*
- * utils.h
- *
- * Utility stuff
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifndef UTILS_H
-#define UTILS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <math.h>
-#include <complex.h>
-#include <string.h>
-#include <stdlib.h>
-#include <pthread.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-
-
-#include "thread-pool.h"
-
-
-/* -------------------------- Fundamental constants ------------------------ */
-
-/* Electron charge in C */
-#define ELECTRON_CHARGE (1.6021773e-19)
-
-/* Planck's constant (Js) */
-#define PLANCK (6.62606896e-34)
-
-/* Speed of light in vacuo (m/s) */
-#define C_VACUO (299792458)
-
-/* Thomson scattering length (m) */
-#define THOMSON_LENGTH (2.81794e-15)
-
-
-/* ------------------------------ Quaternions ------------------------------- */
-
-/**
- * quaternion:
- *
- * <programlisting>
- * struct quaternion
- * {
- * double w
- * double x
- * double y
- * double z
- * };
- * </programlisting>
- *
- * A structure representing a quaternion.
- *
- **/
-struct quaternion;
-
-struct quaternion {
- double w;
- double x;
- double y;
- double z;
-};
-
-extern struct quaternion normalise_quaternion(struct quaternion q);
-extern double quaternion_modulus(struct quaternion q);
-extern struct quaternion random_quaternion(void);
-extern int quaternion_valid(struct quaternion q);
-extern struct rvec quat_rot(struct rvec q, struct quaternion z);
-
-
-/* --------------------------- Useful functions ----------------------------- */
-
-extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r);
-extern size_t notrail(char *s);
-extern void chomp(char *s);
-
-typedef enum {
- ASSPLODE_NONE = 0,
- ASSPLODE_DUPS = 1<<0
-} AssplodeFlag;
-extern int assplode(const char *a, const char *delims, char ***pbits,
- AssplodeFlag flags);
-
-extern void progress_bar(int val, int total, const char *text);
-extern double random_flat(double max);
-extern double flat_noise(double expected, double width);
-extern double gaussian_noise(double expected, double stddev);
-extern int poisson_noise(double expected);
-
-/* Keep these ones inline, to avoid function call overhead */
-static inline struct quaternion invalid_quaternion(void)
-{
- struct quaternion quat;
- quat.w = 0.0;
- quat.x = 0.0;
- quat.y = 0.0;
- quat.z = 0.0;
- return quat;
-}
-
-static inline double distance(double x1, double y1, double x2, double y2)
-{
- return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
-}
-
-static inline double modulus(double x, double y, double z)
-{
- return sqrt(x*x + y*y + z*z);
-}
-
-static inline double modulus_squared(double x, double y, double z) {
- return x*x + y*y + z*z;
-}
-
-static inline double distance3d(double x1, double y1, double z1,
- double x2, double y2, double z2)
-{
- return modulus(x1-x2, y1-y2, z1-z2);
-}
-
-/* Answer in radians */
-static inline double angle_between(double x1, double y1, double z1,
- double x2, double y2, double z2)
-{
- double mod1 = modulus(x1, y1, z1);
- double mod2 = modulus(x2, y2, z2);
- double cosine = (x1*x2 + y1*y2 + z1*z2) / (mod1*mod2);
-
- /* Fix domain if outside due to rounding */
- if ( cosine > 1.0 ) cosine = 1.0;
- if ( cosine < -1.0 ) cosine = -1.0;
-
- return acos(cosine);
-}
-
-static inline int within_tolerance(double a, double b, double percent)
-{
- double tol = fabs(a) * (percent/100.0);
- if ( fabs(b-a) < tol ) return 1;
- return 0;
-}
-
-
-/* ----------------------------- Useful macros ------------------------------ */
-
-#define rad2deg(a) ((a)*180/M_PI)
-#define deg2rad(a) ((a)*M_PI/180)
-
-#define is_odd(a) ((a)%2==1)
-
-#define biggest(a,b) ((a>b) ? (a) : (b))
-#define smallest(a,b) ((a<b) ? (a) : (b))
-
-
-/* Photon energy (J) to wavelength (m) */
-#define ph_en_to_lambda(a) ((PLANCK*C_VACUO)/(a))
-
-/* Photon wavelength (m) to energy (J) */
-#define ph_lambda_to_en(a) ((PLANCK*C_VACUO)/(a))
-
-/* eV to Joules */
-#define eV_to_J(a) ((a)*ELECTRON_CHARGE)
-
-/* Joules to eV */
-#define J_to_eV(a) ((a)/ELECTRON_CHARGE)
-
-#define UNUSED __attribute__((unused))
-
-/* -------------------- Indexed lists for specified types ------------------- */
-
-#include "defs.h"
-
-#define LIST_SIZE (IDIM*IDIM*IDIM)
-
-/* Create functions for storing reflection intensities indexed as h,k,l */
-#define LABEL(x) x##_intensity
-#define TYPE double
-#include "list_tmp.h"
-
-/* CAs above, but for phase values */
-#define LABEL(x) x##_phase
-#define TYPE double
-#include "list_tmp.h"
-
-/* As above, but for (unsigned) integer counts */
-#define LABEL(x) x##_count
-#define TYPE unsigned int
-#include "list_tmp.h"
-
-/* As above, but for simple flags */
-#define LABEL(x) x##_flag
-#define TYPE unsigned char
-#include "list_tmp.h"
-
-
-/* ------------------------------ Message macros ---------------------------- */
-
-extern pthread_mutex_t stderr_lock;
-
-#define ERROR(...) { \
- int error_print_val = get_status_label(); \
- pthread_mutex_lock(&stderr_lock); \
- if ( error_print_val >= 0 ) { \
- fprintf(stderr, "%3i: ", error_print_val); \
- } \
- fprintf(stderr, __VA_ARGS__); \
- pthread_mutex_unlock(&stderr_lock); \
- }
-
-#define STATUS(...) { \
- int status_print_val = get_status_label(); \
- pthread_mutex_lock(&stderr_lock); \
- if ( status_print_val >= 0 ) { \
- fprintf(stderr, "%3i: ", status_print_val); \
- } \
- fprintf(stderr, __VA_ARGS__); \
- pthread_mutex_unlock(&stderr_lock); \
- }
-
-
-/* ------------------------------ File handling ----------------------------- */
-
-extern char *check_prefix(char *prefix);
-extern char *safe_basename(const char *in);
-
-
-#endif /* UTILS_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-06 02:49:38 +0000