% Programa para la construccion calculo de corrientes de cortocircuito
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%             Z B U S              I C C          %%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Este programa emplea el algoritmo de construccion de la amtriz impedancia
% de barra de secuencia positiva para el calculo de falla por cortocircito
% trifasico.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ZbusIcc, Version 1.02
% Disenado y Programado por: Francisco M. Gonzalez-Longatt
% fglongatt@cantv.net  fglongatt@ieee.org
% http://www.giaelec.org
% Derechos de autor reservados
% Solo para ser usado para fines de investigacion o academicos
% Version inicial: Abril, 2003.
% Revision: Francisco M. Gonzalez-Longatt, Enero 2008
%------------------------------------------------------------------------
clear
Z=[]
op=10;
while (op ~= 7)
   clc
   disp(' Programa de matriz Impedancia de Barra')
   disp('        ')
   disp('----------------------------------------------------')
   disp(' [1] Construir Matriz Impedancia de Barra')
   disp(' [2] Modificar Matriz')
   disp(' [3] Calcular corrientes de cortocircuito')
   disp(' [4] Perfil de Tensiones')
   disp(' [5] Contribuciones de Corrientes')
   disp(' [6] Eliminar una linea')
   disp(' [7] Salir')
   disp('----------------------------------------------------')
   op= input(' Indique la operacion :');
   if (op== 1) % Contruccion de la matriz Impedancia de Barra
     clc
     disp(' CONSTRUCCION DE LA MATRIZ IMPEDANCIA DE BARRA')
     disp('----------------------------------------------------')
     disp(' ')
     fg=input('[1] Entrada por Teclado  [2] Entrada por archivo   ');
     if (fg==2)
        disp(' ')
        loadfile=input(' Indique nombre del archivo :','s');
        %load(loadfile,'operacion','acople','ACOPLE');
        load(loadfile,'operacion');
        disp('  ')
        disp(' Archivo Cargado...')
        disp(' ')
        sss=input(' Presione <ENTER>');
        [m,n]=size(operacion);
        %[n_acople,n_acople]=size(acople);
     end
     disp(' ')
     flag=input(' Mostrar Construccion de la matriz [0] No [1] Si  ');
     if (fg==1)
        disp('  ')
        n=input(' Numero de Ramas: ');
        operacion=struct('Barra_Inicio',{},'Barra_Final',{},'Zpq',{});
        disp('   ')
        disp ( ' Introduccion de elementos')
        disp(' ')
        % Entradas de los elementos
        for k=1:n,
           nx=[' Elemento ' int2str(k) ];    
           disp('----------------------------------------------------')
           disp(nx)
           disp('----------------------------------------------------')
           operacion(k).Barra_Inicio = input(' Barra Inicio :');
           operacion(k).Barra_Final = input(' Barra Final :');
           nx=[' Zrama[' int2str(operacion(k).Barra_Inicio) ',' int2str(operacion(k).Barra_Final) ']=' ];
           operacion(k).Zpq=input(nx);
        end
        %disp('   ')
        %fl3=input(' Hay Lineas acopladas [0] No [1] Si  ');
        %if (fl3==1)
        %    disp('  ')
        %    n_acople=input(' Numero de Lineas acopladas : ');    
        %    disp('  ')
        %    acople=struct('Inicio_P',{},'Final_Q',{},'Inicio_R',{},'Final_S',{} ,'Mpq',{});
        %    for t=1:n_acople,
        %        nx=[' Linea Acoplada ' int2str(t) ];    
        %        disp('----------------------------------------------------')
        %        disp(nx)
        %        disp('----------------------------------------------------')
        %        acople(t).Inicio_P = input(' Barra Inicio :');
        %        acople(t).Final_Q = input(' Barra Final :');
        %        disp(' Acoplada con...')
        %        acople(t).Inicio_R = input(' Barra Inicio :');
        %        acople(t).Final_S = input(' Barra Final :');
        %        nx=[' Zmutua[' int2str(acople(t).Inicio_P) '-' int2str(acople(t).Final_Q) ',' int2str(acople(t).Inicio_R) '-' int2str(acople(t).Final_S) ']= ' ];
        %        acople(t).Mpq=input(nx);
        %        ACOPLE(acople(t).Inicio_P,acople(t).Final_Q,acople(t).Inicio_R,acople(t).Final_S)=acople(t).Mpq;
        %    end
        %end
        disp(' Desea Guardar la Lista de Elementos ?')
        fg1=input('[0] No  [1] Si   ');
        if (fg1==1)
           disp(' ')
           savefile=input(' Indique nombre del archivo : ','s')
           %save(savefile,'operacion','acople','ACOPLE');
           save(savefile,'operacion');
           disp('  ')
           disp(' Archivo Guardado...')
           disp(' ')
           sss=input(' Presione <ENTER>');
        end
     end      
     clc
     disp('  Tabla de Construccion')
     disp('----------------------------------------------------')
     
     for k=1:n,
        bus_new=operacion(k).Barra_Final;
        bus_old=operacion(k).Barra_Inicio;
        %Averigua si la linea tiene acople
        %este proecedimiento busca el acople
        %n_bacoples=0;
        %si_acople=0;
        %for t=1:n_acople,
        %  if(acople(t).Inicio_R==bus_old) & (acople(t).Final_S==bus_new)
        %     si_acople=1;
        %     n_bacoples=n_bacoples+1; %dice cuantas barras hay acopladas
        %     hay_acople(n_bacoples)=t;
        % end
        %end
        %for t=1:n_bacoples,
        %   for s=1:n_bacoples;
               
        %   end    
        %end
        
        
        if operacion(k).Barra_Inicio ==0 %No importa la barra final
           % Agregar un elemento entre referencia y barra nueva    
           Z(:,bus_new)=0;
           Z(bus_new,:)=0;
           Z(bus_new,bus_new)=operacion(k).Zpq;
           elemento=[' Elemento ' int2str(0) '-' int2str(bus_new) ' Tipo 1' ];
           disp(elemento)
        end
        [dimension1,dimension2]=size(Z);
        if (operacion(k).Barra_Inicio)>0 
              
            if (operacion(k).Barra_Inicio <= dimension2) 
              % Si se cumple lo anterior se trata de una barra ya existente 
              if (operacion(k).Barra_Final <= dimension2)
                 % Se trata de una Barra que cierra lazo
                 lazo=dimension2+1;
                 Z(:,lazo)=Z(:,bus_new)-Z(:,bus_old);
                 Z(lazo,:)=Z(bus_new,:)-Z(bus_old,:);
                 Z(lazo,lazo)=Z(bus_old,bus_old)+Z(bus_new,bus_new)-2*Z(bus_old,bus_new)+operacion(k).Zpq;
                 nx=[' Elemento ' int2str(bus_old) '-' int2str(bus_new) ' Tipo 3 + Kron'];
                 disp(nx)
                 if (flag==1)
                     disp(Z)
                 end
                 % Se debe aplicar Kron
                 for s = 1:lazo-1,
                  for t=1:lazo-1,
                   Z(s,t)=Z(s,t)-Z(s,lazo)/Z(lazo,lazo)*Z(lazo,t); 
                  end
                 end 
                 Z(:,lazo)=[]; 
                 Z(lazo,:)=[]; 
              end
           end
        end
        if  (operacion(k).Barra_Final > dimension2)
           % La barra inicial es una barra nueva
           if (operacion(k).Barra_Inicio <= dimension2)
              % Se trata de Barra Inicial ya en la matriz y una barra nueva
              %if (si_acople ~=0) 
                  
                  
              %end
              
              
              %if (si_acople==0) 
                Z(:,bus_new)=Z(:,bus_old);
                Z(bus_new,:)=Z(bus_old,:);
                Z(bus_new,bus_new)=Z(bus_old,bus_old)+operacion(k).Zpq;
                nx=[' Elemento ' int2str(bus_old) '-' int2str(bus_new) ' Tipo 2'];
                disp(nx)
                %end
           end  
        end  
     if (flag==1)
        disp(Z)
        disp('-----------------------------------------------------')
     end
     end
     disp('-----------------------------------------------------')
     disp(' Matriz Impedancia de Barra')
     disp('    ')
     Z
     [dimension1,dimension2]=size(Z);
     disp('-----------------------------------------------------')
     nx=[' Orden de la matriz es :' int2str(dimension1) 'x' int2str(dimension2) ];
     disp(nx)
     disp('    ')
     sss=input(' Presione <ENTER>');
   end
   if (op==2)% Modificar matriz
   end
   if (op==3) % Calculo de corrientes de cortocircuito
       clc
       disp('  CALCULO DE CORRIENTES DE CORTOCIRCUITO')
       disp('-----------------------------------------------------')
       [dimension1,dimension2]=size(Z);
       nx=[' Barras del Sistema ' int2str(dimension2) ];
       disp(nx);
       disp(' [1] Todas las Barras')  
       disp(' [2] Solo una Barra')
       disp('-----------------------------------------------------')
       op2=input(' Opcion : ');
       disp(' ')
       E=input(' Introduzca la tension previa a la falla [por unidad] : ');
       disp(' ')
       if (op2==1)
         for k=1:dimension2,
             I(k)=E/Z(k,k);
         end
         disp(' CORRIENTES DE FALLA [EN P.U.]')
         disp('-----------------------------------------------------')
         for k=1:dimension2,
           nx=[' Barra ' int2str(k) ':  I' int2str(k) '= ' num2str(I(k)) 'p.u'];
           disp(nx);
         end
         
       end
       if (op2==2)
         k=input(' Introduzca el Numero de la Barra : ');
         I(k)=E/Z(k,k);
         disp(' CORRIENTES DE FALLA [EN P.U.]')
         disp('-----------------------------------------------------')
         nx=[' Barra ' int2str(k) ':  I' int2str(k) '= ' num2str(I(k)) ' p.u'];
         disp(nx)
       end
       disp(' ')
       sss=input(' Presione <ENTER>');
   end % Final eleccion 3
   if (op==4)
       %calculo de perfil de tensiones    
       clc
       disp('  CALCULO DE CORRIENTES DE CORTOCIRCUITO')
       disp('-----------------------------------------------------')
       disp( ' ')
       k=input(' Introduzca el numero de la Barra Fallada : ');
       disp('  ')
       E=input(' Introduzca la tension previa a la falla [por unidad] : ');
       disp('  ')
       [dimension1,dimension2]=size(Z);
       for p=1:dimension2,
           V(p)=E-(Z(p,k)*E/Z(k,k));   
       end
       disp(' TENSIONES EN LAS BARRAS [EN P.U.]')
       disp('-----------------------------------------------------')
       for p=1:dimension2,
           nx=[' Barra ' int2str(p) ':  E(' int2str(p) ')= ' num2str(V(p)) ' p.u'];
           disp(nx);
       end
   disp(' ')
   sss=input(' Presione <ENTER>');
   end % Final eleccion 4
   
   if (op==5) 
       % Contribuciones de corrientes
       clc
       disp('  CALCULO DE CONTRIBUCIONES DE CORRIENTES DE CORTOCIRCUITO')
       disp('-----------------------------------------------------')
       disp( ' ')
       [dimension1,dimension2]=size(Z);
       k=input(' Introduzca el numero de la Barra Fallada : ');
       disp('  ')
       E=input(' Introduzca la tension previa a la falla [por unidad] : ');
       for p= 1:dimension2,
           for q= 1:dimension2,
               % se debe buscar el Zpq
               if (q ~= p)
                 for t=1:n,
                   pp = operacion(t).Barra_Inicio;
                   qq = operacion(t).Barra_Final;
                   if ((pp==p) & (qq==q))
                     Zpq=operacion(t).Zpq;    
                   end
                 end
                 if (p==2) & (q==3)
                 Zpq   =operacion(t).Zpq;   
                 end    
                 I(p,q)=(Z(q,k)-Z(p,k))*E/(Z(k,k)*Zpq);    
               end
           end
           
       end
       disp(' CONTRIBUCIONES DE CORRIENTES [EN P.U.]')
       disp('-----------------------------------------------------')
       for p=1:dimension2,
           for q=1:dimension2,
              if (q ~= p) & (abs(I(p,q)>0))
              nx=[' Barra ' int2str(p) '- Barra ' int2str(q) ':  I(' int2str(p) ',' int2str(q)  ')= ' num2str(I(p,q)) ' p.u'];
              disp(nx);
              end
           end
       end
       disp(' ')
       sss=input(' Presione <ENTER>');
   end %Final opcion 5
   if (op==6)
     clc
     disp('  SACAR UNA LINEA DE SERVICIO ')
       disp('-----------------------------------------------------')
       disp( ' ')
       p=input(' Barra de Inicio : ');
       disp('  ')
       q=input(' Barra Final: ');
       for t=1:n,
         pp = operacion(t).Barra_Inicio;
         qq = operacion(t).Barra_Final;
         if ((pp==p) & (qq==q))
           Zpq=operacion(t).Zpq;    
         end
       end
       ZZ=Z(p:q,p:q);
       [dime1,dime2]=size(ZZ)
       lazo=dime2+1;
       bus_old=1
       bus_new=dime2;
       ZZ(:,lazo)=ZZ(:,bus_new)-ZZ(:,bus_old);
       ZZ(lazo,:)=ZZ(bus_new,:)-ZZ(bus_old,:);
       ZZ(lazo,lazo)=ZZ(bus_old,bus_old)+ZZ(bus_new,bus_new)-2*ZZ(bus_old,bus_new)+(-1)*Zpq;
       for s = 1:lazo-1,
         for t=1:lazo-1,
           ZZ(s,t)=ZZ(s,t)-ZZ(s,lazo)/ZZ(lazo,lazo)*ZZ(lazo,t); 
         end
       end 
       ZZ(:,lazo)=[]; 
       ZZ(lazo,:)=[]; 
       Z(p:q,p:q)=ZZ;
       disp('-----------------------------------------------------')
       disp(' Matriz Impedancia de Barra')
       disp('    ')
       Z
       [dimension1,dimension2]=size(Z);
       disp('-----------------------------------------------------')
       nx=[' Orden de la matriz es :' int2str(dimension1) 'x' int2str(dimension2) ];
       disp(nx)
       disp('    ')
       sss=input(' Presione <ENTER>');
             
   end %Final de la opcion 6
   
end