aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/geometry.c130
-rw-r--r--libcrystfel/src/geometry.h14
-rw-r--r--libcrystfel/src/integration.c3
-rw-r--r--libcrystfel/src/predict-refine.c8
-rw-r--r--src/partial_sim.c3
-rw-r--r--src/partialator.c3
-rw-r--r--src/post-refinement.c9
-rw-r--r--tests/pr_p_gradient_check.c20
-rw-r--r--tests/prediction_gradient_check.c7
-rw-r--r--tests/prof2d_check.c6
10 files changed, 126 insertions, 77 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index b968dc06..576655de 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -266,15 +266,12 @@ static double partiality(PartialityModel pmodel,
static Reflection *check_reflection(struct image *image, Crystal *cryst,
- PartialityModel pmodel,
signed int h, signed int k, signed int l,
double xl, double yl, double zl,
Reflection *updateme)
{
- const int output = 0;
double tl;
- double rlow, rhigh; /* "Excitation error" */
- double part; /* Partiality */
+ double rlow, rhigh; /* "Excitation error" */
double klow, khigh; /* Wavenumber */
Reflection *refl;
double cet, cez; /* Centre of Ewald sphere */
@@ -312,20 +309,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
&& (fabs(rlow) > pr)
&& (fabs(rhigh) > pr) ) return NULL;
- /* Calculate partiality */
- part = partiality(pmodel, h, k, l, image->serial, rlow, rhigh, pr);
-
- if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) {
- ERROR("Assigning NAN partiality!\n");
- ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
- h, k, l, rlow, rhigh);
- ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
- ERROR("D = %e\n", rlow - rhigh);
- ERROR("klow = %e khigh = %e\n", klow, khigh);
- ERROR("distance(..) = %e\n", distance(cet, cez, tl, zl));
- abort();
- }
-
if ( updateme == NULL ) {
refl = reflection_new(h, k, l);
} else {
@@ -373,15 +356,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
}
- set_partial(refl, rlow, rhigh, part);
+ set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by
+ * calculate_partialities() */
set_lorentz(refl, 1.0);
set_symmetric_indices(refl, h, k, l);
set_redundancy(refl, 1);
- if ( output ) {
- printf("%3i %3i %3i %6f %5.2f\n", h, k, l, 0.0, part);
- }
-
return refl;
}
@@ -472,15 +452,7 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
}
-RefList *find_intersections(struct image *image, Crystal *cryst,
- PartialityModel pmodel)
-{
- return find_intersections_to_res(image, cryst, pmodel, INFINITY);
-}
-
-
-RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
- PartialityModel pmodel, double max_res)
+RefList *predict_to_res(Crystal *cryst, double max_res)
{
double ax, ay, az;
double bx, by, bz;
@@ -509,7 +481,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- mres = largest_q(image);
+ mres = largest_q(crystal_get_image(cryst));
if ( mres > max_res ) mres = max_res;
hmax = mres * modulus(ax, ay, az);
@@ -544,7 +516,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- refl = check_reflection(image, cryst, pmodel,
+ refl = check_reflection(crystal_get_image(cryst), cryst,
h, k, l, xl, yl, zl, NULL);
if ( refl != NULL ) {
@@ -559,12 +531,12 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
}
-static void set_unity_partialities(Crystal *cryst)
+static void set_unity_partialities(RefList *list)
{
Reflection *refl;
RefListIterator *iter;
- for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
+ for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -575,14 +547,85 @@ static void set_unity_partialities(Crystal *cryst)
/**
- * update_partialities:
+ * calculate_partialities:
* @cryst: A %Crystal
* @pmodel: A %PartialityModel
*
- * Updates the partialities of @cryst's reflections according to @pmodel and
+ * Calculates the partialities for the reflections in @cryst, given the current
+ * crystal and image parameters. The crystal's image and reflection lists
+ * must be set. The specified %PartialityModel will be used.
+ *
+ * You must not have changed the crystal or image parameters since you last
+ * called predict_to_res() or update_predictions(), because this function
+ * relies on the limiting wavelength values calculated by those functions.
+ */
+void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
+{
+ RefList *list;
+ Reflection *refl;
+ RefListIterator *iter;
+ double pr;
+ struct image *image;
+
+ list = crystal_get_reflections(cryst);
+ if ( list == NULL ) {
+ ERROR("No reflections for partiality calculation!\n");
+ return;
+ }
+
+ if ( pmodel == PMODEL_UNITY ) {
+ set_unity_partialities(list);
+ return;
+ }
+
+ image = crystal_get_image(cryst);
+ if ( image == NULL ) {
+ ERROR("No image for partiality calculation!\n");
+ return;
+ }
+
+ pr = crystal_get_profile_radius(cryst);
+
+ for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double rlow, rhigh, part;
+ signed int h, k, l;
+
+ get_symmetric_indices(refl, &h, &k, &l);
+ get_partial(refl, &rlow, &rhigh, &part);
+
+ part = partiality(pmodel, h, k, l, image->serial,
+ rlow, rhigh, pr);
+
+ if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) {
+ ERROR("Assigning NAN partiality!\n");
+ ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
+ h, k, l, rlow, rhigh);
+ ERROR("R = %e, bw = %e\n", pr, image->bw);
+ ERROR("D = %e\n", rlow - rhigh);
+ abort();
+ }
+
+ set_partiality(refl, part);
+
+ }
+}
+
+
+/**
+ * update_predictions:
+ * @cryst: A %Crystal
+ *
+ * Updates the predicted reflections (positions and excitation errors, but not
+ * the actual partialities) of @cryst's reflections according to
* the current state of the crystal (e.g. its unit cell parameters).
+ *
+ * If you need to update the partialities as well, call calculate_partialities()
+ * afterwards.
*/
-void update_partialities(Crystal *cryst, PartialityModel pmodel)
+void update_predictions(Crystal *cryst)
{
Reflection *refl;
RefListIterator *iter;
@@ -591,11 +634,6 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel)
double csx, csy, csz;
struct image *image = crystal_get_image(cryst);
- if ( pmodel == PMODEL_UNITY ) {
- set_unity_partialities(cryst);
- return;
- }
-
cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz,
&bsx, &bsy, &bsz, &csx, &csy, &csz);
@@ -613,8 +651,8 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel)
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- check_reflection(image, cryst, pmodel,
- h, k, l, xl, yl, zl, refl);
+ check_reflection(image, cryst, h, k, l,
+ xl, yl, zl, refl);
}
}
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 6b035175..4fbd0bca 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
*
* Geometry of diffraction
*
- * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2016 Thomas White <taw@physics.org>
* 2012 Richard Kirian
*
* This file is part of CrystFEL.
@@ -85,15 +85,13 @@ enum gparam {
};
-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 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_partialities(Crystal *cryst, PartialityModel pmodel);
+extern void update_predictions(Crystal *cryst);
extern void polarisation_correction(RefList *list, UnitCell *cell,
struct image *image);
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index d50f8d6b..96eec92a 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1673,8 +1673,7 @@ void integrate_all_4(struct image *image, IntegrationMethod meth,
image->features);
crystal_set_resolution_limit(image->crystals[i], res);
- list = find_intersections_to_res(image, image->crystals[i],
- pmodel, res+push_res);
+ list = predict_to_res(image->crystals[i], res+push_res);
crystal_set_reflections(image->crystals[i], list);
}
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index e702cf33..3246dbec 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -262,7 +262,7 @@ static int pair_peaks(struct image *image, Crystal *cr,
* doesn't fall on this panel. We're only interested
* in how far away it is from the peak location.
* The predicted position and excitation errors will be
- * filled in by update_partialities(). */
+ * filled in by update_predictions(). */
set_panel(refl, f->p);
rps[n].refl = refl;
@@ -275,7 +275,7 @@ static int pair_peaks(struct image *image, Crystal *cr,
/* Get the excitation errors and detector positions for the candidate
* reflections */
crystal_set_reflections(cr, all_reflist);
- update_partialities(cr, PMODEL_SCSPHERE);
+ update_predictions(cr);
/* Pass over the peaks again, keeping only the ones which look like
* good pairings */
@@ -339,7 +339,7 @@ void refine_radius(Crystal *cr, struct image *image)
return;
}
crystal_set_reflections(cr, reflist);
- update_partialities(cr, PMODEL_SCSPHERE);
+ update_predictions(cr);
crystal_set_reflections(cr, NULL);
qsort(rps, n_acc, sizeof(struct reflpeak), cmpd2);
@@ -635,7 +635,7 @@ int refine_prediction(struct image *image, Crystal *cr)
/* Refine */
for ( i=0; i<MAX_CYCLES; i++ ) {
- update_partialities(cr, PMODEL_SCSPHERE);
+ update_predictions(cr);
if ( iterate(rps, n, crystal_get_cell(cr), image,
&total_x, &total_y, &total_z) ) return 1;
//STATUS("Residual after %i = %e\n", i,
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 2f0a0277..4a348742 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -384,8 +384,9 @@ static void run_job(void *vwargs, int cookie)
snprintf(wargs->image.filename, 255, "dummy.h5");
}
- reflections = find_intersections(&wargs->image, cr, PMODEL_SCSPHERE);
+ reflections = predict_to_res(cr, largest_q(&wargs->image));
crystal_set_reflections(cr, reflections);
+ calculate_partialities(cr, PMODEL_SCSPHERE);
for ( i=0; i<NBINS; i++ ) {
wargs->n_ref[i] = 0;
diff --git a/src/partialator.c b/src/partialator.c
index aded835d..428b0d68 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -1135,7 +1135,8 @@ int main(int argc, char *argv[])
crystal_set_image(cryst, &images[i]);
/* Now it's safe to do the following */
- update_partialities(cryst, pmodel);
+ update_predictions(cryst);
+ calculate_partialities(cryst, pmodel);
}
}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 487ce872..bf30d299 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -606,7 +606,8 @@ static void write_residual_graph(Crystal *cr, const RefList *full)
cell_set_reciprocal(cell, asx, asy, asz,
bsx, bsy, bsz,
csx, csy, csz);
- update_partialities(cr, PMODEL_SCSPHERE);
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_SCSPHERE);
res = residual(cr, full, 0, &n, NULL);
fprintf(fh, "%i %e %e %i\n", i, asx, res, n);
}
@@ -614,7 +615,8 @@ static void write_residual_graph(Crystal *cr, const RefList *full)
cell_set_reciprocal(cell, orig_asx, asy, asz,
bsx, bsy, bsz,
csx, csy, csz);
- update_partialities(cr, PMODEL_SCSPHERE);
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_SCSPHERE);
fclose(fh);
}
@@ -648,7 +650,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
pr_iterate(cr, full, pmodel, NULL, verbose);
- update_partialities(cr, pmodel);
+ update_predictions(cr);
+ calculate_partialities(cr, pmodel);
dev = residual(cr, full, 0, 0, NULL);
if ( fabs(dev - old_dev) < dev*0.0001 ) done = 1;
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 9c2a6633..5322fcca 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -178,7 +178,9 @@ static void calc_either_side(Crystal *cr, double incr_val,
/* Crystal properties */
cr_new = new_shifted_crystal(cr, refine, -incr_val);
- compare = find_intersections(image, cr_new, pmodel);
+ compare = predict_to_res(cr_new, largest_q(image));
+ crystal_set_reflections(cr_new, compare);
+ calculate_partialities(cr_new, pmodel);
scan_partialities(crystal_get_reflections(cr), compare, valid,
vals, 0, pmodel);
cell_free(crystal_get_cell(cr_new));
@@ -186,7 +188,9 @@ static void calc_either_side(Crystal *cr, double incr_val,
reflist_free(compare);
cr_new = new_shifted_crystal(cr, refine, +incr_val);
- compare = find_intersections(image, cr_new, pmodel);
+ compare = predict_to_res(cr_new, largest_q(image));
+ crystal_set_reflections(cr_new, compare);
+ calculate_partialities(cr_new, pmodel);
scan_partialities(crystal_get_reflections(cr), compare, valid,
vals, 2, pmodel);
cell_free(crystal_get_cell(cr_new));
@@ -196,18 +200,24 @@ static void calc_either_side(Crystal *cr, double incr_val,
} else {
struct image im_moved;
+ Crystal *cr_new = crystal_copy(cr);
+ crystal_set_image(cr_new, &im_moved);
/* "Image" properties */
im_moved = *image;
shift_parameter(&im_moved, refine, -incr_val);
- compare = find_intersections(&im_moved, cr, pmodel);
+ compare = predict_to_res(cr_new, largest_q(&im_moved));
+ crystal_set_reflections(cr_new, compare);
+ calculate_partialities(cr_new, pmodel);
scan_partialities(crystal_get_reflections(cr), compare,
valid, vals, 0, pmodel);
reflist_free(compare);
im_moved = *image;
shift_parameter(&im_moved, refine, +incr_val);
- compare = find_intersections(&im_moved, cr, pmodel);
+ compare = predict_to_res(cr_new, largest_q(&im_moved));
+ crystal_set_reflections(cr_new, compare);
+ calculate_partialities(cr_new, pmodel);
scan_partialities(crystal_get_reflections(cr), compare,
valid, vals, 2, pmodel);
reflist_free(compare);
@@ -237,7 +247,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
int n_line;
double cc;
- reflections = find_intersections(crystal_get_image(cr), cr, pmodel);
+ reflections = predict_to_res(cr, largest_q(crystal_get_image(cr)));
crystal_set_reflections(cr, reflections);
nref = num_reflections(reflections);
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index b9f23217..0ffc07ca 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -185,14 +185,14 @@ static void calc_either_side(Crystal *cr, double incr_val,
Crystal *cr_new;
cr_new = new_shifted_crystal(cr, refine, -incr_val);
- compare = find_intersections(image, cr_new, PMODEL_SCSPHERE);
+ compare = predict_to_res(cr_new, largest_q(image));
scan(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, PMODEL_SCSPHERE);
+ compare = predict_to_res(cr_new, largest_q(image));
scan(crystal_get_reflections(cr), compare, valid, vals, 2);
cell_free(crystal_get_cell(cr_new));
crystal_free(cr_new);
@@ -221,8 +221,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
int n_line;
double cc;
- reflections = find_intersections(crystal_get_image(cr), cr,
- PMODEL_SCSPHERE);
+ reflections = predict_to_res(cr, largest_q(crystal_get_image(cr)));
crystal_set_reflections(cr, reflections);
nref = num_reflections(reflections);
diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
index 9fe8b356..8b9de871 100644
--- a/tests/prof2d_check.c
+++ b/tests/prof2d_check.c
@@ -3,11 +3,11 @@
*
* Check 2D profile fitting
*
- * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2013-2014 Thomas White <taw@physics.org>
+ * 2013-2014,2016 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -135,7 +135,7 @@ int main(int argc, char *argv[])
image.n_crystals = 1;
image.crystals = &cr;
- list = find_intersections(&image, cr, PMODEL_SCSPHERE);
+ list = predict_to_res(cr, largest_q(&image));
crystal_set_reflections(cr, list);
for ( fs=0; fs<w; fs++ ) {
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-23 18:32:04 +0000