process_hkl: Take input ("apparent") symmetry on command line

4 files changed, 33 insertions, 52 deletions

diff --git a/src/process_hkl.c b/src/process_hkl.c
index 2f8f3d0a..13a5aa9d 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -66,6 +66,7 @@ static void show_help(const char *s)
 "      --scale               Scale each pattern for best fit with the current\n"
 "                             model.\n"
 "  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
+"  -a, --input-symmetry=<a>  Specify the apparent (input) symmetry.\n"
 "      --reference=<file>    Compare against intensities from <file> when\n"
 "                             scaling or resolving ambiguities.\n"
 "                             The symmetry of the reference list must be the\n"
@@ -552,11 +553,11 @@ int main(int argc, char *argv[])
 	int config_rmerge = 0;
 	unsigned int n_total_patterns;
 	char *sym = NULL;
+	char *in_sym = NULL;
 	char *pdb = NULL;
 	ReflItemList *twins;
 	ReflItemList *observed;
 	int i;
-	const char *holo = NULL;
 	char *histo = NULL;
 	signed int hist_h, hist_k, hist_l;
 	double *hist_vals = NULL;
@@ -580,10 +581,11 @@ int main(int argc, char *argv[])
 		{"sum",                0, &config_sum,         1},
 		{"scale",              0, &config_scale,       1},
 		{"symmetry",           1, NULL,               'y'},
+		{"input-symmetry",     1, NULL,               'a'},
 		{"pdb",                1, NULL,               'p'},
 		{"histogram",          1, NULL,               'g'},
 		{"rmerge",             0, &config_rmerge,      1},
-		{"outstream",          1, NULL,               'a'},
+		{"outstream",          1, NULL,                2},
 		{"reference",          1, NULL,               'r'},
 		{"beam",               1, NULL,               'b'},
 		{0, 0, NULL, 0}
@@ -630,10 +632,14 @@ int main(int argc, char *argv[])
 			reference = strdup(optarg);
 			break;
 
-		case 'a' :
+		case 2 :
 			outstream = strdup(optarg);
 			break;
 
