aboutsummaryrefslogtreecommitdiff
path: root/src/get_hkl.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-11-06 17:05:57 +0100
committerThomas White <taw@physics.org>2013-11-06 17:05:57 +0100
commit41a8d87549ae245d30940943fe814c9434e88b12 (patch)
tree1cb98083f1acb1f934e5510999c910f3fc6d2931 /src/get_hkl.c
parente6508cf8ffb1a686edfe2cda33a0580bb7bd78ec (diff)
get_hkl: Add anisotropic truncation
Diffstat (limited to 'src/get_hkl.c')
-rw-r--r--src/get_hkl.c104
1 files changed, 98 insertions, 6 deletions
diff --git a/src/get_hkl.c b/src/get_hkl.c
index fa5a0d8c..b4025251 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -375,8 +375,11 @@ int main(int argc, char *argv[])
RefList *input;
double adu_per_photon = 0.0;
int have_adu_per_photon = 0;
- int have_cutoff_angstroms = 0;
- double cutoff = 0.0;
+ int have_cutoff_iso = 0;
+ int have_cutoff_aniso = 0;
+ char *cutoff_str = NULL;
+ double cutiso = 0.0;
+ float cutn1, cutn2, cutn3;
char *pdb = NULL;
/* Long options */
@@ -442,8 +445,7 @@ int main(int argc, char *argv[])
break;
case 3 :
- cutoff = strtof(optarg, NULL);
- have_cutoff_angstroms = 1;
+ cutoff_str = strdup(optarg);
break;
case 0 :
@@ -460,6 +462,36 @@ int main(int argc, char *argv[])
}
+ if ( cutoff_str != NULL ) {
+
+ int r;
+
+ r = sscanf(cutoff_str, "%f,%f,%f", &cutn1, &cutn2, &cutn3);
+ if ( r == 3 ) {
+
+ have_cutoff_aniso = 1;
+
+ /* Convert Angstroms -> m */
+ cutn1 /= 1e10; cutn2 /= 1e10; cutn3 /= 1e10;
+
+ } else {
+
+ char *rval;
+
+ errno = 0;
+ cutiso = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --cutoff-angstroms.\n");
+ return 1;
+ }
+
+ have_cutoff_iso = 1;
+
+ }
+ free(cutoff_str);
+
+ }
+
if ( (holo_str != NULL) && (expand_str != NULL) ) {
ERROR("You cannot 'twin' and 'expand' at the same time.\n");
ERROR("Decide which one you want to do first.\n");
@@ -615,7 +647,7 @@ int main(int argc, char *argv[])
}
- if ( have_cutoff_angstroms ) {
+ if ( have_cutoff_iso ) {
RefList *n;
Reflection *refl;
@@ -645,7 +677,7 @@ int main(int argc, char *argv[])
double res;
get_indices(refl, &h, &k, &l);
res = 2.0 * resolution(cell, h, k, l);
- if ( res < 1e10 / cutoff ) {
+ if ( res < 1e10 / cutiso ) {
Reflection *a;
a = add_refl(n, h, k, l);
copy_data(a, refl);
@@ -658,6 +690,66 @@ int main(int argc, char *argv[])
}
+ if ( have_cutoff_aniso ) {
+
+ RefList *n;
+ Reflection *refl;
+ RefListIterator *iter;
+ UnitCell *cell;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double as, bs, cs;
+
+ if ( pdb == NULL ) {
+ ERROR("You must provide a PDB file when using "
+ "--cutoff-angstroms.\n");
+ return 1;
+ }
+
+ cell = load_cell_from_pdb(pdb);
+ if ( cell == NULL ) {
+ ERROR("Failed to load cell from '%s'\n", pdb);
+ return 1;
+ }
+ free(pdb);
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ as = modulus(asx, asy, asz);
+ bs = modulus(bsx, bsy, bsz);
+ cs = modulus(csx, csy, csz);
+
+ n = reflist_new();
+
+ for ( refl = first_refl(input, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+ double sum;
+
+ get_indices(refl, &h, &k, &l);
+
+ sum = pow(h*as*cutn1, 2.0);
+ sum += pow(k*bs*cutn2, 2.0);
+ sum += pow(l*cs*cutn3, 2.0);
+
+ if ( sum < 1.0 ) {
+ Reflection *a;
+ a = add_refl(n, h, k, l);
+ copy_data(a, refl);
+ }
+ }
+
+ cell_free(cell);
+ reflist_free(input);
+ input = n;
+
+ }
+
+
write_reflist(output, input);
reflist_free(input);