clear;clc;
fidin=fopen('電路原始數據.txt','r');%定義輸入，輸出文件
fidout=fopen('輸出結果.txt','w');%初始化
n=fscanf(fidin,"%f",1);
m=fscanf(fidin,"%f",1);
head=zeros(1,m);
tail=zeros(1,m);%記錄回路的首尾信息
R=zeros(1,m);
X=zeros(1,m);
B=zeros(1,m);Y=zeros(n);%節點導納矩陣
P=zeros(1,n+1);
Q=zeros(1,n+1);
U=zeros(1,n);%各點電壓幅值，又潮流計算得出，供精確計算使用，近似計算時認為是1
delta=zeros(1,n);%各點電壓相角，近似計算時認為是0%輸入
for i=1:nU(i)=fscanf(fidin,"%f",1);
end
for i=1:ndelta(i)=fscanf(fidin,"%f",1);delta(i)=delta(i)/180*pi;%化角度為弧度
end
for k=1:m%形成節點導納矩陣i=fscanf(fidin,"%f",1);j=fscanf(fidin,"%f",1);r=fscanf(fidin,"%f",1);x=fscanf(fidin,"%f",1);b=fscanf(fidin,"%f",1);head(k)=i;tail(k)=j;R(k)=r;X(k)=x;B(k)=b;Y(i,i)=Y(i,i)+1/(r+x*1i)+b*1i;Y(j,j)=Y(j,j)+1/(r+x*1i)+b*1i;Y(i,j)=Y(i,j)-1/(r+x*1i);Y(j,i)=Y(j,i)-1/(r+x*1i);
end
%Z=inv(Y);%求導納矩陣的逆矩陣為節點阻抗矩陣
N=fscanf(fidin,"%f",1);
for k=1:Ni=fscanf(fidin,"%f",1);x=fscanf(fidin,"%f",1);U(i)=fscanf(fidin,"%f",1);Y(i,i)=Y(i,i)+1/(x*1i);
end
%Z=inv(Y);%求電源等效后導納矩陣的逆矩陣為節點阻抗矩陣
M=fscanf(fidin,"%f",1);
for k=1:Mi=fscanf(fidin,"%f",1);P(i)=fscanf(fidin,"%f",1);Q(i)=fscanf(fidin,"%f",1);U(i)=fscanf(fidin,"%f",1);Y(i,i)=Y(i,i)+(P(i)-Q(i)*1j)/(U(i)^2);
end
Z=inv(Y);%求電源，負載等效后導納矩陣的逆矩陣為節點阻抗矩陣
K=fscanf(fidin,"%f",1);%獲得短路點編號%短路電流的計算：
Ik_cap=1/abs(Z(4,4));%短路電流幅值的近似值
Ik_theta_cap=-angle(Z(4,4));%短路電流相角的近似值
fprintf(fidout,"短路電流的近似值為：%f，%f°\n",Ik_cap,Ik_theta_cap/pi*180);Ik=abs(U(4)./Z(4,4));%短路電流幅值的精確值
Ik_theta=(delta(4)-angle(Z(4,4)));%短路電流相角的精確值
fprintf(fidout,"短路電流的精確值為：%f，%f°\n",Ik,Ik_theta/pi*180);Ik_erro=Ik-Ik_cap;
Ik_theta_erro=Ik_theta-Ik_theta_cap;
fprintf(fidout,"短路電流的誤差為：%f，%f°\n",Ik_erro,Ik_theta_erro/pi*180);Uk_cap=zeros(1,n);%短路后各個節點的短路電壓近似值：
Uk=zeros(1,n);%短路后各個節點的短路電壓精確值：
for i=1:nUk_cap(i)=1-Z(i,4)*(Ik_cap*cos(Ik_theta_cap)+1i*Ik_cap*sin(Ik_theta_cap));Uk(i)=U(i)*(cos(delta(i))+1i*sin(delta(i)))-Z(i,4)*U(4)*(cos(delta(4))+1i*sin(delta(4)))/Z(4,4);
end
Uk_theta_cap=zeros(1,n);
Uk_theta=zeros(1,n);
for i=1:nfprintf(fidout,"%d節點短路電壓的近似值為：%f，%f°\n",i,abs(Uk_cap(i)),angle(Uk_cap(i))/pi*180);fprintf(fidout,"%d節點短路電壓的精確值為：%f，%f°\n",i,abs(Uk(i)),angle(Uk(i))/pi*180);fprintf(fidout,"%d節點短路電壓的誤差值為：%f\n",i,abs(Uk(i))-abs(Uk_cap(i)));fprintf(fidout,"\n");
end%各支路的短路電流精確值計算（經分析，短路電壓近似值已有較大誤差，故不再計算電路近似值）
Ik_all=zeros(m);
for i=1:mIk_all(head(i),tail(i))=(Uk(head(i))-Uk(tail(i)))/(R(i)+1i*X(i));fprintf(fidout,"%d-%d之間的短路電流的精確值為：%f，%f°\n",head(i),tail(i),abs(Ik_all(head(i),tail(i))),angle(Ik_all(head(i),tail(i)))/pi*180);
endfclose(fidin);
fclose(fidout);%關閉文件
