function perm = order_ngc(startpos, method)
%
% function perm = order_ngc(startpos, method)
%
% N. Taylor, Nov. 2002.
%
% 
% The branch data for the system, as output by read_ngc.m, has ordering 
% based on branch-type and alphabetical order.  The station names are listed 
% alphabetically, but individual buses have no ordering.  To create a bus matrix
% for the PST input, a bus ordering must be assigned.  Rather than using alphbetical
% or type ordering, this function works through the line connectivity so that 
% proximity in network topology/(geography, too) corresponds strongly to proximity
% in the bus matrix.  This also has the effect of making the Ybus matrix quite
% close to being diagonal.
% As there are sometimes several branches from a bus, yet the ordering here has
% to be one-dimensional, the particular ordering will be dependent on the choice
% of starting bus, which is specified by its 6-character name (e.g. 'DIDC42').
%
% IN
% Argument:               startpos (starting bus name, string 1x6),
%                         method (0 for Matlab's SYMRCM function, 1 for home-spun attempt)
% Global variables:       nstat (number of stations, 1x1)
%                         nbus (number of buses, 1x1), 
%                         nbranch (number of branches, 1x1), 
%                         connect (branch to- and from-bus numbers, (nbranch,2) )
%                         elec (matrix of R X B pu for each bus, ordered as in 'connect', (nbranch,3) ),
%
% OUT
% The following variables are going to be generated by the next section...
%                         startpos (starting bus number, replacing the name used before, (1,1) ),
%                         perm   (vector s.t. perm(i) is the index into the original bus order
%                                 for the bus that is now designated as bus number i.
%                                 i.e.  list(perm,1) = newly_ordered_list

jay = sqrt(-1);
global nbranch nstat nbus busnames connect elec



if method==0,   %%%%  SYMRCM  function
    
    % Make a matrix describing pure topological connectivity of the system (1 or 0) 
    % and use the Matlab SYMRCM function to arrange it to be close to diagonal.  
    % This is simple to use, works well, but doesn't allow for sub-ranking by admittance.
    conns = zeros(nbus,nbus); 
    % conns is made to have one for a connection, 0 otherwise
    for i=1:length(connect),
        conns(connect(i,1),connect(i,2))=1; 
        conns(connect(i,2),connect(i,1))=1; 
    end
    % Matlab's symrcm does symmetric reverse Cuthill-McKee
    % permutation, and returns the permutation vector.
    perm = symrcm(conns);
    clear  conns i
    
    
    
elseif method==1,   %%%%  Own method, also accounting for impedances in the ordering
    
    % convert the bus name string in the variable "startpos" into the respective bus-number
    for i=1:nbus,
        if strcmp(busnames{i,1},startpos),
            startpos = i;
            break
        end
    end
    
    % the bus referenced in startpos will be put first in the bus order
    perm = startpos;
    
    done = 0;
    i = 0;
    % loop until the perm vector is full
    while done==0,
        i = i+1;
        
        % find all branches that connect to the current bus,
        froms = find(connect(:,1)==perm(i,1)); % instances of the current bus as a 'from' bus
        tos = find(connect(:,2)==perm(i,1)); % instances of the current bus as a 'to' bus
        farbuses = [];
        branches = [];
        if ~isempty(froms),        
            branches = froms;
            farbuses = connect(froms,2);
        end
        if ~isempty(tos),        
            branches = [branches; tos];
            farbuses = [farbuses; connect(tos,1)];
        end
        
        % remove from the list those branches whose far buses are already in perm
        index = [];
        for k=1:length(perm),
            temp = find(farbuses(:,1)==perm(k,1));
            index = [index; temp]; % indices in the vectors of connected branches 
            % and buses of all already-included buses 
        end
        if ~isempty(index),
            temp = [];
            for k=1:length(farbuses),
                if isempty(find(index(:,1)==k)),
                    temp = [temp; k];
                end
            end
            branches = branches(temp,1);
            farbuses = farbuses(temp,1);
        end
        
        % now farbuses gives numbers of the buses at remote ends of all branches from bus i,
        % and branches gives branch numbers (indices into the 'connect' matrix) for these branches
        % in the same order.  Any parallel branches must now be 'joined' into one, and the 
        % branches should be ranked to make electrically close ones come first in the list...
        
        if ~isempty(farbuses),  % if any more buses are to be added to 'perm'
            
            if length(farbuses)>1, % if any parallelling/ordering is needed
                
                % find any parallel branches (i.e. a repeat occurence of a number in farbuses)
                for k=1:length(farbuses),
                    temp = find(farbuses==farbuses(k,1));
                    
                    % if the farbuses vector contains multiple instances of a bus, not already in newbuses...
                    if sum(temp)>1 & isempty(find(farbuses(k,1)==newbuses)),
                        % find the resultant impedance of all the parallel branches
                        admt = 0;
                        for m==1:length(temp),
                            admt = admt + 1/( elec(branches(temp(m)),1) + jay*elec(branches(temp(m)),2) );
                        end
                        impd = 1/admt;
                        impedance =  [ impedance ]
                    end
                    
                end 
                if
                    % index sorts buses in ascending impedance -> descending admittance
                    [dummy, index] = sort(abs(elec(branches,1)+jay*elec(branches,2)));
                    farbuses = farbuses(index);
                    branches = branches(index);
                end
                
            end
            
            % add the new bus numbers, ranked, to the bus permutation vector 'perm'
            perm = [ perm; farbuses ];
            % check whether all buses have been included yet
            if i==nbus,
                done = 1;
            end
                        
        end % end if ~isempty(farbuses)
        
        
    end % end 'while'
    
    % do a basic check on there being all the expected numbers present..
    for i=1:length(perm),
        temp = find(perm==i);
        if sum(temp)==0,
            disp([sprintf('\n'),' Bus number ',num2str(i),' was not found in the ''perm'' vector! (order_ngc)',sprintf('\n'), ...
                    ' The ''perm'' vector should contain every number from 1 to length(perm).' ])
        elseif sum(temp)<1,
            disp([sprintf('\n'),' Bus number ',num2str(i),' was found ',num2str(sum(temp)), ...
                    ' times in the ''perm'' vector! (order_ngc)',sprintf('\n'),...
                    ' The positions of these entries in ''perm'' are ',num2str(temp'),sprintf('\n'),...
                    ' The  ''perm'' vector should contain one instance of every number from 1 to length(perm).' ])
        end
    end
    
    % clean up defunct variables...
    clear   done i froms tos farbuses branches k temp index dummy
    
end % end method-type 'if'