%図12.11
%
clear all; close all;

%正弦波A sin ωtのパラメータを与える
A = 1; %正弦波の振幅
omega1 = 0.01; %ω=0.01の場合の正弦波の角周波数
omega2 = 1.0; %ω=1.0の場合の正弦波の角周波数
omega3 = 10; %ω=10の場合の正弦波の角周波数
omega4 = 100.0; %ω=100.0の場合の正弦波の角周波数

%G(s)の分子・分母多項式を与える
num = [ 0 10 10 ]; %G(s)の分子多項式
den = [ 1 10.1 1 ]; %G(s)の分母多項式

%伝達関数表現を与える
sys = tf( num, den ); %伝達関数表現

%図12.11(a)

%時間変数の定義
t1 = 0:0.1:1000; %0から1000まで0.1刻み

%正弦波(入力)の計算
u1 = A*sin( omega1*t1 ); %ω=0.1の正弦波

%応答の初期値を与える
x0 = [ 0 ]; %x(0)=0と与える

%周波数応答(出力)を求める
[ y1, t1, x1 ] = lsim( sys, u1, t1, x0 ); %周波数応答

%図12.11(a)のプロット
subplot(2,2,1) %複数の図を並べるためのコマンド.2行2列の1行1列目という意味
plot(t1,u1,'-b',t1,y1,'-r'); %周波数応答をプロット
xlim([0,1000]) %横軸(時間軸)の範囲の指定
ylim([-10 10]) %縦軸の範囲の指定
xticks([0 200 400 600 800 1000]) %横軸の目盛りの値の設定
yticks([-10 -5 0.0 5 10]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('input and output') %縦軸のラベル表示

%図12.11(b)

%時間変数の定義
t2 = 0:0.001:50; %0から50まで0.001刻み

%正弦波(入力)の計算
u2 = A*sin( omega2*t2 ); %ω=1.0の正弦波

%応答の初期値を与える
x0 = [ 0 ]; %x(0)=0と与える

%周波数応答(出力)を求める
[ y2, t2, x1 ] = lsim( sys, u2, t2, x0 ); %周波数応答

%図12.11(b)のプロット
subplot(2,2,2) %複数の図を並べるためのコマンド.2行2列の1行2列目という意味
plot(t2,u2,'-b',t2,y2,'-r'); %周波数応答をプロット
xlim([0,50]) %横軸(時間軸)の範囲の指定
ylim([-1.5 2.5]) %縦軸の範囲の指定
xticks([0 10 20 30 40 50]) %横軸の目盛りの値の設定
yticks([-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('input and output') %縦軸のラベル表示

%図12.11(c)

%正弦波(入力)の計算
u3 = A*sin( omega3*t2 ); %ω=10の正弦波を与える

%応答の初期値を与える
x0 = [ 0 ]; %x(0)=0と与える

%周波数応答(出力)を求める
[ y3, t3, x3 ] = lsim( sys, u3, t2, x0 ); %周波数応答を求める

%図12.11(c)のプロット
subplot(2,2,3) %複数の図を並べるためのコマンド.2行2列の2行1列目という意味
plot(t3,y3); %周波数応答をプロット
xlim([40,50]) %横軸(時間軸)の範囲の指定
ylim([-1.0 1.0]) %縦軸の範囲の指定
xticks([40 42 44 46 48 50]) %横軸の目盛りの値の設定
yticks([-1.0 -0.5 0.0 0.5 1.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('output') %縦軸のラベル表示

%図12.11(d)

%正弦波(入力)の計算
u4 = A*sin( omega4*t2 ); %ω=100.0の正弦波を与える

%応答の初期値を与える
x0 = [ 0 ]; %x(0)=0と与える

%周波数応答(出力)を求める
[ y4, t4, x4 ] = lsim( sys, u4, t2, x0 ); %周波数応答を求める

%図12.11(d)のプロット
subplot(2,2,4) %複数の図を並べるためのコマンド.2行2列の2行目2列目という意味
plot(t4,y4); %周波数応答をプロット
xlim([49,50]) %横軸(時間軸)の範囲の指定
ylim([-0.10 0.10]) %縦軸の範囲の指定
xticks([49.0 49.2 49.4 49.6 49.8 50]) %横軸の目盛りの値の設定
yticks([-0.10 -0.05 0 0.05 0.10]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('output') %縦軸のラベル表示