%図12.12
%
clear all; close all;

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

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

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

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

%図12.12(a)

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

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

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

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

%図12.12(b)

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

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

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

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

%図12.12(c)

%正弦波(入力)の計算
u3 = A*sin( omega3*t ); %ω=1.0の正弦波

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

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

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

%図12.12(d)

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

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

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

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