function ANew=RefineRC(AOld,SplitArray,code)
% ANew=RefineRC(AOld,SplitArray,code)
% refines rows or cols of a FDM model array according to SplitArray
% SplitArray has the form [i1 i2 i3 i4; N1 N2 N3 N4] where
% i refers to the old column and N the number by which the column must be
% subdivided.
% code is either 'x', 'y', 'z' 'X' 'Y  or 'Z' to indicated the direction
% along which the SplitArray is to be applied and to indicate whether the
% array to be split is a grid array or a cell array. A grid array concerns
% the cell face coordinates and non-grid array the cell center. The grid
% array divide the distance betwee adjacent lines according to SplitArray.
% A second arbitraty letter in code indicates that the cell properties are
% to be subdivided and not taken as they are. This is done for instance
% with the cell size arrays.
% Repeat the procedure for each direction sepaprately. To split an already
% splitted array, just run the procedure again with the new array and the
% correct spitarray for the new array. This may get complex though.
% see explanation at the end of this file.
% TO 100602

code=code(1);
if length(code)>1, mustdivide=1; else mustdivide=0; end
    
% Rotate so that all directions are treated equally
switch lower(code)
    case 'x', fprintf('Refining array in x-direction\n');
        AOld=permute(AOld,[2,1,3]);
    case 'y', fprintf('Refining array in y-direction\n');
        AOld=permute(AOld,[1,2,3]);
    case 'z', fprintf('Refining array in z-direction\n');
        AOld=permute(AOld,[3,2,1]);        
    otherwise
        error('Unknown code, see help for this function');
end

% Grid Arrays are treated separately here
if code=='X' || code=='Y' || code=='Z'
        Store=AOld(1,:,:);   % Keep the first coordinates
        AOld=diff(AOld,1,1); % Compete cell sizes in between
        mustdivide=1;
end

if mustdivide
    for i=1:size(SplitArray,2),
        AOld(SplitArray(1,i),:,:)=AOld(SplitArray(1,i),:,:)/SplitArray(2,i);
    end
end

NOld=size(AOld,1);
IOld=ones(2,NOld);       % store split info
IOld(1,:)=1:NOld;        % first line of splitter holds cold indices 1..NOld
IOld(2,SplitArray(1,:))=SplitArray(2,:); % second line holds sudivision numbers

NNew=sum(IOld(2,:));  % Total number of new columns by summing second line
INew=ones(2,NNew);      % Holds indeces of new RC's
INew(1,:)=1:NNew;     % First line is a dummy
for iOld=NOld:-1:1
    INew(2,1:sum(IOld(2,1:iOld)))=iOld; % Reference to old RC
end

ANew=AOld(INew(2,:),:,:);
if code=='X' || code=='Y' || code=='Z'
    AOld=zeros(size(ANew,1)+1,size(ANew,2),size(ANew,3)); % Reuse AOld
    AOld(2:end,:,:)=ANew;  % Put ANew in AOld starting on second row
    AOld(1,:,:)=Store;     % Add head line (initial coordinates
    ANew=cumsum(AOld,1);   % Compute new gridline coordinates
end

%% Rotate back
switch lower(code)
    case 'x', ANew=permute(ANew,[2,1,3]);
    case 'y', ANew=permute(ANew,[1,2,3]);
    case 'z', ANew=permute(ANew,[3,2,1]);
end

%% Procedure
% Split layer according to SplitArray
% The size of the array to split is obtained from the size of A
% The dimension along which to split is obtained from the code a captical
% code indicates that the array is a grid (cell face coordinates instead of
% a cell array). A second letter indicates that the array values must be
% divided according to the new splitted arrays. This is implied by the grid
% arrays.
% All arrays are turned so that the plitted diretion is the first
% dimension. This simplifies the code, all directions are then treated
% equally.
% First complete the SplitArray with all cells. The second line gets ones
% for the cells that are not splitted.
% Then we make the array of the new cells: % NNew=1:sum(IOld)
% INew is a two line array of size 2,NNew. The first line are the cell
% indices fo the new array. The seond line gets the corresponding indices
% of the old array. It is then trivial to generate the new array of the
% correct size filled with the values of the cells of the old array.
% This procedure should always work.
% To treat the grid arrays equally, we divide them into the coordinates of
% first layer (1,:,:) and the differences, i.e. the cell widths. This array
% is then treated like all others with subdivision implied.
% Once the new arrays are generated, we add the old coordinates to the
% front and apply a cumsum to regain the new coordinates.
% In a last step, all arrays are turned into their original positions.

