function [Ucal, Snudo, Srama, Irama] = flujo(fdatos)
% [Ucal, Snudo, Srama] = flujo(fdatos);
% Flujo de Cargas Newton Raphson.
%  
% Entrada:
%    fdatos: Nombre del fichero de datos (e.j. 'red.txt')
% Salida:
%    Ucal : Tensiones nodales complejas (en p.u.)
%    Snudo: Potencia aparente compleja inyectada en los nudos (en p.u.)
%    Srama: Potencia aparente compleja en ramas del nudo n1 al nudo n2 (en p.u.)
%    Irama: Intensidad compleja en ramas del nudo n1 al nudo n2 (en p.u.)
%
%   <Ejemplo de formato de fichero de datos>
%   Copiar desde <INICIO> a <FIN> en un fichero nuevo (p.e. 'red.txt') y
%   quitarle un '%' del incio de cada linea.
%
%# <INICIO>
%#  DATOS DE SISTEMA
%#  ( Lo que va precedido de # o % es un comentario)
%#  ¡¡¡¡ TODOS los datos estan en p.u !!!!
%#  TIPO DE NUDO:  1 Balance,  2 PQ,  3  PV
%# -----------------------------------------------
%
%# Potencia base en MVA
%SB  100
%
%# Datos de las lineas
%LI
%# ---------------------------------------------
%# nº  n1   n2     R       X        B/2      t
%# ---------------------------------------------
%  1   1    2    0.0020   0.0250  0.0150    1.0
%
%# Datos de las nudos
%NU
%# --------------------------------------------------------------
%#     n   Pd     Qd     Pg       Qg      U       G     B    tipo
%# --------------------------------------------------------------
%      1  00.00  0.000  0.000   0.0000  1.0125   0.0   0.0    1
%      2   1.00  1.000  94.06   0.0000  1.0125   0.0   0.0    3
%#<FIN>
%
%  (c)  Prácticas de Lineas y Redes      1999
%       Dpto. de Ingeniería Eléctrica.
%       Universidad de Vigo.
%
% Revisado: Camilo, 2004

disp('---');
disp('--- Flujo de cargas NEWTON-RAPHSON ---');
disp('---');
MaxIter = 5000; %Nº maximo de iteraciones
i = sqrt(-1);
leoli=0;, leonodo=0;
kramas=0;, knudos=0;

% Lectura de datos -------------------------------------
fid = fopen(fdatos);
while 1
   line = fgetl(fid);
   if ~isstr(line), break, end %finaliza el bucle
   if ~all(isspace(line))
       if  length(line) > 1 
           %Eliminamos espacios al principio y al final de la linea
           [r,c] = find( (line~=0) & ~isspace(line) ); 
           line = line(min(c):max(c));
           line = upper(line);
           %Interpretamos el fichero de datos
           if  (line(1) ~= '#') & (line(1) ~= '%') %Es un comentario
               if min(findstr(line,'LI')) == 1 %Lineas
                   leoli=1;
                   leonodo=0;
               elseif min(findstr(line,'NU')) == 1 %Nudos
                   leonodo=1;
                   leoli=0;
               elseif min(findstr(line,'SB')) == 1 %Potencia Base
                   D=sscanf(line,'%s %f',2);
                   SBase = D(1,2);
               else 
                   if leoli == 1
                       kramas=kramas+1;
                       D=sscanf(line,'%d %d %d %f %f %f %f',7);
                       for k=1:7, DRAMAS(kramas,k)=D(k); end; 
                   end
                   if leonodo == 1
                       knudos=knudos+1;
                       D=sscanf(line,'%d %f %f %f %f %f %f %f %d',9);
                       for k=1:9, DNUDOS(knudos,k)=D(k); end
                   end
               end
           end
       end
   end
end %while
fclose(fid);

% Calculamos el número máximo de nudos ( nmax ) --------
nmax=0;
for k=1:kramas
    if DRAMAS(k,2) > nmax,  nmax=DRAMAS(k,2); end
    if DRAMAS(k,3) > nmax,  nmax=DRAMAS(k,3); end
end

