function ngcpst
%
% function ngcpst.m
%
% This function provides an interface for working with system data taken from the
% U.K. National Grid Company's  appendices to its "Seven year Statements".
%
% Some function files are called, too, to do the processing of the initial text 
% files, and to do ordering of the bus numbers and creation of PST data matrices.
%
% As usual, there are near the beginning of the script some variables that can be
% set as constant by the user, and doing this will prevent tese being required as
% input on every run.
%
% N. Taylor,  Nov, 2002., 
% modified (add  display of Z0, prop veloc, etc) Sept 2003,
% corrected display of % as pu, Sept 2004

disp(' ')
disp('************************  NGC data processing  *******************************')
disp(' ')
disp(' This uses system data from text-files, as described by "help read_ngc".')
disp(' ')

% prepare fresh workspace
jay = sqrt(-1);
Ybus = [];
% these need to exist...
branch_data = [];
name_data = [];


%%%%%%%%%%%%%%  SET OR COMMENT OUT THE FOLLOWING, TO AVOID/ALLOW INPUT EACH TIME  %%%%%%%%%%%%%%%%%%%
branch_data = 'ngcdata.txt'; % text file of branch data
name_data = 'ngcnames.txt'; % text file of four-letter station codes and corresponding names
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global statnames   % {nstat,2} first column of cells holds 4-letter string code of a substation,
%           second column holds a variable length string of the corresponding full name.

global busnames    % {nbus,2} same format as statnames.  The difference is that there may be several distinct
%          buses at a station;  their code names (6-character) are the station code + 2 more.

global nbranch     % scalar - number of branches
global nstat       % scalar - number of stations
global nbus        % scalar - number of buses

global connect     % (nbranch,2) assigned reference numbers of to and from buses for each line
%             (these numbers are also the indices of the bus names in "busnames")

global elec        % (nbranch,3) R, X, B  in pu for each branch
global leng        % (nbranch,1) length in km for each branch, with exactly 0 if branch not OHL or UGC
global type        % {nbranch,1} variable length strings to describe the branch type
global rate        % (nbranch,4) columns are Spring, Summer, Autumn, Winter  max MVA loadings for each branch



% Call the read_ngc function to extract structures and 
% matrices by interpreting the text input files.
if isempty(branch_data),
    disp(' ')
    branch_data = input(' Enter the name of the branch data file   > ','s');
end
if isempty(name_data),
    disp(' ')
    name_data = input(' Enter the name of the bus-name data file > ','s');
end
disp(' ')
disp(' Calling read_ngc...')
try
    read_ngc(branch_data,name_data);
catch
end
disp(' ... read_ngc finished.');


% initialise 'perm' in case plotting of permutation is called
% before the ordering function has created perm.
perm = (1:nbus)';

