% GRID_REFINE Refine a grid by adding additional points
%
% This function refines a grid by adding additional points in such a way
% that each interval in the new grid is smaller than the desired value
% delta. All points of the original grid are retained. A typical
% application of this is to create a finer vertical grid for atmospheric
% profiles.
%
% The original grid og must be strictly monotonic (either increasing or
% decreasing, no duplicate values).
%
% See also PROFILE_REFINE, which uses this function.
%
% FORMAT ng = grid_refine(og, delta)
%
% OUT ng New (denser) grid.
% IN og Original grid.
% delta Desired maximum grid spacing.
% 2010-07-02 Created by Stefan Buehler
function ng = grid_refine(og, delta)
% Original grid must have at least two elements.
if length(og) < 2
error( 'Original grid og must have at least two elements.' );
end
% Find out direction of grid and, if necessary, revert it so that it is
% always increasing. (Just to simplify things)
inverted = false;
dir = og(2)-og(1);
if dir < 0
og = og(end:-1:1);
inverted = true;
end
% Check that original grid is monotoneously increasing (even duplicate
% values are not allowed):
d = diff(og);
if min(d) <= 0
error('Original grid og must be strictly monotonic.');
end
% Find out if the original grid is a row or column vector. We will turn
% round the result to match the original grid.
transpose = false;
[rows,cols] = size(og);
if (rows > cols)
transpose = true;
end
% Now the actual work starts...
% Find out how many subdivisions are needed for each interval of the
% original grid.
nsub = ceil(d/delta);
% Calculate subdivision intervals:
sub = d./nsub;
% Create a cell array with the new grid points for each interval in the
% original grid.
ng = arrayfun(@linspace, og(1:end-1), og(2:end)-sub, nsub,...
'UniformOutput', false);
% Get rid of the cell array by concatenating all elements, and add the
% missing last element of the original grid.
ng = [ng{:}, og(end)];
% Turn round the result if the original grid was decreasing:
if inverted==true
ng = ng(end:-1:1);
end
% Transpose the result, if necessary, to match the orientation (row or
% column vector) of the original grid.
if transpose==true
ng = ng';
end