From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= <debian@jff-webhosting.net>
Date: Mon, 1 Sep 2014 13:56:46 +0200
Subject: Initial import of argyll version 1.5.1-8

---
 rspl/sm1.c | 88 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 88 insertions(+)
 create mode 100644 rspl/sm1.c

(limited to 'rspl/sm1.c')

diff --git a/rspl/sm1.c b/rspl/sm1.c
new file mode 100644
index 0000000..5bcf26e
--- /dev/null
+++ b/rspl/sm1.c
@@ -0,0 +1,88 @@
+
+/* Test smoothness scaling behaviour for 1D */
+
+#include <stdio.h>
+#include <math.h>
+
+double trans(double *v, int luord, double vv);
+
+double f(double x) {
+	double fp[5] = { +1.0, 0.7, -0.3, 0.0, 0.0 };
+	double y;
+
+	y = trans(fp, 5, x);
+
+	return y;
+}
+
+int main() {
+	double min = 0.0;
+	double max = 1.0;
+	int i, res = 4;
+	int di = 1;
+
+#define LOC(xx) (min + (max-min) * (xx)/(res-1.0))
+
+	/* For each resolution */
+	for (i = 0; i < 10; i++, res *= 2) {
+		double tse = 0.0;		/* Total squared error */
+		int j;
+
+		/* For each grid point with neigbors */
+		for (j = 1; j < (res-1); j++) {
+			double err;
+			double y1, y2, y3;
+
+			y1 = f(LOC(j-1));
+			y2 = f(LOC(j));
+			y3 = f(LOC(j+1));
+			err = 0.5 * (y3 + y1) - y2;
+			tse += err * err;
+//			tse += fabs(err);
+		}
+		/* Apply adjustments and corrections to error squared */
+		tse *= pow((res-1.0), 4.0);					/* Aprox. geometric resolution factor */
+		tse /= pow((res-2.0),(double)di);			/* Average squared non-smoothness */
+
+//		tse /= (di * pow((res-2.0),(double)di));	/* Average squared non-smoothness */
+		printf("Res %d, tse = %f\n",res,tse);
+	}
+
+
+	return 0;
+}
+
+/* Transfer function */
+double trans(
+double *v,			/* Pointer to first parameter */
+int    luord,		/* Number of parameters */
+double vv			/* Source of value */
+) {
+	double g;
+	int ord;
+
+	for (ord = 0; ord < luord; ord++) {
+		int nsec;			/* Number of sections */
+		double sec;			/* Section */
+
+		g = v[ord];			/* Parameter */
+
+		nsec = ord + 1;		/* Increase sections for each order */
+
+		vv *= (double)nsec;
+
+		sec = floor(vv);
+		if (((int)sec) & 1)
+			g = -g;				/* Alternate action in each section */
+		vv -= sec;
+		if (g >= 0.0) {
+			vv = vv/(g - g * vv + 1.0);
+		} else {
+			vv = (vv - g * vv)/(1.0 - g * vv);
+		}
+		vv += sec;
+		vv /= (double)nsec;
+	}
+
+	return vv;
+} 
-- 
cgit v1.2.3