end_main=0;
while end_main == 0,
    
    disp(' ')
    disp(' Enter:')
    disp('       1 for interactive display of system connection and name data')
    disp('       2 for display of line and cable fundamental parameters')
    disp('       3 to use a command prompt to study data, then return to this menu')    
    disp('       4 to make the Ybus matrix for this system')
    disp('       5 to ''spy'' plot the system connectivity matrix')    
    disp('       6 to calculate admittance between two buses (Ybus already formed)')
    disp('       0 to EXIT')
    disp(' ')
    
    theswitch = input('       ? >> ','s');
    theswitch = str2num(theswitch);
    if isempty(theswitch),
        theswitch = 100;
    end
    
    switch theswitch
        
        
    case 0,
        disp(' ')
        if 'y' == input(' Confirm exit y/n [n] > ','s'),
            disp(' '); disp(' ');
            end_main=0;
            break
        end
        
        
        
    case 1,
        
        % call function to display information 
        % on names and connectivity
        % (function is at bottom of this file)
        try
            disp_ngc  
        catch
            % if disp_ngc fails, just go back to seletion loop...
            disp(' '); disp('  Note: disp_ngc exited abnormally!'); disp(' ');
        end
        
        
    case 2,
        
        % this function (below) displays properties such as surge-impedance,
        % propagation velocity, per-km SI and pu R,X,B.
        branch_props;
        
        
    case 3,
        
        disp(' ')
        disp('  Basic data in the workspace are:')
        disp('     Scalars:  nbranch, nstat, nbus  - numbers of branches, bus-groups (stations), buses')
        disp('     Matrix (nbranch,2):    connect - assigned numbers of to- and from- buses for each line') 
        disp('     Matrix (nbranch,3):    elec  - R,X,B pu')
        disp('     Matrix (nbranch,4):    rate  - Spring, Summer, Autumn, Winter   MVA max')
        disp('     Vector (nbranch,1):    leng  - branch length in km (0 if reactor/transformer, etc.')
        disp('     Structure {nbranch,1}: type  -  description of type of branch, free text')
        disp('     Structure {nstat,2}:   statnames  - first column, 4-letter code; second column full name')
        disp('     Structure {nbus,2}:    busnames  - first column, 6-letter code; second column full description')
        disp(' ')
        disp('  Also, (if suitable options already have been taken to generate them) :')
        disp('     Vector (nbus,1):       perm  - permutation to order buses, i.e. newbus=oldbus(perm);')
        disp('     Matrix (nbus,nbus):    Ybus  - system admittance matrix, with permutation if requested')
        disp(' ')
        disp('  Keyboard mode (K>>) will now be used.  ')
        disp('  To return to the menu, type the word ''return'' and press Enter..')
        disp(' ')
        keyboard
        
        
    case 4,
        
        Ybus = zeros(nbus,nbus);
        halfB = zeros(nbus,1);
        
        for  i=1:nbranch,
            Ybus(connect(i,1),connect(i,2)) = Ybus(connect(i,1),connect(i,2))  - 1 / (elec(i,1) + jay*elec(i,2));
            Ybus(connect(i,2),connect(i,1)) = Ybus(connect(i,2),connect(i,1))  - 1 / (elec(i,1) + jay*elec(i,2));
            halfB(connect(i,1)) = halfB(connect(i,1)) + jay*elec(i,3)/2;
            halfB(connect(i,2)) = halfB(connect(i,2)) + jay*elec(i,3)/2;
        end
        
        temp = -sum(Ybus);
        Ybus = Ybus + diag(temp) + diag(halfB);
        clear i temp halfB
        Ybus = Ybus(perm,perm);
        disp(' ')
        disp(' Ybus calculation is finished.')
        disp(' ')
        
        
    case 5,   
        
        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
        
        figure
        spy(conns(perm,perm))
        title('SPY on connection data - a point if buses x,y are connected')
        xlabel('bus number')
        ylabel('bus nuber')
        
        
    case 6,
        
        if isempty(Ybus),
            disp(' ')
            disp(' The Ybus matrix is empty.  Run Ybus calculation, then retry this option.')
        else
            [Zmut, Vdist] = ztrans_ngc(Ybus);
        end
        
    
    end  % end switch, case
    
end  % end main while


disp(' ')
disp(' The following global variables may be accessed by declaring them as global in the workspace:')

disp(['  statnames   % {nstat,2} first column of cells holds 4-letter string code of a substation,',sprintf('\n'),sprintf('\n'),...
        '                second column holds a variable length string of the corresponding full name.',sprintf('\n'),...
        '  busnames    % {nbus,2} same format as statnames.  The difference is that there may be several distinct',sprintf('\n'),...
        '                 buses at a station;  their code names (6-character) are the station code + 2 more.',sprintf('\n'),...
        '  nbranch     % scalar - number of branches',sprintf('\n'),...
        '  nstat       % scalar - number of stations',sprintf('\n'),...
        '  nbus        % scalar - number of buses',sprintf('\n'),...
        '  connect     % (nbranch,2) assigned reference numbers of to and from buses for each line',sprintf('\n'),...
        '                (these numbers are also the indices of the bus names in "busnames")',sprintf('\n'),...
        '  elec        % (nbranch,3) R, X, B  in pu for each branch',sprintf('\n'),...
        '  leng        % (nbranch,1) length in km for each branch, with exactly 0 if branch not OHL or UGC',sprintf('\n'),...
        '  type        % {nbranch,1} variable length strings to describe the branch type',sprintf('\n'),...
        '  rate        % (nbranch,4) columns are Spring, Summer, Autumn, Winter  max MVA loadings for each branch',sprintf('\n') ])
disp(' ')
disp('  global  statnames busnames nbranch nstat nbus connect elec leng type rate')
disp(' ')
disp(' ')




