--- a/extra/nurbs/inst/private/onebasisfun__.m	Mon Feb 06 13:48:17 2012 +0000
+++ b/extra/nurbs/inst/private/onebasisfun__.m	Mon Feb 06 14:02:18 2012 +0000
@@ -3,28 +3,32 @@
 %  __ONEBASISFUN__: Undocumented internal function
 %
 %   Copyright (C) 2009 Carlo de Falco
+%   Copyright (C) 2012 Rafael Vazquez
 %   This software comes with ABSOLUTELY NO WARRANTY; see the file
 %   COPYING for details.  This is free software, and you are welcome
 %   to distribute it under the conditions laid out in COPYING.
 
-  N = 0;
-  if (~ any (U <= u)) || (~ any (U > u))
-    return;
-  elseif (p == 0)
-      N = 1;
-      return;
-  end
+  N = zeros (size (u));
 
- ln = u - U(1);
- ld = U(end-1) - U(1);
- if (ld ~= 0)
-   N = N + ln * onebasisfun__ (u, p-1, U(1:end-1))/ ld; 
- end
+  for ii = 1:numel (u)
+    if (~ any (U <= u(ii))) || (~ any (U > u(ii)))
+      continue;
+    elseif (p == 0)
+      N(ii) = 1;
+      continue;
+    else
+      ln = u(ii) - U(1);
+      ld = U(end-1) - U(1);
+      if (ld ~= 0)
+        N(ii) = N(ii) + ln * onebasisfun__ (u(ii), p-1, U(1:end-1))/ ld; 
+      end
 
- dn = U(end) - u;
- dd = U(end) - U(2);
- if (dd ~= 0)
-   N = N + dn * onebasisfun__ (u, p-1, U(2:end))/ dd;
- end
+      dn = U(end) - u(ii);
+      dd = U(end) - U(2);
+      if (dd ~= 0)
+        N(ii) = N(ii) + dn * onebasisfun__ (u(ii), p-1, U(2:end))/ dd;
+      end
+    end
+  end
   
 end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/nurbs/inst/private/onebasisfunder__.m	Mon Feb 06 14:02:18 2012 +0000
@@ -0,0 +1,31 @@
+function N = onebasisfunder__ (u, p, U)
+
+%  __ONEBASISFUNDER__: Undocumented internal function
+%
+%   Copyright (C) 2012 Rafael Vazquez
+%   This software comes with ABSOLUTELY NO WARRANTY; see the file
+%   COPYING for details.  This is free software, and you are welcome
+%   to distribute it under the conditions laid out in COPYING.
+
+  N = zeros (size (u));
+
+  for ii = 1:numel (u)
+    if (~ any (U <= u(ii))) || (~ any (U > u(ii)))
+      continue;
+    elseif (p == 0)
+      N(ii) = 0;
+      continue;
+    else
+      ld = U(end-1) - U(1);
+      if (ld ~= 0)
+        N(ii) = N(ii) + p * onebasisfun__ (u(ii), p-1, U(1:end-1))/ ld;
+      end
+
+      dd = U(end) - U(2);
+      if (dd ~= 0)
+        N(ii) = N(ii) - p * onebasisfun__ (u(ii), p-1, U(2:end))/ dd;
+      end
+    end
+  end
+  
+end

--- a/extra/nurbs/inst/tbasisfun.m	Mon Feb 06 13:48:17 2012 +0000
+++ b/extra/nurbs/inst/tbasisfun.m	Mon Feb 06 14:02:18 2012 +0000
@@ -1,12 +1,12 @@
-function N = tbasisfun (u, p, U)
+function [N, Nder] = tbasisfunder (u, p, U)
 %
-% TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.
+% TBASISFUN: Compute a B- or T-Spline basis function, and its derivatives, from its local knot vector.
 %
 % usage:
 %
