線形システムのボードとナイキスト周波数応答

• 1行，19行，33行: 関数の始まり。
• 3行，5行，23行，24行，37行，38行: 局所変数の宣言。
• 6行: 関数の宣言。
• 8行: 多項式変数の定義。
• 9行: 3次実有理関数を入力。
• 11行: 対数空間ベクトル。
• 12行，13行: 複数の戻り値を返す関数Bode()とNyquist()を呼ぶ。
• 20行，21行，34行，35行: 関数の引数の宣言。
• 26行，40行: 虚数単位の定義。
• 27行，41行: 有理関数配列の評価。(配列の各成分毎に評価する)
• 28行，29行: 配列演算。
• 30行，42行: 2つの成分を持つリストを返す。

   Func void main()
   {
     Polynomial s;
     Rational G;
05:  Array w, ga, ph, re, im;
     List Bode(...), Nyquist(...);
   
     s = $;                      // 多項式変数
     G = (s^2 + 3*s + 4)/(s^3 + 5*s^2 + 6*s + 7);
10:                              // 伝達関数
     w = logspace(0.01, 100.0);  // 周波数の範囲
     {ga, ph} = Bode(G, w);      // ゲインと位相
     {re, im} = Nyquist(G, w);   // 実部と虚部

15:  print [[w][re][im]] >> "nyquist.mat"; // データ保存
     print [[w][ga][ph]] >> "bode.mat";    // データ保存
   }
   
   Func List Bode(G, w)
20:  Rational G;
     Array w;
   {
     Complex j;
     Array data, ga, ph;
25:   
     j = (0,1);                    // 虚数単位
     data = eval(G, j*w);          // G(j*w)
     ga = 20.0 * log10(abs(data)); // ゲイン [dB]
     ph = arg(data) * 180.0 / PI;  // 位相  [degree]
30:  return {ga, ph};
   }
   
   Func List Nyquist(G, w)
     Rational G;
35:  Array w;
   {
     Complex j;
     Array data;

40:  j = (0,1);                    // 虚数単位
     data = eval(G, j*w);          // G(j*w)
     return {Re(data), Im(data)};  // 実部と虚部
   }