aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--doc/man/check_hkl.16
-rw-r--r--doc/man/indexamajig.112
-rw-r--r--doc/man/pattern_sim.1106
-rw-r--r--libcrystfel/src/geometry.c163
-rw-r--r--libcrystfel/src/peaks.c34
-rw-r--r--src/dw-hdfsee.c8
-rw-r--r--src/partialator.c2
-rw-r--r--src/post-refinement.c32
-rw-r--r--src/process_hkl.c3
9 files changed, 224 insertions, 142 deletions
diff --git a/doc/man/check_hkl.1 b/doc/man/check_hkl.1
index df1f857e..39dd97ac 100644
--- a/doc/man/check_hkl.1
+++ b/doc/man/check_hkl.1
@@ -37,12 +37,6 @@ Specify the symmetry of the reflections.
Discard reflections with I/sigma(I) < \fIn\fR. Default: -infinity (no cutoff).
.PD 0
-.IP \fB-p\fR \fIunitcell.pdb\fR
-.IP \fB--pdb=\fR\fIunitcell.pdb\fR
-.PD
-Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell.
-
-.PD 0
.IP \fB--rmin=\fR\fI1/d\fR
.PD
Fix the lower resolution limit for the resolutions shells, as 1/d in m^-1.
diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index ac443906..5bc3f7c0 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -216,7 +216,6 @@ these. The defaults are probably not appropriate for your situation.
.PD
The default is \fB--int-radius=4,5,7\fR.
-
.PD 0
.IP \fB--basename\fR
.PD
@@ -296,9 +295,18 @@ Run \fIn\fR analyses in parallel. Default: 1.
Don't attempt to correct the prefix (see \fB--prefix\fR) if it doesn't look correct.
.PD 0
+.IP \fB--closer-peak\fR
+.PD
+If you use this option, indexamajig will integrate around the location of a detected peak instead of the predicted peak location if one is found close to the predicted position, within ten pixels. \fBDon't use this option\fR, because
+there is currently no way to set the definition of 'nearby' to be appropriate
+for your data.
+
+.PD 0
.IP \fB--no-closer-peak\fR
.PD
-Normally, indexamajig will integrate around the location of a detected peak instead of the predicted peak location if one is found close to the predicted position. This option disables this behaviour. Normally it doesn't make much difference either way.
+This is the opposite of \fB--closer-peak\fR, and is provided for compatibility
+with old scripts because this option selects the behaviour which is now the
+default.
.PD 0
.IP \fB--insane\fR
diff --git a/doc/man/pattern_sim.1 b/doc/man/pattern_sim.1
index b0d2fa06..bd2b35e4 100644
--- a/doc/man/pattern_sim.1
+++ b/doc/man/pattern_sim.1
@@ -30,6 +30,96 @@ The result will be written to an HDF5 file in the current directory with the nam
.SH OPTIONS
+.PD 0
+.IP "\fB-p\fR \fIunitcell.pdb\fR"
+.IP \fB--pdb=\fR\fIunitcell.pdb\fR
+.PD
+Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell.
+
+.PD 0
+.IP \fB--gpu\fR
+.PD
+Use the GPU to speed up the calculation. Requires that OpenCL libraries and drivers are available, and that CrystFEL was compiled with OpenCL enabled.
+
+.PD 0
+.IP \fB--gpu-dev=\fRIn\fR
+.PD
+Use GPU device number \fIn\fR. If you omit this option, the list of GPU devices will be shown when you run pattern_sim.
+
+.PD 0
+.IP "\fB-g\fR \fIfilename\fR"
+.IP \fB--geometry=\fR\fIfilename\fR
+.PD
+Read the detector geometry description from \fIfilename\fR. See \fBman crystfel_geometry\fR for more information.
+
+.PD 0
+.IP "\fB-b\fR \fIfilename\fR"
+.IP \fB--beam=\fR\fIfilename\fR
+.PD
+Read the beam description from \fIfilename\fR. See \fBman crystfel_geometry\fR for more information.
+
+.PD 0
+.IP "\fB-n\fR \fn\fR"
+.IP \fB--number=\fR\fIn\fR
+.PD
+Simulate \fIn\fR patterns. Default: \fB-n 1\fR.
+
+.PD 0
+.IP \fB--no-images\fR
+.PD
+Do not save any HDF5 files apart from the powder pattern (if requested).
+
+.PD 0
+.IP "\fB-o\fR \fIfilename\fR"
+.IP \fB--output=\fR\fIfilename\fR
+.PD
+Write the pattern to \fIfilename\fR. The default is \fB--output=sim.5\fR. If more than one pattern is to be simulated (see \fB--number\fR), the filename will be postfixed with a hyphen, the image number and then '.h5'.
+
+.PD 0
+.IP \fB-r\fR
+.IP \fB--random-orientation\fR
+.PD
+Make up a random orientation for each pattern simulated.
+
+.PD 0
+.IP \fB--powder=\fR\fIfilename\fR
+.PD
+Write the sum of all patterns to \fIfilename\fR.
+
+.PD 0
+.IP "\fB-i\fR \ffilename\fR"
+.IP \fB--intensities=\fR\fIfilename\fR
+.PD
+Get the intensities and phases at the reciprocal lattice points from \fIfilename\fR.
+
+.PD 0
+.IP "\fB-y\fR \fIpointgroup\fR"
+.IP \fB--symmetry=\fR\fIpointgroup\fR
+.PD
+Use \fIpointgroup\fR as the symmetry of the intensity list (see \fB--intensities\fR).
+
+.PD 0
+.IP "\fB-t\fR \fImethod\fR"
+.IP \fB--gradients=\fR\fImethod\fR
+.PD
+Use \fImethod\fR as way of calculating the molecular transform between reciprocal lattice points. See the section \fBGRADIENT METHODS\fR below.
+
+.PD 0
+.IP \fB--really-random\fR
+.PD
+Seed the random number generator using the kernel random number generator (/dev/urandom). This means that truly random numbers for the orientation and crystal size, instead of the same sequence being used for each new run.
+
+.PD 0
+.IP \fB--min-size=\fR\fImin\fR
+.IP \fB--min-size=\fR\fImax\fR
+.PD
+Generate random crystal sizes between \fImin\fR and \fImax\fR nanometres. These options must be used together.
+
+.PD 0
+.IP \fB--no-noise\fR
+.PD
+Do not calculate Poisson noise.
+
.SH REFLECTION LISTS
@@ -76,6 +166,22 @@ algorithm. When the intensity is sufficiently high that Knuth's algorithm
would result in machine precision problems, a normal distribution with
standard deviation sqrt(I) is used instead.
+.SH GRADIENT METHODS
+
+The available options for \fB--gradients\fR as as follows:
+
+.IP \fBmosaic\fR
+.PD
+Take the intensity of the nearest Bragg position. This is the fastest method and the only one supported on the GPU, but the least accurate.
+
+.IP \fBinterpolate\fR
+.PD
+Interpolate trilinearly between six adjacent Bragg intensities. This method has intermediate accuracy.
+
+.IP \fBphased\fR
+.PD
+As 'interpolate', but take phase values into account. This is the most accurate method, but the slowest.
+
.SH AUTHOR
This page was written by Thomas White.
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 218c0fee..24669aef 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -99,112 +99,70 @@ static signed int locate_peak(double x, double y, double z, double k,
}
-static double excitation_error(double xl, double yl, double zl,
- double ds, double k, double divergence,
- double tt)
+static double partiality(double rlow, double rhigh, double r)
{
- double al;
- double r;
- double delta;
-
- al = M_PI_2 - asin(-zl/ds);
-
- r = ( ds * sin(al) / sin(tt) ) - k;
-
- delta = sqrt(2.0 * pow(ds, 2.0) * (1.0-cos(divergence)));
- if ( divergence > 0.0 ) {
- r += delta;
- } else {
- r -= delta;
- }
-
- return r;
-}
-
-
-static double partiality(double r1, double r2, double r)
-{
- double q1, q2;
- double p1, p2;
+ double qlow, qhigh;
+ double plow, phigh;
/* Calculate degrees of penetration */
- q1 = (r1 + r)/(2.0*r);
- q2 = (r2 + r)/(2.0*r);
+ qlow = (rlow + r)/(2.0*r);
+ qhigh = (rhigh + r)/(2.0*r);
/* Convert to partiality */
- p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
- p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
+ plow = 3.0*pow(qlow,2.0) - 2.0*pow(qlow,3.0);
+ phigh = 3.0*pow(qhigh,2.0) - 2.0*pow(qhigh,3.0);
- return p2 - p1;
+ return plow - phigh;
}
static Reflection *check_reflection(struct image *image,
signed int h, signed int k, signed int l,
- double asx, double asy, double asz,
- double bsx, double bsy, double bsz,
- double csx, double csy, double csz)
+ double xl, double yl, double zl)
{
const int output = 0;
- double xl, yl, zl;
- double ds, ds_sq;
+ double tl;
double rlow, rhigh; /* "Excitation error" */
signed int p; /* Panel number */
double xda, yda; /* Position on detector */
- int close, inside;
double part; /* Partiality */
- int clamp_low = 0;
- int clamp_high = 0;
- double bandwidth = image->bw;
- double divergence = image->div;
- double lambda = image->lambda;
- double klow, kcen, khigh; /* Wavenumber */
+ int clamp_low, clamp_high;
+ double klow, khigh; /* Wavenumber */
Reflection *refl;
- double tt, psi;
-
- /* "low" gives the largest Ewald sphere,
- * "high" gives the smallest Ewald sphere. */
- klow = 1.0/(lambda - lambda*bandwidth/2.0);
- kcen = 1.0/lambda;
- khigh = 1.0/(lambda + lambda*bandwidth/2.0);
-
- /* Get the coordinates of the reciprocal lattice point */
- xl = h*asx + k*bsx + l*csx;
- yl = h*asy + k*bsy + l*csy;
- zl = h*asz + k*bsz + l*csz;
-
- ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
- ds = sqrt(ds_sq);
+ double cet, cez;
- /* Simple (fast) check to reject reflection if it's "in front" */
- psi = atan2(zl, sqrt(xl*xl + yl*yl));
- if ( psi - atan2(image->profile_radius, ds) > divergence ) return NULL;
+ /* "low" gives the largest Ewald sphere (wavelength short => k large)
+ * "high" gives the smallest Ewald sphere (wavelength long => k small)
+ */
+ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
+ khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
- tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+kcen);
- if ( tt > deg2rad(90.0) ) return NULL;
+ /* If the point is looking "backscattery", reject it straight away */
+ if ( zl < -khigh/2.0 ) return NULL;
- /* Calculate excitation errors */
- rlow = excitation_error(xl, yl, zl, ds, klow, -divergence/2.0, tt);
- rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence/2.0, tt);
+ tl = sqrt(xl*xl + yl*yl);
- /* Is the reciprocal lattice point close to either extreme of
- * the sphere, maybe just outside the "Ewald volume"? */
- close = (fabs(rlow) < image->profile_radius)
- || (fabs(rhigh) < image->profile_radius);
+ cet = -sin(image->div/2.0) * khigh;
+ cez = -cos(image->div/2.0) * khigh;
+ rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */
- /* Is the reciprocal lattice point somewhere between the
- * extremes of the sphere, i.e. inside the "Ewald volume"? */
- inside = signbit(rlow) ^ signbit(rhigh);
+ cet = sin(image->div/2.0) * klow;
+ cez = -cos(image->div/2.0) * klow;
+ rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */
- /* Can't be both inside and close */
- if ( inside ) close = 0;
-
- /* Neither? Skip it. */
- if ( !(close || inside) ) return NULL;
+ if ( (signbit(rlow) == signbit(rhigh))
+ && (fabs(rlow) > image->profile_radius)
+ && (fabs(rhigh) > image->profile_radius) ) return NULL;
/* If the "lower" Ewald sphere is a long way away, use the
* position at which the Ewald sphere would just touch the
- * reflection. */
+ * reflection.
+ *
+ * The six possible combinations of clamp_{low,high} (including
+ * zero) correspond to the six situations in Table 3 of Rossmann
+ * et al. (1979).
+ */
+ clamp_low = 0; clamp_high = 0;
if ( rlow < -image->profile_radius ) {
rlow = -image->profile_radius;
clamp_low = -1;
@@ -213,7 +171,6 @@ static Reflection *check_reflection(struct image *image,
rlow = +image->profile_radius;
clamp_low = +1;
}
- /* Likewise the "higher" Ewald sphere */
if ( rhigh < -image->profile_radius ) {
rhigh = -image->profile_radius;
clamp_high = -1;
@@ -222,16 +179,13 @@ static Reflection *check_reflection(struct image *image,
rhigh = +image->profile_radius;
clamp_high = +1;
}
- assert(clamp_low <= clamp_high);
- /* The six possible combinations of clamp_{low,high} (including
- * zero) correspond to the six situations in Table 3 of Rossmann
- * et al. (1979). */
+ assert(clamp_low >= clamp_high);
/* Calculate partiality */
part = partiality(rlow, rhigh, image->profile_radius);
/* Locate peak on detector. */
- p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
+ p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det, &xda, &yda);
if ( p == -1 ) return NULL;
/* Add peak to list */
@@ -252,6 +206,9 @@ static Reflection *check_reflection(struct image *image,
RefList *find_intersections(struct image *image, UnitCell *cell)
{
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
@@ -270,17 +227,13 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
return NULL;
}
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- /* We add a horrific 20% fudge factor because bandwidth, divergence
- * and so on mean reflections appear beyond the largest q */
- mres = 1.2 * largest_q(image);
+ mres = largest_q(image);
- hmax = mres / modulus(asx, asy, asz);
- kmax = mres / modulus(bsx, bsy, bsz);
- lmax = mres / modulus(csx, csy, csz);
+ hmax = mres * modulus(ax, ay, az);
+ kmax = mres * modulus(bx, by, bz);
+ lmax = mres * modulus(cx, cy, cz);
if ( (hmax >= 256) || (kmax >= 256) || (lmax >= 256) ) {
ERROR("Unit cell is stupidly large.\n");
@@ -290,14 +243,23 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
if ( lmax >= 256 ) lmax = 255;
}
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
for ( h=-hmax; h<=hmax; h++ ) {
for ( k=-kmax; k<=kmax; k++ ) {
for ( l=-lmax; l<=lmax; l++ ) {
Reflection *refl;
+ double xl, yl, zl;
+
+ /* Get the coordinates of the reciprocal lattice point */
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+ zl = h*asz + k*bsz + l*csz;
- refl = check_reflection(image, h, k, l,
- asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+ refl = check_reflection(image, h, k, l, xl, yl, zl);
if ( refl != NULL ) {
add_refl_to_list(refl, reflections);
@@ -333,13 +295,18 @@ void update_partialities(struct image *image)
{
Reflection *vals;
double r1, r2, p, x, y;
+ double xl, yl, zl;
signed int h, k, l;
int clamp1, clamp2;
get_symmetric_indices(refl, &h, &k, &l);
- vals = check_reflection(image, h, k, l,
- asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+ /* Get the coordinates of the reciprocal lattice point */
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+ zl = h*asz + k*bsz + l*csz;
+
+ vals = check_reflection(image, h, k, l, xl, yl, zl);
if ( vals == NULL ) {
set_redundancy(refl, 0);
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 0d4ce64b..dac1a96e 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -52,14 +52,6 @@
#include "beam-parameters.h"
-/* How close a peak must be to an indexed position to be considered "close"
- * for the purposes of double hit detection and sanity checking. */
-#define PEAK_CLOSE (30.0)
-
-/* How close a peak must be to an indexed position to be considered "close"
- * for the purposes of integration. */
-#define PEAK_REALLY_CLOSE (10.0)
-
/* Degree of polarisation of X-ray beam */
#define POL (1.0)
@@ -184,8 +176,8 @@ int integrate_peak(struct image *image, int cfs, int css,
out_lim_sq = pow(ir_out, 2.0);
/* Estimate the background */
- for ( fs=-ir_out; fs<+ir_out; fs++ ) {
- for ( ss=-ir_out; ss<+ir_out; ss++ ) {
+ for ( fs=-ir_out; fs<=+ir_out; fs++ ) {
+ for ( ss=-ir_out; ss<=+ir_out; ss++ ) {
double val;
uint16_t flags;
@@ -233,8 +225,8 @@ int integrate_peak(struct image *image, int cfs, int css,
pk_total = 0.0;
pk_counts = 0;
fsct = 0.0; ssct = 0.0;
- for ( fs=-ir_inn; fs<+ir_inn; fs++ ) {
- for ( ss=-ir_inn; ss<+ir_inn; ss++ ) {
+ for ( fs=-ir_inn; fs<=+ir_inn; fs++ ) {
+ for ( ss=-ir_inn; ss<=+ir_inn; ss++ ) {
double val;
uint16_t flags;
@@ -541,9 +533,6 @@ int peak_sanity_check(struct image *image)
struct integr_ind
{
- signed int h;
- signed int k;
- signed int l;
double res;
Reflection *refl;
};
@@ -585,9 +574,6 @@ static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell,
get_indices(refl, &h, &k, &l);
res = resolution(cell, h, k, l);
- il[i].h = h;
- il[i].k = k;
- il[i].l = l;
il[i].res = res;
il[i].refl = refl;
@@ -627,10 +613,12 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
double pfs, pss;
int r;
Reflection *refl;
+ signed int h, k, l;
refl = il[i].refl;
get_detector_pos(refl, &pfs, &pss);
+ get_indices(refl, &h, &k, &l);
/* Is there a really close feature which was detected? */
if ( use_closer ) {
@@ -643,11 +631,19 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
} else {
f = NULL;
}
- if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) {
+
+ /* FIXME: Horrible hardcoded value */
+ if ( (f != NULL) && (d < 10.0) ) {
+
+ double exe;
+
+ exe = get_excitation_error(refl);
pfs = f->fs;
pss = f->ss;
+ set_detector_pos(refl, exe, pfs, pss);
+
}
}
diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c
index 3881d162..de03c62b 100644
--- a/src/dw-hdfsee.c
+++ b/src/dw-hdfsee.c
@@ -1340,7 +1340,6 @@ static void numbers_update(DisplayWindow *dw)
for ( py=0; py<17; py++ ) {
char s[32];
- float val;
GtkWidget *l;
int x, y;
int invalid;
@@ -1355,11 +1354,13 @@ static void numbers_update(DisplayWindow *dw)
/* Map from unbinned mapped pixel coordinates to a panel */
invalid = reverse_2d_mapping(x, y, &dfs, &dss, dw->image->det);
fs = dfs; ss = dss;
+
if ( !invalid ) {
+
+ float val;
+
val = dw->image->data[fs+ss*dw->image->width];
- }
- if ( !invalid ) {
if ( val > 0.0 ) {
if ( log(val)/log(10.0) < 5 ) {
snprintf(s, 31, "%.0f", val);
@@ -1373,6 +1374,7 @@ static void numbers_update(DisplayWindow *dw)
snprintf(s, 31, "-HUGE");
}
}
+
} else {
strcpy(s, "-");
}
diff --git a/src/partialator.c b/src/partialator.c
index fa0698b5..91f2ff99 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -320,7 +320,7 @@ int main(int argc, char *argv[])
}
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:r:",
+ while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:",
longopts, NULL)) != -1)
{
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 2ff408eb..cf7b2576 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -147,7 +147,7 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
* of excitation error wrt whatever. */
switch ( k ) {
- case REF_DIV :
+ case REF_DIV :
gr = 0.0;
if ( clamp_low == 0 ) {
nom = sqrt(2.0) * ds * sin(image->div/2.0);
@@ -162,7 +162,7 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
if ( isnan(gr) ) gr = 0.0; /* FIXME: This isn't true (?) */
return gr / 4.0; /* FIXME: Shameless fudge factor */
- case REF_R :
+ case REF_R :
g = 0.0;
if ( clamp_low == 0 ) {
g += partiality_rgradient(r1, r);
@@ -172,24 +172,32 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
}
return g;
- /* Cell parameters and orientation */
- case REF_ASX :
+ /* Cell parameters and orientation */
+ case REF_ASX :
return hs * sin(tt) * cos(azix) * g;
- case REF_BSX :
+
+ case REF_BSX :
return ks * sin(tt) * cos(azix) * g;
- case REF_CSX :
+
+ case REF_CSX :
return ls * sin(tt) * cos(azix) * g;
- case REF_ASY :
+
+ case REF_ASY :
return hs * sin(tt) * cos(aziy) * g;
- case REF_BSY :
+
+ case REF_BSY :
return ks * sin(tt) * cos(aziy) * g;
- case REF_CSY :
+
+ case REF_CSY :
return ls * sin(tt) * cos(aziy) * g;
- case REF_ASZ :
+
+ case REF_ASZ :
return hs * cos(tt) * g;
- case REF_BSZ :
+
+ case REF_BSZ :
return ks * cos(tt) * g;
- case REF_CSZ :
+
+ case REF_CSZ :
return ls * cos(tt) * g;
}
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 4164d6f6..26abfd4f 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -502,7 +502,8 @@ int main(int argc, char *argv[])
hist_i = 0;
merge_all(fh, model, NULL, config_startafter, config_stopafter,
- sym, n_total_patterns, NULL, 0, 0, 0, NULL);
+ sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l,
+ &hist_i);
if ( ferror(fh) ) {
ERROR("Stream read error.\n");
return 1;