cubsum1.m

改进等高线拟合代码

function [a,b,c,d,y2,c_flag]=cubsum1(n,x,y,p)

% [a,b,c,d,y2,c_flag]=cubsum1(n,x,y,p)
% Smoothing cubic spline,
% x : x-coordinates monotonic
% y : y-coordinates
% p : weigthing parameters for each point
% a,b,c,d parameters of cubic spline 
% y2 : second derivatives

% dimensions of a,b,c,d,y2, x and y is n
% c_flag = 0  normal excecution
%        = 1  n<5 not allowed
%        = 2  Error in solving system with pentas

 c_flag = 0;
 if n<5,
   c_flag = 1;
   return;
 end
 n1=n-1;
 n2=n-2;
 for k=1:n1;
   d(k)=1.0/(x(k+1)-x(k));
 end
 d(n)=0.0;
 p1=1.0/p(1);
 p2=1.0/p(2);
 h1=d(1);
 h2=d(2);
 r1=(y(2)-y(1))*h1;
 for k=1:n2,
   k1=k+1;
   k2=k+2;
   h3=d(k2);
   p3=1.0/p(k2);
   s=h1+h2;
   a(k)=2.0/h1+2.0/h2+6.0*(h1*h1*p1+s*s*p2+h2*h2*p3);
   b(k)=1.0/h2-6.0*h2*(p2*s+p3*(h2+h3));
   c(k)=6.0*p3*h2*h3;
   r2=(y(k2)-y(k1))*h2;
   y2(k)=6.0*(r2-r1);
   h1=h2;
   h2=h3;
   p1=p2;
   p2=p3;
   r1=r2;
 end
 [f_result,c_flag]=pentas(n2,a,b,c,y2);
 y2=f_result;
 if c_flag ~=0,
    return
 end
 for k=n2:-1:1,
   y2(k+1)=y2(k);
 end
 y2(1)=0.0;
 y2(n)=0.0;
 h1=0.0;
 for k=1:n1,
   k1=k+1;
   b(k)=d(k);
   h2=d(k)*(y2(k1)-y2(k));
   d(k)=h2/6.0;
   a(k)=y(k)-(h2-h1)/p(k);
   c(k)=y2(k)/2.0;
   h1=h2;
 end
 a(n)=y(n)+h1/p(n);
 for k=1:n1,
   k1=k+1;
   h=b(k);
   b(k)=(a(k1)-a(k))*h-(y2(k1)+2.0*y2(k))/(6.0*h);
 end
return