Mercurial > repos > public > sbplib
comparison +grid/old/curve_interp.m @ 0:48b6fb693025
Initial commit.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 17 Sep 2015 10:12:50 +0200 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:48b6fb693025 |
|---|---|
| 1 % Create a cubic spline from the points (x,y) using periodic conditions. | |
| 2 % g = curve_interp(x,y) | |
| 3 function g = curve_interp(x,y) | |
| 4 default_arg('x',[0 2 2 1 1 0]) | |
| 5 default_arg('y',[0 0 2 2 1 1]) | |
| 6 % solve for xp and yp | |
| 7 | |
| 8 % x(t) = at^4 + bt^2+ct+d | |
| 9 | |
| 10 % a = xp1 -2x1 + 2x0 + xp0 | |
| 11 % b = 3x1 -xp1 - 3x0 + 2xp0 | |
| 12 % c = xp0 | |
| 13 % d = x0 | |
| 14 | |
| 15 assert(length(x) == length(y)) | |
| 16 n = length(x); | |
| 17 A = spdiags(ones(n,1)*[2, 8, 2],-1:1,n,n); | |
| 18 A(n,1) = 2; | |
| 19 A(1,n) = 2; | |
| 20 | |
| 21 bx = zeros(n,1); | |
| 22 for i = 2:n-1 | |
| 23 bx(i) = -6*x(i-1)+6*x(i+1); | |
| 24 end | |
| 25 bx(1) = -6*x(n)+6*x(2); | |
| 26 bx(n) = -6*x(n-1)+6*x(1); | |
| 27 | |
| 28 by = zeros(n,1); | |
| 29 for i = 2:n-1 | |
| 30 by(i) = -6*y(i-1)+6*y(i+1); | |
| 31 end | |
| 32 by(1) = -6*y(n)+6*y(2); | |
| 33 by(n) = -6*y(n-1)+6*y(1); | |
| 34 | |
| 35 | |
| 36 xp = A\bx; | |
| 37 yp = A\by; | |
| 38 | |
| 39 x(end+1) = x(1); | |
| 40 y(end+1) = y(1); | |
| 41 | |
| 42 xp(end+1) = xp(1); | |
| 43 yp(end+1) = yp(1); | |
| 44 | |
| 45 function v = g_fun(t) | |
| 46 t = mod(t,1); | |
| 47 i = mod(floor(t*n),n) + 1; | |
| 48 t = t * n -(i-1); | |
| 49 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); | |
| 50 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); | |
| 51 v = [X;Y]; | |
| 52 end | |
| 53 | |
| 54 g = @g_fun; | |
| 55 end | |
| 56 | |
| 57 | |
| 58 |