%%%%%%%%%%%%%%% Function for displaying (long and loopy - hence kept from main program) %%%%%%%%%%%%%%%%%%

function [] = disp_ngc
%
% function [] = disp_ngc
%
% N. Taylor, Nov. 2002.
%
% Interaction with the NGC data matrices and structures, to give such details as
% all other buses directly connected to a bus; impedances; translate between
% full station or bus names and the assigned numbers, etc.

global  nstat nbus nbranch statnames busnames connect type elec leng rate

end_while = 0;
give_list = 1;
while end_while == 0,
    
    disp(' ')
    if give_list==1,
        disp(' Enter the bus or station details:')
        disp('        *  a number is taken as the number of a bus')
        disp('        *  a four-letter upper-case string is taken as the code of a station')
        disp('        *  a six-character string starting with four upper-case letters ')
        disp('                                     is taken as a bus description code')
        disp('        *  any other string is taken as part of the full station name (case sensitive)')
        disp(' ')
        disp('        Type  0 to EXIT, or press Return to REPEAT this list ...')
    end
    give_list = 0;
    request = [];
    disp(' ')
    request = input('        ? >> ','s');
    
    
    if isempty(request),
        % go into while-loop again, ready to display options
        give_list = 1; 
        
        
    else         
        
        % else if string contains no letters
        if sum(isletter(request))==0,
            
            request = str2num(request);
            
            % exit this function if 0 has been entered
            if request==0,
                end_while = 1;
                break
            end
            
            % warn if number not available 
            % (function will continue to input prompt)
            if request<1 | request>nbus,
                disp(' ')
                disp('  The number given is outside the range of bus numbers.')
                disp(' ')
                
                % if all OK, give output
            else
                disp(' ')
                disp(['  Bus ',num2str(request),' has code  ',busnames{request,1},', and name  ',busnames{request,2} ]);
                % find remote buses directly connected to this one (by a single branch)
                rembuses = [];
                for i=1:nbranch,
                    temp1 = find(connect(:,1)==request);                            
                    temp2 = find(connect(:,2)==request);
                    branchindx = [temp1; temp2];
                    rembuses = [connect(temp1,2); connect(temp2,1)];
                    clear  temp1 temp2
                end
                disp(' ')
                disp(['  There are ',num2str(length(rembuses)),' branches from this bus.'])
                disp(' ')
                disp('  The remote buses'' details are:')
                for i=1:length(rembuses),
                    disp(['   #',num2str(rembuses(i)),sprintf('\t'),busnames{rembuses(i),1}, ...
                            ' ',sprintf('\t'),busnames{rembuses(i),2}   ])
                end
                disp(' ')
                disp('  The branch details are:')
                for i=1:length(rembuses),
                    branchno = branchindx(i);
                    disp(['   #',num2str(branchno),'   ',sprintf('\t'),busnames{request,1},' <-> ',busnames{rembuses(i),1}, ...
                            '     ',num2str(100*elec(branchno,1)),'%R  ', ...
                            num2str(100*elec(branchno,2)),'%X  ',num2str(100*elec(branchno,3)),'%B  ' ...
                            '    ',num2str(round(leng(branchno,1))), 'km   ', ...
                            num2str(round(sum(rate(branchno,:))/4)),' MVAmax  ',sprintf('\t'),type{branchno,1} ])
                end
                clear  i rembuses branchindx request branchno
            end
            
            % else if a useable string has been input (now test _what_ sort of letters...)
        else
            
            nchar = length(request);
            
            % string is four upper case characters - treat as station code
            if nchar==4 & strcmp(request,upper(request)),
                
                % find station position in statnames structure
                for i=1:nstat,
                    if strcmp(statnames{i,1},request),
                        station = i;
                        break;
                    end
                end
                
                % complain, or else display details
                if isempty(station),
                    disp(' ')
                    disp('  The search string was not found.')
                    disp(' ')
                else
                    % find buses at the selected station...
                    itsbuses = [];
                    for i=1:nbus,
                        if strcmp(statnames{station,1},busnames{i,1}(1,1:4)),
                            itsbuses = [itsbuses; i ];
                        end
                    end
                    disp(' ')
                    disp(['  Station code ',statnames{station,1},' stands for ',statnames{station,2} ])
                    disp(' ')
                    disp(['  There are ',num2str(length(itsbuses)),' buses at this station:'])
                    disp(' ')
                    for i=1:length(itsbuses),
                        disp(['   ',busnames{itsbuses(i),1},'       ',busnames{itsbuses(i),2}])
                    end
                    disp(' ')
                end
                
                clear i itsbuses
                
                % string is 6 characters, the first four upper case - treat as bus code
            elseif nchar==6 & strcmp(request(1:4),upper(request(1:4))),
                
                % find the (first..) bus matching the 6-character code
                busno = [];
                for i=1:nbus,
                    if strcmp(busnames{i,1},request),
                        busno = i;
                        break;
                    end
                end
                if isempty(busno),
                    disp(' ')
                    disp('  The search string was not found.')
                    disp(' ')
                else
                    % find remote buses directly connected to this one (by a single branch)
                    rembuses = [];
                    for i=1:nbranch,
                        temp1 = find(connect(:,1)==busno);
                        temp2 = find(connect(:,2)==busno);
                        branchindx = [temp1; temp2];
                        rembuses = [connect(temp1,2); connect(temp2,1)];
                        clear temp1 temp2
                    end
                    disp(' ')
                    disp(['  Bus code ',busnames{busno,1},' is bus #',num2str(busno),' ;',sprintf('\t'),busnames{busno,2} ])
                    disp(' ')
                    disp(['  There are ',num2str(length(rembuses)),' branches from this bus.'])
                    disp(' ')
                    disp('  The remote buses'' details are:')
                    for i=1:length(rembuses),
                        disp(['   #',num2str(rembuses(i)),sprintf('\t'),busnames{rembuses(i),1},' ',sprintf('\t'), ...
                                busnames{rembuses(i),2} ])
                    end
                    disp(' ')
                    disp('  The branch details are:')
                    for i=1:length(branchindx),
                        branchno = branchindx(i);
                        disp(['   #',num2str(branchno),'   ',sprintf('\t'),busnames{busno,1},' <-> ',busnames{rembuses(i),1}, ...
                                '     ',num2str(100*elec(branchno,1)),'%R  ', ...
                                num2str(100*elec(branchno,2)),'%X  ',num2str(100*elec(branchno,3)),'%B  ' ...
                                '    ',num2str(round(leng(branchno,1))), 'km ',sprintf('\t'), ...
                                num2str(round(sum(rate(branchno,:))/4)),' MVAmax  ',sprintf('\t'),type{branchno,1} ])
                    end
                    clear  i rembuses busno branchno
                end
                
                % else, try to find a match within all the full bus names...
            else  % if not 4 uppercase, 6 with first 4 upper case, or non-letter (number?)
                
                % find stations whose names contain the input string
                matches = [];
                for i=1:nstat,
                    if strfind(statnames{i,2},request);
                        matches =  [ matches; i ];
                    end
                end
                
                % complain, variously, if the number of matches ~= 1
                if isempty(matches),
                    disp(' ')
                    disp('  There was no match for your search string in the full station names.')
                    disp(' ')
                elseif length(matches) > 1,
                    disp(' ')
                    disp('  More than one station name was found that contains the search string: ')
                    disp(' ')
                    for i=1:length(matches),
                        disp(['   ',statnames{matches(i),1},'     ',statnames{matches(i),2} ])
                    end
                    disp(' ')
                    
                    % if there is but one match, give details    
                else
                    disp(' ')
                    disp('  The following station name contains the given search string: ')
                    disp(' ')
                    disp(['   ',statnames{matches,1},'     ',statnames{matches,2} ])
                    % find buses at the selected station...
                    itsbuses = [];
                    for i=1:nbus,
                        if strcmp(statnames{matches,1},busnames{i,1}(1,1:4)),
                            itsbuses = [itsbuses; i ];
                        end
                    end
                    disp(' ')
                    disp(['  There are ',num2str(length(itsbuses)),' buses at this station:'])
                    disp(' ')
                    for i=1:length(itsbuses),
                        disp(['   ',busnames{itsbuses(i),1},'       ',busnames{itsbuses(i),2}])
                    end
                    disp(' ')
                end
                
                clear i matches itsbuses
                
            end  % end if 4/6/etc letters
            
        end  % end if  numbers/spaces/letters
        
    end  % end if ~isempty(request)
    
