From 629934d82e202ea04b334c49efffe09aaa0f1c4e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 14 Aug 2014 15:07:21 +0200 Subject: Remove "sphere", "thin" and "gaussian" partiality models, add "scgaussian" --- Makefile.am | 13 +- libcrystfel/src/geometry.c | 68 +----- libcrystfel/src/geometry.h | 9 +- libcrystfel/src/integration.c | 3 +- src/partial_sim.c | 2 +- src/partialator.c | 12 +- src/post-refinement.c | 67 ++---- src/post-refinement.h | 12 +- src/process_image.c | 2 +- tests/pr_l_gradient_check.c | 379 ----------------------------- tests/pr_p_gradient_check.c | 31 +-- tests/pr_pl_gradient_check.c | 540 ------------------------------------------ tests/prof2d_check.c | 2 +- 13 files changed, 45 insertions(+), 1095 deletions(-) delete mode 100644 tests/pr_l_gradient_check.c delete mode 100644 tests/pr_pl_gradient_check.c diff --git a/Makefile.am b/Makefile.am index 033ea746..78f3e277 100644 --- a/Makefile.am +++ b/Makefile.am @@ -10,8 +10,7 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_p_gradient_check tests/symmetry_check \ tests/centering_check tests/transformation_check \ - tests/cell_check tests/pr_l_gradient_check \ - tests/pr_pl_gradient_check tests/ring_check \ + tests/cell_check tests/ring_check \ tests/prof2d_check tests/ambi_check MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \ @@ -20,9 +19,7 @@ MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \ PARTIAL_CHECKS = tests/partialator_merge_check_1 \ tests/partialator_merge_check_2 \ tests/partialator_merge_check_3 \ - tests/pr_p_gradient_check \ - tests/pr_l_gradient_check \ - tests/pr_pl_gradient_check + tests/pr_p_gradient_check TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \ tests/integration_check \ @@ -117,12 +114,6 @@ tests_ambi_check_SOURCES = tests/ambi_check.c tests_pr_p_gradient_check_SOURCES = tests/pr_p_gradient_check.c \ src/post-refinement.c -tests_pr_pl_gradient_check_SOURCES = tests/pr_pl_gradient_check.c \ - src/post-refinement.c - -tests_pr_l_gradient_check_SOURCES = tests/pr_l_gradient_check.c \ - src/post-refinement.c - tests_centering_check_SOURCES = tests/centering_check.c tests_transformation_check_SOURCES = tests/transformation_check.c diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index aff19198..8952486d 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -121,18 +121,10 @@ static double partiality(PartialityModel pmodel, case PMODEL_UNITY: return 1.0; - case PMODEL_SPHERE: - plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; - phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh; - return plow - phigh; - - case PMODEL_GAUSSIAN: + case PMODEL_SCGAUSSIAN: plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r)); phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r)); - return plow - phigh; - - case PMODEL_THIN: - return 1.0 - (rmid*rmid)/(r*r); + return 4.0*(plow-phigh)*r / (3.0*D); case PMODEL_SCSPHERE: plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; @@ -157,7 +149,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, Reflection *refl; double cet, cez; /* Centre of Ewald sphere */ double pr; - double L, D; + double D; double del; /* Don't predict 000 */ @@ -198,44 +190,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, return NULL; } - /* Conditions for reflection to be excited at all */ - switch ( pmodel ) { - - case PMODEL_UNITY: /* PMODEL_UNITY shouldn't end up here */ - case PMODEL_SPHERE: - case PMODEL_GAUSSIAN: - case PMODEL_SCSPHERE: - if ( (signbit(rlow) == signbit(rhigh)) - && (fabs(rlow) > pr) - && (fabs(rhigh) > pr) ) return NULL; - break; - - case PMODEL_THIN: - if ( fabs(rmid) > pr ) return NULL; - break; - - } + /* Condition for reflection to be excited at all */ + if ( (signbit(rlow) == signbit(rhigh)) + && (fabs(rlow) > pr) + && (fabs(rhigh) > pr) ) return NULL; D = rlow - rhigh; - /* Lorentz factor is determined direction from the r values, before - * clamping. The multiplication by 0.01e9 to make the - * correction factor vaguely near 1. */ - switch ( pmodel ) { - - case PMODEL_SPHERE: - case PMODEL_GAUSSIAN: - L = LORENTZ_SCALE / D; - break; - - case PMODEL_UNITY: /* PMODEL_UNITY shouldn't end up here */ - case PMODEL_THIN: - case PMODEL_SCSPHERE: - L = 1.0; - break; - - } - /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the * reflection. @@ -282,19 +243,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, set_detector_pos(refl, 0.0, xda, yda); } - if ( pmodel != PMODEL_THIN ) { - set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); - } else { - /* If we are using the TES (Thin Ewald Sphere) model, we abuse - * the fields as follows: - * rlow = the r value for the middle (only) Ewald sphere - * rhigh = 0.0 - * clamp_low = 0 - * clamp_high = +1 - */ - set_partial(refl, rmid, 0.0, part, 0, +1); - } - set_lorentz(refl, L); + set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + set_lorentz(refl, 1.0); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index d8d226f0..ccdea1c2 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -46,23 +46,18 @@ extern "C" { /** * PartialityModel: - * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal - * space. * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. - * @PMODEL_GAUSSIAN : Gaussian profiles in 3D - * @PMODEL_THIN : Thin Ewald sphere intersecting sphere * @PMODEL_SCSPHERE : Sphere model with source coverage factor included + * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included * * A %PartialityModel describes a geometrical model which can be used to * calculate spot partialities and Lorentz correction factors. **/ typedef enum { - PMODEL_SPHERE, PMODEL_UNITY, - PMODEL_GAUSSIAN, - PMODEL_THIN, PMODEL_SCSPHERE, + PMODEL_SCGAUSSIAN, } PartialityModel; diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index 635a7c74..34bde83a 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -1914,7 +1914,8 @@ void integrate_all_2(struct image *image, IntegrationMethod meth, IntDiag int_diag, signed int idh, signed int idk, signed int idl) { - integrate_all_3(image, meth, PMODEL_SPHERE, 0.0, ir_inn, ir_mid, ir_out, + integrate_all_3(image, meth, PMODEL_SCSPHERE, 0.0, + ir_inn, ir_mid, ir_out, int_diag, idh, idk, idl); } diff --git a/src/partial_sim.c b/src/partial_sim.c index 10f5e58d..77b5a73b 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -357,7 +357,7 @@ static void run_job(void *vwargs, int cookie) snprintf(wargs->image.filename, 255, "dummy.h5"); } - reflections = find_intersections(&wargs->image, cr, PMODEL_SPHERE); + reflections = find_intersections(&wargs->image, cr, PMODEL_SCSPHERE); crystal_set_reflections(cr, reflections); for ( i=0; i + * 2010-2014 Thomas White * * This file is part of CrystFEL. * @@ -71,12 +71,8 @@ struct prdata extern struct prdata pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel); -/* Exported so it can be poked by tests/pr_gradient_check */ +/* Exported so it can be poked by tests/pr_p_gradient_check */ extern double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel); - -extern double l_gradient(Crystal *cr, int k, Reflection *refl, - PartialityModel pmodel); - #endif /* POST_REFINEMENT_H */ diff --git a/src/process_image.c b/src/process_image.c index 83973e50..1d1dac82 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -178,7 +178,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, /* Integrate all the crystals at once - need all the crystals so that * overlaps can be detected. */ - integrate_all_4(&image, iargs->int_meth, PMODEL_SPHERE, iargs->push_res, + integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, iargs->push_res, iargs->ir_inn, iargs->ir_mid, iargs->ir_out, iargs->int_diag, iargs->int_diag_h, iargs->int_diag_k, iargs->int_diag_l, results_pipe); diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c deleted file mode 100644 index b18196c9..00000000 --- a/tests/pr_l_gradient_check.c +++ /dev/null @@ -1,379 +0,0 @@ -/* - * pr_l_gradient_check.c - * - * Check Lorentz factor gradients for post refinement - * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2012-2014 Thomas White - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see . - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include "../src/post-refinement.h" - - -static void scan_partialities(RefList *reflections, RefList *compare, - int *valid, long double *vals[3], int idx) -{ - int i; - Reflection *refl; - RefListIterator *iter; - - i = 0; - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - Reflection *refl2; - - get_indices(refl, &h, &k, &l); - refl2 = find_refl(compare, h, k, l); - if ( refl2 == NULL ) { - valid[i] = 0; - i++; - continue; - } - - vals[idx][i] = get_lorentz(refl2); - i++; - } -} - - -static void shift_parameter(struct image *image, int k, double shift) -{ - switch ( k ) - { - case REF_DIV : - image->div += shift; - break; - - default: - ERROR("This test can't check parameter %i\n", k); - exit(1); - } -} - - -static void calc_either_side(Crystal *cr, double incr_val, - int *valid, long double *vals[3], int refine, - PartialityModel pmodel) -{ - RefList *compare; - struct image *image = crystal_get_image(cr); - struct image im_moved; - - im_moved = *image; - shift_parameter(&im_moved, refine, -incr_val); - compare = find_intersections(&im_moved, cr, pmodel); - scan_partialities(crystal_get_reflections(cr), compare, - valid, vals, 0); - reflist_free(compare); - - im_moved = *image; - shift_parameter(&im_moved, refine, +incr_val); - compare = find_intersections(&im_moved, cr, pmodel); - scan_partialities(crystal_get_reflections(cr), compare, - valid, vals, 2); - reflist_free(compare); -} - - - -static double test_gradients(Crystal *cr, double incr_val, int refine, - const char *str, const char *file, - PartialityModel pmodel, int quiet, int plot) -{ - Reflection *refl; - RefListIterator *iter; - long double *vals[3]; - int i; - int *valid; - int nref; - int n_good, n_invalid, n_small, n_nan, n_bad; - RefList *reflections; - FILE *fh = NULL; - int ntot = 0; - double total = 0.0; - char tmp[32]; - double *vec1; - double *vec2; - int n_line; - double cc; - - reflections = find_intersections(crystal_get_image(cr), cr, pmodel); - crystal_set_reflections(cr, reflections); - - nref = num_reflections(reflections); - if ( nref < 10 ) { - ERROR("Too few reflections found. Failing test by default.\n"); - return 0.0; - } - - vals[0] = malloc(nref*sizeof(long double)); - vals[1] = malloc(nref*sizeof(long double)); - vals[2] = malloc(nref*sizeof(long double)); - if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) { - ERROR("Couldn't allocate memory.\n"); - return 0.0; - } - - valid = malloc(nref*sizeof(int)); - if ( valid == NULL ) { - ERROR("Couldn't allocate memory.\n"); - return 0.0; - } - for ( i=0; ipanels[0].res = 13333.3; - image.det->panels[0].clen = 80e-3; - image.det->panels[0].coffset = 0.0; - - image.lambda = ph_en_to_lambda(eV_to_J(8000.0)); - image.div = 1e-3; - image.bw = 0.01; - image.filename = malloc(256); - - cr = crystal_new(); - if ( cr == NULL ) { - ERROR("Failed to allocate crystal.\n"); - return 1; - } - crystal_set_mosaicity(cr, 0.0); - crystal_set_profile_radius(cr, 0.005e9); - crystal_set_image(cr, &image); - - cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9, - deg2rad(90.0), - deg2rad(90.0), - deg2rad(90.0)); - - rng = gsl_rng_alloc(gsl_rng_mt19937); - - for ( i=0; i<3; i++ ) { - - UnitCell *rot; - double val; - PartialityModel pmodel; - - if ( i == 0 ) { - pmodel = PMODEL_SPHERE; - STATUS("Testing flat sphere model:\n"); - } else if ( i == 1 ) { - pmodel = PMODEL_GAUSSIAN; - STATUS("Testing Gaussian model:\n"); - } else if ( i == 2 ) { - pmodel = PMODEL_THIN; - STATUS("No need to test thin Ewald sphere model.\n"); - continue; - } else { - ERROR("WTF?\n"); - return 1; - } - - orientation = random_quaternion(rng); - rot = cell_rotate(cell, orientation); - crystal_set_cell(cr, rot); - - cell_get_reciprocal(rot, - &ax, &ay, &az, &bx, &by, - &bz, &cx, &cy, &cz); - - incr_val = incr_frac * image.div; - val = test_gradients(cr, incr_val, REF_DIV, "div", "div", - pmodel, quiet, plot); - if ( val < 0.99 ) fail = 1; - - } - - gsl_rng_free(rng); - - return fail; -} diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index 82ed2b02..b224f9b4 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -71,21 +71,6 @@ static void scan_partialities(RefList *reflections, RefList *compare, } get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high); - if ( clamp_low && clamp_high && (pmodel != PMODEL_SCSPHERE) ) { - if ( !within_tolerance(p, 1.0, 0.001) ) { - - signed int h, k, l; - - get_indices(refl, &h, &k, &l); - - ERROR("%3i %3i %3i - double clamped but" - " partiality not close to 1.0 (%5.2f)\n", - h, k, l, p); - - } - valid[i] = 0; - } - vals[idx][i] = p; i++; } @@ -451,26 +436,18 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i<4; i++ ) { + for ( i=0; i<2; i++ ) { UnitCell *rot; double val; PartialityModel pmodel; if ( i == 0 ) { - pmodel = PMODEL_SPHERE; - STATUS("Testing flat sphere model:\n"); - } else if ( i == 1 ) { - pmodel = PMODEL_GAUSSIAN; - /* FIXME: Gradients for Gaussian model are not good */ - STATUS("NOT testing Gaussian model.\n"); - continue; - } else if ( i == 2 ) { - pmodel = PMODEL_THIN; - STATUS("Testing Thin Ewald Sphere model:\n"); - } else if ( i == 3 ) { pmodel = PMODEL_SCSPHERE; STATUS("Testing SCSphere model:\n"); + } else if ( i == 1 ) { + pmodel = PMODEL_SCGAUSSIAN; + STATUS("Testing SCGaussian model.\n"); } else { ERROR("WTF?\n"); return 1; diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c deleted file mode 100644 index cb68e0e3..00000000 --- a/tests/pr_pl_gradient_check.c +++ /dev/null @@ -1,540 +0,0 @@ -/* - * pr_p_gradient_check.c - * - * Check partiality gradients for post refinement - * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2012-2014 Thomas White - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see . - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include "../src/post-refinement.h" - - -static void scan_partialities(RefList *reflections, RefList *compare, - int *valid, long double *vals[3], int idx) -{ - int i; - Reflection *refl; - RefListIterator *iter; - - i = 0; - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - Reflection *refl2; - double r1, r2, p; - int clamp_low, clamp_high; - - get_indices(refl, &h, &k, &l); - refl2 = find_refl(compare, h, k, l); - if ( refl2 == NULL ) { - valid[i] = 0; - i++; - continue; - } - - get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high); - if ( clamp_low && clamp_high ) { - if ( !within_tolerance(p, 1.0, 0.001) ) { - - signed int h, k, l; - - get_indices(refl, &h, &k, &l); - - ERROR("%3i %3i %3i - double clamped but" - " partiality not close to 1.0 (%5.2f)\n", - h, k, l, p); - - } - valid[i] = 0; - } - - vals[idx][i] = p * get_lorentz(refl2); - i++; - } -} - - -static UnitCell *new_shifted_cell(UnitCell *input, int k, double shift) -{ - UnitCell *cell; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - - cell = cell_new(); - cell_get_reciprocal(input, &asx, &asy, &asz, &bsx, &bsy, &bsz, - &csx, &csy, &csz); - switch ( k ) - { - case REF_ASX : asx += shift; break; - case REF_ASY : asy += shift; break; - case REF_ASZ : asz += shift; break; - case REF_BSX : bsx += shift; break; - case REF_BSY : bsy += shift; break; - case REF_BSZ : bsz += shift; break; - case REF_CSX : csx += shift; break; - case REF_CSY : csy += shift; break; - case REF_CSZ : csz += shift; break; - } - cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - - return cell; -} - - -static void shift_parameter(struct image *image, int k, double shift) -{ - switch ( k ) - { - case REF_DIV : image->div += shift; break; - } -} - - -static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val) -{ - Crystal *cr_new; - double r; - UnitCell *cell; - - cr_new = crystal_copy(cr); - if ( cr_new == NULL ) { - ERROR("Failed to allocate crystal.\n"); - return NULL; - } - - crystal_set_image(cr_new, crystal_get_image(cr)); - r = crystal_get_profile_radius(cr_new); - - switch ( refine ) { - - case REF_ASX : - case REF_ASY : - case REF_ASZ : - case REF_BSX : - case REF_BSY : - case REF_BSZ : - case REF_CSX : - case REF_CSY : - case REF_CSZ : - cell = new_shifted_cell(crystal_get_cell(cr), refine, - incr_val); - crystal_set_cell(cr_new, cell); - break; - - case REF_R : - cell = cell_new_from_cell(crystal_get_cell(cr)); - crystal_set_cell(cr_new, cell); - crystal_set_profile_radius(cr_new, r + incr_val); - break; - - default : - ERROR("Can't shift %i\n", refine); - break; - - } - - return cr_new; -} - - -static void calc_either_side(Crystal *cr, double incr_val, - int *valid, long double *vals[3], int refine, - PartialityModel pmodel) -{ - RefList *compare; - struct image *image = crystal_get_image(cr); - - if ( (refine != REF_DIV) ) { - - Crystal *cr_new; - - /* Crystal properties */ - cr_new = new_shifted_crystal(cr, refine, -incr_val); - compare = find_intersections(image, cr_new, pmodel); - scan_partialities(crystal_get_reflections(cr), compare, valid, - vals, 0); - cell_free(crystal_get_cell(cr_new)); - crystal_free(cr_new); - reflist_free(compare); - - cr_new = new_shifted_crystal(cr, refine, +incr_val); - compare = find_intersections(image, cr_new, pmodel); - scan_partialities(crystal_get_reflections(cr), compare, valid, - vals, 2); - cell_free(crystal_get_cell(cr_new)); - crystal_free(cr_new); - reflist_free(compare); - - } else { - - struct image im_moved; - - /* "Image" properties */ - im_moved = *image; - shift_parameter(&im_moved, refine, -incr_val); - compare = find_intersections(&im_moved, cr, pmodel); - scan_partialities(crystal_get_reflections(cr), compare, - valid, vals, 0); - reflist_free(compare); - - im_moved = *image; - shift_parameter(&im_moved, refine, +incr_val); - compare = find_intersections(&im_moved, cr, pmodel); - scan_partialities(crystal_get_reflections(cr), compare, - valid, vals, 2); - reflist_free(compare); - - } -} - - -static double test_gradients(Crystal *cr, double incr_val, int refine, - const char *str, const char *file, - PartialityModel pmodel, int quiet, int plot) -{ - Reflection *refl; - RefListIterator *iter; - long double *vals[3]; - int i; - int *valid; - int nref; - int n_good, n_invalid, n_small, n_nan, n_bad; - RefList *reflections; - FILE *fh = NULL; - int ntot = 0; - double total = 0.0; - char tmp[32]; - double *vec1; - double *vec2; - int n_line; - double cc; - - reflections = find_intersections(crystal_get_image(cr), cr, pmodel); - crystal_set_reflections(cr, reflections); - - nref = num_reflections(reflections); - if ( nref < 10 ) { - ERROR("Too few reflections found. Failing test by default.\n"); - return 0.0; - } - - vals[0] = malloc(nref*sizeof(long double)); - vals[1] = malloc(nref*sizeof(long double)); - vals[2] = malloc(nref*sizeof(long double)); - if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) { - ERROR("Couldn't allocate memory.\n"); - return 0.0; - } - - valid = malloc(nref*sizeof(int)); - if ( valid == NULL ) { - ERROR("Couldn't allocate memory.\n"); - return 0.0; - } - for ( i=0; ipanels[0].res = 13333.3; - image.det->panels[0].clen = 80e-3; - image.det->panels[0].coffset = 0.0; - - image.lambda = ph_en_to_lambda(eV_to_J(8000.0)); - image.div = 1e-3; - image.bw = 0.01; - image.filename = malloc(256); - - cr = crystal_new(); - if ( cr == NULL ) { - ERROR("Failed to allocate crystal.\n"); - return 1; - } - crystal_set_mosaicity(cr, 0.0); - crystal_set_profile_radius(cr, 0.005e9); - crystal_set_image(cr, &image); - - cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9, - deg2rad(90.0), - deg2rad(90.0), - deg2rad(90.0)); - - rng = gsl_rng_alloc(gsl_rng_mt19937); - - for ( i=0; i<3; i++ ) { - - UnitCell *rot; - double val; - PartialityModel pmodel; - - if ( i == 0 ) { - pmodel = PMODEL_SPHERE; - STATUS("Testing flat sphere model:\n"); - } else if ( i == 1 ) { - pmodel = PMODEL_GAUSSIAN; - /* FIXME: Gradients for Gaussian model are not good */ - STATUS("NOT testing Gaussian model.\n"); - continue; - } else if ( i == 2 ) { - pmodel = PMODEL_THIN; - STATUS("No need to test thin Ewald sphere model.\n"); - continue; - } else { - ERROR("WTF?\n"); - return 1; - } - - orientation = random_quaternion(rng); - rot = cell_rotate(cell, orientation); - crystal_set_cell(cr, rot); - - cell_get_reciprocal(rot, - &ax, &ay, &az, &bx, &by, - &bz, &cx, &cy, &cz); - - incr_val = incr_frac * image.div; - val = test_gradients(cr, incr_val, REF_DIV, "div", "div", - pmodel, quiet, plot); - if ( val < 0.99 ) fail = 1; - - incr_val = incr_frac * crystal_get_profile_radius(cr); - val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - - incr_val = incr_frac * ax; - val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - incr_val = incr_frac * bx; - val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - incr_val = incr_frac * cx; - val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - - incr_val = incr_frac * ay; - val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - incr_val = incr_frac * by; - val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - incr_val = incr_frac * cy; - val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - - incr_val = incr_frac * az; - val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - incr_val = incr_frac * bz; - val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - incr_val = incr_frac * cz; - val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel, - quiet, plot); - if ( val < 0.99 ) fail = 1; - - } - - gsl_rng_free(rng); - - return fail; -} diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c index 8563c7d9..c7ebd8b2 100644 --- a/tests/prof2d_check.c +++ b/tests/prof2d_check.c @@ -138,7 +138,7 @@ int main(int argc, char *argv[]) image.n_crystals = 1; image.crystals = &cr; - list = find_intersections(&image, cr, PMODEL_SPHERE); + list = find_intersections(&image, cr, PMODEL_SCSPHERE); crystal_set_reflections(cr, list); for ( fs=0; fs