aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/geometry.c405
-rw-r--r--libcrystfel/src/geometry.h12
-rw-r--r--libcrystfel/src/integration.c2
-rw-r--r--libcrystfel/src/predict-refine.c8
-rw-r--r--libcrystfel/src/reflist.c63
-rw-r--r--libcrystfel/src/reflist.h7
6 files changed, 293 insertions, 204 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 576655de..beae9a16 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -238,76 +238,37 @@ static double random_partiality(signed int h, signed int k, signed int l,
}
-static double partiality(PartialityModel pmodel,
- signed int h, signed int k, signed int l,
- int serial,
- double rlow, double rhigh, double pr)
-{
- double D = rlow - rhigh;
-
- /* Convert to partiality */
- switch ( pmodel ) {
-
- default:
- case PMODEL_UNITY:
- return 1.0;
-
- case PMODEL_SCSPHERE:
- return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D);
-
- case PMODEL_SCGAUSSIAN:
- return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D);
-
- case PMODEL_RANDOM:
- return random_partiality(h, k, l, serial);
-
- }
-}
-
-
static Reflection *check_reflection(struct image *image, Crystal *cryst,
signed int h, signed int k, signed int l,
double xl, double yl, double zl,
Reflection *updateme)
{
- double tl;
- double rlow, rhigh; /* "Excitation error" */
- double klow, khigh; /* Wavenumber */
Reflection *refl;
- double cet, cez; /* Centre of Ewald sphere */
- double pr;
- double del;
+ double R, top;
+ double kmin, kmax, k0, knom, k1;
+ double dcs, exerr;
/* Don't predict 000 */
if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
- pr = crystal_get_profile_radius(cryst);
- del = image->div + crystal_get_mosaicity(cryst);
-
- /* "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);
-
- /* If the point is looking "backscattery", reject it straight away */
- if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL;
-
- tl = sqrt(xl*xl + yl*yl);
-
- cet = -sin(del/2.0) * khigh;
- cez = -cos(del/2.0) * khigh;
- rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */
-
- cet = sin(del/2.0) * klow;
- cez = -cos(del/2.0) * klow;
- rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */
-
- /* Condition for reflection to be excited at all */
- if ( (updateme == NULL)
- && (signbit(rlow) == signbit(rhigh))
- && (fabs(rlow) > pr)
- && (fabs(rhigh) > pr) ) return NULL;
+ /* Calculate the limiting wavelengths, lambda0 and lambda1
+ * = 1/k0 and 1/k1 respectively */
+ R = crystal_get_profile_radius(cryst);
+ top = R*R - xl*xl - yl*yl - zl*zl;
+ k0 = top/(2.0*(zl+R));
+ k1 = top/(2.0*(zl-R));
+
+ /* The reflection is excited if any of the reflection is within 2sigma
+ * of the nominal * wavelength of the X-ray beam
+ * (NB image->bw is full width) */
+ kmin = 1.0/(image->lambda + image->lambda*image->bw);
+ knom = 1.0/image->lambda;
+ kmax = 1.0/(image->lambda - image->lambda*image->bw);
+ if ( (k1>kmax) || (k0<kmin) ) return NULL;
+
+ /* Calculate excitation error */
+ dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl);
+ exerr = 1.0/image->lambda - dcs; /* Loss of precision */
if ( updateme == NULL ) {
refl = reflection_new(h, k, l);
@@ -321,7 +282,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
if ( (image->det != NULL) && (updateme != NULL) ) {
double fs, ss;
- locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda,
+ locate_peak_on_panel(xl, yl, zl, knom,
get_panel(updateme), &fs, &ss);
set_detector_pos(refl, fs, ss);
@@ -333,7 +294,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
double fs, ss; /* Position on detector */
signed int p; /* Panel number */
- p = locate_peak(xl, yl, zl, 1.0/image->lambda,
+ p = locate_peak(xl, yl, zl, knom,
image->det, &fs, &ss);
if ( p == -1 ) {
reflection_free(refl);
@@ -344,20 +305,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
- if ( unlikely(rlow < rhigh) ) {
- ERROR("Reflection with rlow < rhigh!\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);
- /* If we are updating, this is (kind of) OK */
- if ( updateme == NULL ) {
- reflection_free(refl);
- return NULL;
- }
- }
-
- set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by
- * calculate_partialities() */
+ set_kpred(refl, knom);
+ set_exerr(refl, exerr);
set_lorentz(refl, 1.0);
set_symmetric_indices(refl, h, k, l);
set_redundancy(refl, 1);
@@ -368,18 +317,12 @@ 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);
+ double tl, phi, azi;
get_symmetric_indices(refl, &hs, &ks, &ls);
@@ -390,26 +333,9 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
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);
+ phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */
+ azi = atan2(yl, xl); /* azimuth */
switch ( k ) {
@@ -531,11 +457,17 @@ RefList *predict_to_res(Crystal *cryst, double max_res)
}
-static void set_unity_partialities(RefList *list)
+static void set_unity_partialities(Crystal *cryst)
{
+ RefList *list;
Reflection *refl;
RefListIterator *iter;
+ list = crystal_get_reflections(cryst);
+ if ( list == NULL ) {
+ ERROR("No reflections for partiality calculation!\n");
+ return;
+ }
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
@@ -546,25 +478,11 @@ static void set_unity_partialities(RefList *list)
}
-/**
- * calculate_partialities:
- * @cryst: A %Crystal
- * @pmodel: A %PartialityModel
- *
- * 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)
+static void set_random_partialities(Crystal *cryst)
{
RefList *list;
Reflection *refl;
RefListIterator *iter;
- double pr;
struct image *image;
list = crystal_get_reflections(cryst);
@@ -573,8 +491,109 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
return;
}
- if ( pmodel == PMODEL_UNITY ) {
- set_unity_partialities(list);
+ image = crystal_get_image(cryst);
+ if ( image == NULL ) {
+ ERROR("No image structure for partiality calculation!\n");
+ return;
+ }
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+ get_symmetric_indices(refl, &h, &k, &l);
+ set_partiality(refl, random_partiality(h, k, l, image->serial));
+ set_lorentz(refl, 1.0);
+ }
+}
+
+
+static double do_integral(double q2, double zl, double R,
+ double lambda, double sig, int verbose)
+{
+ int i;
+ double kmin, kmax, kstart, kfinis;
+ double inc;
+ double total = 0.0;
+ double k0, k1;
+ const int SAMPLES = 50; /* Number of samples for integration */
+ const double N = 1.5; /* Pointiness of spectrum */
+ FILE *fh = NULL;
+
+ k0 = (R*R - q2)/(2.0*(zl+R));
+ k1 = (R*R - q2)/(2.0*(zl-R));
+
+ /* Range over which E is significantly different from zero */
+ kmin = 1.0 / (lambda + 5.0*sig);
+ kmax = 1.0 / (lambda - 5.0*sig);
+
+ kstart = kmin > k1 ? kmin : k1;
+ kfinis = (k0 < 0.0) || (kmax < k0) ? kmax : k0;
+ inc = (kfinis - kstart) / SAMPLES;
+
+ if ( verbose ) {
+ char fn[64];
+ snprintf(fn, 63, "partial%i.graph", verbose);
+ fh = fopen(fn, "w");
+ fprintf(fh, " n p wavelength E P\n");
+ STATUS("Nominal k = %e m^-1\n", 1.0/lambda);
+ STATUS(" (wavelength %e m)\n", lambda);
+ STATUS("Bandwidth %e m\n", sig);
+ STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl));
+ STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl)));
+ STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/k1, 1.0/k0);
+ STATUS("Beam goes from %e to %e m^-1\n", kmin, kmax);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/kmin, 1.0/kmax);
+ STATUS("Integration goes from %e to %e m^-1\n", kstart, kfinis);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/kstart, 1.0/kfinis);
+ }
+
+ for ( i=0; i<SAMPLES; i++ ) {
+
+ double p, kp, lrel;
+ double E, P;
+
+ kp = kstart + i*inc;
+ double pref = sqrt(q2 + kp*kp + 2.0*zl*kp)/(2.0*R);
+ p = pref + 0.5 - kp/(2.0*R);
+
+ /* Spectral energy term */
+ lrel = fabs(1.0/kp - lambda);
+ E = exp(-0.5 * pow(lrel / sig, N));
+ E /= sqrt(2.0 * M_PI * sig);
+
+ /* RLP profile term */
+ P = 4.0*p * (1.0 - p);
+
+ total += E*P*inc;
+
+ if ( fh != NULL ) {
+ fprintf(fh, "%3i %f %e %e %e\n", i, p, 1.0/kp, E, P);
+ }
+ }
+
+ if ( fh != NULL ) fclose(fh);
+
+ return total;
+}
+
+
+
+static void ginn_spectrum_partialities(Crystal *cryst)
+{
+ RefList *list;
+ Reflection *refl;
+ RefListIterator *iter;
+ double r0, m, lambda, sig;
+ struct image *image;
+ UnitCell *cell;
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+
+ list = crystal_get_reflections(cryst);
+ if ( list == NULL ) {
+ ERROR("No reflections for partiality calculation!\n");
return;
}
@@ -584,31 +603,89 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
return;
}
- pr = crystal_get_profile_radius(cryst);
+ cell = crystal_get_cell(cryst);
+ if ( cell == NULL ) {
+ ERROR("No unit cell for partiality calculation!\n");
+ return;
+ }
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ r0 = crystal_get_profile_radius(cryst);
+ m = crystal_get_mosaicity(cryst);
+ lambda = image->lambda;
+ sig = image->bw * lambda;
for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- double rlow, rhigh, part;
+ double R;
signed int h, k, l;
+ double xl, yl, zl;
+ double q2;
+ double total, norm;
get_symmetric_indices(refl, &h, &k, &l);
- get_partial(refl, &rlow, &rhigh, &part);
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+ zl = h*asz + k*bsz + l*csz;
- part = partiality(pmodel, h, k, l, image->serial,
- rlow, rhigh, pr);
+ /* Radius of rlp profile */
+ q2 = xl*xl + yl*yl + zl*zl;
- 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);
+ R = r0 + m * sqrt(q2);
+
+ total = do_integral(q2, zl, R, lambda, sig, 0);
+ norm = do_integral(q2, -0.5*q2*lambda, R, lambda, sig, 0);
+
+ set_partiality(refl, total/norm);
+ set_lorentz(refl, 1.0);
+
+ if ( total > 2.0*norm ) {
+ /* Error! */
+ do_integral(q2, zl, R, lambda, sig, 1);
+ do_integral(q2, -0.5*q2*lambda, R, lambda, sig, 2);
abort();
- }
+ }
- set_partiality(refl, part);
+ }
+}
+
+
+/**
+ * calculate_partialities:
+ * @cryst: A %Crystal
+ * @pmodel: A %PartialityModel
+ *
+ * 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)
+{
+ switch ( pmodel ) {
+
+ case PMODEL_UNITY :
+ set_unity_partialities(cryst);
+ break;
+
+ case PMODEL_XSPHERE :
+ ginn_spectrum_partialities(cryst);
+ break;
+
+ case PMODEL_RANDOM :
+ set_random_partialities(cryst);
+ break;
+
+ default :
+ ERROR("Unknown partiality model %i\n", pmodel);
+ break;
}
}
@@ -694,29 +771,30 @@ void polarisation_correction(RefList *list, UnitCell *cell, struct image *image)
/* Returns dx_h/dP, where P = any parameter */
-double x_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda)
+double x_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p)
{
signed int h, k, l;
- double x, z, wn;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
+ double xl, zl, kpred;
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
get_indices(refl, &h, &k, &l);
- wn = 1.0 / lambda;
- cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- x = h*ax + k*bx + l*cx;
- z = h*az + k*bz + l*cz;
+ 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->clen / (wn+z);
+ return h * p->clen / (kpred + zl);
case GPARAM_BSX :
- return k * p->clen / (wn+z);
+ return k * p->clen / (kpred + zl);
case GPARAM_CSX :
- return l * p->clen / (wn+z);
+ return l * p->clen / (kpred + zl);
case GPARAM_ASY :
return 0.0;
@@ -728,13 +806,13 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
return 0.0;
case GPARAM_ASZ :
- return -h * x * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -h * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_BSZ :
- return -k * x * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -k * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_CSZ :
- return -l * x * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -l * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_DETX :
return -1;
@@ -743,7 +821,7 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
return 0;
case GPARAM_CLEN :
- return x / (wn+z);
+ return xl / (kpred+zl);
}
@@ -753,18 +831,19 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
/* Returns dy_h/dP, where P = any parameter */
-double y_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda)
+double y_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p)
{
signed int h, k, l;
- double y, z, wn;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
+ double yl, zl, kpred;
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
get_indices(refl, &h, &k, &l);
- wn = 1.0 / lambda;
- cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- y = h*ay + k*by + l*cy;
- z = h*az + k*bz + l*cz;
+ 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 ) {
@@ -778,22 +857,22 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
return 0.0;
case GPARAM_ASY :
- return h * p->clen / (wn+z);
+ return h * p->clen / (kpred + zl);
case GPARAM_BSY :
- return k * p->clen / (wn+z);
+ return k * p->clen / (kpred + zl);
case GPARAM_CSY :
- return l * p->clen / (wn+z);
+ return l * p->clen / (kpred + zl);
case GPARAM_ASZ :
- return -h * y * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -h * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_BSZ :
- return -k * y * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -k * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_CSZ :
- return -l * y * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -l * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_DETX :
return 0;
@@ -802,7 +881,7 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
return -1;
case GPARAM_CLEN :
- return y / (wn+z);
+ return yl / (kpred+zl);
}
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 4fbd0bca..c2f2899d 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,7 +3,7 @@
*
* Geometry of diffraction
*
- * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
@@ -47,8 +47,7 @@ extern "C" {
/**
* PartialityModel:
* @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
- * @PMODEL_SCSPHERE : Sphere model with source coverage factor included
- * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included
+ * @PMODEL_XSPHERE : Flat sphere model with super-Gaussian spectrum
* @PMODEL_RANDOM : Randomly assigned partialities
*
* A %PartialityModel describes a geometrical model which can be used to
@@ -57,8 +56,7 @@ extern "C" {
typedef enum {
PMODEL_UNITY,
- PMODEL_SCSPHERE,
- PMODEL_SCGAUSSIAN,
+ PMODEL_XSPHERE,
PMODEL_RANDOM,
} PartialityModel;
@@ -99,9 +97,9 @@ extern double sphere_fraction(double rlow, double rhigh, double pr);
extern double gaussian_fraction(double rlow, double rhigh, double pr);
extern double x_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda);
+ struct panel *p);
extern double y_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda);
+ struct panel *p);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index c18ae110..fcbdf658 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1760,7 +1760,7 @@ void integrate_all_2(struct image *image, IntegrationMethod meth,
IntDiag int_diag,
signed int idh, signed int idk, signed int idl)
{
- integrate_all_3(image, meth, PMODEL_SCSPHERE, 0.0,
+ integrate_all_3(image, meth, PMODEL_XSPHERE, 0.0,
ir_inn, ir_mid, ir_out,
int_diag, idh, idk, idl);
}
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index 3246dbec..14a60bc7 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -91,9 +91,7 @@ static void twod_mapping(double fs, double ss, double *px, double *py,
static double r_dev(struct reflpeak *rp)
{
/* Excitation error term */
- double rlow, rhigh, p;
- get_partial(rp->refl, &rlow, &rhigh, &p);
- return (rlow+rhigh)/2.0;
+ return get_exerr(rp->refl);
}
@@ -425,7 +423,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
/* Positional x terms */
for ( k=0; k<num_params; k++ ) {
gradients[k] = x_gradient(rv[k], rps[i].refl, cell,
- rps[i].panel, image->lambda);
+ rps[i].panel);
}
for ( k=0; k<num_params; k++ ) {
@@ -457,7 +455,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
/* Positional y terms */
for ( k=0; k<num_params; k++ ) {
gradients[k] = y_gradient(rv[k], rps[i].refl, cell,
- rps[i].panel, image->lambda);
+ rps[i].panel);
}
for ( k=0; k<num_params; k++ ) {
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 707c802e..69fb29be 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -67,10 +67,10 @@ struct _refldata {
signed int ls;
/* Partiality and related geometrical stuff */
- double rlow; /* Low excitation error */
- double rhigh; /* High excitation error */
- double p; /* Partiality */
- double L; /* Lorentz factor */
+ double kpred; /* Wavenumber of middle of reflection */
+ double exerr; /* Excitation error */
+ double p; /* Partiality */
+ double L; /* Lorentz factor */
/* Location in image */
double fs;
@@ -422,22 +422,28 @@ double get_intensity(const Reflection *refl)
/**
- * get_partial:
+ * get_kpred:
* @refl: A %Reflection
- * @rlow: Location at which to store the "low" excitation error
- * @rhigh: Location at which to store the "high" excitation error
- * @p: Location at which to store the partiality
*
- * This function is used during post refinement (in conjunction with
- * set_partial()) to get access to the details of the partiality calculation.
+ * Returns: the wavenumber at the centre of the reflection
*
**/
-void get_partial(const Reflection *refl, double *rlow, double *rhigh,
- double *p)
+double get_kpred(const Reflection *refl)
{
- *rlow = refl->data.rlow;
- *rhigh = refl->data.rhigh;
- *p = get_partiality(refl);
+ return refl->data.kpred;
+}
+
+
+/**
+ * get_exerr:
+ * @refl: A %Reflection
+ *
+ * Returns: the excitation error (in m^-1) for this reflection
+ *
+ **/
+double get_exerr(const Reflection *refl)
+{
+ return refl->data.exerr;
}
@@ -614,21 +620,28 @@ void set_panel(Reflection *refl, struct panel *p)
/**
- * set_partial:
+ * set_kpred:
* @refl: A %Reflection
- * @rlow: The "low" excitation error
- * @rhigh: The "high" excitation error
- * @p: The partiality
+ * @kpred: The wavenumber at which the reflection should be predicted
*
- * This function is used during post refinement (in conjunction with
- * get_partial()) to get access to the details of the partiality calculation.
+ * Sets the wavenumber at the centre of the reflection. Used by
+ * predict_to_res() and update_predictions()
+ **/
+void set_kpred(Reflection *refl, double kpred)
+{
+ refl->data.kpred = kpred;
+}
+
+
+/**
+ * set_exerr:
+ * @refl: A %Reflection
+ * @exerr: The excitation error for the reflection
*
**/
-void set_partial(Reflection *refl, double rlow, double rhigh, double p)
+void set_exerr(Reflection *refl, double exerr)
{
- refl->data.rlow = rlow;
- refl->data.rhigh = rhigh;
- refl->data.p = p;
+ refl->data.exerr = exerr;
}
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index f2126ac9..4dd1e56f 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -87,6 +87,8 @@ extern Reflection *next_found_refl(Reflection *refl);
extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
extern struct panel *get_panel(const Reflection *refl);
extern double get_partiality(const Reflection *refl);
+extern double get_kpred(const Reflection *refl);
+extern double get_exerr(const Reflection *refl);
extern double get_lorentz(const Reflection *refl);
extern void get_indices(const Reflection *refl,
signed int *h, signed int *k, signed int *l);
@@ -94,8 +96,6 @@ extern void get_symmetric_indices(const Reflection *refl,
signed int *hs, signed int *ks,
signed int *ls);
extern double get_intensity(const Reflection *refl);
-extern void get_partial(const Reflection *refl, double *rlow, double *rhigh,
- double *p);
extern int get_redundancy(const Reflection *refl);
extern double get_temp1(const Reflection *refl);
extern double get_temp2(const Reflection *refl);
@@ -109,7 +109,8 @@ extern int get_flag(const Reflection *refl);
extern void copy_data(Reflection *to, const Reflection *from);
extern void set_detector_pos(Reflection *refl, double fs, double ss);
extern void set_panel(Reflection *refl, struct panel *p);
-extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
+extern void set_kpred(Reflection *refl, double kpred);
+extern void set_exerr(Reflection *refl, double exerr);
extern void set_partiality(Reflection *refl, double p);
extern void set_lorentz(Reflection *refl, double L);
extern void set_intensity(Reflection *refl, double intensity);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-19 00:58:07 +0000