function out=pointonsurface(p,q,U,V,CP,t,s)

% POINTONSURFACE computes a point on a nonuniform B-spline surface
% of degree (p,q) at given parameter values (t,s). The surface is
% defined by the control net CP and the sequences of knots U and V.

if (nargin ~= 7)
  error('Incorrect number of arguments');
end

m1 = size(U);
m2 = size(V);

n1 = size(CP,2);
n2 = size(CP,3);

if (m1 ~= (n1 + p + 1))
  error('Inconsistency between number of knots, control points and surface degree');
end

if (m2 ~= (n2 + q + 1))
  error('Inconsistency between number of knots, control points and surface degree');
end

% Find knot intervals where parameters t and s are located
i = findinterval(n1,p,t,U);
j = findinterval(n2,q,s,V);

% Evaluate basis functions at t and s
Nt = basisfunc(i,p,t,U);
Ns = basisfunc(j,q,s,V);
tind = i - p - 1;
pt1 = zeros(size(CP,1),1);
for l = 1 : (q + 1)
  pt2 = zeros(size(CP,1),1);
  sind = j - q - 1 + l;
  for k = 1 : (p + 1)
    pt2 = pt2 + (Nt(k) * CP(:,tind + k,sind)); 
  end
  pt1 = pt1 + (Ns(l) * pt2);
end

out = pt1;

return