clear

fid=fopen（‘n.txt‘‘r‘）;
[Ncount]=fscanf（fid‘%d‘inf）;
fclose（fid）;%讀取節點數

fid=fopen（‘t.txt‘‘r‘）;
[TTcount]=fscanf（fid‘%f‘inf）;
fclose（fid）;%讀取時間常數
tt=TT（1）;tmin=TT（2）;tmax=TT（3）;

fid=fopen（‘et.txt‘‘r‘）;
[ETcount]=fscanf（fid‘%f‘inf）;
fclose（fid）;%讀取入侵電壓

fid=fopen（‘ruqin.txt‘‘r‘）;
[Z0count]=fscanf（fid‘%d‘inf）;
fclose（fid）;%讀取入侵波節點及波阻抗
fid=fopen（‘jd_c.txt‘‘r‘）;

[Ccount]=fscanf（fid‘%d?%f‘[2inf]）;%C第一行表示節點，第二行表示該節點接地電容的大小（pf）
[mcnc]=size（C）;
fclose（fid）;%讀取節點對地電容

fid=fopen（‘zhl_z_l.txt‘‘r‘）;
[Zcount]=fscanf（fid‘%d?%d?%f?%f‘[4inf]）;%Z第一二行表示支路的節點，第三行阻抗（歐姆），第四行長度（m）
[mznz]=size（Z）;
fclose（fid）;%讀取支路波阻抗

fid=fopen（‘blq.txt‘‘r‘）;
if?fid==1
????[BLQcount]=fscanf（fid‘%d?%f?%f?%f‘[4inf]）;%第一行：節點；第二行：放電電壓；第三行：C；第四行：a
????[mbnb]=size（BLQ）;
????fclose（fid）;%讀取避雷器
else
????nb=0;
end

T=Z（4:）/300;
Z=[Z;T];%Z第五行為每條支路波每一次傳播的時間us

Y=zeros（NN）;
Y（Z0（1）Z0（1））=Y（Z0（1）Z0（1））+1/Z0（2）;%加入入侵支路阻抗
for?ii=1:nz
????Y（Z（1ii）Z（1ii））=Y（Z（1ii）Z（1ii））+1/Z（3ii）;
????Y（Z（2ii）Z（2ii））=Y（Z（2ii）Z（2ii））+1/Z（3ii）;
end%加入支路阻抗
for?ii=1:nc
????Y（C（1ii）C（1ii））=Y（C（1ii）C（1ii））+2*C（2ii）*（1e-6）/tt;
end%加入接地電容阻抗
Y;
y=inv（Y）;
%…………………………………………………………以上求節點導納矩陣…………………………………………………………%

k0=-1/tt*tmin+1;kmax=1/tt*（tmax-tmin）+1;
%設定t最大小值，為使數組下標為正整數，在t和k間建立映射關系，k=1/tt*（t-tmin）+1即t=（k-1）*tt+tmin，求t=0時k0及tmax時kmax

Iz=zeros（NNkmax）;Ic=zeros（Nkm