|
0
|
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
|
