%図12.5
%
clear all; close all;

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

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

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

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

%図12.5(a)

%正弦波(入力)の計算
u1 = A*sin( omega1*t ); %ω=0.1の正弦波を与える

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

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

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

%図12.5(b)

%正弦波(入力)の計算
u2 = A*sin( omega2*t ); %ω=0.5の正弦波を与える

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

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

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

%図12.5(c)

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

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

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

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

%図12.5(d)

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

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

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

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