Remove old reflections.h

21 files changed, 979 insertions, 925 deletions

diff --git a/Makefile.am b/Makefile.am
index 50631ced..2173bfe4 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -26,7 +26,7 @@ LDADD = $(top_builddir)/lib/libgnu.a @IGNORE_UNUSED_LIBRARIES_CFLAGS@
 src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c src/utils.c \
                           src/image.c src/cell.c src/hdf5-file.c \
                           src/detector.c src/peaks.c \
-                          src/reflections.c src/beam-parameters.c \
+                          src/reflist-utils.c src/beam-parameters.c \
                           src/symmetry.c src/thread-pool.c src/reflist.c
 if HAVE_OPENCL
 src_pattern_sim_SOURCES += src/diffraction-gpu.c src/cl-utils.c
@@ -51,10 +51,10 @@ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
                           src/cell.c src/image.c src/peaks.c src/index.c \
                           src/filters.c src/diffraction.c src/detector.c \
                           src/dirax.c src/mosflm.c \
-                          src/reflections.c src/symmetry.c \
+                          src/reflist-utils.c src/symmetry.c \
                           src/geometry.c src/thread-pool.c \
-                          src/beam-parameters.c src/reflist.c src/stream.c \
-                          src/reflist-utils.c
+                          src/beam-parameters.c src/reflist.c src/stream.c
+
 if HAVE_OPENCL
 src_indexamajig_SOURCES += src/diffraction-gpu.c src/cl-utils.c
 endif
@@ -67,25 +67,25 @@ src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \
 endif
 
 src_get_hkl_SOURCES = src/get_hkl.c src/cell.c src/utils.c \
-                      src/reflections.c src/symmetry.c src/beam-parameters.c \
-                      src/thread-pool.c
+                      src/reflist-utils.c src/symmetry.c src/beam-parameters.c \
+                      src/thread-pool.c src/reflist.c
 
 src_compare_hkl_SOURCES = src/compare_hkl.c src/cell.c src/utils.c \
-                          src/reflections.c src/statistics.c src/symmetry.c \
-                          src/thread-pool.c
+                          src/reflist-utils.c src/statistics.c src/symmetry.c \
+                          src/thread-pool.c src/reflist.c
 
 src_check_hkl_SOURCES = src/check_hkl.c src/cell.c src/utils.c \
-                        src/reflections.c src/statistics.c src/symmetry.c \
-                        src/thread-pool.c
+                        src/reflist-utils.c src/statistics.c src/symmetry.c \
+                        src/thread-pool.c src/reflist.c
 
 src_powder_plot_SOURCES = src/powder_plot.c src/cell.c src/utils.c src/image.c \
                           src/hdf5-file.c src/detector.c src/thread-pool.c \
                           src/beam-parameters.c
 
-src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \
-                         src/utils.c  src/povray.c src/symmetry.c src/render.c \
+src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflist-utils.c \
+                         src/utils.c src/povray.c src/symmetry.c src/render.c \
                          src/hdf5-file.c src/image.c src/filters.c \
-                         src/thread-pool.c
+                         src/thread-pool.c src/reflist.c
 
 src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \
                         src/image.c src/filters.c src/peaks.c src/detector.c \
@@ -97,10 +97,10 @@ src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \
 
 src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
                           src/utils.c src/detector.c src/peaks.c src/image.c \
-                          src/geometry.c src/reflections.c src/stream.c \
+                          src/geometry.c src/reflist-utils.c src/stream.c \
                           src/thread-pool.c src/beam-parameters.c \
                           src/symmetry.c src/post-refinement.c \
-                          src/hrs-scaling.c src/reflist.c src/reflist-utils.c
+                          src/hrs-scaling.c src/reflist.c
 
 if BUILD_CUBEIT
 src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \
@@ -115,7 +115,7 @@ tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \
 INCLUDES = "-I$(top_srcdir)/data"
 
 EXTRA_DIST += src/cell.h src/hdf5-file.h src/image.h src/utils.h \
-              src/diffraction.h src/detector.h src/reflections.h \
+              src/diffraction.h src/detector.h src/reflist-utils.h \
               src/list_tmp.h src/statistics.h src/dw-hdfsee.h \
               src/render.h src/hdfsee.h src/dirax.h src/mosflm.h src/peaks.h \
               src/index.h src/filters.h src/diffraction-gpu.h src/cl-utils.h \
diff --git a/Makefile.in b/Makefile.in
index f8d4bf23..226389f9 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -90,16 +90,17 @@ src_calibrate_detector_OBJECTS = $(am_src_calibrate_detector_OBJECTS)
 src_calibrate_detector_LDADD = $(LDADD)
 src_calibrate_detector_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_src_check_hkl_OBJECTS = src/check_hkl.$(OBJEXT) src/cell.$(OBJEXT) \
-	src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \
+	src/utils.$(OBJEXT) src/reflist-utils.$(OBJEXT) \
 	src/statistics.$(OBJEXT) src/symmetry.$(OBJEXT) \
-	src/thread-pool.$(OBJEXT)
+	src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT)
 src_check_hkl_OBJECTS = $(am_src_check_hkl_OBJECTS)
 src_check_hkl_LDADD = $(LDADD)
 src_check_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_src_compare_hkl_OBJECTS = src/compare_hkl.$(OBJEXT) \
 	src/cell.$(OBJEXT) src/utils.$(OBJEXT) \
-	src/reflections.$(OBJEXT) src/statistics.$(OBJEXT) \
-	src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT)
+	src/reflist-utils.$(OBJEXT) src/statistics.$(OBJEXT) \
+	src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT) \
+	src/reflist.$(OBJEXT)
 src_compare_hkl_OBJECTS = $(am_src_compare_hkl_OBJECTS)
 src_compare_hkl_LDADD = $(LDADD)
 src_compare_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -120,9 +121,9 @@ src_cubeit_OBJECTS = $(am_src_cubeit_OBJECTS)
 src_cubeit_LDADD = $(LDADD)
 src_cubeit_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_src_get_hkl_OBJECTS = src/get_hkl.$(OBJEXT) src/cell.$(OBJEXT) \
-	src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \
+	src/utils.$(OBJEXT) src/reflist-utils.$(OBJEXT) \
 	src/symmetry.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
-	src/thread-pool.$(OBJEXT)
+	src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT)
 src_get_hkl_OBJECTS = $(am_src_get_hkl_OBJECTS)
 src_get_hkl_LDADD = $(LDADD)
 src_get_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -142,10 +143,9 @@ src_hdfsee_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am__src_indexamajig_SOURCES_DIST = src/indexamajig.c src/hdf5-file.c \
 	src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \
 	src/filters.c src/diffraction.c src/detector.c src/dirax.c \
-	src/mosflm.c src/reflections.c src/symmetry.c src/geometry.c \
+	src/mosflm.c src/reflist-utils.c src/symmetry.c src/geometry.c \
 	src/thread-pool.c src/beam-parameters.c src/reflist.c \
-	src/stream.c src/reflist-utils.c src/diffraction-gpu.c \
-	src/cl-utils.c
+	src/stream.c src/diffraction-gpu.c src/cl-utils.c
 @HAVE_OPENCL_TRUE@am__objects_1 = src/diffraction-gpu.$(OBJEXT) \
 @HAVE_OPENCL_TRUE@	src/cl-utils.$(OBJEXT)
 am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
@@ -153,35 +153,34 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
 	src/image.$(OBJEXT) src/peaks.$(OBJEXT) src/index.$(OBJEXT) \
 	src/filters.$(OBJEXT) src/diffraction.$(OBJEXT) \
 	src/detector.$(OBJEXT) src/dirax.$(OBJEXT) \
-	src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \
+	src/mosflm.$(OBJEXT) src/reflist-utils.$(OBJEXT) \
 	src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \
 	src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
-	src/reflist.$(OBJEXT) src/stream.$(OBJEXT) \
-	src/reflist-utils.$(OBJEXT) $(am__objects_1)
+	src/reflist.$(OBJEXT) src/stream.$(OBJEXT) $(am__objects_1)
 src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS)
 src_indexamajig_LDADD = $(LDADD)
 src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_src_partialator_OBJECTS = src/partialator.$(OBJEXT) \
 	src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) src/utils.$(OBJEXT) \
 	src/detector.$(OBJEXT) src/peaks.$(OBJEXT) src/image.$(OBJEXT) \
-	src/geometry.$(OBJEXT) src/reflections.$(OBJEXT) \
+	src/geometry.$(OBJEXT) src/reflist-utils.$(OBJEXT) \
 	src/stream.$(OBJEXT) src/thread-pool.$(OBJEXT) \
 	src/beam-parameters.$(OBJEXT) src/symmetry.$(OBJEXT) \
 	src/post-refinement.$(OBJEXT) src/hrs-scaling.$(OBJEXT) \
-	src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT)
+	src/reflist.$(OBJEXT)
 src_partialator_OBJECTS = $(am_src_partialator_OBJECTS)
 src_partialator_LDADD = $(LDADD)
 src_partialator_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am__src_pattern_sim_SOURCES_DIST = src/pattern_sim.c src/diffraction.c \
 	src/utils.c src/image.c src/cell.c src/hdf5-file.c \
-	src/detector.c src/peaks.c src/reflections.c \
+	src/detector.c src/peaks.c src/reflist-utils.c \
 	src/beam-parameters.c src/symmetry.c src/thread-pool.c \
 	src/reflist.c src/diffraction-gpu.c src/cl-utils.c
 am_src_pattern_sim_OBJECTS = src/pattern_sim.$(OBJEXT) \
 	src/diffraction.$(OBJEXT) src/utils.$(OBJEXT) \
 	src/image.$(OBJEXT) src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
 	src/detector.$(OBJEXT) src/peaks.$(OBJEXT) \
-	src/reflections.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
+	src/reflist-utils.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
 	src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT) \
 	src/reflist.$(OBJEXT) $(am__objects_1)
 src_pattern_sim_OBJECTS = $(am_src_pattern_sim_OBJECTS)
@@ -205,11 +204,12 @@ src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS)
 src_process_hkl_LDADD = $(LDADD)
 src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_src_render_hkl_OBJECTS = src/render_hkl.$(OBJEXT) \
-	src/cell.$(OBJEXT) src/reflections.$(OBJEXT) \
+	src/cell.$(OBJEXT) src/reflist-utils.$(OBJEXT) \
 	src/utils.$(OBJEXT) src/povray.$(OBJEXT) \
 	src/symmetry.$(OBJEXT) src/render.$(OBJEXT) \
 	src/hdf5-file.$(OBJEXT) src/image.$(OBJEXT) \
-	src/filters.$(OBJEXT) src/thread-pool.$(OBJEXT)
+	src/filters.$(OBJEXT) src/thread-pool.$(OBJEXT) \
+	src/reflist.$(OBJEXT)
 src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS)
 src_render_hkl_LDADD = $(LDADD)
 src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -558,7 +558,7 @@ top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
 EXTRA_DIST = configure m4/gnulib-cache.m4 src/cell.h src/hdf5-file.h \
 	src/image.h src/utils.h src/diffraction.h src/detector.h \
-	src/reflections.h src/list_tmp.h src/statistics.h \
+	src/reflist-utils.h src/list_tmp.h src/statistics.h \
 	src/dw-hdfsee.h src/render.h src/hdfsee.h src/dirax.h \
 	src/mosflm.h src/peaks.h src/index.h src/filters.h \
 	src/diffraction-gpu.h src/cl-utils.h src/symmetry.h \
@@ -588,7 +588,7 @@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -I$(top_srcdir)/l
 LDADD = $(top_builddir)/lib/libgnu.a @IGNORE_UNUSED_LIBRARIES_CFLAGS@
 src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \
 	src/utils.c src/image.c src/cell.c src/hdf5-file.c \
-	src/detector.c src/peaks.c src/reflections.c \
+	src/detector.c src/peaks.c src/reflist-utils.c \
 	src/beam-parameters.c src/symmetry.c src/thread-pool.c \
 	src/reflist.c $(am__append_3)
 @HAVE_OPENCL_TRUE@tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src/utils.c \
@@ -606,33 +606,33 @@ src_process_hkl_SOURCES = src/process_hkl.c src/statistics.c \
 src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c \
 	src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \
 	src/filters.c src/diffraction.c src/detector.c src/dirax.c \
-	src/mosflm.c src/reflections.c src/symmetry.c src/geometry.c \
+	src/mosflm.c src/reflist-utils.c src/symmetry.c src/geometry.c \
 	src/thread-pool.c src/beam-parameters.c src/reflist.c \
