Index: src/de-jong.c
===================================================================
--- src/de-jong.c	(revision 11229)
+++ src/de-jong.c	(working copy)
@@ -26,6 +26,7 @@
 #include "math-util.h"
 #include <stdlib.h>
 #include <math.h>
+#include <assert.h>
 
 static void de_jong_class_init(DeJongClass *klass);
 static void de_jong_init(DeJong *self);
@@ -610,6 +611,57 @@
 /********************************************************************** Calculation */
 /************************************************************************************/
 
+/*
+ * Generate histogram data from a flat text file. Each line is of
+ * the form 'x y w', where 'x' and 'y' are integer coordinates, and
+ * 'w' is a weight.
+ */
+static inline void histogram_from_file(double *x, double *y) {
+    static const double scale = 0.01;
+    static FILE *f;
+    static struct {
+	int x, y, w;
+    } *points;
+    static int num_points;
+    static int current_point;
+    static int remaining;
+
+    if (!f) {
+	char buffer[256];
+
+	f = fopen("flights.hist", "r");
+	assert(f);
+
+	while (fgets(buffer, sizeof buffer, f)) {
+	    num_points++;
+	}
+	rewind(f);
+
+	points = malloc(num_points * sizeof *points);
+	while (fgets(buffer, sizeof buffer, f)) {
+	    sscanf(buffer, "%d %d %d",
+		   &points[current_point].x,
+		   &points[current_point].y, 
+		   &points[current_point].w);
+	    current_point++;
+	}
+	assert(current_point == num_points);
+    }
+
+    if (remaining == 0) {
+	current_point++;
+	if (current_point >= num_points) {
+	    current_point = 0;
+	}
+	remaining = points[current_point].w;
+    }
+
+    remaining--;
+    *x = scale * (points[current_point].x); // + uniform_variate() - 0.5);
+    *y = scale * (points[current_point].y); // + uniform_variate() - 0.5);
+}
+
+
 void de_jong_calculate(IterativeMap *map, guint iterations) {
     /* Copy frequently used parameters to local variables */
     DeJong *self = DE_JONG(map);
@@ -748,9 +800,12 @@
 	/* These are the actual Peter de Jong map equations. The new point value
 	 * gets stored into 'point', then we go on and mess with x and y before plotting.
 	 */
-	x = sin(param.a * point_y) - cos(param.b * point_x);
-	y = sin(param.c * point_x) - cos(param.d * point_y);
+	histogram_from_file(&x, &y);
 	/*
+	  x = sin(param.a * point_y) - cos(param.b * point_x);
+	  y = sin(param.c * point_x) - cos(param.d * point_y);
+	*/
+	/*
 	  x = sin(param.a * point_y) + param.c * cos(param.a * point_x);
 	  y = sin(param.b * point_x) + param.d * cos(param.b * point_y);
 	*/