% Formation of Ybus
% This is called by the main program 

Y=sparse(zeros(nb,nb));

   for i=1:nline
incr=1/(nt(i,3)+j*nt(i,4));
Y((nt(i,1)),(nt(i,2)))=Y((nt(i,1)),(nt(i,2)))-incr;
Y((nt(i,2)),(nt(i,1)))=Y((nt(i,2)),(nt(i,1)))-incr;
Y((nt(i,1)),(nt(i,1)))=Y((nt(i,1)),(nt(i,1)))+incr+j*nt(i,5)/2;
Y((nt(i,2)),(nt(i,2)))=Y((nt(i,2)),(nt(i,2)))+incr+j*nt(i,5)/2;
   end

 for i=nline+1:(ntrans+nline)
incr=1/(nt(i,3)+j*nt(i,4));  
incr1=incr/nt(i,5);
incr2=(1-nt(i,5))*incr/(nt(i,5)*nt(i,5));
incr3=(nt(i,5)-1)*incr/nt(i,5);
Y((nt(i,1)),(nt(i,2)))=Y((nt(i,1)),(nt(i,2)))-incr1;
Y((nt(i,2)),(nt(i,1)))=Y((nt(i,2)),(nt(i,1)))-incr1;
Y((nt(i,1)),(nt(i,1)))=Y((nt(i,1)),(nt(i,1)))+incr1+incr2;
Y((nt(i,2)),(nt(i,2)))=Y((nt(i,2)),(nt(i,2)))+incr1+incr3;   
end

   for i=1:nshunt
incr=shunt(i,2)+j*shunt(i,3);   
Y((shunt(i,1)),(shunt(i,1)))=Y((shunt(i,1)),(shunt(i,1)))+incr;
   end

  
YDQmod=[imag(Y) real(Y);real(Y) -imag(Y)];

l=[1:nb;(nb+1):2*nb];
ll=l(:);
YDQ=YDQmod(ll,ll);

clear YDQmod;
clear incr;

% The following portion is for running transient stability programme

kcir1=QL0./VL0.^2;
kci1=PL0./VL0.^2;

for i=1:nload
incr=kci1(i)-j*kcir1(i);  
Y((ld(i,1)),(ld(i,1)))=Y((ld(i,1)),(ld(i,1)))+incr;
end

for i=1:ngen
incr=1/(j*gen(i,4));
Y((gen(i,1)),(gen(i,1)))=Y((gen(i,1)),(gen(i,1)))+incr;
end
YES=input('Enter 1 if you want to run transtability programme for network disturbances, otherwise 0:  ');
if(YES)
display('If NO action to be taken, PRESS ENTER for any/every prompt.')

Tfault=input('Fault initiation time (s), Tfault= ');
Tclear=input('Fault Duration,(s) Tclear= ');
fbus=input('Faulted Bus: ');
Line=input('Line(s) to be tripped, [ , ]= ');

if isempty(Tclear)
      Tclear=0;
   end
   if isempty(Tfault)
      Tfault=1000;
   end
   
   Yf=Y;
   if ~isempty(fbus)
      Yf(fbus,fbus)=Yf(fbus,fbus)+100000;
   end

Ypf=Y;

% Tripping of lines to clear fault
if ~isempty(Line)
for tline = Line
   if (tline>nline)
% For transformer tripping
  incr=1/(nt(tline,3)+j*nt(tline,4));  
  incr1=incr/nt(tline,5);
  incr2=(1-nt(tline,5))*incr/(nt(tline,5)*nt(tline,5));
  incr3=(nt(tline,5)-1)*incr/nt(tline,5);
  Ypf((nt(tline,1)),(nt(tline,2)))=Ypf((nt(tline,1)),(nt(tline,2)))+incr1;
  Ypf((nt(tline,2)),(nt(tline,1)))=Ypf((nt(tline,2)),(nt(tline,1)))+incr1;
  Ypf((nt(tline,1)),(nt(tline,1)))=Ypf((nt(tline,1)),(nt(tline,1)))-incr1-incr2;
  Ypf((nt(tline,2)),(nt(tline,2)))=Ypf((nt(tline,2)),(nt(tline,2)))-incr1-incr3;   