end  % end while

clear   end_while give_list request nchar i
disp(' ')

%%%%%%%%  end of  disp_ngc  function  %%%%%%%%%%%%%%



%%%%%%%%%%  Function for displaying line and cable length-independent properties  %%%%%%%%%%

function [] = branch_props();

% import the data to this function
global  statnames busnames nbranch nstat nbus connect elec leng type rate

mvabase=100;
angf_base=2*pi*50;

% initialise matrices to hold %/km  and Z0 and v for each branch.
line_400_par=[];  cable_400_par=[];
line_275_par=[];  cable_275_par=[];
line_132_par=[];  cable_132_par=[];

for n=1:1:nbranch,
    if ( strcmp(type{n,1},'OHL') | strcmp(type{n,1},'CABLE') ),
        
        % crude way to ascertain the branch voltage from 5th char of bus codes...
        if ( busnames{connect(n,1),1}(1,5) == '4' ),  % 400kV
            
            base_z = (400e3 ^ 2) / ( 1e6 * mvabase );
            % provide factors for conversion of pu impedance to mOhm, mH and uF
            pu2si = [ base_z*1e3  base_z*1e3/angf_base  1e6/(angf_base*base_z) ];
            
            % per km R,X,B in pu (MVA base = mvabase, Vbase = component-rating)
            pc_per_km = elec(n,:)/leng(n,1) * 100;
            % per km R,X,B in SI units (mOhm, mH, mF)
            si_per_km = elec(n,:)/leng(n,1) .* pu2si;
            % surge impedance, pu
            Z0pu = sqrt(elec(n,2)/elec(n,3));
            % surge impedance, Ohms.
            Z0 = base_z * Z0pu;
            % propagation velocity, in pu of speed of light in a vacuum...
            vp = ( angf_base / sqrt( elec(n,2)*elec(n,3) / (1e3*leng(n,1))^2 ) )   / 299792458 ;
            
            temp = [ Z0pu Z0 vp pc_per_km si_per_km ];
            
            if strcmp(type{n,1},'OHL'),
                line_400_par = [ line_400_par; temp ];
            else
                cable_400_par = [ cable_400_par; temp ];
            end
            
            
        elseif ( busnames{connect(n,1),1}(1,5) == '2'),  % 275kV
            
                        base_z = (275e3 ^ 2) / ( 1e6 * mvabase );
            % provide factors for conversion of pu to mOhm, mH and uF
            pu2si = [ base_z*1e3  base_z*1e3/angf_base  1e6/(angf_base*base_z) ];
            
            % per km R,X,B in pu (MVA base = mvabase, Vbase = component-rating)
            pc_per_km = elec(n,:)/leng(n,1) * 100;
            % per km R,X,B in SI units (mOhm, mH, mF)
            si_per_km = elec(n,:)/leng(n,1) .* pu2si;
            % surge impedance, pu
            Z0pu = sqrt(elec(n,2)/elec(n,3));
            % surge impedance, Ohms.
            Z0 = base_z * Z0pu;
            % propagation velocity, in pu of speed of light in a vacuum...
            vp = ( angf_base / sqrt( elec(n,2)*elec(n,3) / (1e3*leng(n,1))^2 ) )   / 299792458 ;
            
            temp = [ Z0pu Z0 vp pc_per_km si_per_km ];
            
            if strcmp(type{n,1},'OHL'),
                line_275_par = [ line_275_par; temp ];
            else
                cable_275_par = [ cable_275_par; temp ];
            end
            
            
        elseif ( busnames{connect(n,1),1}(1,5) == '1'),  % 132kV
            
            base_z = (132e3 ^ 2) / ( 1e6 * mvabase );
            % provide factors for conversion of pu to mOhm, mH and uF
            pu2si = [ base_z*1e3  base_z*1e3/angf_base  1e6/(angf_base*base_z) ];
            
            % per km R,X,B in pu (MVA base = mvabase, Vbase = component-rating)
            pc_per_km = elec(n,:)/leng(n,1) * 100;
            % per km R,X,B in SI units (mOhm, mH, mF)
            si_per_km = elec(n,:)/leng(n,1) .* pu2si;
            % surge impedance, pu
            Z0pu = sqrt(elec(n,2)/elec(n,3));
            % surge impedance, Ohms.
            Z0 = base_z * Z0pu;
            % propagation velocity, in pu of speed of light in a vacuum...
            vp = ( angf_base / sqrt( elec(n,2)*elec(n,3) / (1e3*leng(n,1))^2 ) )   / 299792458 ;
            
            temp = [ Z0pu Z0 vp pc_per_km si_per_km ];
            
            if strcmp(type{n,1},'OHL'),
                line_132_par = [ line_132_par; temp ];
            else
                cable_132_par = [ cable_132_par; temp ];
            end
            
        end
    end
