1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
function v = pchip(x,y,xx)
//this function returns piecewise cubic hermite interpolating polynomial.
// Calling Sequence
// d=pchip(x,y)
// d= pchip(X,,y,xx)
// Parameters
// x: a vector
// y: if Y is vector then it must have the same length as x and Y is matrix then the last dimension of Y must equal length(X).
// xx: Points for interpolation
// v: vector of interpolantant at xx
//// Examples
// x=[0 1 2 3 4 5]
// y=[1 0 1 0 1 0]
// xx=linspace(0,5,800)
// v=pchip(x, y)
// v=pchip(x,y,xx)
//plot(x,y,xx,v,'o')
// Authors
// Jitendra Singh
// NOTE:execute function "pchips" prior executing this function
if argn(2)==3 & ~isreal(xx)
error('Points for interpolation must be real.')
end
nn=size(y,1);
h = diff(x); m = prod(nn);
delta = diff(y,1,2)./repmat(h,m,1);
slopes = zeros(size(y,1),size(y,2));
for r = 1:m
if isreal(delta)
slopes(r,:) = pchips(x,y(r,:),delta(r,:));
else
realslopes = pchips(x,y(r,:),real(delta(r,:)));
imagslopes = pchips(x,y(r,:),imag(delta(r,:)));
slopes(r,:) = complex(realslopes, imagslopes);
end
end
s=slopes;
divdif=delta
d = size(y,1);
dx = diff(x(:).');
dxd = repmat(dx,d,1);
divdif = diff(y,1,2)./dxd;
n = length(x);
dzzdx = (divdif-s(:,1:n-1))./dxd; dzdxdx = (s(:,2:n)-divdif)./dxd;
dnm1 = d*(n-1);
c1=matrix((dzdxdx-dzzdx)./dxd,dnm1,1)
c2=matrix(2*dzzdx-dzdxdx,dnm1,1)
c3=matrix(s(:,1:n-1),dnm1,1)
c4=matrix(y(:,1:n-1),dnm1,1)
v=[c1,c2,c3,c4]
if argn(2)==3
//v = ppval(v,xx);
b=x;
c=v;
l=length(b)-1;
dlk=length(c);
d = size(y,1)
dl=prod(d)*l;
eps=2.2204e-16;
k=fix(dlk/dl+100*eps);
dd=d;
lx = length(xx);
xs = matrix(xx,1,length(xx));
if lx, [cf,idx] = histc([-%inf,b(2:l),%inf], xs);
else idx=ones(1, length(xx));
end
infxs = find(xs==%inf);
if ~isempty(infxs)
index(infxs) = l;
end
nogoodxs = find(idx==0);
if ~isempty (nogoodxs)
xs(nogoodxs) = %nan;
idx(nogoodxs) = 1;
end
xs = xs-b(idx);
d = prod(dd);
sizexx = size(xx)
if d>1
xs = reshape(xs(ones(d,1),:),1,d*lx);
idx = d*idx; temp = (-d:-1).';
idx = reshape(1+idx(ones(d,1),:)+temp(:,ones(1,lx)), d*lx, 1 );
else
if length(sizexx)>1, dd = []; else dd = 1; end
end
v = c(idx,1);
for i=2:k
v = xs(:).*v + c(idx,i);
end
if ~isempty(nogoodxs) & k==1 & l>1
v = matrix(v,d,lx); v(:,nogoodxs) = NaN;
end
v = matrix(v,[dd,sizexx]);
end
endfunction
|