+		case 'a' :
+			in_sym = strdup(optarg);
+			break;
+
 		case 'b' :
 			beam = get_beam_parameters(optarg);
 			if ( beam == NULL ) {
@@ -670,33 +676,35 @@ int main(int argc, char *argv[])
 	cell = load_cell_from_pdb(pdb);
 	free(pdb);
 
+	if ( sym == NULL ) sym = strdup("1");
+
 	/* Show useful symmetry information */
-	if ( sym != NULL ) {
-		holo = get_holohedral(sym);
-		int np = num_general_equivs(holo) / num_general_equivs(sym);
+	if ( (sym != NULL) && (in_sym != NULL) ) {
+
+		int np = num_general_equivs(in_sym) / num_general_equivs(sym);
 		if ( np > 1 ) {
 
 			STATUS("Resolving point group %s into %s "
 			       "(%i possibilities)\n",
-			       holo, sym, np);
+			       in_sym, sym, np);
 			/* Get the list of twin/Bijvoet possibilities */
-			twins = get_twin_possibilities(holo, sym);
-			STATUS("Twin/inversion operation indices from %s are:",
-			       holo);
+			twins = get_twin_possibilities(in_sym, sym);
+			STATUS("Twin operation indices from %s are:", in_sym);
 			for ( i=0; i<num_items(twins); i++ ) {
 				STATUS(" %i", get_item(twins, i)->op);
 			}
 			STATUS("\n");
 
 		} else {
-			STATUS("No twin/inversion resolution necessary.\n");
+			STATUS("No resolution required to get from %s to %s\n",
+			       in_sym, sym);
 			twins = NULL;
 		}
+
 	} else {
-		STATUS("Not performing any twin/inversion resolution.\n");
+		STATUS("No twin resolution requested (use -a otherwise).\n");
 		twins = NULL;
-		sym = strdup("1");
-		holo = strdup("1");
+		in_sym = strdup(sym);
 	}
 
 	if ( histo != NULL ) {
@@ -758,7 +766,7 @@ int main(int argc, char *argv[])
 	merge_all(fh, &model, &observed, &counts,
 	          config_maxonly, config_scale, config_sum,
 	          config_startafter, config_stopafter,
-                  twins, holo, sym, n_total_patterns,
+                  twins, in_sym, sym, n_total_patterns,
                   reference_items, reference_i,
                   hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL,
                   outfh);
@@ -801,7 +809,8 @@ int main(int argc, char *argv[])
 		rewind(fh);
 		merge_all(fh, &model, &observed, &counts,
 		          config_maxonly, config_scale, 0,
-		          config_startafter, config_stopafter, twins, holo, sym,
+		          config_startafter, config_stopafter,
+		          twins, in_sym, sym,
 		          n_total_patterns, reference_items, reference_i,
 		          NULL, 0, 0, 0, NULL, devs, tots, model, NULL);
 
diff --git a/src/symmetry.c b/src/symmetry.c
index 3e9cf5ed..48194c8c 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -328,31 +328,6 @@ void get_asymm(signed int h, signed int k, signed int l,
 }
 
 
-const char *get_holohedral(const char *sym)
-{
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 ) return "-1";
-	if ( strcmp(sym, "-1") == 0 ) return "-1";
-
-	if ( strcmp(sym, "222") == 0 ) return "mmm";
-	if ( strcmp(sym, "mmm") == 0 ) return "mmm";
-
-	/* Tetragonal */
-	if ( strcmp(sym, "422") == 0 ) return "4/mmm";
-	if ( strcmp(sym, "4/mmm") == 0 ) return "4/mmm";
-
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 ) return "6/mmm";
-	if ( strcmp(sym, "6/m") == 0 ) return "6/mmm";
-	if ( strcmp(sym, "6/mmm") == 0 ) return "6/mmm";
-
-	/* TODO: Add more groups here */
-
-	ERROR("Couldn't find holohedral point group for '%s'\n", sym);
-	abort();
-}
-
-
 /* This is kind of like a "numerical" left coset decomposition.
  * Given a reflection index and a point group, it returns the list of twinning
  * possibilities.
diff --git a/src/symmetry.h b/src/symmetry.h
index c338da49..fc1dcaae 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -35,8 +35,6 @@ extern void get_general_equiv(signed int h, signed int k, signed int l,
                               signed int *he, signed int *ke, signed int *le,
                               const char *sym, int idx);
 
-extern const char *get_holohedral(const char *sym);
-
 extern ReflItemList *get_twins(ReflItemList *items,
                                const char *holo, const char *mero);
 
diff --git a/src/templates.c b/src/templates.c
index 4f47c572..4bd85a76 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -113,7 +113,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
                                     double nominal_photon_energy)
 {
 	struct _indexingprivate_template *priv;
-	const char *holo;
+	const char *sym;
 	double omega_max, phi_max;
 	int n_templates;
 	const double omega_step = deg2rad(0.5);
@@ -137,24 +137,23 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
 	priv = calloc(1, sizeof(struct _indexingprivate_template));
 	priv->base.indm = INDEXING_TEMPLATE;
 
-	/* We can only distinguish orientations within the holohedral cell */
-	holo = get_holohedral(cell_get_pointgroup(cell));
+	sym = cell_get_pointgroup(cell);
 
 	/* These define the orientation in space */
-	if ( is_polyhedral(holo) ) {
-		ERROR("WARNING: Holohedral point group is polyhedral.\n");
+	if ( is_polyhedral(sym) ) {
+		ERROR("WARNING: Point group is polyhedral.\n");
 		ERROR("This means I can't properly determine the orientation");
 		ERROR(" ranges for template matching.  Expect trouble.\n");
 	}
-	omega_max = 2.0*M_PI / rotational_order(holo);
-	if ( has_bisecting_mirror_or_diad(holo) ) omega_max /= 2.0;
+	omega_max = 2.0*M_PI / rotational_order(sym);
+	if ( has_bisecting_mirror_or_diad(sym) ) omega_max /= 2.0;
 	phi_max = M_PI;
-	if ( has_perpendicular_mirror(holo) ) phi_max /= 2.0;
+	if ( has_perpendicular_mirror(sym) ) phi_max /= 2.0;
 
 	/* One more axis would define the rotation in the plane of the image */
 
 	STATUS("Orientation ranges in %s: %.0f-%.0f, %.0f-%.0f deg.\n",
-	       holo, 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max));
+	       sym, 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max));
 
 	n_templates = (omega_max * phi_max)/(omega_step * phi_step);
 	STATUS("%i templates to be calculated.\n", n_templates);