aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/crystal.c77
-rw-r--r--libcrystfel/src/crystal.h3
-rw-r--r--libcrystfel/src/dirax.c5
-rw-r--r--libcrystfel/src/image.c18
-rw-r--r--libcrystfel/src/mosflm.c2
-rw-r--r--libcrystfel/src/reax.c66
-rw-r--r--libcrystfel/src/stream.c64
-rw-r--r--libcrystfel/src/stream.h5
-rw-r--r--src/process_hkl.c179
9 files changed, 256 insertions, 163 deletions
diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c
index 3c0649cc..8ba955f6 100644
--- a/libcrystfel/src/crystal.c
+++ b/libcrystfel/src/crystal.c
@@ -55,7 +55,7 @@ struct _crystal
double osf;
double profile_radius;
int pr_dud;
- double diffracting_resolution;
+ double resolution_limit;
/* Integrated (or about-to-be-integrated) reflections */
RefList *reflections;
@@ -63,7 +63,7 @@ struct _crystal
};
-/************************** Setters and Constructors **************************/
+/************************** Constructor / Destructor **************************/
/**
@@ -83,7 +83,7 @@ Crystal *crystal_new()
cryst->cell = NULL;
cryst->reflections = NULL;
- cryst->diffracting_resolution = 0.0;
+ cryst->resolution_limit = 0.0;
cryst->n_saturated = 0;
return cryst;
@@ -97,11 +97,78 @@ Crystal *crystal_new()
* Frees a %Crystal, and all internal resources concerning that crystal.
*
*/
-void crystal_free(UnitCell *cryst)
+void crystal_free(Crystal *cryst)
{
if ( cryst == NULL ) return;
- free(crysta);
+ if ( cryst->cell != NULL ) cell_free(cryst->cell);
+ if ( cryst->reflections != NULL ) reflist_free(cryst->reflections);
+ free(cryst);
}
/********************************** Getters ***********************************/
+
+
+UnitCell *crystal_get_cell(Crystal *cryst)
+{
+ return cryst->cell;
+}
+
+
+double crystal_get_profile_radius(Crystal *cryst)
+{
+ return cryst->profile_radius;
+}
+
+
+RefList *crystal_get_reflections(Crystal *cryst)
+{
+ return cryst->reflections;
+}
+
+
+double crystal_get_resolution_limit(Crystal *cryst)
+{
+ return cryst->resolution_limit;
+}
+
+
+long long int crystal_get_num_saturated_reflections(Crystal *cryst)
+{
+ return cryst->n_saturated;
+}
+
+
+/********************************** Setters ***********************************/
+
+
+void crystal_set_cell(Crystal *cryst, UnitCell *cell)
+{
+ if ( cryst->cell != NULL ) cell_free(cryst->cell);
+ cryst->cell = cell_new_from_cell(cell);
+}
+
+
+void crystal_set_profile_radius(Crystal *cryst, double r)
+{
+ cryst->profile_radius = r;
+}
+
+
+void crystal_set_reflections(Crystal *cryst, RefList *reflist)
+{
+ if ( cryst->reflections != NULL ) reflist_free(reflist);
+ cryst->reflections = reflist;
+}
+
+
+void crystal_set_resolution_limit(Crystal *cryst, double res)
+{
+ cryst->resolution_limit = res;
+}
+
+
+void crystal_set_num_saturated_reflections(Crystal *cryst, long long int n)
+{
+ cryst->n_saturated = n;
+}
diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h
index e5e40589..ca215399 100644
--- a/libcrystfel/src/crystal.h
+++ b/libcrystfel/src/crystal.h
@@ -60,7 +60,8 @@ extern void crystal_set_cell(Crystal *cryst, UnitCell *cell);
extern void crystal_set_profile_radius(Crystal *cryst, double r);
extern void crystal_set_reflections(Crystal *cryst, RefList *reflist);
extern void crystal_set_resolution_limit(Crystal *cryst, double res);
-extern void crystal_set_num_saturated_reflections(Crystal *cryst, int n);
+extern void crystal_set_num_saturated_reflections(Crystal *cryst,
+ long long int n);
#endif /* CRYSTAL_H */
diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c
index 0b2debda..6fbb9663 100644
--- a/libcrystfel/src/dirax.c
+++ b/libcrystfel/src/dirax.c
@@ -128,17 +128,20 @@ static int check_cell(struct dirax_private *dp, struct image *image,
return 0;
}
- crystal_set_cell(cr, cell);
+ crystal_set_cell(cr, out);
if ( dp->indm & INDEXING_CHECK_PEAKS ) {
if ( !peak_sanity_check(image, &cr, 1) ) {
crystal_free(cr); /* Frees the cell as well */
+ cell_free(out);
return 0;
}
}
image_add_crystal(image, cr);
+ cell_free(out); /* Crystal makes its own copy */
+
return 1;
}
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 093e20f2..6e4fd40a 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -159,3 +159,21 @@ void image_remove_feature(ImageFeatureList *flist, int idx)
{
flist->features[idx].valid = 0;
}
+
+
+void image_add_crystal(struct image *image, Crystal *cryst)
+{
+ Crystal **crs;
+ int n;
+
+ n = image->n_crystals;
+ crs = realloc(image->crystals, (n+1)*sizeof(Crystal *));
+ if ( crs == NULL ) {
+ ERROR("Failed to allocate memory for crystals.\n");
+ return;
+ }
+
+ crs[n] = cryst;
+ image->crystals = crs;
+ image->n_crystals = n+1;
+}
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index 8b300012..d1d3c134 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -149,7 +149,7 @@ static int check_cell(struct mosflm_private *mp, struct image *image,
return 0;
}
- crystal_set_cell(cr, cell);
+ crystal_set_cell(cr, out);
if ( mp->indm & INDEXING_CHECK_PEAKS ) {
if ( !peak_sanity_check(image, &cr, 1) ) {
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 7529d12f..b144e266 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -103,7 +103,7 @@ struct reax_search
struct reax_private
{
- IndexingPrivate base;
+ IndexingMethod indm;
struct dvec *directions;
int n_dir;
double angular_inc;
@@ -915,9 +915,39 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew,
}
+static int check_cell(struct reax_private *dp, struct image *image,
+ UnitCell *cell)
+{
+ UnitCell *out;
+ Crystal *cr;
+
+ out = cell_new_from_cell(cell);
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+
+ crystal_set_cell(cr, out);
+
+ if ( dp->indm & INDEXING_CHECK_PEAKS ) {
+ if ( !peak_sanity_check(image, &cr, 1) ) {
+ crystal_free(cr); /* Frees the cell as well */
+ return 0;
+ }
+ }
+
+ image_add_crystal(image, cr);
+
+ return 1;
+}
+
+
static void assemble_cells_from_candidates(struct image *image,
struct reax_search *s,
- UnitCell *cell)
+ UnitCell *cell,
+ struct reax_private *p)
{
int i, j, k;
signed int ti, tj, tk;
@@ -967,7 +997,9 @@ static void assemble_cells_from_candidates(struct image *image,
continue;
}
- peak_lattice_agreement(image, cnew, &fom);
+ /* FIXME! */
+ //peak_lattice_agreement(image, cnew, &fom);
+ fom = 1.0;
add_cell_candidate(&cl, cnew, fom);
}
@@ -985,22 +1017,20 @@ static void assemble_cells_from_candidates(struct image *image,
cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga);
cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA,
&alA, &beA, &gaA);
- if ( (a - aA) > aA/10.0 ) w = 1;
- if ( (b - bA) > bA/10.0 ) w = 1;
- if ( (c - cA) > cA/10.0 ) w = 1;
- if ( (al - alA) > deg2rad(5.0) ) w = 1;
- if ( (be - beA) > deg2rad(5.0) ) w = 1;
- if ( (ga - gaA) > deg2rad(5.0) ) w = 1;
+ if ( fabs(a - aA) > aA/10.0 ) w = 1;
+ if ( fabs(b - bA) > bA/10.0 ) w = 1;
+ if ( fabs(c - cA) > cA/10.0 ) w = 1;
+ if ( fabs(al - alA) > deg2rad(5.0) ) w = 1;
+ if ( fabs(be - beA) > deg2rad(5.0) ) w = 1;
+ if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1;
if ( w ) {
STATUS("This cell is a long way from that sought:\n");
cell_print(cl.cand[i].cell);
}
}
- image->ncells = cl.n_cand;
- assert(image->ncells <= MAX_CELL_CANDIDATES);
for ( i=0; i<cl.n_cand; i++ ) {
- image->candidate_cells[i] = cl.cand[i].cell;
+ if ( check_cell(p, image, cl.cand[i].cell) ) break;
}
free(cl.cand);
@@ -1016,7 +1046,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
struct reax_search *s;
int i;
- assert(pp->indm == INDEXING_REAX);
p = (struct reax_private *)pp;
fft_in = fftw_malloc(p->nel*sizeof(double));
@@ -1039,7 +1068,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
// fft_in, fft_out, p->plan, smin, smax,
// image->det, p);
- assemble_cells_from_candidates(image, s, cell);
+ assemble_cells_from_candidates(image, s, cell, p);
for ( i=0; i<s->n_search; i++ ) {
free(s->search[i].cand);
@@ -1051,7 +1080,9 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
-IndexingPrivate *reax_prepare()
+IndexingPrivate *reax_prepare(IndexingMethod indm, UnitCell *cell,
+ const char *filename, struct detector *det,
+ struct beam_params *beam, float *ltl)
{
struct reax_private *p;
int samp;
@@ -1060,8 +1091,6 @@ IndexingPrivate *reax_prepare()
p = calloc(1, sizeof(*p));
if ( p == NULL ) return NULL;
- p->base.indm = INDEXING_REAX;
-
p->angular_inc = deg2rad(1.0);
/* Reserve memory, over-estimating the number of directions */
@@ -1122,6 +1151,8 @@ IndexingPrivate *reax_prepare()
p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out,
1, FFTW_MEASURE);
+ p->indm = indm;
+
return (IndexingPrivate *)p;
}
@@ -1130,7 +1161,6 @@ void reax_cleanup(IndexingPrivate *pp)
{
struct reax_private *p;
- assert(pp->indm == INDEXING_REAX);
p = (struct reax_private *)pp;
free(p->directions);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 943ee656..31fb8527 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -69,30 +69,6 @@ struct _stream
};
-int count_patterns(FILE *fh)
-{
- char *rval;
-
- int n_total_patterns = 0;
- do {
- char line[1024];
-
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
- chomp(line);
- if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++;
-
- } while ( rval != NULL );
-
- if ( ferror(fh) ) {
- ERROR("Read error while counting patterns.\n");
- return 0;
- }
-
- return n_total_patterns;
-}
-
-
static int read_peaks(FILE *fh, struct image *image)
{
char *rval = NULL;
@@ -482,28 +458,6 @@ void write_stream_header(FILE *ofh, int argc, char *argv[])
}
-int skip_some_files(FILE *fh, int n)
-{
- char *rval = NULL;
- int n_patterns = 0;
-
- do {
-
- char line[1024];
-
- if ( n_patterns == n ) return 0;
-
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
- chomp(line);
- if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
-
- } while ( rval != NULL );
-
- return 1;
-}
-
-
Stream *open_stream_for_read(const char *filename)
{
Stream *st;
@@ -591,3 +545,21 @@ int is_stream(const char *filename)
return 0;
}
+
+
+/**
+ * rewind_stream:
+ * @st: A %Stream
+ *
+ * Attempts to set the file pointer for @st to the start of the stream, so that
+ * later calls to read_chunk() will repeat the sequence of chunks from the
+ * start.
+ *
+ * Programs must not assume that this operation always succeeds!
+ *
+ * Returns: non-zero if the stream could not be rewound.
+ */
+int rewind_stream(Stream *st)
+{
+ return fseek(st->fh, 0, SEEK_SET);
+}
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index db74583d..9340ca0d 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -49,9 +49,8 @@ extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile,
extern int read_chunk(Stream *st, struct image *image);
-/* Nasty functions that should be avoided */
-extern int count_patterns(FILE *fh);
-extern int skip_some_files(FILE *fh, int n);
+extern int rewind_stream(Stream *st);
+
extern int is_stream(const char *filename);
#endif /* STREAM_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 8bfbea2b..e9a8daa5 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -3,12 +3,12 @@
*
* Assemble and process FEL Bragg intensities
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2013 Thomas White <taw@physics.org>
* 2011 Andrew Martin <andrew.martin@desy.de>
* 2012 Lorenzo Galli <lorenzo.galli@desy.de>
*
@@ -48,6 +48,8 @@
#include "stream.h"
#include "reflist.h"
#include "image.h"
+#include "crystal.h"
+#include "thread-pool.h"
static void show_help(const char *s)
@@ -170,11 +172,11 @@ static double scale_intensities(RefList *reference, RefList *new,
}
-static int merge_pattern(RefList *model, struct image *new, RefList *reference,
- const SymOpList *sym,
- double *hist_vals, signed int hist_h,
- signed int hist_k, signed int hist_l, int *hist_n,
- int config_nopolar)
+static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
+ RefList *reference, const SymOpList *sym,
+ double *hist_vals, signed int hist_h,
+ signed int hist_k, signed int hist_l, int *hist_n,
+ int config_nopolar)
{
Reflection *refl;
RefListIterator *iter;
@@ -184,17 +186,18 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
double csx, csy, csz;
if ( reference != NULL ) {
- scale = scale_intensities(reference, new->reflections, sym);
+ scale = scale_intensities(reference,
+ crystal_get_reflections(cr), sym);
} else {
scale = 1.0;
}
if ( isnan(scale) ) return 1;
- cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
+ cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
- for ( refl = first_refl(new->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -226,7 +229,7 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- ool = 1.0 / new->lambda;
+ ool = 1.0 / image->lambda;
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool);
phi = atan2(yl, xl);
pa = pow(sin(phi)*sin(tt), 2.0);
@@ -260,63 +263,72 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
}
-static void merge_all(FILE *fh, RefList *model, RefList *reference,
- int config_startafter, int config_stopafter,
- const SymOpList *sym,
- int n_total_patterns,
- double *hist_vals, signed int hist_h,
- signed int hist_k, signed int hist_l,
- int *hist_i, int config_nopolar, int min_measurements)
+static void display_progress(int n_images, int n_crystals, int n_crystals_used)
+{
+ if ( !isatty(STDERR_FILENO) ) return;
+ if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
+
+ pthread_mutex_lock(&stderr_lock);
+ fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used",
+ n_images, n_crystals, n_crystals_used);
+ pthread_mutex_unlock(&stderr_lock);
+
+ fflush(stdout);
+}
+
+
+
+static int merge_all(Stream *st, RefList *model, RefList *reference,
+ const SymOpList *sym,
+ double *hist_vals, signed int hist_h,
+ signed int hist_k, signed int hist_l,
+ int *hist_i, int config_nopolar, int min_measurements)
{
int rval;
- int n_patterns = 0;
- int n_used = 0;
+ int n_images = 0;
+ int n_crystals = 0;
+ int n_crystals_used = 0;
Reflection *refl;
RefListIterator *iter;
- if ( skip_some_files(fh, config_startafter) ) {
- ERROR("Failed to skip first %i files.\n", config_startafter);
- return;
- }
-
do {
struct image image;
+ int i;
image.det = NULL;
/* Get data from next chunk */
- rval = read_chunk(fh, &image);
+ rval = read_chunk(st, &image);
if ( rval ) break;
- n_patterns++;
+ n_images++;
- if ( (image.reflections != NULL) && (image.indexed_cell) ) {
+ for ( i=0; i<image.n_crystals; i++ ) {
int r;
+ Crystal *cr = image.crystals[i];
- r = merge_pattern(model, &image, reference, sym,
+ n_crystals++;
+ r = merge_crystal(model, &image, cr, reference, sym,
hist_vals, hist_h, hist_k, hist_l,
hist_i, config_nopolar);
- if ( r == 0 ) n_used++;
+ if ( r == 0 ) n_crystals_used++;
+
+ crystal_free(cr);
}
free(image.filename);
- reflist_free(image.reflections);
image_feature_list_free(image.features);
- cell_free(image.indexed_cell);
-
- progress_bar(n_patterns, n_total_patterns-config_startafter,
- "Merging");
- if ( config_stopafter ) {
- if ( n_patterns == config_stopafter ) break;
- }
+ display_progress(n_images, n_crystals, n_crystals_used);
} while ( rval == 0 );
+ if ( rval ) return 1;
+
for ( refl = first_refl(model, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
@@ -339,7 +351,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference,
}
- STATUS("%i of the patterns could be used.\n", n_used);
+ return 0;
}
@@ -348,14 +360,11 @@ int main(int argc, char *argv[])
int c;
char *filename = NULL;
char *output = NULL;
- FILE *fh;
+ Stream *st;
RefList *model;
int config_maxonly = 0;
- int config_startafter = 0;
- int config_stopafter = 0;
int config_sum = 0;
int config_scale = 0;
- unsigned int n_total_patterns;
char *sym_str = NULL;
SymOpList *sym;
char *histo = NULL;
@@ -369,6 +378,7 @@ int main(int argc, char *argv[])
int config_nopolar = 0;
char *rval;
int min_measurements = 2;
+ int r;
/* Long options */
const struct option longopts[] = {
@@ -409,11 +419,11 @@ int main(int argc, char *argv[])
break;
case 's' :
- config_stopafter = atoi(optarg);
+ ERROR("The option '--stop-after' no longer works.\n");
break;
case 'f' :
- config_startafter = atoi(optarg);
+ ERROR("The option '--start-after' no longer works.\n");
break;
case 'y' :
@@ -465,25 +475,11 @@ int main(int argc, char *argv[])
free(sym_str);
/* Open the data stream */
- if ( strcmp(filename, "-") == 0 ) {
- fh = stdin;
- } else {
- fh = fopen(filename, "r");
- }
- free(filename);
- if ( fh == NULL ) {
- ERROR("Failed to open input file\n");
- return 1;
- }
-
- /* Count the number of patterns in the file */
- n_total_patterns = count_patterns(fh);
- if ( n_total_patterns == 0 ) {
- ERROR("No patterns to process.\n");
+ st = open_stream_for_read(filename);
+ if ( st == NULL ) {
+ ERROR("Failed to open stream.\n");
return 1;
}
- STATUS("There are %i patterns to process\n", n_total_patterns);
- rewind(fh);
model = reflist_new();
@@ -497,7 +493,8 @@ int main(int argc, char *argv[])
return 1;
}
- space_for_hist = n_total_patterns * num_equivs(sym, NULL);
+ /* FIXME: This array must grow as necessary. */
+ space_for_hist = 0 * num_equivs(sym, NULL);
hist_vals = malloc(space_for_hist * sizeof(double));
free(histo);
STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
@@ -529,40 +526,46 @@ int main(int argc, char *argv[])
}
hist_i = 0;
- merge_all(fh, model, NULL, config_startafter, config_stopafter,
- sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l,
- &hist_i, config_nopolar, min_measurements);
- if ( ferror(fh) ) {
- ERROR("Stream read error.\n");
+ r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l,
+ &hist_i, config_nopolar, min_measurements);
+ fprintf(stderr, "\n");
+ if ( r ) {
+ ERROR("Error while reading stream.\n");
return 1;
}
- rewind(fh);
if ( config_scale ) {
RefList *reference;
- STATUS("Extra pass for scaling...\n");
+ if ( rewind_stream(st) ) {
- reference = copy_reflist(model);
+ ERROR("Couldn't rewind stream - scaling cannot be "
+ "performed.\n");
- reflist_free(model);
- model = reflist_new();
+ } else {
- rewind(fh);
+ int r;
- merge_all(fh, model, reference,
- config_startafter, config_stopafter, sym,
- n_total_patterns,
- hist_vals, hist_h, hist_k, hist_l, &hist_i,
- config_nopolar, min_measurements);
+ STATUS("Extra pass for scaling...\n");
- if ( ferror(fh) ) {
- ERROR("Stream read error.\n");
- return 1;
- }
+ reference = copy_reflist(model);
+
+ reflist_free(model);
+ model = reflist_new();
+
+ r = merge_all(st, model, reference, sym,
+ hist_vals, hist_h, hist_k, hist_l, &hist_i,
+ config_nopolar, min_measurements);
+ fprintf(stderr, "\n");
+ if ( r ) {
+ ERROR("Error while reading stream.\n");
+ return 1;
+ }
+
+ reflist_free(reference);
- reflist_free(reference);
+ }
}
@@ -579,7 +582,7 @@ int main(int argc, char *argv[])
write_reflist(output, model);
- fclose(fh);
+ close_stream(st);
free(sym);
reflist_free(model);