aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-03-19 17:56:37 +0100
committerThomas White <taw@physics.org>2015-04-20 15:50:39 +0200
commit52fc1532ba4a93f63b82ed4f5382d637dccac3c5 (patch)
tree49656dc8a4787732f076294a69cab7961f398aca /libcrystfel
parent254c019a28b7be2b39674cb84e293be16f68ff4c (diff)
Factorise dr/da part of gradient calculation for PR and prediction
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/geometry.c81
-rw-r--r--libcrystfel/src/geometry.h19
2 files changed, 100 insertions, 0 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 35492fde..66786cdb 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -297,6 +297,87 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
+double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
+{
+ double azi;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xl, yl, zl;
+ signed int hs, ks, ls;
+ double rlow, rhigh, p;
+ double philow, phihigh, phi;
+ double khigh, klow;
+ double tl, cet, cez;
+
+ get_partial(refl, &rlow, &rhigh, &p);
+
+ 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;
+
+ /* "low" gives the largest Ewald sphere (wavelength short => k large)
+ * "high" gives the smallest Ewald sphere (wavelength long => k small)
+ */
+ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
+ khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
+
+ tl = sqrt(xl*xl + yl*yl);
+
+ cet = -sin(image->div/2.0) * klow;
+ cez = -cos(image->div/2.0) * klow;
+ philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
+
+ cet = -sin(image->div/2.0) * khigh;
+ cez = -cos(image->div/2.0) * khigh;
+ phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
+
+ /* Approximation: philow and phihigh are very similar */
+ phi = (philow + phihigh) / 2.0;
+
+ azi = atan2(yl, xl);
+
+ 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);
+
+ }
+
+ ERROR("No r gradient defined for parameter %i\n", k);
+ abort();
+}
+
+
RefList *find_intersections(struct image *image, Crystal *cryst,
PartialityModel pmodel)
{
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 162c0355..e0d21dda 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -61,12 +61,31 @@ typedef enum {
} PartialityModel;
+
+/* Enumeration of parameters which may want to be refined */
+enum {
+ GPARAM_ASX,
+ GPARAM_ASY,
+ GPARAM_ASZ,
+ GPARAM_BSX,
+ GPARAM_BSY,
+ GPARAM_BSZ,
+ GPARAM_CSX,
+ GPARAM_CSY,
+ GPARAM_CSZ,
+ GPARAM_R,
+ GPARAM_DIV,
+};
+
+
extern RefList *find_intersections(struct image *image, Crystal *cryst,
PartialityModel pmodel);
extern RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
PartialityModel pmodel,
double max_res);
+extern double r_gradient(UnitCell *cell, int k, Reflection *refl,
+ struct image *image);
extern void update_partialities(Crystal *cryst, PartialityModel pmodel);
extern void polarisation_correction(RefList *list, UnitCell *cell,
struct image *image);