--- a/main/splines/inst/csaps.m	Sat Nov 24 22:16:49 2012 +0000
+++ b/main/splines/inst/csaps.m	Sun Nov 25 22:44:06 2012 +0000
@@ -85,10 +85,6 @@
   u = (6*(1-p)*QT*diag(1 ./ w)*QT' + p*R) \ (QT*y);
   a = y - 6*(1-p)*diag(1 ./ w)*QT'*u;
 
-## note: add knots to either end of spline pp-form to ensure linear extrapolation
-  xminus = x(1) - eps(x(1));
-  xplus = x(end) + eps(x(end));
-  x = [xminus; x; xplus];
   
 ## derivatives for the piecewise cubic spline  
   aa = bb = cc = dd = zeros (n+1, m);
@@ -101,13 +97,22 @@
   bb(2:n, :) = diff(a) ./ h - (cc(2:n, :)/2).*h - (dd(2:n, :)/6).*(h.^2);
  unwind_protect_cleanup
  warning (warn_state, "Octave:broadcast");
- end_unwind_protect  
-  bb(1, :) = bb(2, :); #linear extension of splines
-  bb(n + 1, :) = bb(n, :);
+ end_unwind_protect
+
+## note: add knots to either end of spline pp-form to ensure linear extrapolation
+  xminus = x(1) - eps(x(1));
+  xplus = x(end) + eps(x(end));
+  x = [xminus; x; xplus];  
+  slope_minus = bb(2, :);
+  slope_plus = bb(n, :) + cc(n, :)*h(n-1) + (dd(n, :)/2)*h(n-1)^2;
+  bb(1, :) = slope_minus; #linear extension of splines
+  bb(n + 1, :) = slope_plus;
   aa(1, :) = a(1, :) - eps(x(1))*bb(1, :);
   
   ret = mkpp (x, cat (2, dd'(:)/6, cc'(:)/2, bb'(:), aa'(:)), size(y, 2));
 
+
+
   if ~isempty (xi)
     ret = ppval (ret, xi);
   endif
@@ -122,4 +127,10 @@
 %!assert (csaps(x',y',1,x'), y', 10*eps);
 %!assert (csaps(x',y',1,x), y, 10*eps);
 %!assert (csaps(x,[y 2*y],1,x)', [y 2*y], 10*eps);
+%{
+x = ([1:10 10.5 11.3])'; y = sin(x);
+xx = 0:0.1:12;
+[yy, p] = csaps(x,y,1,xx);
+plot(x, y, 's', xx, yy)
 
+%}

--- a/main/splines/inst/csaps_sel.m	Sat Nov 24 22:16:49 2012 +0000
+++ b/main/splines/inst/csaps_sel.m	Sun Nov 25 22:44:06 2012 +0000
@@ -121,23 +121,8 @@
   sigma2 = mean(MSR(:)) * (n / (n-Ht)); #estimated data error variance (wahba83)
   unc_y = sqrt(sigma2 * Hd ./ w); #uncertainty (SD) of fitted curve at each input x-value (hutchinson86)
 
-## solve for the scaled second derivatives u and for the function values a at the knots (if p = 1, a = y) 
-  u = (6*(1-p)*QT*diag(1 ./ w)*QT' + p*R) \ (QT*y);
-  a = y - 6*(1-p)*diag(1 ./ w)*QT'*u; #where y is calculated
-
-## derivatives at all but the last knot for the piecewise cubic spline
-  aa = a(1:(end-1), :);
-  cc = zeros(size(y)); 
-  cc(2:(n-1), :) = 6*p*u; #cc([1 n], :) = 0 [natural spline]
-  dd = diff(cc) ./ h;
-  cc = cc(1:(end-1), :);
-  bb = diff(a) ./ h - (cc/2).*h - (dd/6).*(h.^2);
-
-  ret = mkpp (x, cat (2, dd'(:)/6, cc'(:)/2, bb'(:), aa'(:)), size(y, 2));
-
-  if ~isempty(xi)
-    ret = ppval (ret, xi);
-  endif
+## construct the fitted smoothing spline 
+  [ret,p]=csaps(x,y,p,xi,w);
 
 endfunction