% Inicializamos la Ynod --------------------------------
Ynod = zeros(nmax);
for k=1: kramas
    n1=DRAMAS(k,2);, n2=DRAMAS(k,3);
    y=1/(DRAMAS(k,4)+DRAMAS(k,5)*i); % 1/(R+i·X)
    Ynod(n1,n1)=Ynod(n1,n1)+y+DRAMAS(k,6)*i; % ... + i·B
    Ynod(n2,n2)=Ynod(n2,n2)+y*(1/(DRAMAS(k,7)*DRAMAS(k,7)))+DRAMAS(k,6)*i;
    Ynod(n1,n2)=Ynod(n1,n2)-y*( 1.0/DRAMAS(k,7) );
    Ynod(n2,n1)=Ynod(n2,n1)-y*( 1.0/DRAMAS(k,7) );
end
for k=1: knudos
    n=DNUDOS(k,1);
    Ynod(n,n)=Ynod(n,n)+(DNUDOS(k,7)+DNUDOS(k,8)*i);
end

disp('===== Matriz admitancias nodales =====')
for k=1:nmax
    fprintf('[%d]> ',k);
    for m=1:nmax
       fprintf('%7.3f%+7.3fj  ',real(Ynod (k,m)),imag(Ynod (k,m)));
    end
    fprintf('\n');
end
fprintf('\n');

%  Vector de Potencias y Tensiones Especificadas -------
Pesp = zeros(nmax,1);
Qesp = zeros(nmax,1);
Uesp = ones(nmax,1);
Fesp = zeros(nmax,1);
Tesp = ones(nmax,1)*2;

for k=1:knudos
    n=DNUDOS(k,1);
    Tesp(n)= DNUDOS(k,9);
    if Tesp(n) == 2 | Tesp(n) == 3, Pesp(n)= (DNUDOS(k,4)-DNUDOS(k,2)); end
    if Tesp(n) == 2,                Qesp(n)= (DNUDOS(k,5)-DNUDOS(k,3)); end
    if Tesp(n) == 1 | Tesp(n) == 3, Uesp(n)= (DNUDOS(k,6)); end
end;
Sesp = Pesp + Qesp * i;


% Vector de Tensiones Inicial --------------------------
disp('===== Bucle de iteraciones =====')
Ucal = Uesp .* ( cos(Fesp) + i*sin(Fesp) );
Uread=Ucal;

%-------------------------------------------------------
% INICIO BUCLE DE ITERACIONES --------------------------
t0 = cputime;
itera=0; %Contador de iteraciones
while  1
    Ical  = Ynod * Ucal;
    Scal  = Ucal .* conj(Ical);
    DScal = Sesp - Scal;

    % Calculamos Delta P y Deltas Q --------------------
    ks=0;
    for k=1: nmax
        if Tesp(k) ~= 1,  ks=ks+1;, DS(ks)=real( DScal(k) ); end
    end
    for k=1: nmax
        if Tesp(k) ~= 1 & Tesp(k) ~=3,  ks=ks+1;, DS(ks)=imag( DScal(k) ); end
    end
    % --------------------------------------------------

    ERROR = DS*DS';
    fprintf('*** Iteración: %d -> error: %G\n', itera, ERROR );

    % CONDICIONES DEL FIN DE BUCLE
    if isnan(ERROR)>0.5 | isinf(ERROR)>0.5
        error('Hay problemas de convergencia');
        break
    end
    if itera > MaxIter
        error(['Se ha excedido el numero maximo de iteraciones: ' MaxIter]);
        break
    end
    if ERROR < 1E-10,  break, end   %% Finaliza el bucle

    itera=itera+1;

    % Calculo del Jacobiano ----------------------------
    Jacob = zeros(2*nmax, 2*nmax);
    for k=1: nmax
        for m=1: nmax
            Gkm = real( Ynod(k,m) );, Bkm = imag( Ynod(k,m) );
            Pk= real( Scal(k) );
            Qk = imag( Scal(k) );
            Pm= real( Scal(m) );
            Qm = imag( Scal(m) );
            Uk= abs( Ucal(k) );
            Um = abs( Ucal(m) );
            Fkm = angle( Ucal(k) ) - angle( Ucal(m) );
            if k == m
                Jacob(k,m)            = - Qk - Bkm * Uk*Uk;
                Jacob(k,m+nmax)       =   Gkm*Uk + Pk/Uk;
                Jacob(k+nmax, m)      =   Pk - Gkm * Uk*Uk;
                Jacob(k+nmax, m+nmax) = - Bkm*Uk + Qk/Uk;
            else
                Jacob(k,m)            =   Uk * Um * ( Gkm*sin(Fkm) - Bkm*cos(Fkm) );
                Jacob(k,m+nmax)       =        Uk * ( Gkm*cos(Fkm) + Bkm*sin(Fkm) );
                Jacob(k+nmax, m)      = - Uk * Um * ( Gkm*cos(Fkm) + Bkm*sin(Fkm) );
                Jacob(k+nmax, m+nmax) =        Uk * ( Gkm*sin(Fkm) - Bkm*cos(Fkm) );
            end       
        end
    end

    for k=nmax: -1 : 1
        if Tesp(k)== 1 | Tesp(k) == 3, Jacob(k+nmax,:)=[];, Jacob(:,k+nmax)=[];, end
    end
    for k=nmax: -1 : 1
        if Tesp(k)==1, Jacob(k,:)=[];, Jacob(:,k)=[];, end
    end
    % --------------------------------------------------

    DU  = inv(Jacob) * DS';
    
    %DPQ = DS; display(DPQ);
    %J=Jacob;  display(J);
    %DTU = DU; display(DTU);

    modUcal=abs(Ucal);
    fasUcal=angle(Ucal);

    ku=0;
    for k=1: nmax
        if Tesp(k) ~= 1, ku=ku+1; fasUcal(k)= fasUcal(k) + DU(ku);, end
    end
    for k=1: nmax
        if Tesp(k) ~= 1 & Tesp(k) ~= 3, ku=ku+1;, modUcal(k)= modUcal(k) + DU(ku); end
    end

    Ucal = modUcal.*(cos(fasUcal)+i*sin(fasUcal));

