diff options
Diffstat (limited to 'libcrystfel/src/predict-refine.c')
-rw-r--r-- | libcrystfel/src/predict-refine.c | 194 |
1 files changed, 193 insertions, 1 deletions
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; |