end

% for non-empty matrices, set an index to show non-emptiness, and calc. mean values
l400=0; c400=0; l275=0; c275=0; l132=0; c132=0;
if ~isempty(line_400_par)
    line_400_mean  = sum(line_400_par,1)  / length(line_400_par(:,1));
    l400=1;
end
if ~isempty(cable_400_par)
    cable_400_mean = sum(cable_400_par,1) / length(cable_400_par(:,1));
    c400=1;
end
if ~isempty(line_275_par)
    line_275_mean  = sum(line_275_par,1)  / length(line_275_par(:,1));
    l275=1;
end
if ~isempty(cable_275_par)
    cable_275_mean = sum(cable_275_par,1) / length(cable_275_par(:,1));
    c275=1;
end
if ~isempty(line_132_par)
    line_132_mean  = sum(line_132_par,1)  / length(line_132_par(:,1));
    l132=1;
end
if ~isempty(cable_132_par)
    cable_132_mean = sum(cable_132_par,1) / length(cable_132_par(:,1));
    c132=1;
end

end_while = 0;
while end_while == 0,
    
    disp(' ')
    disp(' Select:')
    disp('        1  -  for mean values of each voltage and branch-type')
    if l400; disp('        2  -  for 400kV OHL data'); end
    if c400; disp('        3  -  for 400kV UGC data'); end
    if l275; disp('        4  -  for 275kV OHL data'); end
    if c275; disp('        5  -  for 275kV UGC data'); end
    if l132; disp('        6  -  for 132kV OHL data'); end
    if c132; disp('        7  -  for 132kV UGC data'); end
    disp(' ')
    disp('        Type  0 to EXIT, or press Return to REPEAT this list ...')
    request = [];
    disp(' ')
    request = input('        ? >> ');
    
    if isempty(request),
        % nothing -- repeat menu
        
    elseif request==0,
        end_while = 1;
        
    elseif request==1,
        disp(' ')
        disp('====== Mean values of parameters (per kilometre) ======  ');  
        disp(' ')
        
        if l400,
            disp(' 400kV, line')
            fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                line_400_mean );
            disp(' ')
        end
        if c400,
            disp(' 400kV, cable')        
            fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                cable_400_mean );
            disp(' ')
        end
        if l275,
            disp(' 275kV, line')
            fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                line_275_mean );
            disp(' ')
        end
        if c275,
            disp(' 275kV, cable')
            fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                cable_275_mean );
            disp(' ')
        end
        if l132,
            disp(' 132kV, line')
            fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                line_132_mean );
            disp(' ')
        end
        if c132,
            disp(' 132kV, cable')
            fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                cable_132_mean );
            disp(' ')
        end
        
    elseif ( request>1 & request<8 ),
        all_arrays={ line_400_par, cable_400_par, line_275_par, cable_275_par, line_132_par, cable_132_par };
        disp(' ')
        fprintf(1,'  Z0%%   Z0       vp \t R%%     X%%     B%% \t R mOhm  L mH    C uF \n');
        branch_data = all_arrays{1,request-1};
        for pos=1:1:length(branch_data(:,1)),
            fprintf(1,'  %3.2f   %3.2f   %1.2f \t %2.3f  %2.3f  %2.3f \t %3.2f   %3.2f   %3.2f  \n', ...
                branch_data(pos,:)  );
        end
    end
