Add --force-bandwidth and --force-radius

2 files changed, 39 insertions, 1 deletions

diff --git a/src/partialator.c b/src/partialator.c
index e9053841..0d5a8a84 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -329,7 +329,9 @@ static void show_help(const char *s)
 "      --max-rel-B            Maximum allowable relative |B| factor.\n"
 "      --no-logs              Do not write extensive log files.\n"
 "  -w <pg>                    Apparent point group for resolving ambiguities.\n"
-"      --operator=<op>        Indexing ambiguity operator for resolving.\n");
+"      --operator=<op>        Indexing ambiguity operator for resolving.\n"
+"      --force-bandwidth=<n>  Set all bandwidths to <n> (fraction).\n"
+"      --force-radius=<n>     Set all profile radii to <n> nm^-1.\n");
 }
 
 
@@ -858,6 +860,8 @@ int main(int argc, char *argv[])
 	int no_logs = 0;
 	char *w_sym_str = NULL;
 	char *operator = NULL;
+	double force_bandwidth = -1.0;
+	double force_radius = -1.0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -882,6 +886,8 @@ int main(int argc, char *argv[])
 		{"max-rel-b",          1, NULL,                7}, /* compat */
 		{"reference",          1, NULL,                8}, /* ssshhh! */
 		{"operator",           1, NULL,                9},
+		{"force-bandwidth",    1, NULL,               10},
+		{"force-radius",       1, NULL,               11},
 
 		{"no-scale",           0, &no_scale,           1},
 		{"no-pr",              0, &no_pr,              1},
@@ -1020,6 +1026,25 @@ int main(int argc, char *argv[])
 			operator = strdup(optarg);
 			break;
 
+			case 10 :
+			errno = 0;
+			force_bandwidth = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --force-bandwidth.\n");
+				return 1;
+			}
+			break;
+
+			case 11 :
+			errno = 0;
+			force_radius = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --force-radius.\n");
+				return 1;
+			}
+			force_radius *= 1e9;
+			break;
+
 			case 0 :
 			break;
 
@@ -1278,6 +1303,12 @@ int main(int argc, char *argv[])
 	STATUS("Initial partiality calculation...\n");
 	for ( i=0; i<n_crystals; i++ ) {
 		Crystal *cr = crystals[i];
+		if ( force_bandwidth > 0.0 ) {
+			crystal_get_image(cr)->bw = force_bandwidth;
+		}
+		if ( force_radius > 0.0 ) {
+			crystal_set_profile_radius(cr, force_radius);
+		}
 		update_predictions(cr);
 		calculate_partialities(cr, pmodel);
 	}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 60810349..acafb8aa 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -330,6 +330,13 @@ static double residual_f(const gsl_vector *v, void *pp)
 		return GSL_NAN;
 	}
 
+	/* Can happen with grid scans and certain --force-radius values */
+	if ( fabs(crystal_get_profile_radius(cr)) < 0.0000001e9 ) {
+		crystal_free(cr);
+		if ( pv->verbose ) STATUS("radius very small\n");
+		return GSL_NAN;
+	}
+
 	if ( im.lambda <= 0.0 ) {
 		crystal_free(cr);
 		if ( pv->verbose ) STATUS("lambda < 0\n");