diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/check_hkl.c | 17 | ||||
| -rw-r--r-- | src/get_hkl.c | 104 |
2 files changed, 112 insertions, 9 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c index 71f107da..c2ab0788 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -59,12 +59,14 @@ static void show_help(const char *s) " --rmax=<res> Fix upper resolution limit for resolution shells. (m^-1).\n" " --sigma-cutoff=<n> Discard reflections with I/sigma(I) < n.\n" " --nshells=<n> Use <n> resolution shells.\n" +" --shell-file=<file> Write resolution shells to <file>.\n" "\n"); } static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, - double rmin_fix, double rmax_fix, int nshells) + double rmin_fix, double rmax_fix, int nshells, + const char *shell_file) { int *possible; unsigned int *measurements; @@ -127,7 +129,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, return; } - fh = fopen("shells.dat", "w"); + fh = fopen(shell_file, "w"); if ( fh == NULL ) { ERROR("Couldn't open 'shells.dat'\n"); return; @@ -371,6 +373,7 @@ int main(int argc, char *argv[]) float rmax_fix = -1.0; float sigma_cutoff = -INFINITY; int nshells = 10; + char *shell_file = NULL; /* Long options */ const struct option longopts[] = { @@ -381,6 +384,7 @@ int main(int argc, char *argv[]) {"rmax", 1, NULL, 3}, {"sigma-cutoff", 1, NULL, 4}, {"nshells", 1, NULL, 5}, + {"shell-file", 1, NULL, 6}, {0, 0, NULL, 0} }; @@ -432,6 +436,10 @@ int main(int argc, char *argv[]) } break; + case 6 : + shell_file = strdup(optarg); + break; + case '?' : break; @@ -470,6 +478,8 @@ int main(int argc, char *argv[]) } free(file); + if ( shell_file == NULL ) shell_file = strdup("shells.dat"); + /* Check that the intensities have the correct symmetry */ if ( check_list_symmetry(raw_list, sym) ) { ERROR("The input reflection list does not appear to" @@ -508,11 +518,12 @@ int main(int argc, char *argv[]) rej, num_reflections(raw_list), sigma_cutoff); reflist_free(raw_list); - plot_shells(list, cell, sym, rmin_fix, rmax_fix, nshells); + plot_shells(list, cell, sym, rmin_fix, rmax_fix, nshells, shell_file); free_symoplist(sym); reflist_free(list); cell_free(cell); + free(shell_file); return 0; } 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); |