オープン-ショート法による,通信ケーブル(不平衡線路)の一次定数測定用計算プログラム

通信ケーブル(不平衡線路)の一次定数を求める場合,適当に切ったケーブルの終端を開放,短絡して入力インピーダンスを測定し,その値から一次定数を計算により求めることができます.下記はそのためのプログラムです.参考文献は,下記のものです.
データは,周波数,振幅,位相の順で空白区切りあるいはタブ区切りのデータとします.

(参考文献)R.Croze,L.Simon著,林 憲一訳,有線電話伝送工学-線路理論-,学献社,pp.98-101.

(プログラム)

Length=1;                                      # ケーブル長の設定
Open=load("open.txt");                         # 終端開放時のデータの読み込み.周波数,振幅,位相の順
Short=load("short.txt");                       # 終端短絡時のデータの読み込み.データは空白区切り
N= size (Open,1);                              # データの行数の読み込み
for i=1:N
  f(i)=Open(i,1);                              # 周波数の読込
  Zf(i)=Open(i, 2).*exp(1i*Open(i,3)*pi/180);  # データを複素数に変換
  Zs(i)=Short(i,2).*exp(1i*Short(i,3)*pi/180); # データを複素数に変換
  Theta_s(i)=Short(i,3)*pi/180;                # 角度をラジアンに変換
  Theta_f(i)=-Open(i,3)*pi/180;                # 角度をラジアンに変換
  #測定したのはZfで参考文献では,1/Zf=Yfの位相をθfとしているので,マイナスが必要
  fai(i)=(Theta_s(i)+Theta_f(i))/2;            # φ=(θs+θf)/2の計算
  cosfai(i)=cos(fai(i));                       # cos(φ)
  sinfai(i)=sin(fai(i));                       # sin(φ)
  guzai(i)=(Theta_s(i)-Theta_f(i))/2;          # ξ=(θs-θf)/2の計算
  cosguzai(i)=cos(guzai(i));                   # cos(ξ)
  singuzai(i)=sin(guzai(i));                   # sin(ξ)
  ZsYf(i)=sqrt(abs(Zs(i)./Zf(i)));             # √|ZsYf|=√|Zs/Zf| の計算
  ZsZf(i)=sqrt(abs(Zs(i).*Zf(i)));             # √|ZsZf| の計算
  Alpha(i)=1/2/Length*atanh((2*ZsYf(i).*cosfai(i)) ./(1+ZsYf(i)));                  # αの計算
  Beta(i)=1/2/Length*atanh((2*ZsYf(i).*sinfai(i)) ./(1-ZsYf(i)));                   # βの計算
  R(i)=ZsZf(i).*(Alpha(i).*cosguzai(i)-Beta(i).*singuzai(i));                       # Rの計算
  L(i)=(1/(2*pi*f(i))).*ZsZf(i).*(Alpha(i).*singuzai(i).+Beta(i).*cosguzai(i));     # Lの計算
  G(i)=(1/ZsZf(i)).*(Alpha(i).*cosguzai(i)+Beta(i).*singuzai(i));                   # Gの計算
  C(i)=(1/(2*pi*f(i))).*(1/ZsZf(i)).*(-Alpha(i).*singuzai(i)+Beta(i).*cosguzai(i)); # Cの計算
endfor
#
fid=fopen('kekka.txt','w');                                                         # ファイルのオープン
for i=1:N
  fprintf(fid,'%8.2f %12.3e %12.3e %12.3e %12.3e\n',f(i),R(i),L(i),G(i),C(i));      # ファイル出力
endfor
fclose(fid);                                                                        # ファイルのクローズ