partialator: Add --start-params

1 files changed, 83 insertions, 5 deletions

diff --git a/src/partialator.c b/src/partialator.c
index a756ae3d..4b1d5b18 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -244,6 +244,45 @@ static RefList *apply_max_adu(RefList *list, double max_adu)
 }
 
 
+static void skip_to_end(FILE *fh)
+{
+	int c;
+	do {
+		c = fgetc(fh);
+	} while ( (c != '\n') && (c != EOF) );
+}
+
+
+static int set_initial_params(Crystal *cr, FILE *fh)
+{
+	if ( fh != NULL ) {
+
+		int err;
+		int n;
+		float osf, B;
+
+		err = fscanf(fh, "%i %f %f", &n, &osf, &B);
+		if ( err != 3 ) {
+			ERROR("Failed to read parameters.\n");
+			return 1;
+		}
+
+		crystal_set_osf(cr, osf);
+		crystal_set_Bfac(cr, B*1e-20);
+
+		skip_to_end(fh);
+
+	} else {
+
+		crystal_set_osf(cr, 1.0);
+		crystal_set_Bfac(cr, 0.0);
+
+	}
+
+	return 0;
+}
+
+
 static void show_duds(Crystal **crystals, int n_crystals)
 {
 	int j;
@@ -326,6 +365,8 @@ int main(int argc, char *argv[])
 	struct srdata srdata;
 	int polarisation = 1;
 	double max_adu = +INFINITY;
+	char *sparams_fn = NULL;
+	FILE *sparams_fh;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -341,6 +382,7 @@ int main(int argc, char *argv[])
 
 		{"min-measurements",   1, NULL,                2},
 		{"max-adu",            1, NULL,                3},
+		{"start-params",       1, NULL,                4},
 
 		{"no-scale",           0, &noscale,            1},
 		{"no-polarisation",    0, &polarisation,       0},
@@ -415,6 +457,10 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 4 :
+			sparams_fn = strdup(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -473,6 +519,20 @@ int main(int argc, char *argv[])
 	n_crystals = 0;
 	images = NULL;
 	crystals = NULL;
+	if ( sparams_fn != NULL ) {
+		char line[1024];
+		sparams_fh = fopen(sparams_fn, "r");
+		if ( sparams_fh == NULL ) {
+			ERROR("Failed to open '%s'\n", sparams_fn);
+			return 1;
+		}
+		fgets(line, 1024, sparams_fh);
+		STATUS("Reading initial scaling factors (G,B) from '%s'\n",
+		       sparams_fn);
+		free(sparams_fn);
+	} else {
+		sparams_fh = NULL;
+	}
 
 	do {
 
@@ -540,6 +600,11 @@ int main(int argc, char *argv[])
 			crystal_set_user_flag(cr, 0);
 			reflist_free(cr_refl);
 
+			if ( set_initial_params(cr, sparams_fh) ) {
+				ERROR("Failed to set initial parameters\n");
+				return 1;
+			}
+
 			n_crystals++;
 
 		}
@@ -551,6 +616,7 @@ int main(int argc, char *argv[])
 	} while ( 1 );
 	display_progress(n_images, n_crystals);
 	fprintf(stderr, "\n");
+	if ( sparams_fh != NULL ) fclose(sparams_fh);
 
 	close_stream(st);
 
@@ -581,9 +647,14 @@ int main(int argc, char *argv[])
 
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");
-	if ( noscale ) STATUS("Scale factors fixed at 1.\n");
-	full = scale_intensities(crystals, n_crystals,
-	                         nthreads, noscale, pmodel, min_measurements);
+	if ( noscale ) {
+		STATUS("Skipping scaling step (--no-scale).\n");
+		full = lsq_intensities(crystals, n_crystals, nthreads, pmodel);
+	} else {
+		full = scale_intensities(crystals, n_crystals, nthreads, pmodel,
+		                         min_measurements);
+	}
+
 	check_rejection(crystals, n_crystals);
 
 	srdata.crystals = crystals;
@@ -614,8 +685,15 @@ int main(int argc, char *argv[])
 
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
-		full = scale_intensities(crystals, n_crystals, nthreads,
-		                         noscale, pmodel, min_measurements);
+		if ( noscale ) {
+			STATUS("Skipping scaling step (--no-scale).\n");
+			full = lsq_intensities(crystals, n_crystals, nthreads,
+			                       pmodel);
+		} else {
+			full = scale_intensities(crystals, n_crystals, nthreads,
+			                         pmodel, min_measurements);
+		}
+
 		check_rejection(crystals, n_crystals);
 
 		srdata.full = full;