% % branch              %R      %X       %B        km   MVAmean  type
% 
% % 275kV
% WARL20 <-> TILB21     0.061   0.552    3.768B    14   1110     OHL (4 cond.)
% CHSI20 <-> WWEY21     0.069   0.623    4.25      15    871	 OHL (2 cond.)
% COVE20 <-> HAMH20     0.046   0.708    6.46      20   1280     OHL (4 cond.)
% COVE20 <-> BESW20     0.084   0.76     5.186     19    871     OHL (2 cond.)
% SJOW20 <-> TOTW21     0.029   0.36   177.92      16    786	 CABLE
% 
% % 400kV
% ECLA40 <-> COWL40     0.039   0.604   24.66      36 	2533	 OHL (4 cond.)
% ECLA40 <-> AMEM4A     0.062   0.665   21.887     36 	1895	 OHL (2 cond.)
% SJOW4A <-> CITR40     0.006   0.056  123.894      7 	1294	 CABLE


% Use the above example data (from UK National Grid Co, Annual Report, 2001) to 
% calculate Surge Imedance loadings etc. for the line and cable types.
function example_branch_data();

S_base = 100e6;
omega = 2*pi*50;
disp(' '); disp(' per km values:');


% 275 kV
V_bases = 1e3 * [ 275; 275; 275; 400; 400; 400 ];

pu_kms  =   [  ([0.069 0.623 4.25 ]/15+[0.084 0.76 5.186]/19)/2
               ([0.061 0.552 3.768]/14+[0.046 0.708 6.46 ]/20)/2
                [0.029 0.36  177.92 ]/16
                [0.062 0.665 21.887 ]/36
                [0.039 0.604 24.66 ]/36
                [0.006 0.056 123.894 ]/7   ]   /  100 ;
    
    
mva_means  =  [  (871+871)/2;
                 (1110+1290)/2;
                 786;
                 1895;
                 2533;
                 1294 ];
             
set_names = { '275kV, OHL, twin-conductor';
              '275kV, OHL, quad-conductor';
              '275kV, Cable';
              '400kV, OHL, twin-conductor';
              '400kV, OHL, quad-conductor';
              '400kV, Cable' };
            
for ii=1:1:6,
    calc_branch(V_bases(ii,:),  pu_kms(ii,:), mva_means(ii,:), omega, S_base, set_names{ii,1} );
end


%%%%%%%%%%%%%%%% caclulate
function [] = calc_branch(V_base, pu_km, mva_mean, omega, S_base, set_name)

Z_base = V_base^2 / S_base;
Ohm_km = [  pu_km(1,1:2)*Z_base  ,  pu_km(1,3)/Z_base  ];
RLC_km = [  Ohm_km(1,1),  Ohm_km(1,2)/omega,  Ohm_km(1,3)/omega  ];
Z_0    = sqrt(RLC_km(1,2) / RLC_km(1,3));
SIL    = 3 * V_base^2 / Z_0 / 1e6;
vprop  = 1 / sqrt(RLC_km(1,2)/1000 * RLC_km(1,3)/1000) / 299792458;

disp(' ')
fprintf(' %s,   Z_base = %4.0f \n', set_name, Z_base);
fprintf('  OhmR   OhmX   mSiemB\tR(Ohm) L(mH)  C(mF) \tZ_0  SIL(MVA) limit(MVA) prop_vel \n');
fprintf('  %1.3f  %1.3f  %1.3f \t%1.3f  %1.3f  %1.3f \t%2.1f  %4.0f  %4.0f  %1.3f \n', ...
    [ Ohm_km.*[1,1,1e3], RLC_km.*[1,1e3,1e3], Z_0,  SIL,  mva_mean, vprop ] );

