function [o1,o2,o3] = seval(i1,i2,i3,i4,i5,i6,i7,i8);
% function p = seval(a,u,x); (function or curve)
%
% function p = seval(a,rho,u,x); (rational)
%
% function p = seval(a,u,v,x,y); (tensor product function)
%
% function p = seval(a,rho,u,v,x,y); (rational)
%
% function [p1,p2,p3] = seval(a1,a2,a3,u,v,x,y); (tensor product surface)
%
% function [p1,p2,p3] = seval(a1,a2,a3,rho,u,v,x,y); (rational)
%
% evaluation
%
% input:  a,a1,a2,a3 ... coefficient matrix
%         rho        ... weights
%         u,v        ... knot vector
%         x,y        ... points to be evaluated
%
% output: p,p1,p2,p3 ... values at points x (x,y)
                                                                            
% written by B. Harzer and Kai Holthaus in 1993
%
% changed by Juergen Koch, June 1994

if nargin == 3   % spline function or spline curve
  
  [ma,la] = size(i1); lu = length(i2); lx = length(i3);
  d = lu - la - 1;
  
  o1 = zeros(ma,lx); 			% just to get better performance
  
  for i = 1:lx
  
    xx = i3(i);
    if ((i2(d+1) <= xx) & (xx <= i2(lu-d))) 
      
      a = i1;
      
      % find index l with u_l <= x < u_{l+1} 
      ll = find(xx == i2); m = length(ll);
      if m == 0
	ll = find(xx < i2); l = ll(1) - 1;
      elseif m <= d
	l = ll(m);
      else 				% fix problem with first knot 
	l = ll(m) + (xx == i2(d+1)); 	% having multiplicity > d
      end

      for k = d-1:-1:m
	J1 = [l-k:l-m];
	w = ones(ma,1)*((xx - i2(J1))./(i2(J1+k+1) - i2(J1)));
	a(:,J1) = w.*a(:,J1) + (1-w).*a(:,J1-1);
      end
    
      o1(:,i) = a(:,l-m);
      
    else
      warning('seval: x not in [u_0,u_n] !');
    end
  end

elseif nargin == 4 			% rat. spline function or spline curve

  d = size(i1,1);
  r = seval([(ones(d,1)*i2).*i1; i2],i3,i4); % use seval for rho.*a and rho
  o1 = r(1:d,:)./(ones(d,1)*r(d+1,:));
  
elseif nargin == 5 			% tensor product spline function

  o1 = seval(seval(i1,i3,i5)',i2,i4)';
  
elseif nargin == 6 			% rat. tensor product spline function

  o1 = seval(i1.*i2,i3,i4,i5,i6)./seval(i2,i3,i4,i5,i6);
    
elseif nargin == 7 			% tensor product spline surface
  
  o1 = seval(i1,i4,i5,i6,i7);
  o2 = seval(i2,i4,i5,i6,i7);
  o3 = seval(i3,i4,i5,i6,i7);
  
elseif nargin == 8 			% rat. tensor product spline surface
    
  o1 = seval(i1,i4,i5,i6,i7,i8);
  o2 = seval(i2,i4,i5,i6,i7,i8);
  o3 = seval(i3,i4,i5,i6,i7,i8);
    
else
  disp('seval: nargin <> 3,4,5,6,7,8 !');
end