diff options
-rw-r--r-- | libcrystfel/src/crystfel-mille.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 37 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.h | 1 | ||||
-rw-r--r-- | tests/gradient_check.c | 51 | ||||
-rw-r--r-- | tests/meson.build | 22 |
5 files changed, 95 insertions, 20 deletions
diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c index 1bc126ed..93c03028 100644 --- a/libcrystfel/src/crystfel-mille.c +++ b/libcrystfel/src/crystfel-mille.c @@ -93,7 +93,7 @@ void write_mille(Mille *mille, int n, UnitCell *cell, for ( j=0; j<9; j++ ) { fs_ss_gradient(rv[j], rps[i].refl, cell, &image->detgeom->panels[rps[i].peak->pn], - Minvs[rps[i].peak->pn], + Minvs[rps[i].peak->pn], 0, 0, 0, &local_gradients_fs[j], &local_gradients_ss[j]); } @@ -105,7 +105,7 @@ void write_mille(Mille *mille, int n, UnitCell *cell, fs_ss_gradient(GPARAM_DET_TX, rps[i].refl, cell, &image->detgeom->panels[rps[i].peak->pn], - Minvs[rps[i].peak->pn], + Minvs[rps[i].peak->pn], 0, 0, 0, &global_gradients_fs[j], &global_gradients_ss[j]); diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 0e8754ed..d8dbbca0 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -239,7 +239,7 @@ int fs_ss_gradient_physics(int param, Reflection *refl, UnitCell *cell, int fs_ss_gradient_panel(int param, Reflection *refl, UnitCell *cell, struct detgeom_panel *p, gsl_matrix *Minv, double fs, double ss, double mu, - gsl_vector *t, + gsl_vector *t, double cx, double cy, double cz, float *fsg, float *ssg) { gsl_vector *v; @@ -261,19 +261,31 @@ int fs_ss_gradient_panel(int param, Reflection *refl, UnitCell *cell, break; case GPARAM_DET_RX : - *fsg = 0.0; - *ssg = 0.0; - return 0; + gsl_matrix_set(dMdp, 1, 0, p->cnz-cz); + gsl_matrix_set(dMdp, 2, 0, cy-p->cny); + gsl_matrix_set(dMdp, 1, 1, p->fsz); + gsl_matrix_set(dMdp, 2, 1, -p->fsy); + gsl_matrix_set(dMdp, 1, 2, p->ssz); + gsl_matrix_set(dMdp, 2, 2, -p->ssy); + break; case GPARAM_DET_RY : - *fsg = 0.0; - *ssg = 0.0; - return 0; + gsl_matrix_set(dMdp, 0, 0, cz-p->cnz); + gsl_matrix_set(dMdp, 2, 0, p->cnx-cx); + gsl_matrix_set(dMdp, 0, 1, -p->fsz); + gsl_matrix_set(dMdp, 2, 1, p->fsx); + gsl_matrix_set(dMdp, 0, 2, -p->ssz); + gsl_matrix_set(dMdp, 2, 2, p->ssx); + break; case GPARAM_DET_RZ : - *fsg = 0.0; - *ssg = 0.0; - return 0; + gsl_matrix_set(dMdp, 0, 0, cy-p->cny); + gsl_matrix_set(dMdp, 1, 0, p->cnx-cx); + gsl_matrix_set(dMdp, 0, 1, -p->fsy); + gsl_matrix_set(dMdp, 1, 1, p->fsx); + gsl_matrix_set(dMdp, 0, 2, -p->ssy); + gsl_matrix_set(dMdp, 1, 2, p->ssx); + break; default: ERROR("Invalid panel gradient %i\n", param); @@ -303,6 +315,7 @@ int fs_ss_gradient_panel(int param, Reflection *refl, UnitCell *cell, * in metres (not pixels) */ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, struct detgeom_panel *p, gsl_matrix *Minv, + double cx, double cy, double cz, float *fsg, float *ssg) { signed int h, k, l; @@ -371,7 +384,7 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, } else { return fs_ss_gradient_panel(param, refl, cell, p, Minv, fs, ss, mu, t, - fsg, ssg); + cx, cy, cz, fsg, ssg); } } @@ -648,7 +661,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, image->lambda); fs_ss_gradient(rv[k], rps[i].refl, cell, &image->detgeom->panels[rps[i].peak->pn], - Minvs[rps[i].peak->pn], + Minvs[rps[i].peak->pn], 0, 0, 0, &fs_gradients[k], &ss_gradients[k]); } diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h index e2f912d0..604799c0 100644 --- a/libcrystfel/src/predict-refine.h +++ b/libcrystfel/src/predict-refine.h @@ -82,6 +82,7 @@ extern double r_gradient(int param, Reflection *refl, UnitCell *cell, extern int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, struct detgeom_panel *p, gsl_matrix *panel_Minv, + double cx, double cy, double cz, float *fsg, float *ssg); #endif /* PREDICT_REFINE_H */ diff --git a/tests/gradient_check.c b/tests/gradient_check.c index 46d366d6..7396b8ad 100644 --- a/tests/gradient_check.c +++ b/tests/gradient_check.c @@ -52,15 +52,47 @@ int main(int argc, char *argv[]) int fail = 0; double step; gsl_matrix **panel_matrices; + int didsomething = 0; + const double cx = -200.0; + const double cy = +100.0; + const double cz = -50.0; rps = make_test_image(&n_refls, &image); panel_matrices = make_panel_minvs(image.detgeom); before = make_dev_list(rps, n_refls, image.detgeom); - #ifdef MOVE_PANEL + #ifdef TRANSLATE_PANEL step = 0.1; /* Pixels */ image.detgeom->panels[0].THING_TO_MOVE += step; + didsomething = 1; + #endif + + #ifdef ROTATE_PANEL_X + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = deg2rad(0.01); + rotate2d(&p->cnz, &p->cny, cz, cy, step); + rotate2d(&p->fsz, &p->fsy, 0, 0, step); + rotate2d(&p->ssz, &p->ssy, 0, 0, step); + didsomething = 1; + #endif + + #ifdef ROTATE_PANEL_Y + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = deg2rad(0.01); + rotate2d(&p->cnx, &p->cnz, cx, cz, step); + rotate2d(&p->fsx, &p->fsz, 0, 0, step); + rotate2d(&p->ssx, &p->ssz, 0, 0, step); + didsomething = 1; + #endif + + #ifdef ROTATE_PANEL_Z + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = deg2rad(0.01); + rotate2d(&p->cnx, &p->cny, cx, cy, step); + rotate2d(&p->fsx, &p->fsy, 0, 0, step); + rotate2d(&p->ssx, &p->ssy, 0, 0, step); + didsomething = 1; #endif #ifdef CHANGE_CELL @@ -74,8 +106,14 @@ int main(int argc, char *argv[]) cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + didsomething = 1; #endif + if ( !didsomething ) { + fprintf(stderr, "Nothing changed. Check the build system.\n"); + return 1; + } + update_predictions(image.crystals[0]); after = make_dev_list(rps, n_refls, image.detgeom); @@ -91,20 +129,27 @@ int main(int argc, char *argv[]) fs_ss_gradient(TEST_GPARAM, rps[i].refl, crystal_get_cell(image.crystals[0]), &image.detgeom->panels[rps[i].peak->pn], - panel_matrices[rps[i].peak->pn], + panel_matrices[rps[i].peak->pn], cx, cy, cz, &calc[1], &calc[2]); obs[0] = (after[0][i] - before[0][i]) / step; obs[1] = (after[1][i] - before[1][i]) / step; obs[2] = (after[2][i] - before[2][i]) / step; - #ifdef MOVE_PANEL + #ifdef TRANSLATE_PANEL if ( fabs(calc[0]) > 1e-12 ) n_wrong_r++; /* Should be zero */ if ( fabs(obs[0]) > 1e-12 ) n_wrong_obsr++; /* Should also be zero */ if ( fabs(obs[1] - calc[1]) > 1e-3 ) n_wrong_fs++; if ( fabs(obs[2] - calc[2]) > 1e-3 ) n_wrong_ss++; #endif + #if defined(ROTATE_PANEL_X) || defined(ROTATE_PANEL_Y) || defined(ROTATE_PANEL_Z) + if ( fabs(calc[0]) > 1e-12 ) n_wrong_r++; /* Should be zero */ + if ( fabs(obs[0]) > 1e-12 ) n_wrong_obsr++; /* Should also be zero */ + if ( fabs(obs[1] - calc[1]) > 1.0 ) n_wrong_fs++; /* Units are pixels/rad */ + if ( fabs(obs[2] - calc[2]) > 1.0 ) n_wrong_ss++; /* (numbers are big) */ + #endif + #ifdef CHANGE_CELL if ( fabs(obs[0] - calc[0]) > 1e-2 ) n_wrong_r++; if ( fabs(obs[1] - calc[1]) > 1e-8 ) n_wrong_fs++; diff --git a/tests/meson.build b/tests/meson.build index 424310f4..9dc245a9 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -96,7 +96,7 @@ exe = executable('prof2d_check', test('prof2d_check', exe) -# Refinement gradient checks, part 1: panel positions +# Refinement gradient checks, part 1: panel translations panel_gradient_tests = [['gradient_panel_x', 'cnx', 'GPARAM_DET_TX'], ['gradient_panel_y', 'cny', 'GPARAM_DET_TY'], ['gradient_panel_z', 'cnz', 'GPARAM_DET_TZ']] @@ -107,13 +107,29 @@ foreach name : panel_gradient_tests 'gradient_check_utils.c'], c_args : ['-DTHING_TO_MOVE='+name[1], '-DTEST_GPARAM='+name[2], - '-DMOVE_PANEL=1'], + '-DTRANSLATE_PANEL=1'], dependencies : [libcrystfeldep, mdep, gsldep]) test(name[0], exe) endforeach -# Refinement gradient checks, part 2: diffraction physics +# Refinement gradient checks, part 2: panel rotations +panel_gradient_tests = [['gradient_panel_rx', '-DROTATE_PANEL_X', 'GPARAM_DET_RX'], + ['gradient_panel_ry', '-DROTATE_PANEL_Y', 'GPARAM_DET_RY'], + ['gradient_panel_rz', '-DROTATE_PANEL_Z', 'GPARAM_DET_RZ']] + +foreach name : panel_gradient_tests + exe = executable(name[0], + ['gradient_check.c', + 'gradient_check_utils.c'], + c_args : [name[1], + '-DTEST_GPARAM='+name[2]], + dependencies : [libcrystfeldep, mdep, gsldep]) + test(name[0], exe) +endforeach + + +# Refinement gradient checks, part 3: diffraction physics panel_gradient_tests = [['gradient_cell_asx', 'asx', 'GPARAM_ASX'], ['gradient_cell_asy', 'asy', 'GPARAM_ASY'], ['gradient_cell_asz', 'asz', 'GPARAM_ASZ'], |