clc

v_rendimiento=zeros(1,6)
v_reg=zeros(1,6)
fp=[0.75,0.80,0.85,0.90,0.95,0.99]

for i=1:1:6

v_factorPotencia=fp(i)
seccion=725
resistencia=0.0424
separacion=0.4
largoLinea=450

angulo=acos(v_factorPotencia)*180/%pi
radio=(sqrt(seccion/%pi))/10^3
largo=largoLinea
e0=8.8542*(10^(-12))
R=resistencia/4
w=2*%pi*50
subconductores=4
G=0
n=4
mT1=[1,complex(0.00266,0.007999);
          0,1]
mT2=[1,complex(0.03,0.12);
          0,1]


if n==1 then reql=0.7788007831*radio                                            //// Condicion para REQL
   elseif n==2 then reql=sqrt(0.778800783*radio*separacion)
   elseif n==3 then reql=(0.778800783*radio*separacion^2)^(1/3)
   elseif n==4 then reql=1.09*((0.778800783*radio*separacion^3)^(1/4))
   else reql=radio
end

deg=(13*13*26)^(1/3)
L=2*(10^(-7))*log(deg/reql)
L=L*(10^3)

if n==1 then reqc=radio*separacion               //// Condicion para REQC
elseif n==2 then reqc=sqrt(radio*separacion)
elseif n==3 then reqc= (radio*separacion^(2))^(1/3)
elseif n==4 then reqc=1.09*((radio*(separacion^(3)))^(1/4))
else reqc=radio
end

C=((2*+%pi)*e0)/(log(deg/reqc))
C=C*(10^(3))
gama=((complex(R,w*L))*(complex(G,w*C)))^(0.5)
zc=((complex(R,w*L))/(complex(G,w*C)))^(0.5)
zb=((500)^2)/(10)
        // Calculo Modelo exacto matriz
a=cosh(gama*largo)
b=(zc*sinh(gama*largo))/zb
ce=((sinh(gama*largo))/zc)*zb
de=cosh(gama*largo)

V2=complex(0.875*cos(0*%pi/180),0.875*sin(0*%pi/180))
I2=complex(1.1428*cos((-angulo*%pi)/180),1.1428*sin((-angulo*%pi)/180))
mL=[a,b;ce,de]
mC=[V2;I2]
multiplicacion=mT1*mL*mT2*mC
V1=multiplicacion(1,1)
I1=multiplicacion(2,1)
disp(V1,'Matriz Transfo 1')
disp(I1,'Matriz Transfo 1')
/*
disp(mT1,'Matriz Transfo 1')
disp(mT2,'Matriz Transfo 2')
disp(factorPotencia,'Phi')
disp(mT1*mL*mT2,'Multiplicacion de Matrices T1,ABCD,T2,V2 e I2')
disp(mL,'Matriz ABCD: ')
disp(det(mL),'Det ABCD: ')
*/

//disp(V1,'Voltaje   1: ')
//disp(I1,'Corriente 1: ')
S1=V1*(complex(real(I1),-imag(I1))) // sin falla
//disp(S1,'Potencia aparente entrada')
aI2=real(I2)
bI2=imag(I2)
aI1=real(I1)
bI1=imag(I1)
reg=((sqrt(real(V1)^2+imag(V1)^2))-(sqrt(real(V2)^2+imag(V2)^2)))/(sqrt(real(V1)^2+imag(V1)^2))
rendimiento=real(V2*complex(aI2,-bI2))/real(V1*complex(aI1,-bI1))*100
v_reg(i)=reg




end


plot2d(fp,v_reg)
xlabel('Regulacion')
xlabel('Factor de potencia')
title('Regulacion vs Factor de potencia')