interp.m

MDPSAS工具箱是马里兰大学开发的

function B = interp(A,codir,locat)

% interp.m  Evaluate an scalarfield object at a chosen point in
%           one of the coordinate directions; for example
%
%           B = interp(A,'z',0.2)
%
%           where A and B are scalarfield objects, 'z'
%           is one of the coordinate directions in the
%           physical domain, the 0.2 the numerical value
%           at which interpolation is to be performed.
%
%           NOTE: This is a linear interpolation algorithm

% Use recursively if locat is a double array

warning('Linear interpolation; please switch to linterp')

if prod(size(locat)) > 1
   sizelocat = size(locat);
   for i = 1:prod(sizelocat)
      interparray(i) = interp(A,codir,locat(i));
      end
   B = reshape(interparray,sizelocat);
   return
   end

dir = -1;

X = get(A,'pd');
name = get(X,'name');

for i = 1:length(name)
   if strcmp(name{i},codir), dir = i; end
   end

if dir == -1
   error('Incorrect coordinate direction name')
   return
   end

qp = get(X,'qp');
qpDirVal =locat;
for i = 2:size(qp{dir},1)
   if qp{dir}(i-1) <= qpDirVal & qpDirVal <= qp{dir}(i)
      indexLo = i-1;
      indexHi  = i;
      interpFactorHi = (qpDirVal - qp{dir}(i-1))/(qp{dir}(i) - qp{dir}(i-1)); % =1 for Hi pt
      interpFactorLo = (qp{dir}(i) - qpDirVal)  /(qp{dir}(i) - qp{dir}(i-1)); % =1 for Lo pt
      end
   end

pd = normal(X,codir); % Remove the coordinate along which interpolation was performed

switch dir
   case 1
      val =  interpFactorLo*A.val(indexLo,:,:,:) ...
           + interpFactorHi*A.val(indexHi,:,:,:);
   case 2
      val =  interpFactorLo*A.val(:,indexLo,:,:) ...
           + interpFactorHi*A.val(:,indexHi,:,:);
   case 3
      val =  interpFactorLo*A.val(:,:,indexLo,:) ...
           + interpFactorHi*A.val(:,:,indexHi,:);
   case 4
      val =  interpFactorLo*A.val(:,:,:,indexLo) ...
           + interpFactorHi*A.val(:,:,:,indexHi);
   otherwise
      error('Not set up for higher-dim spaces')
      end

% Remove the direction along which interpolation was performed
% from the scalar field

newqp = get(pd,'qp');

newshape = [length(newqp{1}) 1];

if length(newqp) == 2
   newshape = [length(newqp{1}) length(newqp{2})];
elseif length(newqp) == 3
   newshape = [length(newqp{1}) length(newqp{2}) length(newqp{3})];
elseif length(newqp) == 4
   newshape = [length(newqp{1}) length(newqp{2}) length(newqp{3}) length(newqp{4})];
elseif length(newqp) > 4
   error('Not set up for higher-dim spaces')
   end

if isempty(newqp{1}) % return a scalar, if no pd
   B = val;
else
   val = reshape(val,newshape);
   B = scalarfield(pd,val);
   end

% ------------------------------------------------------ %

function B = ebd(A,dir,qpDirVal)

% ebd.m Evaluate quadgrid object at one of its boundaries
%       Called as:
%                  B = ebd(A,dir,lindex)
%
%     INPUT/OUTPUT:
%      A,B : quadgrid objects
%

geom = get(A,'geom');
qp   = get(A,'qp');
w    = get(A,'w');
d    = get(A,'d');
dd   = get(A,'dd');
Q    = get(A,'Q');
name = get(A,'name');

qp{dir} = qpDirVal;
w{dir}  = 1;

B = quadgrid(geom,qp,w,d,dd,Q,name);