%図14.3
%
clear all; close all;

%伝達関数の分子・分母多項式を与える
num = [1]; %分子多項式
den1 = [1 1]; %P_1(s)の分母多項式
den2 = [1 0]; %P_2(s)の分母多項式

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

%制御対象とコントローラの伝達関数表現を与える
sys1 = tf( num, den1 ); %P_1(s)の伝達関数表現
sys2 = tf( num, den2 ); %P_2(s)の伝達関数表現
c = tf( numc, denc); %コントローラの伝達関数表現

%開ループ伝達関数を求める
sysL1 = c*sys1; %P_1(s)の場合の開ループ伝達関数
sysL2 = c*sys2; %P_2(s)の場合の開ループ伝達関数

%閉ループ伝達関数を求める
sysC1 = feedback( sysL1, 1 ); %P_1(s)の場合の閉ループ伝達関数
sysC2 = feedback( sysL2, 1 ); %P_2(s)の場合の閉ループ伝達関数

%角周波数の範囲を指定
w = logspace(-1, 2, 1000); %対数的に等間隔なベクトルの生成(10^{-1}から10^{2}で1000点)

%図14.3のプロット

%図14.3(a),(b)のプロット

figure(1) %図のウィンドウを開く

[ gainL1 phaseL1 wL1] = bode( sysL1, w ); %ゲインと位相
gainL1 = gainL1(:); %P_1(s)の場合のゲインの配列を1次元ベクトルに変換
gainL1dB = 20*log10(gainL1); %P_1(s)の場合のゲインをデシベル値に変換
phaseL1 = phaseL1(:); %P_1(s)の場合の位相の配列を1次元ベクトルに変換

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

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

%図14.3(b)のプロット
[ gainL2 phaseL2 wL2 ] = bode( sysL2, w ); %ゲインと位相
gainL2 = gainL2(:); %P_2(s)の場合のゲインの配列を1次元ベクトルに変換
gainL2dB = 20*log10(gainL2); %P_2(s)の場合のゲインをデシベル値に変換
phaseL2 = phaseL2(:); %P_2(s)の場合の位相の配列を1次元ベクトルに変換

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

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

%図14.3(c),(d)のプロット

figure(2) %図のウィンドウを開く

%時間変数の定義
t = 0:0.01:5; %0から5まで0.01刻み

%ステップ応答の計算
[ y1, t1 ] = step( sysC1, t ); %ステップ応答
[ y2, t2 ] = step( sysC2, t ); %ステップ応答

%図14.2(c)のプロット
subplot(1,2,1) %複数の図を並べるためのコマンド.1行2列の1列目という意味
plot(t1,y1); %ステップ応答をプロット
xlim([0 5]) %横軸(時間軸)の範囲の指定
ylim([0.0 1.2]) %縦軸の範囲の指定
xticks([0 1 2 3 4 5]) %横軸の目盛りの値の設定
yticks([0.0 0.5 1.0]) %横軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('y(t)') %縦軸のラベル表示
title('P_{1}(s), C(s)') %タイトルの表示

%図14.2(d)のプロット
subplot(1,2,2) %複数の図を並べるためのコマンド.1行2列の2列目という意味
plot(t2,y2); %ステップ応答をプロット
xlim([0 5]) %横軸(時間軸)の範囲の指定
ylim([0.0 1.2]) %縦軸の範囲の指定
xticks([0 1 2 3 4 5]) %横軸の目盛りの値の設定
yticks([0.0 0.5 1.0]) %横軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸(時間軸)のラベル表示
ylabel('y(t)') %縦軸のラベル表示
title('P_{2}(s), C(s)') %タイトルの表示