1 files changed, 110 insertions, 0 deletions
diff --git a/rspl/sm2.c b/rspl/sm2.c
new file mode 100644
index 0000000..84bc8cf
--- /dev/null
+++ b/rspl/sm2.c
@@ -0,0 +1,110 @@
+
+/* Test smoothness scaling behaviour for 2D */
+
+#include <stdio.h>
+#include <math.h>
+
+double trans(double *v, int luord, double vv);
+
+double f(double x, double y) {
+	double fp[2][5] = {
+		{ 1.0, 0.7, -0.3, 0.0, 0.0 },
+		{ 1.0, 0.7, -0.3, 0.0, 0.0 }
+	};
+	double v;
+
+#ifdef NEVER
+	/* 1D function */
+	v = trans(fp[0], 5, x);
+#else
+	/* 1D on angle */
+	v = trans(fp[0], 5, 0.5 * (x+y));
+	v *= 2.0 * sqrt(2.0);		// ?????
+#endif
+
+//	v = trans(fp[0], 5, x)
+//	  + trans(fp[0], 5, y);
+
+	return v;
+}
+
+int main() {
+	double min = 0.0;
+	double max = 1.0;
+	int di = 2;
+	int i, res = 4;
+
+#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, k;
+
+		/* For each grid point with neigbors */
+		for (j = 1; j < (res-1); j++) {
+			for (k = 1; k < (res-1); k++) {
+				double err;
+				double y1, y2, y3;
+	
+				y1 = f(LOC(j-1), LOC(k));
+				y2 = f(LOC(j+0), LOC(k));
+				y3 = f(LOC(j+1), LOC(k));
+				err = 0.5 * (y3 + y1) - y2;
+				tse += err * err;
+//				tse += fabs(err);
+	
+				y1 = f(LOC(j), LOC(k-1));
+				y2 = f(LOC(j), LOC(k+0));
+				y3 = f(LOC(j), LOC(k+1));
+				err = 0.5 * (y3 + y1) - y2;
+				tse += err * err;
+//				tse += fabs(err);
+			}
+		}
+		/* Apply adjustments and corrections */
+		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;
+}