grad7.m

等高线拟合源代码（matlab源码）

% GRAD7   : Determination of values for the first derivatives        [f]
%
% [B] = GRAD7 (X,Y)    or  [B] = GRAD7 (X,Y,EPSS)
%
% Returns in B the estimates for the first derivatives for a
% spline fit trough the points given in X and Y.
% If the value EPSS is not given this routine uses the matlab
% precision value eps.

% uses a as a workspace.

function [b]=grad7(x,y,epss)

if  nargin ==2,
 epss = eps;
end

n= length(x);

if n <3,
  error('N moet groter zijn dan 2');
end
n1 = n-1;
n2 = n-2;
for k=1:n1,
  k1=k+1;
  a(k)=(y(k1)-y(k))/(x(k1)-x(k));
end

for k=2:n1,
  b(k)=(a(k)+a(k-1))/2.0;
end
b(1)=a(1);
b(n)=a(n1);

for k=1:n2,
  k1=k+1;
  if (abs(a(k))<= epss),
    b(k)=0.0;
    b(k1)=0.0;
  elseif ((a(k)*a(k1))<0.0) 
    b(k1)=0.0;
  end
end

if (abs(a(n1))<=epss),
  b(n1)=0.0;
  b(n)=0.0;
end

for k=1:n1,
  k1=k+1;
  if ((b(k)*a(k)) < 0.0),
    b(k)=0.0
  end
  if((b(k1)*a(k)) < 0.0),
    b(k1) = 0.0;
  end
end
if ((b(n)*a(n1)) < 0.0),
  b(n)=0.0;
end

% return % toegevoedgd, schalering die volgt lijkt niet goed te werken


for k=1:n1,
  pk = 0.0;
  qk = 0.0;
  k1 = k+1;
  if ((abs(b(k))+abs(b(k1)))> epss),
    h=3.0 * a(k);
    rk=(2.0 * b(k) + b(k1))/h;
    sk=(2.0 * b(k1) + b(k))/h;
    if (1.0 - rk)*(sk-1.0) < 0.0,
      pk = min([rk,sk]);
    end
  end
  if (abs(a(k))> eps),
    h=3.0*a(k);
    qk=(b(k)+b(k1))/h - sqrt(b(k)*b(k1))/abs(h);
  end
  a(k) = max([qk pk 1.0]);
end

b(1)=b(1)/a(1);
b(n)=b(n)/a(n1);
for k=1:n2;
  k1=k+1;
  b(k1)=b(k1)/max([a(k) a(k1)]);
end
return