end

%%%%%%%% end of branch_props function %%%%%%%%%%



%%%%%%%%%%  Function for calculating network impedance between two buses  %%%%%%%%%%%%%%%

function [Zmut, Vbus] = ztrans_ngc(Ybus)
%
% function [Zmut, Vbus] = ztrans_ngc(Ybus)
%
% The Ybus matrix is partitioned into four submatrices. The corresponding parts of the current 
% and voltage vectors are: 
%      the two buses whose transfer impedance is sought:   known voltage (1pu and 0pu), unknown current
%                                       all other buses:   unknown voltage, known current (0pu)
% Unknown voltage buses correspond to known currents, and vice versa    

global  nbus busnames

jay = sqrt(-1);


doneflag = 0;
while doneflag == 0,
    disp(' ')
    bus1 = input(' Reference number of the first bus  [1] > ');
    if isempty(bus1)
        bus1 = 1;
    end
    if bus1 < 1  |  bus1 > nbus,
        disp(' ')
        disp(' The given number does not have a corresponding bus.  ')
    else
        doneflag = 1;
    end
end

doneflag = 0;
while doneflag == 0,
    disp(' ')
    bus2 = input([' Reference number of the other bus [',num2str(nbus),'] > ']);
    if isempty(bus2),
        bus2 = nbus;
    end
    if bus2 < 1  |  bus2 > nbus,
        disp(' ')
        disp(' The given number does not have a corresponding bus')
    else
        doneflag = 1;
    end