else
%For line tripping   
  incr=1/(nt(tline,3)+j*nt(tline,4));
  Ypf(nt(tline,1),nt(tline,1))=Ypf(nt(tline,1),nt(tline,1))-incr-j*nt(tline,5)/2;
  Ypf(nt(tline,2),nt(tline,2))=Ypf(nt(tline,2),nt(tline,2))-incr-j*nt(tline,5)/2;
  Ypf(nt(tline,1),nt(tline,2))=Ypf(nt(tline,1),nt(tline,2))+incr;
  Ypf(nt(tline,2),nt(tline,1))=Ypf(nt(tline,2),nt(tline,1))+incr;  
 end
end   
end
purt_vref=0;
else 
purt_vref= input('Enter 1 if you want to run transtability programme for perturbation of VREF, otherwise 0:  ');
if(purt_vref)
 Tclear = 0;
 Tfault = 1000;
 Yf=Y;
 Ypf=Y;
Gvref = input('Enter the generator number whose Vref needs to be perturbed: ');

if(isempty(ng_static))
   ng_static=0;
end

if (sum(ng_static==Gvref)*(~AVR(Gvref))~=0)
    for gg=1:1:size(exc_static,1)
     if(exc_static(gg,1)==Gvref)
        Gn_vref=gg;
     end
    end
  Vref_static(Gn_vref) = Vref_static(Gn_vref) + 0.01;
end

if(isempty(ng_DC1A))
   ng_DC1A=0;
end

if (sum(ng_DC1A==Gvref)*(~AVR(Gvref))~=0)
    for gg=1:1:size(exc_DC1A,1)
     if(exc_DC1A(gg,1)==Gvref)
        Gn_vref=gg;
     end
    end
  Vref_DC1A(Gn_vref) = Vref_DC1A(Gn_vref) + 0.01;
end

if(isempty(ng_AC4A))
   ng_AC4A=0;
end

if (sum(ng_AC4A==Gvref)*(~AVR(Gvref))~=0)
    for gg=1:1:size(exc_AC4A,1)
     if(exc_AC4A(gg,1)==Gvref)
        Gn_vref=gg;
     end
    end
  Vref_AC4A(Gn_vref) = Vref_AC4A(Gn_vref) + 0.01;
end
end %end of purt_verf
end  %end of YES


if((YES+purt_vref)==0)
purt_Tm= input('Enter 1 if you want to run transtability programme for ramping of Tm, otherwise 0:  ');
if(purt_Tm)
 Tclear = 0;
 Tfault = 1000;
 Yf=Y;
 Ypf=Y;
 gtm = input('Enter the generator number whose Tm needs to be ramped-up/down: ');
if(TURB(gtm)==0) 
   disp('The Turbine-governor on the chosen generator must be disabled using TURB selector')
end
else
   gtm=gen(1,1); % adummy number
   Tclear = 0;
   Tfault = 1000;
   Yf=Y;
   Ypf=Y;
end %end of purt_Tm

else
   gtm=gen(1,1); % adummy number
   purt_Tm=0;
end %end 

%------------------------------------------------------------
%For ramping of Tm with a rate of Tm/s.
   [mm order_gen_tm]=sort([setxor(1:nb,gtm),gtm]);
   t_low=1; %time instant of ramping low
   t_up=11; %time instant of ramping up
 %----------------------------------------------------------  
%Frequency Dependancy of Load parameters:: Not accounted.
%please do not tamper this.
kpf=1.5*0;
kqf=2*0;
if((YES+purt_vref+purt_Tm)==0)
   disp('----You can run the transient stability programme without any disturbance----')
end

