From 4ab308554c1f170545bf1dd9d9990f27fe1caf4e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 1 Aug 2013 16:06:58 +0200 Subject: Add test for gradient of Lorentz factor (as well as partiality) --- Makefile.am | 9 +- src/post-refinement.h | 3 + tests/.gitignore | 3 +- tests/pr_gradient_check.c | 511 ------------------------------------------- tests/pr_l_gradient_check.c | 358 ++++++++++++++++++++++++++++++ tests/pr_p_gradient_check.c | 515 ++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 884 insertions(+), 515 deletions(-) delete mode 100644 tests/pr_gradient_check.c create mode 100644 tests/pr_l_gradient_check.c create mode 100644 tests/pr_p_gradient_check.c diff --git a/Makefile.am b/Makefile.am index 8a17d502..04eb9629 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,9 +7,9 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ src/compare_hkl src/partialator src/check_hkl src/partial_sim noinst_PROGRAMS = tests/list_check tests/integration_check \ - tests/pr_gradient_check tests/symmetry_check \ + tests/pr_p_gradient_check tests/symmetry_check \ tests/centering_check tests/transformation_check \ - tests/cell_check + tests/cell_check tests/pr_l_gradient_check MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check @@ -86,7 +86,10 @@ tests_integration_check_SOURCES = tests/integration_check.c tests_symmetry_check_SOURCES = tests/symmetry_check.c -tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \ +tests_pr_p_gradient_check_SOURCES = tests/pr_p_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 diff --git a/src/post-refinement.h b/src/post-refinement.h index bdd1eb10..7e13090b 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -68,4 +68,7 @@ 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/tests/.gitignore b/tests/.gitignore index dbe97b34..1a5ece46 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -3,7 +3,8 @@ list_check gpu_sim_check integration_check -pr_gradient_check +pr_l_gradient_check +pr_p_gradient_check symmetry_check centering_check transformation_check diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c deleted file mode 100644 index 4600db9b..00000000 --- a/tests/pr_gradient_check.c +++ /dev/null @@ -1,511 +0,0 @@ -/* - * pr_gradient_check.c - * - * Check gradients for post refinement - * - * Copyright © 2012-2013 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; - 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) -{ - 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); - 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); - 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); - 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); - 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; - 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); - 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)); - - for ( i=0; i<1; i++ ) { - - UnitCell *rot; - double val; - - orientation = random_quaternion(); - 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; - - } - - return fail; -} diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c new file mode 100644 index 00000000..a4fd4526 --- /dev/null +++ b/tests/pr_l_gradient_check.c @@ -0,0 +1,358 @@ +/* + * pr_l_gradient_check.c + * + * Check Lorentz factor gradients for post refinement + * + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2013 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) +{ + 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); + 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); + 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; + 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); + 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)); + + for ( i=0; i<1; i++ ) { + + UnitCell *rot; + double val; + + orientation = random_quaternion(); + 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; + + } + + return fail; +} diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c new file mode 100644 index 00000000..9adf4780 --- /dev/null +++ b/tests/pr_p_gradient_check.c @@ -0,0 +1,515 @@ +/* + * pr_p_gradient_check.c + * + * Check partiality gradients for post refinement + * + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2013 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; + 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) +{ + 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); + 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); + 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); + 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); + 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; + 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); + 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)); + + for ( i=0; i<1; i++ ) { + + UnitCell *rot; + double val; + + orientation = random_quaternion(); + 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; + + } + + return fail; +} -- cgit v1.2.3