end

disp(' ')
disp(' The first bus will be set to 1pu voltage, the other bus to 0pu.')
disp(' Others will be left floating, and the currents from the 0pu bus')
disp(' will be taken to indicate the approximate impedance between buses')
disp(' The accuracy of this method is diminished when there are large')
disp(' shunt susceptances in branches, especially if in a branch directly')
disp(' between the buses (the full Ybus, including B, is used).  ')


% set up indices of known and unknown quantities ...
% know V, do not know I,
indx1 = [bus1; bus2];  
% know I, do not know V,
temp = ones(nbus,1);
temp(indx1) = [0; 0];
indx2 = find(temp);
clear  temp

% Now, with I and V partitioned into type 1 and 2 buses,
% have that:
%               [I1]  =  [ A   B ]  [V1] 
%               [I2]     [ C   D ]  [V2]
%               
%   therefore:
%               I1  =  [ A  -  B * inv(D) * C ] * V1
%               
%   and:
%               V2  =  - inv(D) * C * V1
% 
%   which can be used to give I1 more succinctly,
%   as:
%               I1  =  [ A  B ] * [V1]              
%                                 [V2]

% set up Volt and Current vectors with all known values and with 
% arbitrary ones at unknown buses
Vbus = zeros(nbus,1);
Ibus = zeros(nbus,1);
Vbus([bus1,bus2]) = [1; 0];

% use partitioned Ybus matrix to calculate the unknown quantities
Vbus(indx2) = -inv(Ybus(indx2,indx2)) * Ybus(indx2,indx1) * Vbus(indx1);
Ibus(indx1) = Ybus(indx1,:) * Vbus;        

% calculate admittance & impedance between the given buses
Ymut = Ibus(bus2);  
Zmut = 1/Ymut;

% display the single impedance
disp(' ')
disp(' The calculated impedance between ')
disp(' ')
disp(['      #',num2str(bus1),sprintf('\t'),busnames{bus1,1},'     ',busnames{bus1,2}])
disp(' and')
disp(['      #',num2str(bus2),sprintf('\t'),busnames{bus2,1},'     ',busnames{bus2,2}])
disp(' ')
disp([' is  ',num2str(Zmut),' pu, '])
disp([' which has absolute value  ',num2str(abs(Zmut)),' pu'])

disp(' ')
disp([' To see a plot of all system bus-voltage magnitudes with bus ',num2str(bus1),' at 1pu'])
doplot = input([' and bus ',num2str(bus2),' at 0pu voltage, enter ''y''  y/n  [n] > '],'s');
disp(' ')
if isempty(doplot),
    doplot = 'n';
end
if doplot == 'y',
    bar(abs(Vbus))
    title(['Vbus vs. ordered bus number, resulting from bus ',num2str(bus1),sprintf('\n'), ...
            'at 1pu and bus ',num2str(bus2),' at 0pu, all others free'])
    xlabel('bus number')
    ylabel('pu voltage magnitude')
end

%%%%%%%%%%%%%%% end of ztrans_ngc function %%%%%%%%%%%%%%%%%