function [GV,RangeFactor] = evalSurf(blk,varargin)
%EVALSURF  Evaluates gains at specific operating points.
%
%   EVALSURF evaluates the gain surface formula over a set or grid of
%   scheduling variable values.
%
%   GV = EVALSURF(GS,X) evaluates the gain surface GS at the points X.
%   X is an N-by-M array where N is the number of points and M is the
%   number of scheduling variables (fields in GS.SamplingGrid). The output
%   GV is an array of gain values of size [N,Ny,Nu] if each gain is of size
%   [Ny,Nu].
%
%   GV = EVALSURF(GS,X1,...,XM) evaluates the gain surface GS over the
%   rectangular grid generated by the vectors X1,...,XM (one per scheduling
%   variable). The result gain array is of size [N1 ... NM Ny Nu] where Nj
%   is the length of Xj (number of points along the j-th grid dimension).
%
%   GV = EVALSURF(GS,...,GRIDFLAG) controls the layout of GV. When
%   GridFlag = 'gridfirst', GV is of size [N1 ... NM Ny Nu] with the grid 
%   dimensions first and the gain dimensions last. This is the natural 
%   format for scalar gains. When GRIDFLAG = 'gridlast', GV is of size 
%   [Ny Nu N1 ... NM] with the gain dimensions first. This format is more
%   readable for matrix-valued gains.
%
%   This function is useful to turn gain surfaces into lookup tables. Set 
%   X1,...,XM to the desired table breakpoints and use GV as table data.
%
%   See also GETDATA, SETDATA, VIEWSURF, tunableSurface.

%   Copyright 1986-2015 The MathWorks, Inc.

% Check inputs
narginchk(2,Inf)
if nargin>2 && ischar(varargin{end})
   GridFirst = strcmpi(varargin{end},'gridfirst');
   varargin = varargin(1:end-1);
else
   GridFirst = true;  % default
end
if ~all(cellfun(@(x) isnumeric(x) && isreal(x),varargin))
   error(message('Control:lftmodel:evalSurf1'))
end
nx = numel(varargin);
m = blk.nVar_;
if nx==1
   % Scattered data or 1D grid
   X = varargin{1};
   mDim = (size(X)==m);
   if ~(ismatrix(X) && any(mDim))
      error(message('Control:lftmodel:evalSurf2'))
   elseif mDim(1) && ~mDim(2)
      X = X.';
   end
   CG = num2cell(X,1);
elseif nx==m
   if ~all(cellfun(@(x) isvector(x),varargin))
      error(message('Control:lftmodel:evalSurf3'))
   end
   % Construct ND grid
   [varargin{1:nx}] = ndgrid(varargin{1:nx});
   CG = varargin;
else
   error(message('Control:lftmodel:evalSurf4'))
end

% Evaluate gains on specified grid
nf = blk.nFun_;
gridSize = size(CG{1});
if nf==0
   % Constant gain
   GV = repmat(getScaledCoefficients(blk),[1 1 gridSize]);
else
   % Warn if some points are far outside the tuning range
   SG = struct2cell(getData(blk.SamplingGrid_));
   RangeFactor = 0;
   for ct=1:m
      % Center and scale with respect to SamplingGrid range
      xData = SG{ct}(:);
      xmin = min(xData);
      xmax = max(xData);
      normRange = max(abs((CG{ct}(:)-(xmin+xmax)/2)));
      if xmin~=xmax
         normRange = 2*normRange/(xmax-xmin);
      end
      RangeFactor = max(RangeFactor,normRange);
   end
   if nargout<2 && RangeFactor>1.25   % 25% slack
      warning(message('Control:lftmodel:evalSurf5'))
   end
   % Compute gain values
   CG = tunableSurface.applyNormalization(CG,blk.Normalization);
   try
      fVals = tunableSurface.evalBF(CG,blk.BasisFunctions_,nf);
   catch ME
      throw(ME)
   end
   G = getScaledCoefficients(blk);
   [ny,nu] = size(blk);
   GV = zeros([ny nu gridSize]);
   for ct1=1:prod(gridSize)
      M = G(:,1:nu);
      jG = nu;
      for ct2=2:nf+1
         M = M + fVals(ct2,ct1) * G(:,jG+1:jG+nu);
         jG = jG+nu;
      end
      GV(:,:,ct1) = M;
   end
end

% Honor GRIDLAST flag
if GridFirst
   GV = permute(GV,[3:nx+2 1 2]);
end
