aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/crystal.c18
-rw-r--r--libcrystfel/src/crystal.h2
-rw-r--r--libcrystfel/src/dirax.c2
-rw-r--r--src/partial_sim.c5
-rw-r--r--src/post-refinement.c6
-rw-r--r--src/post-refinement.h2
-rw-r--r--tests/pr_gradient_check.c144
7 files changed, 126 insertions, 53 deletions
diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c
index 6ff251e0..664d933e 100644
--- a/libcrystfel/src/crystal.c
+++ b/libcrystfel/src/crystal.c
@@ -106,8 +106,6 @@ Crystal *crystal_new()
void crystal_free(Crystal *cryst)
{
if ( cryst == NULL ) return;
- if ( cryst->cell != NULL ) cell_free(cryst->cell);
- if ( cryst->reflections != NULL ) reflist_free(cryst->reflections);
free(cryst);
}
@@ -163,13 +161,18 @@ int crystal_get_user_flag(Crystal *cryst)
}
+double crystal_get_mosaicity(Crystal *cryst)
+{
+ return cryst->m;
+}
+
+
/********************************** Setters ***********************************/
void crystal_set_cell(Crystal *cryst, UnitCell *cell)
{
- if ( cryst->cell != NULL ) cell_free(cryst->cell);
- cryst->cell = cell_new_from_cell(cell);
+ cryst->cell = cell;
}
@@ -181,7 +184,6 @@ void crystal_set_profile_radius(Crystal *cryst, double r)
void crystal_set_reflections(Crystal *cryst, RefList *reflist)
{
- if ( cryst->reflections != NULL ) reflist_free(reflist);
cryst->reflections = reflist;
}
@@ -214,3 +216,9 @@ void crystal_set_user_flag(Crystal *cryst, int user_flag)
{
cryst->user_flag = user_flag;
}
+
+
+void crystal_set_mosaicity(Crystal *cryst, double m)
+{
+ cryst->m = m;
+}
diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h
index f7c49a81..3d0ad9d1 100644
--- a/libcrystfel/src/crystal.h
+++ b/libcrystfel/src/crystal.h
@@ -58,6 +58,7 @@ extern long long int crystal_get_num_saturated_reflections(Crystal *cryst);
extern int crystal_get_user_flag(Crystal *cryst);
extern double crystal_get_osf(Crystal *cryst);
extern struct image *crystal_get_image(Crystal *cryst);
+extern double crystal_get_mosaicity(Crystal *cryst);
extern void crystal_set_cell(Crystal *cryst, UnitCell *cell);
extern void crystal_set_profile_radius(Crystal *cryst, double r);
@@ -68,5 +69,6 @@ extern void crystal_set_num_saturated_reflections(Crystal *cryst,
extern void crystal_set_user_flag(Crystal *cryst, int flag);
extern void crystal_set_osf(Crystal *cryst, double osf);
extern void crystal_set_image(Crystal *cryst, struct image *image);
+extern void crystal_set_mosaicity(Crystal *cryst, double m);
#endif /* CRYSTAL_H */
diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c
index 6fbb9663..2cbb0ec1 100644
--- a/libcrystfel/src/dirax.c
+++ b/libcrystfel/src/dirax.c
@@ -140,8 +140,6 @@ static int check_cell(struct dirax_private *dp, struct image *image,
image_add_crystal(image, cr);
- cell_free(out); /* Crystal makes its own copy */
-
return 1;
}
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 4e96fbbc..502e7028 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -252,7 +252,6 @@ static void run_job(void *vwargs, int cookie)
struct queue_args *qargs = wargs->qargs;
int i;
Crystal *cr;
- UnitCell *cell;
RefList *reflections;
cr = crystal_new();
@@ -268,9 +267,7 @@ static void run_job(void *vwargs, int cookie)
/* Set up a random orientation */
orientation = random_quaternion();
- cell = cell_rotate(qargs->cell, orientation);
- crystal_set_cell(cr, cell);
- cell_free(cell);
+ crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
snprintf(wargs->image.filename, 255, "dummy.h5");
reflections = find_intersections(&wargs->image, cr);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 31558ec7..1439b148 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius)
/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(Crystal *cr, int k, Reflection *refl, double r)
+double gradient(Crystal *cr, int k, Reflection *refl)
{
double ds, azix, aziy;
double ttlow, tthigh, tt;
@@ -104,6 +104,7 @@ double gradient(Crystal *cr, int k, Reflection *refl, double r)
double klow, khigh;
double gr;
struct image *image = crystal_get_image(cr);
+ double r = crystal_get_profile_radius(cr);
get_symmetric_indices(refl, &hs, &ks, &ls);
@@ -419,8 +420,7 @@ static double pr_iterate(Crystal *cr, const RefList *full)
/* Calculate all gradients for this reflection */
for ( k=0; k<NUM_PARAMS; k++ ) {
double gr;
- gr = gradient(cr, k, refl,
- crystal_get_profile_radius(cr));
+ gr = gradient(cr, k, refl);
gradients[k] = gr;
}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 52e70b2b..fe171882 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -62,7 +62,7 @@ enum {
extern void pr_refine(Crystal *cr, const RefList *full);
/* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(Crystal *cr, int k, Reflection *refl, double r);
+extern double gradient(Crystal *cr, int k, Reflection *refl);
#endif /* POST_REFINEMENT_H */
diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c
index b39dfd05..9d713e0b 100644
--- a/tests/pr_gradient_check.c
+++ b/tests/pr_gradient_check.c
@@ -122,41 +122,98 @@ static void shift_parameter(struct image *image, int k, double shift)
}
-static void calc_either_side(struct image *image, double incr_val,
+static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
+{
+ Crystal *cr_new;
+ double r;
+ UnitCell *cell;
+
+ cr_new = crystal_new();
+ if ( cr_new == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return NULL;
+ }
+
+ crystal_set_image(cr_new, crystal_get_image(cr));
+ r = crystal_get_profile_radius(cr_new);
+
+ switch ( refine ) {
+
+ case REF_ASX :
+ case REF_ASY :
+ case REF_ASZ :
+ case REF_BSX :
+ case REF_BSY :
+ case REF_BSZ :
+ case REF_CSX :
+ case REF_CSY :
+ case REF_CSZ :
+ cell = new_shifted_cell(crystal_get_cell(cr), refine,
+ -incr_val);
+ crystal_set_cell(cr_new, cell);
+ break;
+
+ case REF_R :
+ cell = cell_new_from_cell(crystal_get_cell(cr));
+ crystal_set_profile_radius(cr_new, r + incr_val);
+ break;
+
+ default :
+ ERROR("Can't shift %i\n", refine);
+ break;
+
+ }
+
+ return cr_new;
+}
+
+static void calc_either_side(Crystal *cr, double incr_val,
int *valid, long double *vals[3], int refine)
{
RefList *compare;
- UnitCell *cell;
+ struct image *image = crystal_get_image(cr);
+
+ if ( (refine != REF_DIV) ) {
+
+ /* Crystal properties */
- if ( (refine != REF_DIV) && (refine != REF_R) ) {
+ Crystal *cr_new;
- cell = new_shifted_cell(image->indexed_cell, refine, -incr_val);
- compare = find_intersections(image, cell);
- scan_partialities(image->reflections, compare, valid, vals, 0);
- cell_free(cell);
+ cr_new = new_shifted_crystal(cr, refine, -incr_val);
+
+ compare = find_intersections(image, cr_new);
+ scan_partialities(crystal_get_reflections(cr), compare, valid,
+ vals, 0);
+ cell_free(crystal_get_cell(cr_new));
+ crystal_free(cr_new);
reflist_free(compare);
- cell = new_shifted_cell(image->indexed_cell, refine, +incr_val);
- compare = find_intersections(image, cell);
- scan_partialities(image->reflections, compare, valid, vals, 2);
- cell_free(cell);
+ cr_new = new_shifted_crystal(cr, refine, +incr_val);
+
+ compare = find_intersections(image, cr_new);
+ scan_partialities(crystal_get_reflections(cr), compare, valid,
+ vals, 0);
+ cell_free(crystal_get_cell(cr_new));
+ crystal_free(cr_new);
reflist_free(compare);
} else {
+ /* "Image" properties */
+
struct image im_moved;
im_moved = *image;
shift_parameter(&im_moved, refine, -incr_val);
- compare = find_intersections(&im_moved, im_moved.indexed_cell);
- scan_partialities(im_moved.reflections, compare,
+ compare = find_intersections(&im_moved, cr);
+ scan_partialities(crystal_get_reflections(cr), compare,
valid, vals, 0);
reflist_free(compare);
im_moved = *image;
shift_parameter(&im_moved, refine, +incr_val);
- compare = find_intersections(&im_moved, im_moved.indexed_cell);
- scan_partialities(im_moved.reflections, compare,
+ compare = find_intersections(&im_moved, cr);
+ scan_partialities(crystal_get_reflections(cr), compare,
valid, vals, 2);
reflist_free(compare);
@@ -165,7 +222,7 @@ static void calc_either_side(struct image *image, double incr_val,
-static int test_gradients(struct image *image, double incr_val, int refine,
+static int test_gradients(Crystal *cr, double incr_val, int refine,
const char *str)
{
Reflection *refl;
@@ -175,11 +232,13 @@ static int test_gradients(struct image *image, double incr_val, int refine,
int *valid;
int nref;
int n_acc, n_valid;
+ RefList *reflections;
//FILE *fh;
- image->reflections = find_intersections(image, image->indexed_cell);
+ reflections = find_intersections(crystal_get_image(cr), cr);
+ crystal_set_reflections(cr, reflections);
- nref = num_reflections(image->reflections);
+ nref = num_reflections(reflections);
if ( nref < 10 ) {
ERROR("Too few reflections found. Failing test by default.\n");
return -1;
@@ -200,16 +259,15 @@ static int test_gradients(struct image *image, double incr_val, int refine,
}
for ( i=0; i<nref; i++ ) valid[i] = 1;
- scan_partialities(image->reflections, image->reflections,
- valid, vals, 1);
+ scan_partialities(reflections, reflections, valid, vals, 1);
- calc_either_side(image, incr_val, valid, vals, refine);
+ calc_either_side(cr, incr_val, valid, vals, refine);
//fh = fopen("wrongness.dat", "a");
n_valid = nref; n_acc = 0;
i = 0;
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(reflections, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -231,8 +289,7 @@ static int test_gradients(struct image *image, double incr_val, int refine,
grad2 = (vals[2][i] - vals[1][i]) / incr_val;
grad = (grad1 + grad2) / 2.0;
- cgrad = gradient(image, refine, refl,
- image->profile_radius);
+ cgrad = gradient(cr, refine, refl);
get_partial(refl, &r1, &r2, &p, &cl, &ch);
@@ -289,6 +346,7 @@ int main(int argc, char *argv[])
double bx, by, bz;
double cx, cy, cz;
UnitCell *cell;
+ Crystal *cr;
struct quaternion orientation;
int i;
int val;
@@ -303,10 +361,17 @@ int main(int argc, char *argv[])
image.lambda = ph_en_to_lambda(eV_to_J(8000.0));
image.div = 1e-3;
image.bw = 0.01;
- image.m = 0.0;
- image.profile_radius = 0.005e9;
image.filename = malloc(256);
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 1;
+ }
+ crystal_set_mosaicity(cr, 0.0);
+ crystal_set_profile_radius(cr, 0.005e9);
+ crystal_set_image(cr, &image);
+
cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9,
deg2rad(90.0),
deg2rad(90.0),
@@ -316,36 +381,39 @@ int main(int argc, char *argv[])
for ( i=0; i<1; i++ ) {
+ UnitCell *rot;
+
orientation = random_quaternion();
- image.indexed_cell = cell_rotate(cell, orientation);
+ rot = cell_rotate(cell, orientation);
+ crystal_set_cell(cr, rot);
- cell_get_reciprocal(image.indexed_cell,
+ cell_get_reciprocal(rot,
&ax, &ay, &az, &bx, &by,
&bz, &cx, &cy, &cz);
incr_val = incr_frac * image.div;
- val += test_gradients(&image, incr_val, REF_DIV, "div");
+ val += test_gradients(cr, incr_val, REF_DIV, "div");
incr_val = incr_frac * ax;
- val += test_gradients(&image, incr_val, REF_ASX, "ax*");
+ val += test_gradients(cr, incr_val, REF_ASX, "ax*");
incr_val = incr_frac * ay;
- val += test_gradients(&image, incr_val, REF_ASY, "ay*");
+ val += test_gradients(cr, incr_val, REF_ASY, "ay*");
incr_val = incr_frac * az;
- val += test_gradients(&image, incr_val, REF_ASZ, "az*");
+ val += test_gradients(cr, incr_val, REF_ASZ, "az*");
incr_val = incr_frac * bx;
- val += test_gradients(&image, incr_val, REF_BSX, "bx*");
+ val += test_gradients(cr, incr_val, REF_BSX, "bx*");
incr_val = incr_frac * by;
- val += test_gradients(&image, incr_val, REF_BSY, "by*");
+ val += test_gradients(cr, incr_val, REF_BSY, "by*");
incr_val = incr_frac * bz;
- val += test_gradients(&image, incr_val, REF_BSZ, "bz*");
+ val += test_gradients(cr, incr_val, REF_BSZ, "bz*");
incr_val = incr_frac * cx;
- val += test_gradients(&image, incr_val, REF_CSX, "cx*");
+ val += test_gradients(cr, incr_val, REF_CSX, "cx*");
incr_val = incr_frac * cy;
- val += test_gradients(&image, incr_val, REF_CSY, "cy*");
+ val += test_gradients(cr, incr_val, REF_CSY, "cy*");
incr_val = incr_frac * cz;
- val += test_gradients(&image, incr_val, REF_CSZ, "cz*");
+ val += test_gradients(cr, incr_val, REF_CSZ, "cz*");
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 16:40:40 +0000