--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/old/curve_interp.m	Thu Sep 17 10:12:50 2015 +0200
@@ -0,0 +1,58 @@
+% Create a cubic spline from the points (x,y) using periodic conditions.
+%   g = curve_interp(x,y)
+function g = curve_interp(x,y)
+    default_arg('x',[0 2 2 1 1 0])
+    default_arg('y',[0 0 2 2 1 1])
+    % solve for xp and yp
+
+    % x(t) = at^4 + bt^2+ct+d
+
+    % a = xp1 -2x1 + 2x0 +  xp0
+    % b = 3x1 -xp1 - 3x0 + 2xp0
+    % c = xp0
+    % d = x0
+
+    assert(length(x) == length(y))
+    n = length(x);
+    A = spdiags(ones(n,1)*[2, 8, 2],-1:1,n,n);
+    A(n,1) = 2;
+    A(1,n) = 2;
+
+    bx = zeros(n,1);
+    for i = 2:n-1
+        bx(i) = -6*x(i-1)+6*x(i+1);
+    end
+    bx(1) = -6*x(n)+6*x(2);
+    bx(n) = -6*x(n-1)+6*x(1);
+
+    by = zeros(n,1);
+    for i = 2:n-1
+        by(i) = -6*y(i-1)+6*y(i+1);
+    end
+    by(1) = -6*y(n)+6*y(2);
+    by(n) = -6*y(n-1)+6*y(1);
+
+
+    xp = A\bx;
+    yp = A\by;
+
+    x(end+1) = x(1);
+    y(end+1) = y(1);
+
+    xp(end+1) = xp(1);
+    yp(end+1) = yp(1);
+
+    function v = g_fun(t)
+        t = mod(t,1);
+        i = mod(floor(t*n),n) + 1;
+        t = t * n -(i-1);
+        X = (2*x(i)-2*x(i+1)+xp(i)+xp(i+1))*t.^3 + (-3*x(i)+3*x(i+1)-2*xp(i)-xp(i+1))*t.^2 + (xp(i))*t + x(i);
+        Y = (2*y(i)-2*y(i+1)+yp(i)+yp(i+1))*t.^3 + (-3*y(i)+3*y(i+1)-2*yp(i)-yp(i+1))*t.^2 + (yp(i))*t + y(i);
+        v = [X;Y];
+    end
+
+    g = @g_fun;
+end
+
+
+