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 arraywarning('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   enddir = -1;X = get(A,'pd');name = get(X,'name');for i = 1:length(name)   if strcmp(name{i},codir), dir = i; end   endif dir == -1   error('Incorrect coordinate direction name')   return   endqp = 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   endpd = normal(X,codir); % Remove the coordinate along which interpolation was performedswitch 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 fieldnewqp = 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')   endif 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);