%図6.4
%
clear all; close all;

%パラメータの値を与える
K = 1; %K=1
zeta = 2; %ζ=2
omegan1 = 1.1; %ω_n=1.1
omegan2 = 5.0; %ω_n=5
omegan3 = 10.0; %ω_n=10

%伝達関数の分子・分母多項式を与える
num1 = [ 0 0 K*omegan1^2 ]; %ω_n=1.1の場合の分子多項式
den1 = [ 1 2*zeta*omegan1 omegan1^2 ]; %ω_n=1.1の場合の分母多項式
num2 = [ 0 0 K*omegan2^2 ]; %ω_n=5.0の場合の分子多項式
den2 = [ 1 2*zeta*omegan2 omegan2^2 ]; %ω_n=5.0の場合の分母多項式
num3 = [ 0 0 K*omegan3^2 ]; %ω_n=10.0の場合の分子多項式
den3 = [ 1 2*zeta*omegan3 omegan3^2 ]; %ω_n=10.0の場合の分母多項式

%伝達関数表現を与える
sys1 = tf( num1, den1 ); %ω_n=1.1の場合の伝達関数表現
sys2 = tf( num2, den2 ); %ω_n=5.0の場合の伝達関数表現
sys3 = tf( num3, den3 ); %ω_n=10.0の場合の伝達関数表現

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

%インパルス応答の計算
[ y1, t1 ] = impulse( sys1, t ); %ω_n=1.1の場合のインパルス応答
[ y2, t2 ] = impulse( sys2, t ); %ω_n=5.0の場合のインパルス応答
[ y3, t3 ] = impulse( sys3, t ); %ω_n=10.0の場合のインパルス応答

%図6.4のプロット
plot(t1,y1,'-b',t2,y2,'-g',t3,y3,'-r'); %インパルス応答をプロット
xlim([0,5]) %横軸(時間軸)の範囲の指定
ylim([0,2.5]) %縦軸の範囲の指定
xticks([0 1 2 3 4 5]) %横軸の目盛りの値の設定
yticks([0.0 0.5 1.0 1.5 2.0 2.5]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('y(t)') %縦軸のラベル表示
legend('ω_n = 1.0','ω_n = 5.0','ω_n = 10'); %凡例の表示