-	src/stream.c src/reflist-utils.c $(am__append_6)
+	src/stream.c $(am__append_6)
 @BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \
 @BUILD_HDFSEE_TRUE@                     src/hdf5-file.c src/utils.c src/image.c src/filters.c \
 @BUILD_HDFSEE_TRUE@                     src/thread-pool.c src/detector.c
 
 src_get_hkl_SOURCES = src/get_hkl.c src/cell.c src/utils.c \
-                      src/reflections.c src/symmetry.c src/beam-parameters.c \
-                      src/thread-pool.c
+                      src/reflist-utils.c src/symmetry.c src/beam-parameters.c \
+                      src/thread-pool.c src/reflist.c
 
 src_compare_hkl_SOURCES = src/compare_hkl.c src/cell.c src/utils.c \
-                          src/reflections.c src/statistics.c src/symmetry.c \
-                          src/thread-pool.c
+                          src/reflist-utils.c src/statistics.c src/symmetry.c \
+                          src/thread-pool.c src/reflist.c
 
 src_check_hkl_SOURCES = src/check_hkl.c src/cell.c src/utils.c \
-                        src/reflections.c src/statistics.c src/symmetry.c \
-                        src/thread-pool.c
+                        src/reflist-utils.c src/statistics.c src/symmetry.c \
+                        src/thread-pool.c src/reflist.c
 
 src_powder_plot_SOURCES = src/powder_plot.c src/cell.c src/utils.c src/image.c \
                           src/hdf5-file.c src/detector.c src/thread-pool.c \
                           src/beam-parameters.c
 
-src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \
-                         src/utils.c  src/povray.c src/symmetry.c src/render.c \
+src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflist-utils.c \
+                         src/utils.c src/povray.c src/symmetry.c src/render.c \
                          src/hdf5-file.c src/image.c src/filters.c \
-                         src/thread-pool.c
+                         src/thread-pool.c src/reflist.c
 
 src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \
                         src/image.c src/filters.c src/peaks.c src/detector.c \
@@ -644,10 +644,10 @@ src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \
 
 src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
                           src/utils.c src/detector.c src/peaks.c src/image.c \
-                          src/geometry.c src/reflections.c src/stream.c \
+                          src/geometry.c src/reflist-utils.c src/stream.c \
                           src/thread-pool.c src/beam-parameters.c \
                           src/symmetry.c src/post-refinement.c \
-                          src/hrs-scaling.c src/reflist.c src/reflist-utils.c
+                          src/hrs-scaling.c src/reflist.c
 
 @BUILD_CUBEIT_TRUE@src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \
 @BUILD_CUBEIT_TRUE@                     src/detector.c src/render.c src/filters.c src/image.c \
@@ -803,12 +803,14 @@ src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibra
 src/check_hkl.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/cell.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
-src/reflections.$(OBJEXT): src/$(am__dirstamp) \
+src/reflist-utils.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/statistics.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/symmetry.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
+src/reflist.$(OBJEXT): src/$(am__dirstamp) \
+	src/$(DEPDIR)/$(am__dirstamp)
 src/check_hkl$(EXEEXT): $(src_check_hkl_OBJECTS) $(src_check_hkl_DEPENDENCIES) src/$(am__dirstamp)
 	@rm -f src/check_hkl$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(src_check_hkl_OBJECTS) $(src_check_hkl_LDADD) $(LIBS)
@@ -825,10 +827,6 @@ src/filters.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/stream.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
-src/reflist.$(OBJEXT): src/$(am__dirstamp) \
-	src/$(DEPDIR)/$(am__dirstamp)
-src/reflist-utils.$(OBJEXT): src/$(am__dirstamp) \
-	src/$(DEPDIR)/$(am__dirstamp)
 src/cubeit$(EXEEXT): $(src_cubeit_OBJECTS) $(src_cubeit_DEPENDENCIES) src/$(am__dirstamp)
 	@rm -f src/cubeit$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(src_cubeit_OBJECTS) $(src_cubeit_LDADD) $(LIBS)
@@ -948,7 +946,6 @@ mostlyclean-compile:
 	-rm -f src/povray.$(OBJEXT)
 	-rm -f src/powder_plot.$(OBJEXT)
 	-rm -f src/process_hkl.$(OBJEXT)
-	-rm -f src/reflections.$(OBJEXT)
 	-rm -f src/reflist-utils.$(OBJEXT)
 	-rm -f src/reflist.$(OBJEXT)
 	-rm -f src/render.$(OBJEXT)
@@ -994,7 +991,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/povray.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/powder_plot.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/process_hkl.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflections.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist-utils.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render.Po@am__quote@
diff --git a/src/check_hkl.c b/src/check_hkl.c
index fb6d02d5..903188a1 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -22,9 +22,10 @@
 #include <getopt.h>
 
 #include "utils.h"
-#include "reflections.h"
 #include "statistics.h"
 #include "symmetry.h"
+#include "reflist.h"
+#include "reflist-utils.h"
 
 
 /* Number of bins for plot of resolution shells */
@@ -46,9 +47,8 @@ static void show_help(const char *s)
 }
 
 
-static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
-                        const char *sym, unsigned int *counts,
-                        const double *sigma, double rmin_fix, double rmax_fix)
+static void plot_shells(RefList *list, UnitCell *cell, const char *sym,
+                        double rmin_fix, double rmax_fix)
 {
 	double num[NBINS];
 	int cts[NBINS];
@@ -70,6 +70,8 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 	int nmeas = 0;
 	int nmeastot = 0;
 	int nout = 0;
+	Reflection *refl;
+	RefListIterator *iter;
 
 	if ( cell == NULL ) {
 		ERROR("Need the unit cell to plot resolution shells.\n");
@@ -93,16 +95,17 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 		mean[i] = 0;
 	}
 
-	rmin = +INFINITY;
-	rmax = 0.0;
-	for ( i=0; i<num_items(items); i++ ) {
+	/* Iterate over all common reflections and calculate min and max
+	 * resolution */
+	rmin = +INFINITY;  rmax = 0.0;
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
 		double d;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl, &h, &k, &l);
 
 		d = resolution(cell, h, k, l) * 2.0;
 		if ( d > rmax ) rmax = d;
@@ -179,16 +182,16 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 	free(counted);
 
 	/* Calculate means */
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
 		double d;
 		int bin;
 		int j;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl, &h, &k, &l);
 
 		d = resolution(cell, h, k, l) * 2.0;
 
@@ -205,7 +208,7 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 		}
 
 		measured[bin]++;
-		mean[bin] += lookup_intensity(ref, h, k, l);
+		mean[bin] += get_intensity(refl);
 
 	}
 
@@ -214,21 +217,21 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 	}
 
 	/* Characterise the data set */
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
 		double d;
 		int bin;
 		int j;
 		double val, esd;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl, &h, &k, &l);
 
 		d = resolution(cell, h, k, l) * 2.0;
