function newM = Interpolate2D_v(M)
% Perform 2-d interpolation on the real-valued data in nr-by-nc matrix M.
% The interpolated data are added between existing data points so newM is 
%   (2*nr-1)-by-(2*nc-1). Use the simple average as the interpolated value.
%
% VECTORIZED

[nr,nc]= size(M);

% Interpolate one dimension at a time.
  % New dimensions after interpolating in x-direction: nr-by-(nc*2-1)
  wideM= zeros(nr, nc*2-1);
  % Insert a column of interpolated data between existing columns of data
  for c= 1:nc-1
    wideIm(:,c*2-1)= M(:,c);
    wideIm(:,c*2)= (M(:,c)+M(:,c+1))./2;
  end
  wideIm(:,nc*2-1)= M(:,nc);

  % New dimensions after interpolating in y-direction: (nr*2-c)-by-(nc*2-1)
  newM= zeros(nr*2-1, nc*2-1);
  % Insert a row of interpolated data between existing rows of data
  for r= 1:nr-1
      newM(r*2-1,:)= wideIm(r,:);
      newM(r*2,:)= (wideIm(r,:)+wideIm(r+1,:))./2;
  end
  newM(nr*2-1,:)= wideIm(nr,:);



