diff options
-rw-r--r-- | src/post-refinement.c | 43 |
1 files changed, 29 insertions, 14 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c index a532cd07..e98f7253 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -54,9 +54,6 @@ /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ #define MAX_CYCLES (10) -/* Number of parameters to refine, in the order they appear in the enum */ -#define NUM_PARAMS (9) - /* Returns dp(gauss)/dr at "r" */ static double gaussian_fraction_gradient(double r, double R) { @@ -311,13 +308,31 @@ static double pr_iterate(Crystal *cr, const RefList *full, double max_shift; int nref = 0; const int verbose = 0; + int num_params = 0; + enum gparam rv[32]; *n_filtered = 0; + /* If partiality model is anything other than "unity", refine all the + * geometrical parameters */ + if ( pmodel != PMODEL_UNITY ) { + rv[num_params++] = GPARAM_ASX; + rv[num_params++] = GPARAM_ASY; + rv[num_params++] = GPARAM_ASZ; + rv[num_params++] = GPARAM_BSX; + rv[num_params++] = GPARAM_BSY; + rv[num_params++] = GPARAM_BSZ; + rv[num_params++] = GPARAM_CSX; + rv[num_params++] = GPARAM_CSY; + rv[num_params++] = GPARAM_CSZ; + } + + STATUS("Refining %i parameters.\n", num_params); + reflections = crystal_get_reflections(cr); - M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); - v = gsl_vector_calloc(NUM_PARAMS); + M = gsl_matrix_calloc(num_params, num_params); + v = gsl_vector_calloc(num_params); /* Construct the equations, one per reflection in this image */ for ( refl = first_refl(reflections, &iter); @@ -331,7 +346,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, int k; double p, l; Reflection *match; - double gradients[NUM_PARAMS]; + double gradients[num_params]; /* Find the full version */ get_indices(refl, &ha, &ka, &la); @@ -355,16 +370,16 @@ static double pr_iterate(Crystal *cr, const RefList *full, w = pow(w, -1.0); /* Calculate all gradients for this reflection */ - for ( k=0; k<NUM_PARAMS; k++ ) { - gradients[k] = p_gradient(cr, k, refl, pmodel) * l; + for ( k=0; k<num_params; k++ ) { + gradients[k] = p_gradient(cr, rv[k], refl, pmodel) * l; } - for ( k=0; k<NUM_PARAMS; k++ ) { + for ( k=0; k<num_params; k++ ) { int g; double v_c, v_curr; - for ( g=0; g<NUM_PARAMS; g++ ) { + for ( g=0; g<num_params; g++ ) { double M_c, M_curr; @@ -406,9 +421,9 @@ static double pr_iterate(Crystal *cr, const RefList *full, shifts = solve_svd(v, M, n_filtered, verbose); if ( shifts != NULL ) { - for ( param=0; param<NUM_PARAMS; param++ ) { + for ( param=0; param<num_params; param++ ) { double shift = gsl_vector_get(shifts, param); - apply_shift(cr, param, shift); + apply_shift(cr, rv[param], shift); //STATUS("Shift %i: %e\n", param, shift); if ( fabs(shift) > max_shift ) max_shift = fabs(shift); } |