aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/predict-refine.c62
-rw-r--r--tests/prediction_gradient_check.c38
2 files changed, 45 insertions, 55 deletions
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index 6ef31805..6784aae2 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -349,23 +349,16 @@ void refine_radius(Crystal *cr, struct image *image)
}
-/* Returns d(xh-xpk)/dP, where P = any parameter */
-static double x_gradient(int param, struct reflpeak *rp, struct detector *det,
- double lambda, UnitCell *cell)
+/* Returns dx_h/dP, where P = any parameter */
+static double x_gradient(int param, Reflection *refl, UnitCell *cell,
+ struct panel *p, double lambda)
{
signed int h, k, l;
- double xpk, ypk, xh, yh;
- double fsh, ssh;
double x, z, wn;
double ax, ay, az, bx, by, bz, cx, cy, cz;
- twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
- get_detector_pos(rp->refl, &fsh, &ssh);
- twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
- get_indices(rp->refl, &h, &k, &l);
-
+ 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;
@@ -373,13 +366,13 @@ static double x_gradient(int param, struct reflpeak *rp, struct detector *det,
switch ( param ) {
case GPARAM_ASX :
- return h * rp->panel->clen / (wn+z);
+ return h * p->clen / (wn+z);
case GPARAM_BSX :
- return k * rp->panel->clen / (wn+z);
+ return k * p->clen / (wn+z);
case GPARAM_CSX :
- return l * rp->panel->clen / (wn+z);
+ return l * p->clen / (wn+z);
case GPARAM_ASY :
return 0.0;
@@ -391,13 +384,13 @@ static double x_gradient(int param, struct reflpeak *rp, struct detector *det,
return 0.0;
case GPARAM_ASZ :
- return -h * x * rp->panel->clen / (wn*wn + 2*wn*z + z*z);
+ return -h * x * p->clen / (wn*wn + 2*wn*z + z*z);
case GPARAM_BSZ :
- return -k * x * rp->panel->clen / (wn*wn + 2*wn*z + z*z);
+ return -k * x * p->clen / (wn*wn + 2*wn*z + z*z);
case GPARAM_CSZ :
- return -l * x * rp->panel->clen / (wn*wn + 2*wn*z + z*z);
+ return -l * x * p->clen / (wn*wn + 2*wn*z + z*z);
case GPARAM_DETX :
return -1;
@@ -415,23 +408,16 @@ static double x_gradient(int param, struct reflpeak *rp, struct detector *det,
}
-/* Returns d(yh-ypk)/dP, where P = any parameter */
-static double y_gradient(int param, struct reflpeak *rp, struct detector *det,
- double lambda, UnitCell *cell)
+/* Returns dy_h/dP, where P = any parameter */
+static double y_gradient(int param, Reflection *refl, UnitCell *cell,
+ struct panel *p, double lambda)
{
signed int h, k, l;
- double xpk, ypk, xh, yh;
- double fsh, ssh;
double y, z, wn;
double ax, ay, az, bx, by, bz, cx, cy, cz;
- twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
- get_detector_pos(rp->refl, &fsh, &ssh);
- twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
- get_indices(rp->refl, &h, &k, &l);
-
+ 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;
@@ -448,22 +434,22 @@ static double y_gradient(int param, struct reflpeak *rp, struct detector *det,
return 0.0;
case GPARAM_ASY :
- return h * rp->panel->clen / (wn+z);
+ return h * p->clen / (wn+z);
case GPARAM_BSY :
- return k * rp->panel->clen / (wn+z);
+ return k * p->clen / (wn+z);
case GPARAM_CSY :
- return l * rp->panel->clen / (wn+z);
+ return l * p->clen / (wn+z);
case GPARAM_ASZ :
- return -h * y * rp->panel->clen / (wn*wn + 2*wn*z + z*z);
+ return -h * y * p->clen / (wn*wn + 2*wn*z + z*z);
case GPARAM_BSZ :
- return -k * y * rp->panel->clen / (wn*wn + 2*wn*z + z*z);
+ return -k * y * p->clen / (wn*wn + 2*wn*z + z*z);
case GPARAM_CSZ :
- return -l * y * rp->panel->clen / (wn*wn + 2*wn*z + z*z);
+ return -l * y * p->clen / (wn*wn + 2*wn*z + z*z);
case GPARAM_DETX :
return 0;
@@ -553,8 +539,8 @@ 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], image->det,
- image->lambda, cell);
+ gradients[k] = x_gradient(rv[k], rps[i].refl, cell,
+ rps[i].panel, image->lambda);
}
for ( k=0; k<num_params; k++ ) {
@@ -585,8 +571,8 @@ 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], image->det,
- image->lambda, cell);
+ gradients[k] = y_gradient(rv[k], rps[i].refl, cell,
+ rps[i].panel, image->lambda);
}
for ( k=0; k<num_params; k++ ) {
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index 18c9f023..c5d4206b 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -41,12 +41,24 @@
#include <cell-utils.h>
#include <geometry.h>
#include <reflist.h>
-#include "../src/predict-refine.c"
int checkrxy;
+static void twod_mapping(double fs, double ss, double *px, double *py,
+ struct panel *p)
+{
+ double xs, ys;
+
+ xs = fs*p->fsx + ss*p->ssx;
+ ys = fs*p->fsy + ss*p->ssy;
+
+ *px = (xs + p->cnx) / p->res;
+ *py = (ys + p->cny) / p->res;
+}
+
+
static void scan(RefList *reflections, RefList *compare,
int *valid, long double *vals[3], int idx)
{
@@ -260,7 +272,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
long double grad1, grad2, grad;
double cgrad;
signed int h, k, l;
- struct reflpeak rp;
get_indices(refl, &h, &k, &l);
@@ -283,27 +294,20 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
} else {
- struct imagefeature pk;
struct image *image;
- pk.fs = 0.0;
- pk.ss = 0.0;
-
image = crystal_get_image(cr);
- rp.refl = refl;
- rp.peak = &pk;
- rp.panel = &image->det->panels[0];
if ( checkrxy == 1 ) {
- cgrad = x_gradient(refine, &rp,
- crystal_get_image(cr)->det,
- crystal_get_image(cr)->lambda,
- crystal_get_cell(cr));
+ cgrad = x_gradient(refine, refl,
+ crystal_get_cell(cr),
+ &image->det->panels[0],
+ crystal_get_image(cr)->lambda);
} else {
- cgrad = y_gradient(refine, &rp,
- crystal_get_image(cr)->det,
- crystal_get_image(cr)->lambda,
- crystal_get_cell(cr));
+ cgrad = y_gradient(refine, refl,
+ crystal_get_cell(cr),
+ &image->det->panels[0],
+ crystal_get_image(cr)->lambda);
}
}