end
t1 = cputime;
fprintf('------ CONVERGE en %d iteraciones y en %f segundos -------\n',itera, t1-t0);
% FIN DEL BUCLE DE ITERACIONES -------------------------
%-------------------------------------------------------


% Calculo de la potencia -------------------------------
Icalculada = Ynod * Ucal;

Iconjugada = conj(Icalculada);
Potencia   = Iconjugada .* Ucal;
Snudo = Potencia; %Potencia NUDOS
Qcalculada=imag(Potencia);

% Resultados -------------------------------------------
fprintf('\n===== Resultados de nudos\n');
fprintf(' ---------------------------------------------------------------------- \n');
fprintf(' nudo   U      fase(rad)         Pcal		Qcal			 \n');
fprintf(' ---------------------------------------------------------------------- \n');
for k=1: nmax
    fprintf(' %d  %f  %f   %f   %f\n',k, abs(Ucal(k)), angle(Ucal(k)) ,real(Potencia(k)), imag(Potencia(k)) );
end

fprintf('\n===== Resultados de líneas\n');
fprintf(' --------------------------------------------- \n');
fprintf(' nº  n1  n2   P12   Q12   Pp    Qp    I12(p.u) \n');
fprintf(' ----------------------------------------------\n');

Srama = zeros(kramas,1);
Irama = zeros(kramas,1);
for k=1: kramas
    nram= DRAMAS(k,1);
    n1 =DRAMAS(k,2);, n2=DRAMAS(k,3);
    r  =DRAMAS(k,4);, x=DRAMAS(k,5);, b=DRAMAS(k,6);, t  =DRAMAS(k,7);
    y  = 1/(r+x*i);
    y1 = y/t;,  y2 = (y*(t-1))/t;,  y3 = (y*(1-t)) / (t*t);
    Iser=(Ucal(n1)-Ucal(n2))*y1;
    In12= Iser + (Ucal(n1)*y2) + Ucal(n1)*j*b;
    In21=-Iser + (Ucal(n2)*y3) + Ucal(n2)*j*b;
    S12= Ucal(n1) * conj(In12);
    S21= Ucal(n2) * conj(In21);
    Per = real(S12+S21);
    Qer = imag(S12+S21);
    Srama(k) = S12;
    Irama(k) = In12;
    fprintf('%d  %d  %d   %f   %f   %f   %f   %f\n',k,n1,n2,real(S12), imag(S12), Per, Qer, abs(In12) );
end