From 538e990e4667c45e537755799d82ab8a56fdd1f4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 20 Jun 2023 11:10:25 +0200 Subject: Move gradients to predict-refine Makes sense for the _dev and _gradient functions to be in the same place., --- libcrystfel/src/crystfel-mille.h | 1 - libcrystfel/src/geometry.c | 192 -------------------------------------- libcrystfel/src/geometry.h | 36 -------- libcrystfel/src/predict-refine.c | 194 ++++++++++++++++++++++++++++++++++++++- libcrystfel/src/predict-refine.h | 41 +++++++++ 5 files changed, 234 insertions(+), 230 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/crystfel-mille.h b/libcrystfel/src/crystfel-mille.h index 00d0415d..bdd66f37 100644 --- a/libcrystfel/src/crystfel-mille.h +++ b/libcrystfel/src/crystfel-mille.h @@ -34,7 +34,6 @@ typedef void *Mille; #include "cell.h" #include "image.h" #include "predict-refine.h" -#include "geometry.h" /** * \file crystfel-mille.h diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 5289856b..9dc471ee 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -429,69 +429,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } -double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double xl, yl, zl; - signed int hs, ks, ls; - double tl, phi, azi; - - get_symmetric_indices(refl, &hs, &ks, &ls); - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - xl = hs*asx + ks*bsx + ls*csx; - yl = hs*asy + ks*bsy + ls*csy; - zl = hs*asz + ks*bsz + ls*csz; - - tl = sqrt(xl*xl + yl*yl); - phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */ - azi = atan2(yl, xl); /* azimuth */ - - switch ( k ) { - - case GPARAM_ASX : - return - hs * sin(phi) * cos(azi); - - case GPARAM_BSX : - return - ks * sin(phi) * cos(azi); - - case GPARAM_CSX : - return - ls * sin(phi) * cos(azi); - - case GPARAM_ASY : - return - hs * sin(phi) * sin(azi); - - case GPARAM_BSY : - return - ks * sin(phi) * sin(azi); - - case GPARAM_CSY : - return - ls * sin(phi) * sin(azi); - - case GPARAM_ASZ : - return - hs * cos(phi); - - case GPARAM_BSZ : - return - ks * cos(phi); - - case GPARAM_CSZ : - return - ls * cos(phi); - - case GPARAM_DETX : - case GPARAM_DETY : - case GPARAM_CLEN : - return 0.0; - - } - - ERROR("No r gradient defined for parameter %i\n", k); - abort(); -} - - /** * \param cryst: A \ref Crystal * \param max_res: Maximum resolution to predict to (m^-1) @@ -1084,132 +1021,3 @@ void polarisation_correction(RefList *list, UnitCell *cell, set_esd_intensity(refl, sigma / pol); } } - - -/* Returns dfs_refl/dP, where P = any parameter * - * fs_refl is the fast scan position of refl in panel 'p', - * in metres (not pixels) */ -double fs_gradient(int param, Reflection *refl, UnitCell *cell, - struct detgeom_panel *p) -{ - signed int h, k, l; - double xl, zl, kpred; - double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; - - get_indices(refl, &h, &k, &l); - kpred = get_kpred(refl); - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - xl = h*asx + k*bsx + l*csx; - zl = h*asz + k*bsz + l*csz; - - switch ( param ) { - - case GPARAM_ASX : - return h * p->cnz * p->pixel_pitch / (kpred + zl); - - case GPARAM_BSX : - return k * p->cnz * p->pixel_pitch / (kpred + zl); - - case GPARAM_CSX : - return l * p->cnz * p->pixel_pitch / (kpred + zl); - - case GPARAM_ASY : - return 0.0; - - case GPARAM_BSY : - return 0.0; - - case GPARAM_CSY : - return 0.0; - - case GPARAM_ASZ : - return -h * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); - - case GPARAM_BSZ : - return -k * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); - - case GPARAM_CSZ : - return -l * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); - - case GPARAM_DETX : - return -1; - - case GPARAM_DETY : - return 0; - - case GPARAM_CLEN : - return -xl / (kpred+zl); - - default: - ERROR("Positional gradient requested for parameter %i?\n", param); - abort(); - - } - -} - - -/* Returns dss_refl/dP, where P = any parameter - * ss_refl is the slow scan position of refl in panel 'p', - * in metres (not pixels) */ -double ss_gradient(int param, Reflection *refl, UnitCell *cell, - struct detgeom_panel *p) -{ - signed int h, k, l; - double yl, zl, kpred; - double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; - - get_indices(refl, &h, &k, &l); - kpred = get_kpred(refl); - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - yl = h*asy + k*bsy + l*csy; - zl = h*asz + k*bsz + l*csz; - - switch ( param ) { - - case GPARAM_ASX : - return 0.0; - - case GPARAM_BSX : - return 0.0; - - case GPARAM_CSX : - return 0.0; - - case GPARAM_ASY : - return h * p->cnz * p->pixel_pitch / (kpred + zl); - - case GPARAM_BSY : - return k * p->cnz * p->pixel_pitch / (kpred + zl); - - case GPARAM_CSY : - return l * p->cnz * p->pixel_pitch / (kpred + zl); - - case GPARAM_ASZ : - return -h * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); - - case GPARAM_BSZ : - return -k * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); - - case GPARAM_CSZ : - return -l * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); - - case GPARAM_DETX : - return 0; - - case GPARAM_DETY : - return -1; - - case GPARAM_CLEN : - return -yl / (kpred+zl); - - default : - ERROR("Positional gradient requested for parameter %i?\n", param); - abort(); - - } -} diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index a7057f12..d05c8832 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -62,35 +62,6 @@ typedef enum { } PartialityModel; -/** Enumeration of parameters which may want to be refined */ -enum gparam { - GPARAM_ASX, - GPARAM_ASY, - GPARAM_ASZ, - GPARAM_BSX, - GPARAM_BSY, - GPARAM_BSZ, - GPARAM_CSX, - GPARAM_CSY, - GPARAM_CSZ, - GPARAM_R, - GPARAM_DIV, - GPARAM_DETX, - GPARAM_DETY, - GPARAM_CLEN, - GPARAM_OSF, /* Linear scale factor */ - GPARAM_BFAC, /* D-W scale factor */ - GPARAM_ANG1, /* Out of plane rotation angles of crystal */ - GPARAM_ANG2, - GPARAM_WAVELENGTH, - GPARAM_ROTX, /* Detector panel rotation about +x */ - GPARAM_ROTY, /* Detector panel rotation about +y */ - GPARAM_ROTZ, /* Detector panel rotation about +z */ - - GPARAM_EOL /* End of list */ -}; - - /** * This structure represents the polarisation of the incident radiation */ @@ -110,8 +81,6 @@ extern RefList *predict_to_res(Crystal *cryst, double max_res); extern void calculate_partialities(Crystal *cryst, PartialityModel pmodel); -extern double r_gradient(UnitCell *cell, int k, Reflection *refl, - struct image *image); extern void update_predictions(Crystal *cryst); extern struct polarisation parse_polarisation(const char *text); extern void polarisation_correction(RefList *list, UnitCell *cell, @@ -120,11 +89,6 @@ extern void polarisation_correction(RefList *list, UnitCell *cell, extern double sphere_fraction(double rlow, double rhigh, double pr); extern double gaussian_fraction(double rlow, double rhigh, double pr); -extern double fs_gradient(int param, Reflection *refl, UnitCell *cell, - struct detgeom_panel *p); -extern double ss_gradient(int param, Reflection *refl, UnitCell *cell, - struct detgeom_panel *p); - #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index ad7f5584..27f97700 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -68,7 +68,7 @@ static const enum gparam rv[] = }; -static double r_dev(struct reflpeak *rp) +double r_dev(struct reflpeak *rp) { /* Excitation error term */ return get_exerr(rp->refl); @@ -91,6 +91,198 @@ double ss_dev(struct reflpeak *rp, struct detgeom *det) } +double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double xl, yl, zl; + signed int hs, ks, ls; + double tl, phi, azi; + + get_symmetric_indices(refl, &hs, &ks, &ls); + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + xl = hs*asx + ks*bsx + ls*csx; + yl = hs*asy + ks*bsy + ls*csy; + zl = hs*asz + ks*bsz + ls*csz; + + tl = sqrt(xl*xl + yl*yl); + phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */ + azi = atan2(yl, xl); /* azimuth */ + + switch ( k ) { + + case GPARAM_ASX : + return - hs * sin(phi) * cos(azi); + + case GPARAM_BSX : + return - ks * sin(phi) * cos(azi); + + case GPARAM_CSX : + return - ls * sin(phi) * cos(azi); + + case GPARAM_ASY : + return - hs * sin(phi) * sin(azi); + + case GPARAM_BSY : + return - ks * sin(phi) * sin(azi); + + case GPARAM_CSY : + return - ls * sin(phi) * sin(azi); + + case GPARAM_ASZ : + return - hs * cos(phi); + + case GPARAM_BSZ : + return - ks * cos(phi); + + case GPARAM_CSZ : + return - ls * cos(phi); + + case GPARAM_DETX : + case GPARAM_DETY : + case GPARAM_CLEN : + return 0.0; + + } + + ERROR("No r gradient defined for parameter %i\n", k); + abort(); +} + + +/* Returns dfs_refl/dP, where P = any parameter * + * fs_refl is the fast scan position of refl in panel 'p', + * in metres (not pixels) */ +double fs_gradient(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p) +{ + signed int h, k, l; + double xl, zl, kpred; + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + + get_indices(refl, &h, &k, &l); + kpred = get_kpred(refl); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + xl = h*asx + k*bsx + l*csx; + zl = h*asz + k*bsz + l*csz; + + switch ( param ) { + + case GPARAM_ASX : + return h * p->cnz * p->pixel_pitch / (kpred + zl); + + case GPARAM_BSX : + return k * p->cnz * p->pixel_pitch / (kpred + zl); + + case GPARAM_CSX : + return l * p->cnz * p->pixel_pitch / (kpred + zl); + + case GPARAM_ASY : + return 0.0; + + case GPARAM_BSY : + return 0.0; + + case GPARAM_CSY : + return 0.0; + + case GPARAM_ASZ : + return -h * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); + + case GPARAM_BSZ : + return -k * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); + + case GPARAM_CSZ : + return -l * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); + + case GPARAM_DETX : + return -1; + + case GPARAM_DETY : + return 0; + + case GPARAM_CLEN : + return -xl / (kpred+zl); + + default: + ERROR("Positional gradient requested for parameter %i?\n", param); + abort(); + + } + +} + + +/* Returns dss_refl/dP, where P = any parameter + * ss_refl is the slow scan position of refl in panel 'p', + * in metres (not pixels) */ +double ss_gradient(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p) +{ + signed int h, k, l; + double yl, zl, kpred; + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + + get_indices(refl, &h, &k, &l); + kpred = get_kpred(refl); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; + + switch ( param ) { + + case GPARAM_ASX : + return 0.0; + + case GPARAM_BSX : + return 0.0; + + case GPARAM_CSX : + return 0.0; + + case GPARAM_ASY : + return h * p->cnz * p->pixel_pitch / (kpred + zl); + + case GPARAM_BSY : + return k * p->cnz * p->pixel_pitch / (kpred + zl); + + case GPARAM_CSY : + return l * p->cnz * p->pixel_pitch / (kpred + zl); + + case GPARAM_ASZ : + return -h * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); + + case GPARAM_BSZ : + return -k * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); + + case GPARAM_CSZ : + return -l * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl); + + case GPARAM_DETX : + return 0; + + case GPARAM_DETY : + return -1; + + case GPARAM_CLEN : + return -yl / (kpred+zl); + + default : + ERROR("Positional gradient requested for parameter %i?\n", param); + abort(); + + } +} + + static int cmpd2(const void *av, const void *bv) { struct reflpeak *a, *b; diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h index e74e00a8..c667dc6d 100644 --- a/libcrystfel/src/predict-refine.h +++ b/libcrystfel/src/predict-refine.h @@ -30,6 +30,36 @@ #define PREDICT_REFINE_H struct reflpeak; + +/** Enumeration of parameters which may want to be refined */ +enum gparam { + GPARAM_ASX, + GPARAM_ASY, + GPARAM_ASZ, + GPARAM_BSX, + GPARAM_BSY, + GPARAM_BSZ, + GPARAM_CSX, + GPARAM_CSY, + GPARAM_CSZ, + GPARAM_R, + GPARAM_DIV, + GPARAM_DETX, + GPARAM_DETY, + GPARAM_CLEN, + GPARAM_OSF, /* Linear scale factor */ + GPARAM_BFAC, /* D-W scale factor */ + GPARAM_ANG1, /* Out of plane rotation angles of crystal */ + GPARAM_ANG2, + GPARAM_WAVELENGTH, + GPARAM_ROTX, /* Detector panel rotation about +x */ + GPARAM_ROTY, /* Detector panel rotation about +y */ + GPARAM_ROTZ, /* Detector panel rotation about +z */ + + GPARAM_EOL /* End of list */ +}; + + #include "crystal.h" #include "crystfel-mille.h" @@ -50,8 +80,19 @@ extern int refine_prediction(struct image *image, Crystal *cr, Mille *mille); extern int refine_radius(Crystal *cr, struct image *image); +extern double r_dev(struct reflpeak *rp); + extern double fs_dev(struct reflpeak *rp, struct detgeom *det); extern double ss_dev(struct reflpeak *rp, struct detgeom *det); +extern double r_gradient(UnitCell *cell, int k, Reflection *refl, + struct image *image); + +extern double fs_gradient(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p); + +extern double ss_gradient(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p); + #endif /* PREDICT_REFINE_H */ -- cgit v1.2.3