%図13.11
%
clear all; close all;

%伝達関数の分子・分母多項式を与える
num = [ 1 ]; %分子多項式
den = [ 1 2 1 ]; %分母多項式

%コントローラの分子・分母多項式を与える.
numc = [ 0.5 ]; %K_i=0.5の場合の分子多項式
denc = [ 1 0 ]; %分母多項式

%制御対象とコントローラの伝達関数表現を与える
sys = tf( num, den ); %制御対象の伝達関数表現
c = tf( numc, denc ); %K_i=0.5の場合のコントローラの伝達関数表現

%開ループ伝達関数を求める
sysL = c*sys; %開ループ伝達関数

%図13.11(a)のプロット
figure(1) %図のウィンドウを開く
nyquist(sysL) %ベクトル軌跡のプロット(ナイキスト線図)
xlim([-2,0]) %横軸(実軸)の範囲の指定
ylim([-1 1]) %縦軸(虚軸)の範囲の指定

%図13.11(b)のプロット
figure(2) %図のウィンドウを開く
%角周波数の範囲を指定
w = logspace(-1, 1, 1000); %対数的に等間隔なベクトルの生成(10^{-2}から10^{2}で1000点)

%ゲインと位相の計算
[gain1 phase1 w1] = bode(sysL, w); %ゲインと位相
gain1 = gain1(:); %ゲインの配列を1次元ベクトルに変換
gain1dB = 20*log10(gain1); %K_{i}=0.5の場合のゲインをデシベル値に変換
phase1 = phase1(:); %K_{i}=0.5の場合の位相の配列を1次元ベクトルに変換

%図13.11のプロット
subplot(2,1,1) %複数の図を並べるためのコマンド.2行2列の1行目という意味
semilogx(w1, gain1dB); %ゲイン線図をプロット
xlim([10^(-1),10^1]) %横軸(角周波数)の範囲の指定
ylim([-30 40]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1)]) %横軸の目盛りの値の設定
yticks([-20 0 20 40]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
ylabel('Gain [dB]'); %縦軸のラベル表示

subplot(2,1,2) %複数の図を並べるためのコマンド.2行1列の2行目という意味
semilogx(w1, phase1); %位相線図をプロット
xlim([10^(-1),10^1]) %横軸(角周波数)の範囲の指定
ylim([-270 -90]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2) 10^(3)]) %横軸の目盛りの値の設定
yticks([-270 -225 -180 -135 -90]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
ylabel('Phase [deg]'); %縦軸のラベル表示