function  [bus, line] = write_pst(nbus, nbranch, connect, elec, busnames);
%
% 
% Uses the output of the read_ngc and order_ngc functions to produce matrices
% "bus" and "branch", suitable for use as the input data (bus and line) for
% the Matlab Power System Toolbox, as described below.
% 
% Note that only kV(n), from, to, Rpu, Xpu and Bpu  are specifically given in
% the input data used by the read_ngc function - other parameters are estimates.
%
%
%  Bus column format, with examples:
%
%      #    |V|    <V     Pg    Qg     Pl    Ql     G     B    type   Qgmax   Qgmin  kV(n)  Vpumax  Vpumin
%      12   1.02   -3.4   6     2      1     0      0     1     2      999    -999     400      3      0
%
%  Notes:  number (#) is arbitrary. the order of buses in this matrix is that used in the program, and 
%             the given number is just a name for the bus.  With the NGC data, ordered to keep 
%             electrically close buses close in the numbering too, the bus numbers are simply  1:nbus. 
%          |V| is in pu
%          angle (<V) is in degrees.
%          Pg, Qg, Pl, Ql  are pu generated and load quantities at a bus.
%          type  is 1 for the swing bus, 2 for a PV bus, and 3 for a PQ bus.
%          Qgmax and Qgmin are not used here  - generator limits are not implemented.
%          nominal voltage (kV(n)) is included to allow the output to show real (not pu) values.  
%          Vpumax and Vpumin are used to limit each solution iteration's output voltages. These are 
%             initially set to extreme values so as not to have effect in remotely plausible situations.
% 
%  Line column format, with examples:
%
%   from   to     Rpu     Xpu    Bpu    |r|    <r     |r|max   |r|min    delr  
%   12     120    0.001   0.01   0.5    1       0     2         0        0
%  
%  Notes: from and to use the 'given' bus numbers (the arbitrary user-assigned names, of col. 1 of
%             the bus matrix).  There may be more than one line of the same end numbers - this is OK.
%             With the NGC data, ordered to keep electrically close buses close in the numbering too, 
%             the bus numbers are simply  1:nbus.
%         R and X are pu series values.
%         Bpu is the total shunt susceptance modelled as B/2 at either end of the line  (pi model).  
%             If the branch is a transformer, Bpu can be made negative to represent the transformer 
%             magnetising susceptance.
%         |r|, <r |r|max&min and delr  all relate to tap changing and phase shifting, which are not 
%             inluded in this program.
%
%                 |r|  :  1
%      |              > <       | Vs'              ____________                  |
%   Vs |--------------> <-------|------|----------|___R + jX___|----------|------|  Vr
%      |              > <       |     _|_                                _|_     |
%     FROM bus                        ___ B/2                        B/2 ___   TO bus
%                                      |                                  |

disp(' ')
disp(' write_pst.m ')
disp(['  For the moment, a very basic setup:  only branch impedances and bus numbers are known.',sprintf('\n'), ...
    '  All other details will be filled in to cause tap-changing to be ignored, and bus types to be',sprintf('\n'), ...
    '  to be all PQ apart from the first (swing).  No loads or generation are included.',sprintf('\n'), ...
    '  Liberal limits are put on bus voltage, Q generation, etc.',sprintf('\n'), ...
    '  Further development is needed to make use of any generator and load data that can',sprintf('\n'), ...
    '  be found for the system.',sprintf('\n'),sprintf('\n') ])

% For the moment, a very basic setup:  only branch impedances and bus numbers are known.
% All other details will be filled in to cause tap-changing to be ignored, and bus types 
% all to be PQ apart from the first (swing).  No loads or generation are included.
% Liberal limits are put on bus voltage, Q generation, etc.
% Further development is needed to make use of any generator and load data that can
% be found for the system.


% Make bus matrix

% zeros
bus = zeros(nbus,15);

%                bus#          |V|pu          <V 
bus(:,1:3) = [ (1:nbus)' , ones(nbus,1) , zeros(nbus,1) ];

%             swing at #1, Pq else
bus(:,10) = [ 1; 3*ones(nbus-1,1) ];

% interpret the first of the two descriptor codes (voltage rating) in the bus code-name
for i=1:nbus,
    temp = busnames{i,1}(1,5);
    if     temp(1)=='6';  Vnom =  66.0;
    elseif temp(1)=='4';  Vnom = 400.0;
    elseif temp(1)=='3';  Vnom =  33.0;
    elseif temp(1)=='2';  Vnom = 275.0;
    elseif temp(1)=='1';  Vnom = 132.0;
    end
    bus(i,13) = Vnom;
end

%                   Qgmax               Qgmin
bus(:,11:12) = [ 999*ones(nbus,1), -999*ones(nbus,1) ];

%                   Vmax              Vmin
bus(:,14:15) = [ 2*ones(nbus,1),  zeros(nbus,1) ];



% Make 'line' (branch) matrix

% zeros
line = zeros(nbranch,10);

%             bus numbers, to and from
line(:,1:2) =   connect;

%             pu R, X, B
line(:,3:5) =   elec;

%                |r| (tapratio)    <r (tap angle)    max and min r     tap step
line(:,6:10) = [ ones(nbranch,1), zeros(nbranch,1), ones(nbranch,2), zeros(nbranch,1) ];



% Now write to a file, as a PST system data mfile, if desired.
disp(' The "bus" and "line" matrices have been composed, and are the output of ')
disp(' this function.  It is now possible to write them to a PST-style m-file.')
disp(' ')

% select write or not
write_file = input('  Write file? y/n [y] > ','s');
if isempty(write_file) | write_file~='n',
    write_file = 'y';
end
if write_file=='n',
    return
end

% wait till non-empty string input, for filename
while 1,
    disp(' ')
    filename = input('  Enter the desired filename (it will get a ''.m'' extension) > ','s');
    if ~isempty(filename),
        break
    end
end

% add the  .m  extension,
filename = [ filename, '.m' ];

% find hostname, to go in the written file
if isunix,
    % qeury for hostname, and remove end-of-line character
    [dummy,nameofhost] = unix('hostname');
    nameofhost = nameofhost(1:length(nameofhost)-1)
    clear dummy
else nameofhost = ' ';
end

% open the file for writing, and insert formatted data and comment strings
fileid = fopen(filename,'w');

% fprintf(fileid,['% ',filename,sprintf('\n'),'% ',sprintf('\n'),'% Loadflow data in Power System Toolbox format.', ...
%         sprintf('\n'),'% ',sprintf('\n'),'% Written by ''write_pst.m'',  ',datestr(now), ...
%         sprintf('\n'),'% ',sprintf('\n'),'% ',sprintf('\n'),sprintf('\n'),' bus = [ ...',sprintf('\n') ]); 

fprintf(fileid,['%% %s \n%% \n%% System loadflow data, in Power System Toolbox format.' ...
                '\n%% \n%% Written by write_pst,  %s , on  %s \n%% \n \n bus = [ ... \n'], ...
               filename, datestr(now), nameofhost);


fprintf(fileid,['     %3.0f    %1.4f    %3.2f     %3.4f    %3.4f    %3.4f    %3.4f    '...   
        '%3.4f     %3.4f    %1.0f    %3.2f    %3.2f     %3.1f    %1.3f    %1.3f \n'] ,bus');

fprintf(fileid,']; \n\n line = [ ... \n');

fprintf(fileid,['     %3.0f    %3.0f    %3.4f     %3.4f     %3.4f    '...   
                 '%1.2f     %1.2f    %1.2f    %1.2f    %1.2f \n'] ,line');

fprintf(fileid,']; \n');
         
fclose(fileid);


% display exit status
disp(' ')
disp([' write_pst :  file ',filename,' written.'])
disp(' ')