-		val = lookup_intensity(ref, h, k, l);
-		esd = lookup_sigma(sigma, h, k, l);
+		val = get_intensity(refl);
+		esd = get_esd_intensity(refl);
 
 		bin = -1;
 		for ( j=0; j<NBINS; j++ ) {
@@ -245,11 +248,11 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 		if ( !isfinite(val/esd) ) continue;
 
 		/* measured[bin] was done earlier */
-		measurements[bin] += lookup_count(counts, h, k, l);
+		measurements[bin] += get_redundancy(refl);
 		snr[bin] += val / esd;
 		snr_total += val / esd;
 		nmeas++;
-		nmeastot += lookup_count(counts, h, k, l);
+		nmeastot += get_redundancy(refl);
 
 		var[bin] += pow(val-mean[bin], 2.0);
 
@@ -289,17 +292,15 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
 int main(int argc, char *argv[])
 {
 	int c;
-	double *ref;
 	UnitCell *cell;
 	char *file = NULL;
 	char *sym = NULL;
-	int i;
-	ReflItemList *items;
-	ReflItemList *good_items;
+	RefList *raw_list;
+	RefList *list;
+	Reflection *refl;
+	RefListIterator *iter;
 	char *pdb = NULL;
-	double *esd;
 	int rej = 0;
-	unsigned int *cts;
 	float rmin_fix = -1.0;
 	float rmax_fix = -1.0;
 
@@ -369,37 +370,35 @@ int main(int argc, char *argv[])
 	cell = load_cell_from_pdb(pdb);
 	free(pdb);
 
-	ref = new_list_intensity();
-	esd = new_list_sigma();
-	cts = new_list_count();
-	items = read_reflections(file, ref, NULL, cts, esd);
-	if ( items == NULL ) {
-		ERROR("Couldn't open file '%s'\n", file);
+	raw_list = read_reflections(file);
+	if ( raw_list == NULL ) {
+		ERROR("Couldn't read file '%s'\n", file);
 		return 1;
 	}
 
 	/* Check that the intensities have the correct symmetry */
-	if ( check_symmetry(items, sym) ) {
+	if ( check_list_symmetry(raw_list, sym) ) {
 		ERROR("The input reflection list does not appear to"
 		      " have symmetry %s\n", sym);
 		return 1;
 	}
 
-	/* Reject reflections */
-	good_items = new_items();
-	for ( i=0; i<num_items(items); i++ ) {
+	/* Reject some reflections */
+	list = reflist_new();
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
 		double val, sig;
 		int ig = 0;
 		double d;
+		Reflection *new;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl, &h, &k, &l);
 
-		val = lookup_intensity(ref, h, k, l);
-		sig = lookup_sigma(esd, h, k, l);
+		val = get_intensity(refl);
+		sig = get_esd_intensity(refl);
 
 		if ( val < 3.0 * sig ) {
 			rej++;
@@ -412,11 +411,12 @@ int main(int argc, char *argv[])
 
 		//if ( ig ) continue;
 
-		add_item(good_items, h, k, l);
+		new = add_refl(list, h, k, l);
+		copy_data(new, refl);
 
 	}
 
-	plot_shells(ref, items, cell, sym, cts, esd, rmin_fix, rmax_fix);
+	plot_shells(list, cell, sym, rmin_fix, rmax_fix);
 
 	return 0;
 }
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index abe8bc6a..16a744ea 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -20,11 +20,12 @@
 #include <string.h>
 #include <unistd.h>
 #include <getopt.h>
+#include <assert.h>
 
 #include "utils.h"
-#include "reflections.h"
 #include "statistics.h"
 #include "symmetry.h"
+#include "reflist-utils.h"
 
 
 /* Number of bins for plot of resolution shells */
@@ -38,7 +39,7 @@ static void show_help(const char *s)
 "Compare intensity lists.\n"
 "\n"
 "  -h, --help                 Display this help message.\n"
-"  -o, --output=<filename>    Specify output filename for correction factor.\n"
+"  -o, --ratio=<filename>     Specify output filename for ratios.\n"
 "  -y, --symmetry=<sym>       The symmetry of both the input files.\n"
 "  -p, --pdb=<filename>       PDB file to use (default: molecule.pdb).\n"
 "      --shells               Plot the figures of merit by resolution.\n"
@@ -48,9 +49,9 @@ static void show_help(const char *s)
 }
 
 
-static void plot_shells(const double *ref1, const double *ref2,
-                        ReflItemList *items, double scale, UnitCell *cell,
-                        const char *sym, double rmin_fix, double rmax_fix)
+static void plot_shells(RefList *list1, RefList *list2, double scale,
+                        UnitCell *cell, const char *sym,
+                        double rmin_fix, double rmax_fix)
 {
 	double num[NBINS];
 	int cts[NBINS];
@@ -69,6 +70,12 @@ static void plot_shells(const double *ref1, const double *ref2,
 	int ctot;
 	FILE *fh;
 	int nout = 0;
+	Reflection *refl1;
+	RefListIterator *iter;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	int hmax, kmax, lmax;
 
 	if ( cell == NULL ) {
 		ERROR("Need the unit cell to plot resolution shells.\n");
@@ -90,16 +97,20 @@ static void plot_shells(const double *ref1, const double *ref2,
 		snr[i] = 0;
 	}
 
-	rmin = +INFINITY;
-	rmax = 0.0;
-	for ( i=0; i<num_items(items); i++ ) {
+	/* Iterate over all common reflections and calculate min and max
+	 * resolution */
+	rmin = +INFINITY;  rmax = 0.0;
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
 		double d;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
 		d = resolution(cell, h, k, l) * 2.0;
 		if ( d > rmax ) rmax = d;
@@ -117,6 +128,7 @@ static void plot_shells(const double *ref1, const double *ref2,
 	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
 	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
 
+	/* Calculate the resolution bins */
 	total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
 	vol_per_shell = total_vol / NBINS;
 	rmins[0] = rmin;
@@ -145,9 +157,16 @@ static void plot_shells(const double *ref1, const double *ref2,
 
 	/* Count the number of reflections possible in each shell */
 	counted = new_list_count();
-	for ( h=-50; h<=+50; h++ ) {
-	for ( k=-50; k<=+50; k++ ) {
-	for ( l=-50; l<=+50; l++ ) {
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	hmax = rmax / modulus(asx, asy, asz);
+	kmax = rmax / modulus(bsx, bsy, bsz);
+	lmax = rmax / modulus(csx, csy, csz);
+
+	for ( h=-hmax; h<hmax; h++ ) {
+	for ( k=-kmax; k<kmax; k++ ) {
+	for ( l=-lmax; l<lmax; l++ ) {
 
 		double d;
 		signed int hs, ks, ls;
@@ -178,17 +197,21 @@ static void plot_shells(const double *ref1, const double *ref2,
 	den = 0.0;
 	ctot = 0;
 	nout = 0;
-	for ( i=0; i<num_items(items); i++ ) {
 
-		struct refl_item *it;
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
+
 		signed int h, k, l;
 		double d;
 		int bin;
 		double i1, i2;
 		int j;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
 		d = resolution(cell, h, k, l) * 2.0;
 
@@ -206,8 +229,8 @@ static void plot_shells(const double *ref1, const double *ref2,
 			continue;
 		}
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = lookup_intensity(ref2, h, k, l);
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 		i2 *= scale;
 
 		num[bin] += fabs(i1 - i2);
@@ -238,33 +261,32 @@ static void plot_shells(const double *ref1, const double *ref2,
 int main(int argc, char *argv[])
 {
 	int c;
-	double *ref1;
-	double *ref2;
-	double *ref2_transformed;
-	double *out;
 	UnitCell *cell;
-	char *outfile = NULL;
+	char *ratiofile = NULL;
 	char *afile = NULL;
 	char *bfile = NULL;
 	char *sym = NULL;
 	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
 	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
-	int i, ncom;
-	ReflItemList *i1, *i2, *icommon;
+	int ncom;
+	RefList *list1;
+	RefList *list2;
+	RefList *list2_transformed;
+	RefList *ratio;
+	RefList *deleteme;
 	int config_shells = 0;
 	char *pdb = NULL;
-	double *esd1;
-	double *esd2;
 	int rej1 = 0;
 	int rej2 = 0;
-	unsigned int *cts1;
 	float rmin_fix = -1.0;
 	float rmax_fix = -1.0;
+	Reflection *refl1;
+	RefListIterator *iter;
 
 	/* Long options */
 	const struct option longopts[] = {
 		{"help",               0, NULL,               'h'},
-		{"output",             1, NULL,               'o'},
+		{"ratio" ,             1, NULL,               'o'},
 		{"symmetry",           1, NULL,               'y'},
 		{"shells",             0, &config_shells,     1},
 		{"pdb",                1, NULL,               'p'},
@@ -282,7 +304,7 @@ int main(int argc, char *argv[])
 			return 0;
 
 		case 'o' :
-			outfile = strdup(optarg);
+			ratiofile = strdup(optarg);
 			break;
 
 		case 'y' :
@@ -335,63 +357,62 @@ int main(int argc, char *argv[])
 	cell = load_cell_from_pdb(pdb);
 	free(pdb);
 
-	ref1 = new_list_intensity();
-	esd1 = new_list_sigma();
-	cts1 = new_list_count();
-	i1 = read_reflections(afile, ref1, NULL, cts1, esd1);
-	if ( i1 == NULL ) {
-		ERROR("Couldn't open file '%s'\n", afile);
+	list1 = read_reflections(afile);
+	if ( list1 == NULL ) {
+		ERROR("Couldn't read file '%s'\n", afile);
 		return 1;
 	}
-	ref2 = new_list_intensity();
-	esd2 = new_list_sigma();
-	i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
-	if ( i2 == NULL ) {
-		ERROR("Couldn't open file '%s'\n", bfile);
+
+	list2 = read_reflections(bfile);
+	if ( list2 == NULL ) {
+		ERROR("Couldn't read file '%s'\n", bfile);
 		return 1;
 	}
 
 	/* Check that the intensities have the correct symmetry */
-	if ( check_symmetry(i1, sym) ) {
+	if ( check_list_symmetry(list1, sym) ) {
 		ERROR("The first input reflection list does not appear to"
 		      " have symmetry %s\n", sym);
 		return 1;
 	}
-	if ( check_symmetry(i2, sym) ) {
+	if ( check_list_symmetry(list2, sym) ) {
 		ERROR("The second input reflection list does not appear to"
 		      " have symmetry %s\n", sym);
 		return 1;
 	}
 
-	/* List for output scale factor map */
-	out = new_list_intensity();
-
 	/* Find common reflections (taking symmetry into account) */
-	icommon = new_items();
-	ref2_transformed = new_list_intensity();
-	for ( i=0; i<num_items(i1); i++ ) {
+	list2_transformed = reflist_new();
+	ratio = reflist_new();
+	deleteme = reflist_new();
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
 		signed int he, ke, le;
 		double val1, val2;
 		double sig1, sig2;
 		int ig = 0;
 		double d;
+		Reflection *refl2;
+		Reflection *tr;
 
-		it = get_item(i1, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
 
-		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
-			//STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
-			//       1.0/(2.0*resolution(cell, h, k, l)/1e9));
+		if ( !find_equiv_in_list(list2, h, k, l, sym, &he, &ke, &le) ) {
+			/* No common reflection */
+			add_refl(deleteme, h, k, l);
 			continue;
 		}
 
-		val1 = lookup_intensity(ref1, h, k, l);
-		val2 = lookup_intensity(ref2, he, ke, le);
-		sig1 = lookup_sigma(esd1, h, k, l);
-		sig2 = lookup_sigma(esd2, he, ke, le);
+		refl2 = find_refl(list2, he, ke, le);
+		assert(refl2 != NULL);
+
+		val1 = get_intensity(refl1);
+		val2 = get_intensity(refl2);
+		sig1 = get_esd_intensity(refl1);
+		sig2 = get_esd_intensity(refl2);
 
 		if ( val1 < 3.0 * sig1 ) {
 			rej1++;
@@ -408,67 +429,90 @@ int main(int argc, char *argv[])
 
 		//if ( ig ) continue;
 
-		set_intensity(ref2_transformed, h, k, l, val2);
-		set_intensity(out, h, k, l, val1/val2);
-		add_item(icommon, h, k, l);
+		/* Add the old data from 'refl2' to a new list with the same
+		 * indices as its equivalent in 'list1' */
+		tr = add_refl(list2_transformed, h, k, l);
+		copy_data(tr, refl2);
+
+		/* Add divided version to 'output' list */
+		tr = add_refl(ratio, h, k, l);
+		set_int(tr, val1/val2);
 
 	}
-	ncom = num_items(icommon);
+
 	STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej1, afile);
 	STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej2, bfile);
 
+	ncom = num_reflections(list2_transformed);
 	STATUS("%i,%i reflections: %i in common\n",
-	       num_items(i1), num_items(i2), ncom);
+	       num_reflections(list1), num_reflections(list2), ncom);
+
+	/* Trim reflections from 'list1' which had no equivalents in 'list2' */
+	for ( refl1 = first_refl(deleteme, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-	R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale_r1fi);
-	STATUS("R1(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+		signed int h, k, l;
+		Reflection *del;
+
+		get_indices(refl1, &h, &k, &l);
+		del = find_refl(list1, h, k, l);
+		assert(del != NULL);
+
+		delete_refl(del);
+
+	}
+	reflist_free(deleteme);
+	reflist_free(list2);
+
+	R1 = stat_r1_ignore(list1, list2_transformed, &scale_r1fi);
+	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
+	       " (ignoring negative intensities)\n",
 	       R1*100.0, scale_r1fi);
 
-	R1 = stat_r1_zero(ref1, ref2_transformed, icommon, &scale_r1);
-	STATUS("R1(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+	R1 = stat_r1_zero(list1, list2_transformed, &scale_r1);
+	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
+	       " (zeroing negative intensities)\n",
 	       R1*100.0, scale_r1);
 
-	R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2);
+	R2 = stat_r2(list1, list2_transformed, &scale_r2);
 	STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
 
-	R1i = stat_r1_i(ref1, ref2_transformed, icommon, &scale_r1i);
+	R1i = stat_r1_i(list1, list2_transformed, &scale_r1i);
 	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
 
-	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
-	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+	Rdiff = stat_rdiff_ignore(list1, list2_transformed, &scale_rdig);
+	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
+	       " (ignoring negative intensities)\n",
 	       Rdiff*100.0, scale_rdig);
 
-	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
-	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+	Rdiff = stat_rdiff_zero(list1, list2_transformed, &scale);
+	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
+	       " (zeroing negative intensities)\n",
 	       Rdiff*100.0, scale);
 
-	Rdiff = stat_rdiff_intensity(ref1, ref2_transformed, icommon,
-	                             &scale_rintint);
+	Rdiff = stat_rdiff_intensity(list1, list2_transformed, &scale_rintint);
 	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
 	       Rdiff*100.0, scale_rintint);
 
-	pearson = stat_pearson_i(ref1, ref2_transformed, icommon);
+	pearson = stat_pearson_i(list1, list2_transformed);
 	STATUS("Pearson r(I) = %5.4f\n", pearson);
 
-	pearson = stat_pearson_f_ignore(ref1, ref2_transformed, icommon);
+	pearson = stat_pearson_f_ignore(list1, list2_transformed);
 	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
 	       pearson);
 
-	pearson = stat_pearson_f_zero(ref1, ref2_transformed, icommon);
+	pearson = stat_pearson_f_zero(list1, list2_transformed);
 	STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
 	       pearson);
 
 	if ( config_shells ) {
-		plot_shells(ref1, ref2_transformed, icommon, scale_r1fi,
+		plot_shells(list1, list2_transformed, scale_r1fi,
 		            cell, sym, rmin_fix, rmax_fix);
 	}
 
-	if ( outfile != NULL ) {
-
-		write_reflections(outfile, icommon, out, NULL,
-		                  NULL, NULL, cell);
-		STATUS("Sigma(I) values in output file are not meaningful.\n");
-
+	if ( ratiofile != NULL ) {
+		write_reflist(ratiofile, ratio, cell);
 	}
 
 	return 0;
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 1b7e4406..026198ab 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -22,7 +22,7 @@
 #include <getopt.h>
 
 #include "utils.h"
-#include "reflections.h"
+#include "reflist-utils.h"
 #include "symmetry.h"
 #include "beam-parameters.h"
 
@@ -45,7 +45,8 @@ static void show_help(const char *s)
 "To calculate Poisson samples accurately, you must also give:\n"
 "  -b, --beam=<file>          Get beam parameters from file.\n"
 "\n"
-"You can artificially 'twin' the reflections, or expand them out:\n"
+"You can artificially 'twin' the reflections, or expand them out.  You can also"
+" do both, in which case the 'twinning' will be done first:\n"
 "  -w, --twin=<sym>           Generate twinned data according to the given\n"
 "                              point group.\n"
 "  -e, --expand=<sym>         Expand reflections to this point group.\n"
@@ -59,99 +60,123 @@ static void show_help(const char *s)
 "\n"
 "Don't forget to specify the output filename:\n"
 "  -o, --output=<filename>    Output filename (default: stdout).\n"
+"  -p, --pdb=<filename>       Use unit cell parameters from this PDB file to\n"
+"                              generate resolution values in the output file.\n"
 );
 }
 
 
 /* Apply Poisson noise to all reflections */
-static void poisson_reflections(double *ref, ReflItemList *items,
-                                double adu_per_photon)
+static void poisson_reflections(RefList *list, double adu_per_photon)
 {
-	int i;
-	const int n = num_items(items);
+	Reflection *refl;
+	RefListIterator *iter;
 
-	for ( i=0; i<n; i++ ) {
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
-		double val;
-		int c;
+		double val, c;
 
-		it = get_item(items, i);
+		val = get_intensity(refl);
 
-		val = lookup_intensity(ref, it->h, it->k, it->l);
 		c = adu_per_photon * poisson_noise(val/adu_per_photon);
-		set_intensity(ref, it->h, it->k, it->l, c);
-
-		progress_bar(i, n-1, "Simulating noise");
+		set_int(refl, c);
 
 	}
 }
 
 
 /* Apply 10% uniform noise to all reflections */
-static void noise_reflections(double *ref, ReflItemList *items)
+static void noise_reflections(RefList *list)
 {
-	int i;
-	const int n = num_items(items);
-
-	for ( i=0; i<n; i++ ) {
+	Reflection *refl;
+	RefListIterator *iter;
 
-		struct refl_item *it;
-		double val;
-		double r;
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		it = get_item(items, i);
+		double val, r;
 
-		val = lookup_intensity(ref, it->h, it->k, it->l);
+		val = get_intensity(refl);
 
 		r = (double)random()/RAND_MAX;
 		val += 0.1 * val * r;
 
-		set_intensity(ref, it->h, it->k, it->l, val);
-
-		progress_bar(i, n-1, "Simulating noise");
+		set_int(refl, val);
 
 	}
 }
 
 
-static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
-                                      const char *holo, const char *mero,
-                                      double *esds)
+static RefList *template_reflections(RefList *list, RefList *template)
 {
-	int i;
-	ReflItemList *new;
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *out;
 
-	new = new_items();
+	out = reflist_new();
+
+	for ( refl = first_refl(template, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+		Reflection *new;
+		Reflection *old;
+
+		get_indices(refl, &h, &k, &l);
+
+		old = find_refl(list, h, k, l);
+		if ( old == NULL ) continue;
+
+		new = add_refl(out, h, k, l);
+		copy_data(new, old);
+
+	}
+
+	return out;
+}
+
+
+static RefList *twin_reflections(RefList *in,
+                                 const char *holo, const char *mero)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *out;
 
 	if ( num_general_equivs(holo) < num_general_equivs(mero) ) {
 		ERROR("%s is not a subgroup of %s!\n", mero, holo);
 		return NULL;
 	}
 
-	for ( i=0; i<num_items(items); i++ ) {
+	out = reflist_new();
+
+	for ( refl = first_refl(in, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
 		double total, sigma;
-		struct refl_item *it;
 		signed int h, k, l;
 		int n, j;
 		int skip;
 
-		it = get_item(items, i);
+		get_indices(refl, &h, &k, &l);
 
 		/* There is a many-to-one correspondence between reflections
 		 * in the merohedral and holohedral groups.  Do the calculation
 		 * only once for each reflection in the holohedral group, which
 		 * contains fewer reflections.
 		 */
-		get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
-		if ( find_item(new, h, k, l) ) continue;
-
-		n = num_equivs(h, k, l, holo);
+		get_asymm(h, k, l, &h, &k, &l, holo);
+		if ( find_refl(out, h, k, l) != NULL ) continue;
 
 		total = 0.0;
 		sigma = 0.0;
 		skip = 0;
+		n = num_equivs(h, k, l, holo);
 		for ( j=0; j<n; j++ ) {
 
 			signed int he, ke, le;
@@ -164,8 +189,8 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
 			 * equivalent which belongs to our definition of the
 			 * asymmetric unit cell, so check them all.
 			 */
-			if ( !find_unique_equiv(items, he, ke, le, mero,
-			                        &hu, &ku, &lu) ) {
+			if ( !find_equiv_in_list(in, he, ke, le, mero,
+			                         &hu, &ku, &lu) ) {
 				/* Don't have this reflection, so bail out */
 				ERROR("Twinning %i %i %i requires the %i %i %i "
 				      "reflection (or an equivalent in %s), "
@@ -176,61 +201,57 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
 				break;
 			}
 
-			total += lookup_intensity(ref, hu, ku, lu);
-			sigma += pow(lookup_sigma(esds, hu, ku, lu), 2.0);
+			total += get_intensity(refl);
+			sigma += pow(get_esd_intensity(refl), 2.0);
 
 		}
 
 		if ( !skip ) {
 
-			set_intensity(ref, h, k, l, total);
-			set_sigma(esds, h, k, l, sqrt(sigma));
-			add_item(new, h, k, l);
+			Reflection *new = add_refl(out, h, k, l);
+			set_int(new, total);
+			set_esd_intensity(new, sqrt(sigma));
 
 		}
 
 	}
 
-	return new;
+	return out;
 }
 
 
-static ReflItemList *expand_reflections(double *ref, ReflItemList *items,
-                                        const char *target, const char *initial)
+static RefList *expand_reflections(RefList *in,
+                                   const char *target, const char *initial)
 {
-	int i;
-	ReflItemList *new;
-
-	new = new_items();
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *out;
 
 	if ( num_general_equivs(target) > num_general_equivs(initial) ) {
 		ERROR("%s is not a subgroup of %s!\n", initial, target);
 		return NULL;
 	}
 
-	for ( i=0; i<num_items(items); i++ ) {
+	out = reflist_new();
+
+	for ( refl = first_refl(in, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
 		signed int h, k, l;
-		signed int hd, kd, ld;
 		int n, j;
 		double intensity;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
-
-		/* Actually we don't really care what the equivalent is,
-		 * we just want to be sure that there is nly be one version of
-		 * this reflection. */
-		find_unique_equiv(items, h, k, l, initial, &hd, &kd, &ld);
+		get_indices(refl, &h, &k, &l);
+		intensity = get_intensity(refl);
 
-		/* Now find out how many reflections need to be filled in */
 		n = num_equivs(h, k, l, initial);
-		intensity = lookup_intensity(ref, h, k, l);
 
+		/* For each equivalent in the higher symmetry group */
 		for ( j=0; j<n; j++ ) {
 
 			signed int he, ke, le;
+			Reflection *new;
 
 			/* Get the equivalent */
 			get_equiv(h, k, l, &he, &ke, &le, initial, j);
@@ -239,43 +260,33 @@ static ReflItemList *expand_reflections(double *ref, ReflItemList *items,
 			get_asymm(he, ke, le, &he, &ke, &le, target);
 
 			/* Make sure the intensity is in the right place */
-			set_intensity(ref, he, ke, le, intensity);
-
-			/* Add the reflection, but only once */
-			if ( !find_item(new, he, ke, le) ) {
-				add_item(new, he, ke, le);
-			}
+			new = add_refl(out, he, ke, le);
+			copy_data(new, refl);
 
 		}
 
 	}
 
-	return new;
+	return out;
 }
 
 
 int main(int argc, char *argv[])
 {
 	int c;
-	double *input_ref;
-	double *phases;
-	double *esds;
-	char *template = NULL;
 	int config_noise = 0;
 	int config_poisson = 0;
-	int config_nophase = 0;
 	int config_multi = 0;
 	char *holo = NULL;
 	char *mero = NULL;
 	char *expand = NULL;
+	char *input_file = NULL;
+	char *template = NULL;
 	char *output = NULL;
-	char *input = NULL;
-	char *filename = NULL;
-	ReflItemList *input_items;
-	ReflItemList *write_items;
-	UnitCell *cell = NULL;
 	char *beamfile = NULL;
 	struct beam_params *beam = NULL;
+	RefList *input;
+	UnitCell *cell = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -290,6 +301,7 @@ int main(int argc, char *argv[])
 		{"intensities",        1, NULL,               'i'},
 		{"multiplicity",       0, &config_multi,       1},
 		{"beam",               1, NULL,               'b'},
+		{"pdb",                1, NULL,               'p'},
 		{0, 0, NULL, 0}
 	};
 
@@ -311,7 +323,7 @@ int main(int argc, char *argv[])
 			break;
 
 		case 'i' :
-			input = strdup(optarg);
+			input_file = strdup(optarg);
 			break;
 
 		case 'y' :
@@ -330,6 +342,14 @@ int main(int argc, char *argv[])
 			beamfile = strdup(optarg);
 			break;
 
+		case 'p' :
+			cell = load_cell_from_pdb(optarg);
+			if ( cell == NULL ) {
+				ERROR("Failed to get cell from '%s'\n", optarg);
+				return 1;
+			}
+			break;
+
 		case 0 :
 			break;
 
@@ -339,10 +359,6 @@ int main(int argc, char *argv[])
 
 	}
 
-	if ( filename == NULL ) {
-		filename = strdup("molecule.pdb");
-	}
-
 	if ( (holo != NULL) && (expand != NULL) ) {
 		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
 		ERROR("Decide which one you want to do first.\n");
@@ -358,19 +374,14 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	cell = load_cell_from_pdb(filename);
-	if ( !config_nophase ) {
-		phases = new_list_phase();
-	} else {
-		phases = NULL;
+	if ( cell == NULL ) {
+		ERROR("You need to give a PDB file with the unit cell.\n");
+		return 1;
 	}
-	esds = new_list_sigma();
 
-	input_ref = new_list_intensity();
-	input_items = read_reflections(input, input_ref, phases,
-	                               NULL, esds);
+	input = read_reflections(input_file);
 	free(input);
-	if ( check_symmetry(input_items, mero) ) {
+	if ( check_list_symmetry(input, mero) ) {
 		ERROR("The input reflection list does not appear to"
 		      " have symmetry %s\n", mero);
 		return 1;
@@ -378,8 +389,7 @@ int main(int argc, char *argv[])
 
 	if ( config_poisson ) {
 		if ( beam != NULL ) {
-			poisson_reflections(input_ref, input_items,
-			                    beam->adu_per_photon);
+			poisson_reflections(input, beam->adu_per_photon);
 		} else {
 			ERROR("You must give a beam parameters file in order"
 			      " to calculate Poisson noise.\n");
@@ -387,72 +397,69 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	if ( config_noise ) noise_reflections(input_ref, input_items);
+	if ( config_noise ) noise_reflections(input);
 
 	if ( holo != NULL ) {
 
-		ReflItemList *new;
+		RefList *new;
 		STATUS("Twinning from %s into %s\n", mero, holo);
-		new = twin_reflections(input_ref, input_items,
-		                       holo, mero, esds);
-		delete_items(input_items);
-		input_items = new;
+		new = twin_reflections(input, holo, mero);
+
+		/* Replace old with new */
+		reflist_free(input);
+		input = new;
+
+		/* The symmetry of the list has changed */
+		free(mero);
+		mero = holo;
 
 	}
 
 	if ( expand != NULL ) {
 
-		ReflItemList *new;
+		RefList *new;
 		STATUS("Expanding from %s into %s\n", mero, expand);
-		new = expand_reflections(input_ref, input_items, expand, mero);
-		delete_items(input_items);
-		input_items = new;
+		new = expand_reflections(input, expand, mero);
+
+		/* Replace old with new */
+		reflist_free(input);
+		input = new;
 
 	}
 
 	if ( config_multi ) {
 
-		int i;
+		Reflection *refl;
+		RefListIterator *iter;
 
-		for ( i=0; i<num_items(input_items); i++ ) {
+		for ( refl = first_refl(input, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) ) {
 
-			struct refl_item *it;
 			double inty;
+			signed int h, k, l;
+
+			get_indices(refl, &h, &k, &l);
+			inty = get_intensity(refl);
 
-			it = get_item(input_items, i);
-			inty = lookup_intensity(input_ref, it->h, it->k, it->l);
-			inty *= num_equivs(it->h, it->k, it->l, mero);
-			set_intensity(input_ref, it->h, it->k, it->l, inty);
-			STATUS("%i %i %i %i\n", it->h, it->k, it->l,
-			       num_equivs(it->h, it->k, it->l, mero));
+			inty *= (double)num_equivs(h, k, l, mero);
+			set_int(refl, inty);
 
 		}
 	}
 
 	if ( template ) {
 
-		/* Write out only reflections which are in the template
-		 * (and which we have in the input) */
-		ReflItemList *template_items;
-		template_items = read_reflections(template,
-		                                  NULL, NULL, NULL, NULL);
-		write_items = intersection_items(input_items, template_items);
-		delete_items(template_items);
-
-	} else {
-
-		/* Write out all reflections */
-		write_items = new_items();
-		/* (quick way of copying a list) */
-		union_items(write_items, input_items);
+		RefList *t = read_reflections(template);
+		RefList *new = template_reflections(input, t);
+		reflist_free(input);
+		input = new;
 
 	}
 
-	write_reflections(output, write_items, input_ref, esds, phases,
-	                  NULL, cell);
+	write_reflist(output, input, cell);
 
-	delete_items(input_items);
-	delete_items(write_items);
+	reflist_free(input);
 
 	return 0;
 }
diff --git a/src/indexamajig.c b/src/indexamajig.c
index e1e0e8fd..c0757edd 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -31,11 +31,11 @@
 #include "peaks.h"
 #include "detector.h"
 #include "filters.h"
-#include "reflections.h"
 #include "thread-pool.h"
 #include "beam-parameters.h"
 #include "geometry.h"
 #include "stream.h"
+#include "reflist-utils.h"
 
 
 enum {
diff --git a/src/partialator.c b/src/partialator.c
index 88926582..553a216a 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -26,7 +26,6 @@
 #include "utils.h"
 #include "hdf5-file.h"
 #include "symmetry.h"
-#include "reflections.h"
 #include "stream.h"
 #include "geometry.h"
 #include "peaks.h"
@@ -35,6 +34,7 @@
 #include "post-refinement.h"
 #include "hrs-scaling.h"
 #include "reflist.h"
+#include "reflist-utils.h"
 
 
 static void show_help(const char *s)
@@ -418,7 +418,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Output results */
-	write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL);
+	//write_reflist(outfile, obs, I_full, NULL, NULL, cts, NULL);
 
 	/* Clean up */
 	for ( i=0; i<n_total_patterns; i++ ) {
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 6d3c07ab..ba760e32 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -29,9 +29,10 @@
 #include "hdf5-file.h"
 #include "detector.h"
 #include "peaks.h"
-#include "reflections.h"
 #include "beam-parameters.h"
 #include "symmetry.h"
+#include "reflist.h"
+#include "reflist-utils.h"
 
 
 static void show_help(const char *s)
@@ -412,34 +413,29 @@ int main(int argc, char *argv[])
 
 	} else {
 
-		int i;
-		ReflItemList *items;
+		RefList *reflections;
+
+		reflections = read_reflections(intfile);
+		free(intfile);
 
 		if ( grad == GRADIENT_PHASED ) {
-			phases = new_list_phase();
+			phases = phases_from_list(reflections);
 		} else {
 			phases = NULL;
 		}
-		intensities = new_list_intensity();
-		phases = new_list_phase();
-		flags = new_list_flag();
-		items = read_reflections(intfile, intensities, phases,
-		                         NULL, NULL);
-		free(intfile);
-
-		for ( i=0; i<num_items(items); i++ ) {
-			struct refl_item *it = get_item(items, i);
-			set_flag(flags, it->h, it->k, it->l, 1);
-		}
+		intensities = intensities_from_list(reflections);
+		phases = phases_from_list(reflections);
+		flags = flags_from_list(reflections);
 
 		/* Check that the intensities have the correct symmetry */
-		if ( check_symmetry(items, sym) ) {
+		if ( check_list_symmetry(reflections, sym) ) {
 			ERROR("The input reflection list does not appear to"
 			      " have symmetry %s\n", sym);
 			return 1;
 		}
 
-		delete_items(items);
+		reflist_free(reflections);
+
 	}
 
 	image.det = get_detector_geometry(geometry);
diff --git a/src/povray.c b/src/povray.c
index 808d3e22..8d0ec048 100644
--- a/src/povray.c
+++ b/src/povray.c
@@ -23,13 +23,13 @@
 #include "utils.h"
 #include "symmetry.h"
 #include "render_hkl.h"
+#include "povray.h"
 
 
 #define MAX_PROC (256)
 
 
-int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
-                            ReflItemList *items, unsigned int nproc,
+int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
                             const char *sym, int wght, double boost)
 {
 	FILE *fh;
@@ -39,6 +39,8 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
 	pid_t pids[MAX_PROC];
 	float max;
 	int i;
+	Reflection *refl;
+	RefListIterator *iter;
 
 	if ( (nproc > MAX_PROC) || (nproc < 1) ) {
 		ERROR("Number of processes must be a number between 1 and %i\n",
@@ -165,27 +167,29 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
 	fprintf(fh, "}\n");
 
 	max = 0.0;
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
 		float val;
+		signed int h, k, l;
 
-		it = get_item(items, i);
+		get_indices(refl, &h, &k, &l);
 
 		switch ( wght ) {
 		case WGHT_I :
-			val = lookup_intensity(ref, it->h, it->k, it->l);
+			val = get_intensity(refl);
 			break;
 		case WGHT_SQRTI :
-			val = lookup_intensity(ref, it->h, it->k, it->l);
+			val = get_intensity(refl);
 			val = (val>0.0) ? sqrt(val) : 0.0;
 			break;
 		case WGHT_COUNTS :
-			val = lookup_count(counts, it->h, it->k, it->l);
-			val /= (float)num_equivs(it->h, it->k, it->l, sym);
+			val = get_redundancy(refl);
+			val /= (float)num_equivs(h, k, l, sym);
 			break;
 		case WGHT_RAWCOUNTS :
-			val = lookup_count(counts, it->h, it->k, it->l);
+			val = get_redundancy(refl);
 			break;
 		default :
 			ERROR("Invalid weighting.\n");
@@ -197,30 +201,31 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
 	}
 	max /= boost;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
 
-		struct refl_item *it;
-		float radius;
+		signed int h, k, l;float radius;
 		int s;
 		float val, p, r, g, b, trans;
 		int j;
 
-		it = get_item(items, i);
+		get_indices(refl, &h, &k, &l);
 
 		switch ( wght ) {
 		case WGHT_I :
-			val = lookup_intensity(ref, it->h, it->k, it->l);
+			val = get_intensity(refl);
 			break;
 		case WGHT_SQRTI :
-			val = lookup_intensity(ref, it->h, it->k, it->l);
+			val = get_intensity(refl);
 			val = (val>0.0) ? sqrt(val) : 0.0;
 			break;
 		case WGHT_COUNTS :
-			val = lookup_count(counts, it->h, it->k, it->l);
-			val /= (float)num_equivs(it->h, it->k, it->l, sym);
+			val = get_redundancy(refl);
+			val /= (float)num_equivs(h, k, l, sym);
 			break;
 		case WGHT_RAWCOUNTS :
-			val = lookup_count(counts, it->h, it->k, it->l);
+			val = get_redundancy(refl);
 			break;
 		default :
 			ERROR("Invalid weighting.\n");
@@ -269,12 +274,12 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
 		trans = 1.0-(val/max);
 
 		/* For each equivalent */
-		for ( j=0; j<num_equivs(it->h, it->k, it->l, sym); j++ ) {
+		for ( j=0; j<num_equivs(h, k, l, sym); j++ ) {
 
 			signed int he, ke, le;
 			float x, y, z;
 
-			get_equiv(it->h, it->k, it->l, &he, &ke, &le, sym, j);
+			get_equiv(h, k, l, &he, &ke, &le, sym, j);
 
 			x = asx*he + bsx*ke + csx*le;
 			y = asy*he + bsy*ke + csy*le;
diff --git a/src/povray.h b/src/povray.h
index 69bf657e..0b0b3446 100644
--- a/src/povray.h
+++ b/src/povray.h
@@ -3,7 +3,7 @@
  *
  * Invoke POV-ray
  *
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
  *
  * Part of CrystFEL - crystallography with a FEL
  *
@@ -16,8 +16,10 @@
 #ifndef POVRAY_H
 #define POVRAY_H
 
-extern int povray_render_animation(UnitCell *cell, double *ref,
-                                   unsigned int *counts, ReflItemList *items,
+#include "reflist.h"
+#include "cell.h"
+
+extern int povray_render_animation(UnitCell *cell, RefList *list,
                                    unsigned int nproc, const char *sym,
                                    int wght, double boost);
 
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 6ff40159..f5e705bd 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -468,11 +468,11 @@ int main(int argc, char *argv[])
 		output = strdup("processed.hkl");
 	}
 
-	if ( pdb == NULL ) {
-		pdb = strdup("molecule.pdb");
-	}
-
 	cell = load_cell_from_pdb(pdb);
+	if ( cell == NULL ) {
+		ERROR("Failed to load cell from '%s'\n", pdb);
+		return 1;
+	}
 	free(pdb);
 
 	if ( sym == NULL ) sym = strdup("1");
@@ -537,7 +537,6 @@ int main(int argc, char *argv[])
 
 	write_reflist(output, model, cell);
 
-	cell_free(cell);
 	fclose(fh);
 	free(sym);
 	reflist_free(model);
diff --git a/src/reflections.c b/src/reflections.c
deleted file mode 100644
index 4743aeee..00000000
--- a/src/reflections.c
+++ /dev/null
@@ -1,230 +0,0 @@
-/*
- * reflections.c
- *
- * Utilities for handling reflections
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <complex.h>
-#include <string.h>
-
-#include "utils.h"
-#include "cell.h"
-#include "reflections.h"
-#include "beam-parameters.h"
-
-
-double *poisson_esds(double *intensities, ReflItemList *items,
-                     double adu_per_photon)
-{
-	int i;
-	double *esds = new_list_sigma();
-
-	for ( i=0; i<num_items(items); i++ ) {
-
-		struct refl_item *it;
-		signed int h, k, l;
-		double intensity, sigma;
-
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
-
-		intensity = lookup_intensity(intensities, h, k, l);
-
-		if ( intensity > 0.0 ) {
-			sigma = adu_per_photon * sqrt(intensity/adu_per_photon);
-		} else {
-			sigma = 0.0;
-		}
-
-		set_sigma(esds, h, k, l, sigma);
-
-	}
-
-	return esds;
-}
-
-
-void write_reflections(const char *filename, ReflItemList *items,
-                       double *intensities, double *esds, double *phases,
-                       unsigned int *counts, UnitCell *cell)
-{
-	FILE *fh;
-	int i;
-
-	if ( filename == NULL ) {
-		fh = stdout;
-	} else {
-		fh = fopen(filename, "w");
-	}
-
-	if ( fh == NULL ) {
-		ERROR("Couldn't open output file '%s'.\n", filename);
-		return;
-	}
-
-	/* Write spacings and angle if zone axis pattern */
-	fprintf(fh, "  h   k   l          I    phase   sigma(I) "
-		    " 1/d(nm^-1)  counts\n");
-
-	for ( i=0; i<num_items(items); i++ ) {
-
-		struct refl_item *it;
-		signed int h, k, l;
-
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
-
-		int N;
-		double intensity, s, sigma;
-		char ph[32];
-
-		if ( counts != NULL ) {
-			N = lookup_count(counts, h, k, l);
-		} else {
-			N = 1;
-		}
-
-		intensity = lookup_intensity(intensities, h, k, l);
-
-		if ( phases != NULL ) {
-			double p;
-			p = lookup_phase(phases, h, k, l);
-			snprintf(ph, 31, "%8.6f", p);
-		} else {
-			strncpy(ph, "       -", 31);
-		}
-
-		if ( cell != NULL ) {
-			s = 2.0*resolution(cell, h, k, l);
-		} else {
-			s = 0.0;
-		}
-
-		if ( esds != NULL ) {
-			sigma = lookup_sigma(esds, h, k, l);
-		} else {
-			sigma = 0.0;
-		}
-
-		/* h, k, l, I, sigma(I), s */
-		fprintf(fh, "%3i %3i %3i %10.2f %s %10.2f  %10.2f %7i\n",
-		        h, k, l, intensity, ph, sigma, s/1.0e9, N);
-
-	}
-	fclose(fh);
-}
-
-
-/* Read reflections from file.  Returns the list of reflections successfully
- * read in.  "intensities", "phases" and "counts" are lists which will be
- * populated with the values read from the file.  Existing values in either list
- * will be overwritten if the reflection is read from the file, but other values
- * will be left intact.
- *
- * "intensities", "phases" or "counts" can be NULL, if you don't need them.
- */
-ReflItemList *read_reflections(const char *filename,
-                               double *intensities, double *phases,
-                               unsigned int *counts, double *esds)
-{
-	FILE *fh;
-	char *rval;
-	ReflItemList *items;
-
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-		ERROR("Failed to open input file '%s'\n", filename);
-		return NULL;
-	}
-
-	items = new_items();
-
-	do {
-
-		char line[1024];
-		signed int h, k, l;
-		float intensity, ph, res, sigma;
-		char phs[1024];
-		int r;
-		int cts;
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		r = sscanf(line, "%i %i %i %f %s %f %f %i",
-		           &h, &k, &l, &intensity, phs, &sigma, &res, &cts);
-		if ( r >= 8 ) {
-			/* Woohoo */
-		} else if ( r >= 7 ) {
-			/* No "counts", that's fine.. */
-			cts = 1;
-		} else if ( r >= 6 ) {
-			/* No resolution.  Didn't want it anyway. */
-			res = 0.0;
-		} else if ( r >= 5 ) {
-			/* No sigma.  It's OK today, but one
-			 * day I'll get you... */
-			sigma = 0.0;
-		} else if ( r >= 4 ) {
-			/* No phase.  Better not need it.. */
-			if ( phases != NULL ) {
-				ERROR("Need phases and none were specified!\n");
-				abort();
-			}
-		} else {
-			/* You lose. */
-			continue;
-		}
-
-		add_item(items, h, k, l);
-
-		if ( intensities != NULL ) {
-			set_intensity(intensities, h, k, l, intensity);
-		}
-		if ( phases != NULL ) {
-			ph = atof(phs);
-			set_phase(phases, h, k, l, ph);
-		}
-		if ( counts != NULL ) {
-			set_count(counts, h, k, l, cts);
-		}
-		if ( esds != NULL ) {
-			set_sigma(esds, h, k, l, sigma);
-		}
-
-	} while ( rval != NULL );
-
-	fclose(fh);
-
-	return items;
-}
-
-
-double *ideal_intensities(double complex *sfac)
-{
-	double *ref;
-	signed int h, k, l;
-
-	ref = new_list_intensity();
-
-	/* Generate ideal reflections from complex structure factors */
-	for ( h=-INDMAX; h<=INDMAX; h++ ) {
-	for ( k=-INDMAX; k<=INDMAX; k++ ) {
-	for ( l=-INDMAX; l<=INDMAX; l++ ) {
-		double complex F = lookup_sfac(sfac, h, k, l);
-		double intensity = pow(cabs(F), 2.0);
-		set_intensity(ref, h, k, l, intensity);
-	}
-	}
-	}
-
-	return ref;
-}
diff --git a/src/reflections.h b/src/reflections.h
deleted file mode 100644
index 20060f21..00000000
--- a/src/reflections.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/*
- * reflections.h
- *
- * Utilities for handling reflections
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef REFLECTIONS_H
-#define REFLECTIONS_H
-
-
-#include "cell.h"
-#include "utils.h"
-
-
-extern double *poisson_esds(double *intensities, ReflItemList *items,
-                            double adu_per_photon);
-
-extern void write_reflections(const char *filename, ReflItemList *items,
-                              double *intensities, double *esds, double *phases,
-                              unsigned int *counts, UnitCell *cell);
-
-extern ReflItemList *read_reflections(const char *filename,
-                                      double *intensities, double *phases,
-                                      unsigned int *counts, double *esds);
-
-extern double *ideal_intensities(double complex *sfac);
-
-
-#endif	/* REFLECTIONS_H */
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index 33d373b8..aac94fbf 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -15,6 +15,153 @@
 
 #include "reflist.h"
 #include "cell.h"
+#include "utils.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+
+
+double *intensities_from_list(RefList *list)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	double *out = new_list_intensity();
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+		double intensity = get_intensity(refl);
+
+		get_indices(refl, &h, &k, &l);
+
+		set_intensity(out, h, k, l, intensity);
+
+	}
+
+	return out;
+}
+
+
+double *phases_from_list(RefList *list)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	double *out = new_list_phase();
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+		double phase = get_phase(refl);
+
+		get_indices(refl, &h, &k, &l);
+
+		set_phase(out, h, k, l, phase);
+
+	}
+
+	return out;
+
+}
+
+
+unsigned char *flags_from_list(RefList *list)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	unsigned char *out = new_list_flag();
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+
+		get_indices(refl, &h, &k, &l);
+
+		set_flag(out, h, k, l, 1);
+
+	}
+
+	return out;
+
+}
+
+
+int check_list_symmetry(RefList *list, const char *sym)
+{
+	unsigned char *flags;
+	Reflection *refl;
+	RefListIterator *iter;
+
+	flags = flags_from_list(list);
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		int j;
+		int found = 0;
+		signed int h, k, l;
+
+		get_indices(refl, &h, &k, &l);
+
+		for ( j=0; j<num_equivs(h, k, l, sym); j++ ) {
+
+			signed int he, ke, le;
+			get_equiv(h, k, l, &he, &ke, &le, sym, j);
+
+			if ( abs(he) > INDMAX ) continue;
+			if ( abs(le) > INDMAX ) continue;
+			if ( abs(ke) > INDMAX ) continue;
+
+			found += lookup_flag(flags, he, ke, le);
+
+		}
+
+		if ( found > 1 ) {
+			free(flags);
+			return 1;  /* Symmetry is wrong! */
+		}
+
+	}
+
+	free(flags);
+
+	return 0;
+}
+
+
+int find_equiv_in_list(RefList *list, signed int h, signed int k,
+                       signed int l, const char *sym, signed int *hu,
+                       signed int *ku, signed int *lu)
+{
+	int i;
+	int found = 0;
+
+	for ( i=0; i<num_equivs(h, k, l, sym); i++ ) {
+
+		signed int he, ke, le;
+		Reflection *f;
+		get_equiv(h, k, l, &he, &ke, &le, sym, i);
+		f = find_refl(list, he, ke, le);
+
+		/* There must only be one equivalent.  If there are more, it
+		 * indicates that the user lied about the input symmetry.
+		 * This situation should have been checked for earlier by
+		 * calling check_symmetry() with 'items' and 'mero'. */
+
+		if ( (f != NULL) && !found ) {
+			*hu = he;  *ku = ke;  *lu = le;
+			return 1;
+		}
+
+	}
+
+	return 0;
+}
 
 
 void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
@@ -31,18 +178,26 @@ void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
 
 		signed int h, k, l;
 		double intensity, esd_i, s;
+		int red;
 		double fs, ss;
+		char res[16];
 
 		get_indices(refl, &h, &k, &l);
 		get_detector_pos(refl, &fs, &ss);
 		intensity = get_intensity(refl);
-		esd_i = 0.0; /* FIXME! */
-		s = 0.0; /* FIXME! */
+		esd_i = get_esd_intensity(refl);
+		red = get_redundancy(refl);
+
+		if ( cell != NULL ) {
+			s = resolution(cell, h, k, l);
+			snprintf(res, 16, "%10.2f", s/1e9);
+		} else {
+			strcpy(res, "         -");
+		}
 
-		/* h, k, l, I, sigma(I), s */
 		fprintf(fh,
-		       "%3i %3i %3i %10.2f %s %10.2f  %10.2f %7i %6.1f %6.1f\n",
-		       h, k, l, intensity, "       -", esd_i, s/1.0e9, 1,
+		       "%3i %3i %3i %10.2f %s %10.2f  %s %7i %6.1f %6.1f\n",
+		       h, k, l, intensity, "       -", esd_i, res, red,
 		       fs, ss);
 
 	}
@@ -70,3 +225,85 @@ int write_reflist(const char *filename, RefList *list, UnitCell *cell)
 
 	return 0;
 }
+
+
+RefList *read_reflections_from_file(FILE *fh)
+{
+	char *rval = NULL;
+	int first = 1;
+	RefList *out;
+
+	out = reflist_new();
+
+	do {
+
+		char line[1024];
+		signed int h, k, l;
+		float intensity, sigma, fs, ss;
+		char phs[1024];
+		char ress[1024];
+		int cts;
+		int r;
+		Reflection *refl;
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+		chomp(line);
+
+		if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out;
+
+		r = sscanf(line, "%i %i %i %f %s %f %s %i %f %f",
+		           &h, &k, &l, &intensity, phs, &sigma, ress, &cts,
+		           &fs, &ss);
+		if ( (r != 10) && (!first) ) {
+			reflist_free(out);
+			return NULL;
+		}
+
+		first = 0;
+		if ( r == 10 ) {
+
+			double ph;
+			char *v;
+
+			refl = add_refl(out, h, k, l);
+			set_int(refl, intensity);
+			set_detector_pos(refl, fs, ss, 0.0);
+			set_esd_intensity(refl, sigma);
+
+			ph = strtod(phs, &v);
+			if ( v != NULL ) set_ph(refl, ph);
+
+			/* The 1/d value is actually ignored. */
+
+		}
+
+	} while ( rval != NULL );
+
+	/* Got read error of some kind before finding PEAK_LIST_END_MARKER */
+	return NULL;
+}
+
+
+RefList *read_reflections(const char *filename)
+{
+	FILE *fh;
+	RefList *out;
+
+	if ( filename == NULL ) {
+		fh = stdout;
+	} else {
+		fh = fopen(filename, "r");
+	}
+
+	if ( fh == NULL ) {
+		ERROR("Couldn't open input file '%s'.\n", filename);
+		return NULL;
+	}
+
+	out = read_reflections_from_file(fh);
+
+	fclose(fh);
+
+	return out;
+}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
index d5c1e655..b7ecd3ea 100644
--- a/src/reflist-utils.h
+++ b/src/reflist-utils.h
@@ -21,10 +21,24 @@
 #include "cell.h"
 
 
+#define REFLECTION_END_MARKER "End of reflections"
+
+
 extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell);
 
-extern int write_reflist(const char *filename, RefList *list,
-                             UnitCell *cell);
+extern int write_reflist(const char *filename, RefList *list, UnitCell *cell);
+
+extern RefList *read_reflections_from_file(FILE *fh);
+
+extern RefList *read_reflections(const char *filename);
+
+extern double *intensities_from_list(RefList *list);
+extern double *phases_from_list(RefList *list);
+extern unsigned char *flags_from_list(RefList *list);
 
+extern int check_list_symmetry(RefList *list, const char *sym);
+extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
+                              signed int l, const char *sym, signed int *hu,
+                              signed int *ku, signed int *lu);
 
 #endif	/* REFLIST_UTILS_H */
diff --git a/src/reflist.c b/src/reflist.c
index 3181b418..a81c738e 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -46,6 +46,9 @@ struct _refldata {
 	double intensity;
 	double esd_i;
 
+	/* Phase */
+	double phase;
+
 	/* Redundancy */
 	int redundancy;
 
@@ -250,8 +253,20 @@ double get_esd_intensity(Reflection *refl)
 }
 
 
+double get_phase(Reflection *refl)
+{
+	return refl->data.phase;
+}
+
+
 /********************************** Setters ***********************************/
 
+void copy_data(Reflection *to, Reflection *from)
+{
+	memcpy(&to->data, &from->data, sizeof(struct _refldata));
+}
+
+
 void set_detector_pos(Reflection *refl, double exerr, double x, double y)
 {
 	refl->data.excitation_error = exerr;
@@ -301,6 +316,12 @@ void set_esd_intensity(Reflection *refl, double esd)
 }
 
 
+void set_ph(Reflection *refl, double phase)
+{
+	refl->data.phase = phase;
+}
+
+
 /********************************* Insertion **********************************/
 
 static void insert_node(Reflection *head, Reflection *new)
@@ -689,3 +710,9 @@ void optimise_reflist(RefList *list)
 		       recursive_depth(list->head->child[0]));
 	}
 }
+
+
+int num_reflections(RefList *list)
+{
+	return recursive_count(list->head->child[0]);
+}
diff --git a/src/reflist.h b/src/reflist.h
index f6d01ecf..5b647102 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -44,8 +44,10 @@ extern int get_scalable(Reflection *refl);
 extern int get_redundancy(Reflection *refl);
 extern double get_sum_squared_dev(Reflection *refl);
 extern double get_esd_intensity(Reflection *refl);
+extern double get_phase(Reflection *refl);
 
 /* Set */
+extern void copy_data(Reflection *to, Reflection *from);
 extern void set_detector_pos(Reflection *refl, double exerr,
                              double x, double y);
 extern void set_partial(Reflection *refl, double r1, double r2, double p,
@@ -55,6 +57,7 @@ extern void set_scalable(Reflection *refl, int scalable);
 extern void set_redundancy(Reflection *refl, int red);
 extern void set_sum_squared_dev(Reflection *refl, double dev);
 extern void set_esd_intensity(Reflection *refl, double esd);
+extern void set_ph(Reflection *refl, double phase);
 
 /* Insertion */
 extern Reflection *add_refl(RefList *list, INDICES);
@@ -66,7 +69,8 @@ extern void delete_refl(Reflection *refl);
 extern Reflection *first_refl(RefList *list, RefListIterator **iterator);
 extern Reflection *next_refl(Reflection *refl, RefListIterator *iter);
 
-/* Voodoo */
+/* Misc */
 extern void optimise_reflist(RefList *list);
+extern int num_reflections(RefList *list);
 
 #endif	/* REFLIST_H */
diff --git a/src/render_hkl.c b/src/render_hkl.c
index fc1c5121..aec6d714 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -26,11 +26,12 @@
 #endif
 
 #include "utils.h"
-#include "reflections.h"
 #include "povray.h"
 #include "symmetry.h"
 #include "render.h"
 #include "render_hkl.h"
+#include "reflist.h"
+#include "reflist-utils.h"
 
 
 static void show_help(const char *s)
@@ -81,8 +82,7 @@ static void show_help(const char *s)
 static void draw_circles(signed int xh, signed int xk, signed int xl,
                          signed int yh, signed int yk, signed int yl,
                          signed int zh, signed int zk, signed int zl,
-                         double *ref, unsigned int *counts, ReflItemList *items,
-                         const char *sym,
+                         RefList *list, const char *sym,
                          cairo_t *dctx, int wght, double boost, int colscale,
                          UnitCell *cell, double radius, double theta,
                          double as, double bs, double cx, double cy,
@@ -98,36 +98,36 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
 		*max_res = 0.0;  *max_ux = 0;  *max_uy = 0;
 	}
 
-	/* Loop across the two basis directions */
+	/* Iterate over possible reflections in this zone */
 	for ( xi=-INDMAX; xi<INDMAX; xi++ ) {
 	for ( yi=-INDMAX; yi<INDMAX; yi++ ) {
 
 		double u, v, val, res;
 		signed int h, k, l;
-		signed int he, ke, le;
+		Reflection *refl;
 
 		h = xi*xh + yi*yh;
 		k = xi*xk + yi*yk;
 		l = xi*xl + yi*yl;
 
 		/* Got this reflection? */
-		if ( find_unique_equiv(items, h, k, l, sym,
-		                       &he, &ke, &le) == 0 ) continue;
+		refl = find_refl(list, h, k, l);
+		if ( refl == NULL ) continue;
 
 		switch ( wght ) {
 		case WGHT_I :
-			val = lookup_intensity(ref, he, ke, le);
+			val = get_intensity(refl);
 			break;
 		case WGHT_SQRTI :
-			val = lookup_intensity(ref, he, ke, le);
+			val = get_intensity(refl);
 			val = (val>0.0) ? sqrt(val) : 0.0;
 			break;
 		case WGHT_COUNTS :
-			val = lookup_count(counts, he, ke, le);
-			val /= (float)num_equivs(he, ke, le, sym);
+			val = get_redundancy(refl);
+			val /= (float)num_equivs(h, k, l, sym);
 			break;
 		case WGHT_RAWCOUNTS :
-			val = lookup_count(counts, he, ke, le);
+			val = get_redundancy(refl);
 			break;
 		default :
 			ERROR("Invalid weighting.\n");
@@ -228,8 +228,7 @@ static void render_overlined_indices(cairo_t *dctx,
 }
 
 
-static void render_za(UnitCell *cell, ReflItemList *items,
-                      double *ref, unsigned int *counts,
+static void render_za(UnitCell *cell, RefList *list,
                       double boost, const char *sym, int wght, int colscale,
                       signed int xh, signed int xk, signed int xl,
                       signed int yh, signed int yk, signed int yl,
@@ -284,7 +283,7 @@ static void render_za(UnitCell *cell, ReflItemList *items,
 
 	scale = 1.0;
 	draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
-	             ref, counts, items, sym, NULL, wght, boost, colscale, cell,
+	             list, sym, NULL, wght, boost, colscale, cell,
 	             0.0, theta, as, bs, 0.0, 0.0, scale,
 	             &max_ux, &max_uy, &max_val, &max_u, &max_v, &max_res);
 
@@ -342,7 +341,7 @@ static void render_za(UnitCell *cell, ReflItemList *items,
 	cy = 512.0 - 20.0;
 
 	draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
-	             ref, counts, items, sym, dctx, wght, boost, colscale, cell,
+	             list, sym, dctx, wght, boost, colscale, cell,
 	             max_r, theta, as, bs, cx, cy, scale,
 	             NULL, NULL, &max_val, NULL, NULL, NULL);
 
@@ -462,8 +461,8 @@ int main(int argc, char *argv[])
 {
 	int c;
 	UnitCell *cell;
+	RefList *list;
 	char *infile;
-	double *ref;
 	int config_povray = 0;
 	int config_zoneaxis = 0;
 	int config_sqrt = 0;
@@ -477,7 +476,6 @@ int main(int argc, char *argv[])
 	int wght;
 	int colscale;
 	char *cscale = NULL;
-	unsigned int *cts;
 	signed int dh=1, dk=0, dl=0;
 	signed int rh=0, rk=1, rl=0;
 	char *down = NULL;
@@ -613,7 +611,8 @@ int main(int argc, char *argv[])
 	if ( config_zoneaxis ) {
 		if ( (( down == NULL ) && ( right != NULL ))
 		  || (( down != NULL ) && ( right == NULL )) ) {
-			ERROR("Either specify both 'down' and 'right', or neither.\n");
+			ERROR("Either specify both 'down' and 'right',"
+			      " or neither.\n");
 			return 1;
 		}
 		if ( down != NULL ) {
@@ -641,24 +640,22 @@ int main(int argc, char *argv[])
 		ERROR("Couldn't load unit cell from %s\n", pdb);
 		return 1;
 	}
-	ref = new_list_intensity();
-	cts = new_list_count();
-	ReflItemList *items = read_reflections(infile, ref, NULL, cts, NULL);
-	if ( items == NULL ) {
-		ERROR("Couldn't open file '%s'\n", infile);
+	list = read_reflections(infile);
+	if ( list == NULL ) {
+		ERROR("Couldn't read file '%s'\n", infile);
 		return 1;
 	}
-	if ( check_symmetry(items, sym) ) {
+	if ( check_list_symmetry(list, sym) ) {
 		ERROR("The input reflection list does not appear to"
 		      " have symmetry %s\n", sym);
 		return 1;
 	}
 
 	if ( config_povray ) {
-		r = povray_render_animation(cell, ref, cts, items,
+		r = povray_render_animation(cell, list,
 		                            nproc, sym, wght, boost);
 	} else if ( config_zoneaxis ) {
-		render_za(cell, items, ref, cts, boost, sym, wght, colscale,
+		render_za(cell, list, boost, sym, wght, colscale,
 		          rh, rk, rl, dh, dk, dl, outfile);
 	} else {
 		ERROR("Try again with either --povray or --zone-axis.\n");
@@ -666,7 +663,7 @@ int main(int argc, char *argv[])
 
 	free(pdb);
 	free(sym);
-	delete_items(items);
+	reflist_free(list);
 	if ( outfile != NULL ) free(outfile);
 
 	return r;
diff --git a/src/statistics.c b/src/statistics.c
index d0d3ae85..77e7919f 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -25,9 +25,8 @@
 
 
 struct r_params {
-	const double *ref1;
-	const double *ref2;
-	ReflItemList *items;  /* Which reflections to use */
+	RefList *list1;
+	RefList *list2;
 	int fom;              /* Which FoM to use (see the enum just below) */
 };
 
@@ -43,27 +42,30 @@ enum {
 
 
 /* Return the least squares optimal scaling factor when comparing intensities.
- * ref1,ref2 are the two intensity lists to compare.  "items" is a ReflItemList
+ * list1,list2 are the two intensity lists to compare.  "items" is a ReflItemList
  * containing the reflections which should be taken into account.
  */
-double stat_scale_intensity(const double *ref1, const double *ref2,
-                            ReflItemList *items)
+double stat_scale_intensity(RefList *list1, RefList *list2)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
 		double i1, i2;
-		struct refl_item *it;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = lookup_intensity(ref2, h, k, l);
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
 		top += i1 * i2;
 		bot += i2 * i2;
@@ -77,30 +79,36 @@ double stat_scale_intensity(const double *ref1, const double *ref2,
 /* Return the least squares optimal scaling factor when comparing the square
  * roots of the intensities (i.e. one approximation to the structure factor
  * moduli).
- * ref1,ref2 are the two intensity lists to compare (they contain intensities,
+ * list1,list2 are the two intensity lists to compare (they contain intensities,
  * not square rooted intensities).  "items" is a ReflItemList containing the
  * reflections which should be taken into account.
  */
-static double stat_scale_sqrti(const double *ref1, const double *ref2,
-                               ReflItemList *items)
+static double stat_scale_sqrti(RefList *list1, RefList *list2)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, i2, f1, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
+
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
 		if ( i1 < 0.0 ) continue;
 		f1 = sqrt(i1);
-		i2 = lookup_intensity(ref2, h, k, l);
+
 		if ( i2 < 0.0 ) continue;
 		f2 = sqrt(i2);
 
@@ -113,26 +121,33 @@ static double stat_scale_sqrti(const double *ref1, const double *ref2,
 }
 
 
-static double internal_r1_ignorenegs(const double *ref1, const double *ref2,
-                                     ReflItemList *items, double scale)
+static double internal_r1_ignorenegs(RefList *list1, RefList *list2,
+                                     double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, f1, i2, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
+
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
 		if ( i1 < 0.0 ) continue;
 		f1 = sqrt(i1);
-		i2 = lookup_intensity(ref2, h, k, l);
+
 		if ( i2 < 0.0 ) continue;
 		f2 = sqrt(i2);
 		f2 *= scale;
@@ -146,25 +161,32 @@ static double internal_r1_ignorenegs(const double *ref1, const double *ref2,
 }
 
 
-static double internal_r1_negstozero(const double *ref1, const double *ref2,
-                                     ReflItemList *items, double scale)
+static double internal_r1_negstozero(RefList *list1, RefList *list2,
+                                     double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, f1, i2, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
+
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
 		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
-		i2 = lookup_intensity(ref2, h, k, l);
+
 		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
 		f2 *= scale;
 
@@ -177,24 +199,27 @@ static double internal_r1_negstozero(const double *ref1, const double *ref2,
 }
 
 
-static double internal_r2(const double *ref1, const double *ref2,
-                          ReflItemList *items, double scale)
+static double internal_r2(RefList *list1, RefList *list2, double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
 		double i1, i2;
-		struct refl_item *it;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = scale * lookup_intensity(ref2, h, k, l);
+		i1 = get_intensity(refl1);
+		i2 = scale * get_intensity(refl2);
 
 		top += pow(i1 - i2, 2.0);
 		bot += pow(i1, 2.0);
@@ -205,24 +230,27 @@ static double internal_r2(const double *ref1, const double *ref2,
 }
 
 
-static double internal_r_i(const double *ref1, const double *ref2,
-                          ReflItemList *items, double scale)
+static double internal_r_i(RefList *list1, RefList *list2, double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
 		double i1, i2;
-		struct refl_item *it;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = scale * lookup_intensity(ref2, h, k, l);
+		i1 = get_intensity(refl1);
+		i2 = scale * get_intensity(refl2);
 
 		top += fabs(i1-i2);
 		bot += fabs(i1);
@@ -233,24 +261,29 @@ static double internal_r_i(const double *ref1, const double *ref2,
 }
 
 
-static double internal_rdiff_intensity(const double *ref1, const double *ref2,
-                                       ReflItemList *items, double scale)
+static double internal_rdiff_intensity(RefList *list1, RefList *list2,
+                                       double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
 		double i1, i2;
-		struct refl_item *it;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
+
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = lookup_intensity(ref2, h, k, l);
 		i2 *= scale;
 
 		top += fabs(i1 - i2);
@@ -262,25 +295,32 @@ static double internal_rdiff_intensity(const double *ref1, const double *ref2,
 }
 
 
-static double internal_rdiff_negstozero(const double *ref1, const double *ref2,
-                                        ReflItemList *items, double scale)
+static double internal_rdiff_negstozero(RefList *list1, RefList *list2,
+                                        double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, i2, f1, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
+
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
 		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
-		i2 = lookup_intensity(ref2, h, k, l);
+
 		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
 		f2 *= scale;
 
@@ -293,26 +333,33 @@ static double internal_rdiff_negstozero(const double *ref1, const double *ref2,
 }
 
 
-static double internal_rdiff_ignorenegs(const double *ref1, const double *ref2,
-                                        ReflItemList *items, double scale)
+static double internal_rdiff_ignorenegs(RefList *list1, RefList *list2,
+                                        double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
-	int i;
+	Reflection *refl1;
+	RefListIterator *iter;
 
-	for ( i=0; i<num_items(items); i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, i2, f1, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
+
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
 		if ( i1 < 0.0 ) continue;
 		f1 = sqrt(i1);
-		i2 = lookup_intensity(ref2, h, k, l);
+
 		if ( i2 < 0.0 ) continue;
 		f2 = sqrt(i2);
 		f2 *= scale;
@@ -332,26 +379,21 @@ static double calc_r(double scale, void *params)
 
 	switch ( rp->fom) {
 	case R_1_ZERO :
-		return internal_r1_negstozero(rp->ref1, rp->ref2,
-		                              rp->items, scale);
+		return internal_r1_negstozero(rp->list1, rp->list2, scale);
 	case R_1_IGNORE :
-		return internal_r1_ignorenegs(rp->ref1, rp->ref2,
-		                              rp->items, scale);
+		return internal_r1_ignorenegs(rp->list1, rp->list2, scale);
 	case R_2 :
-		return internal_r2(rp->ref1, rp->ref2, rp->items, scale);
+		return internal_r2(rp->list1, rp->list2, scale);
 
 	case R_1_I :
-		return internal_r_i(rp->ref1, rp->ref2, rp->items, scale);
+		return internal_r_i(rp->list1, rp->list2, scale);
 
 	case R_DIFF_ZERO :
-		return internal_rdiff_negstozero(rp->ref1, rp->ref2,
-		                                 rp->items, scale);
+		return internal_rdiff_negstozero(rp->list1, rp->list2,scale);
 	case R_DIFF_IGNORE :
-		return internal_rdiff_ignorenegs(rp->ref1, rp->ref2,
-		                                 rp->items, scale);
+		return internal_rdiff_ignorenegs(rp->list1, rp->list2, scale);
 	case R_DIFF_INTENSITY :
-		return internal_rdiff_intensity(rp->ref1, rp->ref2,
-		                                rp->items, scale);
+		return internal_rdiff_intensity(rp->list1, rp->list2, scale);
 	}
 
 	ERROR("No such FoM!\n");
@@ -359,8 +401,8 @@ static double calc_r(double scale, void *params)
 }
 
 
-static double r_minimised(const double *ref1, const double *ref2,
-                          ReflItemList *items, double *scalep, int fom)
+static double r_minimised(RefList *list1, RefList *list2,
+                          double *scalep, int fom)
 {
 	gsl_function F;
 	gsl_min_fminimizer *s;
@@ -369,9 +411,8 @@ static double r_minimised(const double *ref1, const double *ref2,
 	struct r_params rp;
 	int iter = 0;
 
-	rp.ref1 = ref1;
-	rp.ref2 = ref2;
-	rp.items = items;
+	rp.list1 = list1;
+	rp.list2 = list2;
 	rp.fom = fom;
 
 	F.function = &calc_r;
@@ -385,12 +426,12 @@ static double r_minimised(const double *ref1, const double *ref2,
 	case R_1_IGNORE :
 	case R_DIFF_ZERO :
 	case R_DIFF_IGNORE :
-		scale = stat_scale_sqrti(ref1, ref2, items);
+		scale = stat_scale_sqrti(list1, list2);
 		break;
 	case R_2 :
 	case R_1_I :
 	case R_DIFF_INTENSITY :
-		scale = stat_scale_intensity(ref1, ref2, items);
+		scale = stat_scale_intensity(list1, list2);
 		break;
 	}
 	//STATUS("Initial scale factor estimate: %5.2e\n", scale);
@@ -429,77 +470,74 @@ static double r_minimised(const double *ref1, const double *ref2,
 }
 
 
-double stat_r1_ignore(const double *ref1, const double *ref2,
-                      ReflItemList *items, double *scalep)
+double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_1_IGNORE);
+	return r_minimised(list1, list2, scalep, R_1_IGNORE);
 }
 
 
-double stat_r1_zero(const double *ref1, const double *ref2,
-                    ReflItemList *items, double *scalep)
+double stat_r1_zero(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_1_ZERO);
+	return r_minimised(list1, list2, scalep, R_1_ZERO);
 }
 
 
-double stat_r2(const double *ref1, const double *ref2,
-               ReflItemList *items, double *scalep)
+double stat_r2(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_2);
+	return r_minimised(list1, list2, scalep, R_2);
 }
 
 
-double stat_r1_i(const double *ref1, const double *ref2,
-                 ReflItemList *items, double *scalep)
+double stat_r1_i(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_1_I);
+	return r_minimised(list1, list2, scalep, R_1_I);
 }
 
 
-double stat_rdiff_zero(const double *ref1, const double *ref2,
-                       ReflItemList *items, double *scalep)
+double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_DIFF_ZERO);
+	return r_minimised(list1, list2, scalep, R_DIFF_ZERO);
 }
 
 
-double stat_rdiff_ignore(const double *ref1, const double *ref2,
-                         ReflItemList *items, double *scalep)
+double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_DIFF_IGNORE);
+	return r_minimised(list1, list2, scalep, R_DIFF_IGNORE);
 }
 
 
-double stat_rdiff_intensity(const double *ref1, const double *ref2,
-                            ReflItemList *items, double *scalep)
+double stat_rdiff_intensity(RefList *list1, RefList *list2, double *scalep)
 {
-	return r_minimised(ref1, ref2, items, scalep, R_DIFF_INTENSITY);
+	return r_minimised(list1, list2,  scalep, R_DIFF_INTENSITY);
 }
 
-double stat_pearson_i(const double *ref1, const double *ref2,
-                      ReflItemList *items)
+
+double stat_pearson_i(RefList *list1, RefList *list2)
 {
 	double *vec1, *vec2;
-	int i = 0;
-	int ni = num_items(items);
+	int ni = num_reflections(list1) + num_reflections(list2);
 	double val;
 	int nacc = 0;
+	Reflection *refl1;
+	RefListIterator *iter;
 
 	vec1 = malloc(ni*sizeof(double));
 	vec2 = malloc(ni*sizeof(double));
 
-	for ( i=0; i<ni; i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
 		double i1, i2;
-		struct refl_item *it;
 		signed int h, k, l;
+		Reflection *refl2;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = lookup_intensity(ref2, h, k, l);
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
 		vec1[nacc] = i1;
 		vec2[nacc] = i2;
@@ -514,40 +552,42 @@ double stat_pearson_i(const double *ref1, const double *ref2,
 }
 
 
-double stat_pearson_f_ignore(const double *ref1, const double *ref2,
-                             ReflItemList *items)
+double stat_pearson_f_ignore(RefList *list1, RefList *list2)
 {
 	double *vec1, *vec2;
-	int i = 0;
-	int ni = num_items(items);
+	int ni = num_reflections(list1) + num_reflections(list2);
 	double val;
 	int nacc = 0;
+	Reflection *refl1;
+	RefListIterator *iter;
 
 	vec1 = malloc(ni*sizeof(double));
 	vec2 = malloc(ni*sizeof(double));
 
-	for ( i=0; i<ni; i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, i2, f1, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
 		int ig = 0;
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
+
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		if ( i1 < 0.0 ) {
-			ig=1;
-		}
+		if ( i1 < 0.0 ) ig = 1;
 		f1 = sqrt(i1);
-		i2 = lookup_intensity(ref2, h, k, l);
-		if ( i2 < 0.0 ) {
-			ig=1;
-		}
+
+		if ( i2 < 0.0 ) ig = 1;
 		f2 = sqrt(i2);
 
-		if ( ig  ) continue;
+		if ( ig ) continue;
 
 		vec1[nacc] = f1;
 		vec2[nacc] = f2;
@@ -563,30 +603,35 @@ double stat_pearson_f_ignore(const double *ref1, const double *ref2,
 }
 
 
-double stat_pearson_f_zero(const double *ref1, const double *ref2,
-                      ReflItemList *items)
+double stat_pearson_f_zero(RefList *list1, RefList *list2)
 {
 	double *vec1, *vec2;
-	int i = 0;
-	int ni = num_items(items);
+	int ni = num_reflections(list1) + num_reflections(list2);
 	double val;
 	int nacc = 0;
+	Reflection *refl1;
+	RefListIterator *iter;
 
 	vec1 = malloc(ni*sizeof(double));
 	vec2 = malloc(ni*sizeof(double));
 
-	for ( i=0; i<ni; i++ ) {
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) ) {
 
-		double i1, i2, f1, f2;
-		struct refl_item *it;
+		double i1, i2;
+		double f1, f2;
 		signed int h, k, l;
+		Reflection *refl2;
+
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) continue;  /* No common reflection */
 
-		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 
-		i1 = lookup_intensity(ref1, h, k, l);
 		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
-		i2 = lookup_intensity(ref2, h, k, l);
 		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
 
 		vec1[nacc] = f1;
diff --git a/src/statistics.h b/src/statistics.h
index 18603357..34ef2406 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -17,35 +17,25 @@
 #define STATISTICS_H
 
 
-#include "utils.h"
-
-extern double stat_scale_intensity(const double *ref1, const double *ref2,
-                                   ReflItemList *items);
-
-extern double stat_r1_zero(const double *ref1, const double *ref2,
-                           ReflItemList *items, double *scalep);
-extern double stat_r1_ignore(const double *ref1, const double *ref2,
-                             ReflItemList *items, double *scalep);
-
-extern double stat_r2(const double *ref1, const double *ref2,
-                      ReflItemList *items, double *scalep);
-
-extern double stat_r1_i(const double *ref1, const double *ref2,
-                        ReflItemList *items, double *scalep);
-
-extern double stat_rdiff_zero(const double *ref1, const double *ref2,
-                         ReflItemList *items, double *scalep);
-extern double stat_rdiff_ignore(const double *ref1, const double *ref2,
-                         ReflItemList *items, double *scalep);
-extern double stat_rdiff_intensity(const double *ref1, const double *ref2,
-                         ReflItemList *items, double *scalep);
-
-extern double stat_pearson_i(const double *ref1, const double *ref2,
-                             ReflItemList *items);
-extern double stat_pearson_f_zero(const double *ref1, const double *ref2,
-                                  ReflItemList *items);
-extern double stat_pearson_f_ignore(const double *ref1, const double *ref2,
-                                    ReflItemList *items);
+#include "reflist.h"
+
+extern double stat_scale_intensity(RefList *list1, RefList *list2);
+
+extern double stat_r1_zero(RefList *list1, RefList *list2, double *scalep);
+extern double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep);
+
+extern double stat_r2(RefList *list1, RefList *list2, double *scalep);
+
+extern double stat_r1_i(RefList *list1, RefList *list2, double *scalep);
+
+extern double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep);
+extern double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep);
+extern double stat_rdiff_intensity(RefList *list1, RefList *list2,
+                                   double *scalep);
+
+extern double stat_pearson_i(RefList *list1, RefList *list2);
+extern double stat_pearson_f_zero(RefList *list1, RefList *list2);
+extern double stat_pearson_f_ignore(RefList *list1, RefList *list2);
 
 
 #endif	/* STATISTICS_H */
diff --git a/src/stream.c b/src/stream.c
index 545bed43..97c38b77 100644
--- a/src/stream.c
+++ b/src/stream.c
@@ -31,7 +31,8 @@
 #define PEAK_LIST_START_MARKER "Peaks from peak search"
 #define PEAK_LIST_END_MARKER "End of peak list"
 #define REFLECTION_START_MARKER "Reflections measured after indexing"
-#define REFLECTION_END_MARKER "End of reflections"
+/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also
+ * used to terminate a standalone list of reflections */
 
 static void exclusive(const char *a, const char *b)
 {
@@ -132,49 +133,6 @@ int count_patterns(FILE *fh)
 }
 
 
-static int read_reflections(FILE *fh, struct image *image)
-{
-	char *rval = NULL;
-	int first = 1;
-
-	image->reflections = reflist_new();
-
-	do {
-
-		char line[1024];
-		signed int h, k, l;
-		float intensity, sigma, res, fs, ss;
-		char phs[1024];
-		int cts;
-		int r;
-		Reflection *refl;
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		chomp(line);
-
-		if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return 0;
-
-		r = sscanf(line, "%i %i %i %f %s %f %f %i %f %f",
-		           &h, &k, &l, &intensity, phs, &sigma, &res, &cts,
-		           &fs, &ss);
-		if ( (r != 10) && (!first) ) return 1;
-
-		first = 0;
-		if ( r == 10 ) {
-			refl = add_refl(image->reflections, h, k, l);
-			set_int(refl, intensity);
-			set_detector_pos(refl, fs, ss, 0.0);
-			set_esd_intensity(refl, sigma);
-		}
-
-	} while ( rval != NULL );
-
-	/* Got read error of some kind before finding PEAK_LIST_END_MARKER */
-	return 1;
-}
-
-
 static int read_peaks(FILE *fh, struct image *image)
 {
 	char *rval = NULL;
@@ -409,7 +367,8 @@ int read_chunk(FILE *fh, struct image *image)
 		}
 
 		if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) {
-			if ( read_reflections(fh, image) ) {
+			image->reflections = read_reflections_from_file(fh);
+			if ( image->reflections == NULL ) {
 				ERROR("Failed while reading reflections\n");
 				return 1;
 			}