-% N = tbasisfun (u, p, U)
-% N = tbasisfun ([u; v], [p q], {U, V})
-% N = tbasisfun ([u; v; w], [p q r], {U, V, W})
+% [N, Nder] = tbasisfun (u, p, U)
+% [N, Nder] = tbasisfun ([u; v], [p q], {U, V})
+% [N, Nder] = tbasisfun ([u; v; w], [p q r], {U, V, W})
 % 
 % INPUT:
 %
@@ -19,9 +19,11 @@
 %
 % OUTPUT:
 %
-%  N : basis function evaluated at the given parametric points 
+%  N    : basis function evaluated at the given parametric points 
+%  Nder : basis function gradient evaluated at the given parametric points 
 %
 %    Copyright (C) 2009 Carlo de Falco
+%    Copyright (C) 2012 Rafael Vazquez
 %
 %    This program is free software: you can redistribute it and/or modify
 %    it under the terms of the GNU General Public License as published by
@@ -40,44 +42,48 @@
     U = sort (U);
     assert (numel (U) == p+2)
     
-    N = zeros(1,numel(u));
-    for ii=1:numel(u)
-      N(ii) = onebasisfun__ (u(ii), p, U);
+    N = onebasisfun__ (u, p, U);
+
+    if (nargout == 2)
+      Nder = onebasisfunder__ (u, p, U);
     end
-
-  elseif size(U,2) == 2
+    
+  elseif (size(U,2) == 2)
     U{1} = sort(U{1}); U{2} = sort(U{2});
     assert (numel(U{1}) == p(1)+2 && numel(U{2}) == p(2)+2)
     
-    Nu = zeros(1,numel(u(1,:))); Nv = zeros(1,numel(u(1,:)));
-    for ii=1:numel(u(1,:))
-      Nu(ii) = onebasisfun__ (u(1,ii), p(1), U{1});
-    end
-
-    for ii=1:numel(u(1,:))
-      Nv(ii) = onebasisfun__ (u(2,ii), p(2), U{2});
-    end
+    Nu = onebasisfun__ (u(1,:), p(1), U{1});
+    Nv = onebasisfun__ (u(2,:), p(2), U{2});
 
     N = Nu.*Nv;
 
-  elseif size(U,2) == 3
+    if (nargout == 2)
+      Ndu = onebasisfunder__ (u(1,:), p(1), U{1});
+      Ndv = onebasisfunder__ (u(2,:), p(2), U{2});
+      
+      Nder(1,:) = Ndu.*Nv;
+      Nder(2,:) = Nu.*Ndv;
+    end
+    
+  elseif (size(U,2) == 3)
     U{1} = sort(U{1}); U{2} = sort(U{2}); U{3} = sort(U{3});
     assert (numel(U{1}) == p(1)+2 && numel(U{2}) == p(2)+2 && numel(U{3}) == p(3)+2)
     
-    Nu = zeros(1,numel(u(1,:))); Nv = zeros(1,numel(u(1,:))); Nw = zeros(1,numel(u(1,:)));
-    for ii=1:numel(u(1,:))
-      Nu(ii) = onebasisfun__ (u(1,ii), p(1), U{1});
-    end
-
-    for ii=1:numel(u(1,:))
-      Nv(ii) = onebasisfun__ (u(2,ii), p(2), U{2});
-    end
-
-    for ii=1:numel(u(1,:))
-      Nw(ii) = onebasisfun__ (u(3,ii), p(3), U{3});
-    end
+    Nu = onebasisfun__ (u(1,:), p(1), U{1});
+    Nv = onebasisfun__ (u(2,:), p(2), U{2});
+    Nw = onebasisfun__ (u(3,:), p(3), U{3});
 
     N = Nu.*Nv.*Nw;
+
+    if (nargout == 2)
+      Ndu = onebasisfunder__ (u(1,:), p(1), U{1});
+      Ndv = onebasisfunder__ (u(2,:), p(2), U{2});
+      Ndw = onebasisfunder__ (u(3,:), p(3), U{3});
+      
+      Nder(1,:) = Ndu.*Nv.*Nw;
+      Nder(2,:) = Nu.*Ndv.*Nw;
+      Nder(3,:) = Nu.*Nv.*Ndw;
+    end
   end
 
 end