echo off;
clear all;
close all;
clc;

echo on;

% --------------------------------------------------------
%
%            Control Engineering
%      Nagaoka Univercity of Technology
%
%  Prof. Kenji HIRATA and Teaching Assistants.
%          
%         -- Computer Exercise 3 --
%
%  ver.1.0  2006
%  ver.1.1  2007
%  ver.1.2  2008 by Keigo NISHIGAMI.
%  ver.1.3  2011 by Masahiro YOSHITAKE.
% --------------------------------------------------------

pause	% press any key to continue.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computer Exercise 3 では，Matlab を用いてゲイン補償および
% 位相遅れ補償による制御系設計を行い，制御系設計に関する理解を
% 深めることを目的とする．
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pause	% press any key to continue.

% --------------------------------------------------------
% [ 課題 ]
% 制御対象が次の 2 次系の伝達関数で与えられる場合を考える．
%
%                               20
%                 p(s) = -----------------
%                         (s + 10)(s + 1)
%
% p(s) によって構成されるフィードバック系に対して，定常（位置）
% 偏差が十分小さく，位相余裕 PM が 40 〜 60 [deg] となるように
% 位相遅れ (PI) 補償器を設計したい．
% 
% 以下に示される各ステップの指示に従って，
% 位相遅れ補償器を設計しなさい．
%
% (わからない場合は講義スライド 13 か，
%  テキスト P.159 を参考にしてください)
% --------------------------------------------------------

pause	% press any key to continue.

% ----------------------- [ 準備 1 ] -----------------------
% はじめに，制御対象の伝達関数 p(s) の定義と p(s) のみで構成
% される閉ループ系の伝達関数 t(s) を計算します．
% --------------------------------------------------------

pause	% press any key to continue.

%%% 制御対象 p(s) は，次で与えられます．
p = tf([20],[1 11 10]);
zpk(p)

pause	% press any key to continue.

%%%%% 閉ループ伝達関数 t(s) は，次で与えられます．
t = feedback(p,1);
zpk(t)

pause	% press any key to continue.

% ----------------------- [ 準備 2 ] -----------------------
% 次に，補償なしの場合の開ループ系のボード線図と閉ループ系の
% ステップ応答を計算し，グラフに表示します．表示されたグラフ
% から，このシステムの特性を把握してください．
% --------------------------------------------------------

pause	% press any key to continue.

%%%%% plant bode 線図のデータ
[magp,phasep,wp]=bode(p,{1e-2,1e2});
[gmp,pmp,wpcp,wgcp] = margin(p);
logmagp = 20*log10(arr2vec(magp,wp));
php = arr2vec(phasep,wp);

%%%%% 閉ループ t のステップ応答のデータ
ts = 0.01;
finalt = 1.6;
time = 0:ts:finalt;
uu = 1*ones(size(time));
[yt,tt] = lsim(t,uu,time);

pause	% press any key to continue.

%%%%%%%%%%%%%%% plant bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',2);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',2);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcp = {sprintf('\\omega_{gc}=%2.2f',wgcp)};
text(0.77,0.85,swgcp,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcpd = {sprintf('uncompensated system')};
text(0.55,0.9,swgcpd,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcpdd = {sprintf('PM=%2.1f',pmp)};
text(0.75,0.4,swgcpdd,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
%%%%% file name
tgain = text(0.01,0.01,'p01dodep.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt(size(time,2)) 1],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット %%%%%
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt(size(time,2)),'-.');
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%% file name
tgain = text(0.01,0.01,'p01yt.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);
swgcpd = {sprintf('uncompensated system')};
text(0.55,0.9,swgcpd,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

pause	% press any key to continue.

close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以下では，ゲイン補償 k を決定するために，各定数ゲイン k に対する開ループ
% 系のボード線図と閉ループ系のステップ応答を表示します．
% ピンク色の線が補償後の応答を表します．
% 各定数ゲイン k に対するゲイン交差周波数 wgc (画面に表示)，
% 位相余裕 PM, 立ち上がり時間 Tr を読み取ってください．
%
% 立ち上がり時間 Tr は，速応性の指標となります．
% ここでは，立ち上がり時間 Tr を，応答が定常値の 10% から 90% に
% 達するまでの時間と定義します．
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pause	% press any key to continue.

% --------------------------------------------------------------------
% ゲイン補償 k = 0.1 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 0.1
sprintf('%s','proportional compensation')
k=0.1;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 yt1(size(time,2)) 0.5 1],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

pause	% press any key to continue.

% --------------------------------------------------------------------
% 定数ゲイン k = 0.5 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 0.5
sprintf('%s','proportional compensation')
k=0.5;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 yt1(size(time,2)) 1],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

pause	% press any key to continue.

% --------------------------------------------------------------------
% 定数ゲイン k = 1.0 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 1.0
sprintf('%s','proportional compensation')
k=1.0;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt1(size(time,2)) 1],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);
pause	% press any key to continue.
     
% --------------------------------------------------------------------
% 定数ゲイン k = 7 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 7
sprintf('%s','proportional compensation')
k=7;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt1(size(time,2)) 1],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);
     
pause	% press any key to continue.
     
% --------------------------------------------------------------------
% 定数ゲイン k = 10 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 10
sprintf('%s','proportional compensation')
k=10;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt1(size(time,2)) 1],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

pause	% press any key to continue.
     
% --------------------------------------------------------------------
% 定数ゲイン k = 20 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 20
sprintf('%s','proportional compensation')
k=20;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt1(size(time,2))],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

pause	% press any key to continue.
     
% --------------------------------------------------------------------
% 定数ゲイン k = 50 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 50
sprintf('%s','proportional compensation')
k = 50;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);


%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt1(size(time,2))],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

pause % press any key to continue.

% --------------------------------------------------------------------
% 定数ゲイン k = 100 に対する開ループ系のボード線図と閉ループ系のステップ
% 応答を表示します．
% --------------------------------------------------------------------

pause	% press any key to continue.

%%%%% proportional control prameter
%%%%% k = 100
sprintf('%s','proportional compensation')
k=100;
c1 = k
l1 = p*c1;

%%%%% l1 bode 線図のデータ
[magl1,phasel1,wl1]=bode(l1,{1e-2,1e2});
[gml1,pml1,wpcl1,wgcl1] = margin(l1);
logmagl1 = 20*log10(arr2vec(magl1,wl1));
phl1 = arr2vec(phasel1,wl1);

%%%%% 閉ループ系 t
t1 = feedback(l1,1);

%%%%% 閉ループ t のステップ応答のデータ
[yt1,tt1] = lsim(t1,uu,time);

%%%%%%%%%%%%%%% bode 線図のプロット
fig01 = figure(1); clf;
set(fig01,'position',[20 140 560 750]);

%%%%% gain 線図の axes
hgain = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','log',...
	     'xlim',[1e-1 1e2],...
	     'xtick',[1e-1 1e0 1e1 1e2],...
	     'xgrid','off',...
	     'ylim',[-30 50],...
	     'ytick',[-20 0 20 40],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(fig01,'currentaxes',hgain); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[100 -100],':k',...
	      [1e1 1e1],[100 -100],':k');
set(get(hgain,'xlabel'),...
    'string','\omega[rad/s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(hgain,'ylabel'),...
    'string','gain [dB]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hgain); %%%%% 0 [dB] の線 
hold on; l0gain = plot([1e-1 1e2],[0 0],'-.b');

%%%%% 位相線図の axes
hphase = axes('position',[0.2 0.15 0.75 0.34],...
	      'box','on',...
	      'fontname','Helvetica',...
	      'fontsize',14,...
	      'xscale','log',...
	      'xlim',[1e-1 1e2],...
	      'xtick',[1e-1 1e0 1e1 1e2],...
	      'xgrid','off',...
	      'ylim',[-200 0],...
	      'ytick',[-180 -135 -90 -45 0],...
	      'ygrid','on',...
	      'gridlinestyle',':');
set(fig01,'currentaxes',hphase); %%%%% x 軸のグリッドを書いておく
hold on; plot([1e-0 1e-0],[-200 0],':k',...
	      [1e1 1e1],[-200 0],':k');
set(get(hphase,'xlabel'),'string','\omega[rad/s]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% x 軸ラベル
set(get(hphase,'ylabel'),'string','phase [deg]',...
		  'fontname','helvetica',...
		  'fontsize',14); %%%%% y 軸ラベル
set(fig01,'currentaxes',hphase); %%%%% -180 [deg] の線 
hold on; lm180phase = plot([1e-1 1e2],[-180 -180],'-.b');

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% gain 線図のプロット %%%%% P
set(fig01,'currentaxes',hgain);
hold on; plogmagp = plot(wp,logmagp,'-c'); set(plogmagp,'linewidth',1);
hold on; pwgp = plot([wgcp wgcp],[0 -30],'-.b');

%%%%% gain 線図のプロット %%%%% L1
set(fig01,'currentaxes',hgain);
hold on; plogmagl1 = plot(wl1,logmagl1,'-m'); set(plogmagl1,'linewidth',2);
hold on; pwgp = plot([wgcl1 wgcl1],[0 -30],'-.b');

%%%%% 位相線図のプロット %%%%% P
set(fig01,'currentaxes',hphase);
hold on; pphp = plot(wp,php,'-c'); set(pphp,'linewidth',1);
hold on; pwgcp = plot([wgcp wgcp],[-180 0],'-.b');
hold on; ppmp = plot([wgcp wgcp],[-180+pmp -180],'-g');
set(ppmp,'linewidth',1);

%%%%% 位相線図のプロット %%%%% L1
set(fig01,'currentaxes',hphase);
hold on; pphl1 = plot(wl1,phl1,'-m'); set(pphl1,'linewidth',2);
hold on; pwgcl1 = plot([wgcl1 wgcl1],[-180 0],'-.b');
hold on; ppml1 = plot([wgcl1 wgcl1],[-180+pml1 -180],'-g');
set(ppml1,'linewidth',2);

%%%%% コメント
set(fig01,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) gain',...
	     'fontname','helvetica',...
	     'fontsize',14);
tgain = text(0.5,0.05,'(b) phase',...
	     'fontname','helvetica',...
	     'fontsize',14);
swgcl1 = {sprintf('\\omega_{gc}=%2.2f',wgcl1)};
text(0.77,0.85,swgcl1,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01dodel1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

%%%%%%%%%% ステップ応答のプロット
fig02 = figure(2); clf;
set(fig02,'position',[600 140 560 750]);

%%%%% y(t) 線図の axes
houtput = axes('position',[0.2 0.65 0.75 0.34],...
	     'box','on',...
	     'fontName','Helvetica',...
	     'fontsize',14,...
	     'xscale','linear',...
	     'xlim',[0 1.5],...
	     'xtick',[0 0.5 1 1.5],...
	     'xgrid','on',...
	     'ylim',[0 1.4],...
	     'ytick',[0 0.5 yt1(size(time,2))],...
	     'ygrid','on',...
	     'gridlinestyle',':');
set(get(houtput,'xlabel'),...
    'string','\itt \rm[s]',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% x 軸ラベル 
set(get(houtput,'ylabel'),...
    'string','\ity\rm(\itt\rm)',...
    'fontname','helvetica',...
    'fontsize',14); %%%%% y 軸ラベル

%%%%% コメント用の axes
hcomment = axes('position',[0 0 1 1],...
		'visible','off',...
		'xlim',[0 1],...
		'ylim',[0 1]);

%%%%% 出力 y のプロット P
set(fig02,'currentaxes',houtput);
hold on; pyt = plot(tt,yt,'-c'); set(pyt,'linewidth',1);

%%%%% 出力 y のプロット L1
set(fig02,'currentaxes',houtput);
hold on; ptr = plot(tt,ones(1,size(time,2))*yt1(size(time,2)),'-.');
hold on; pyt1 = plot(tt1,yt1,'-m'); set(pyt1,'linewidth',2);

%%%%% コメント
set(fig02,'currentaxes',hcomment);
tgain = text(0.5,0.55,'(a) output \ity\rm(\itt\rm)',...
	     'fontname','helvetica',...
	     'fontsize',14);
%%%%%
swgcl1d = {sprintf('k=%2.1f',k)};
text(0.8,0.9,swgcl1d,...
     'FontName','Times',...
     'FontSize',14,...
     'Color','b');

%%%%% file name
tgain = text(0.01,0.01,'p01yt1.eps',...
	     'fontname','helvetica',...
	     'fontsize',8);

pause % press any key to continue.
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 定数ゲイン k が大きくなるにしたがって，位相余裕が小さくなっていることを
% 確認してください．
% また，k が大きくなると，定常特性は改善されますが，定常偏差を 0 にするには
% いくらでも大きなゲイン k が必要なことに注意してください．
% さらに，ゲイン k が大きくなるにしたがって，ステップ応答が振動的になって
% いることを確認してください．
%
% 確認ができたら，読み取った位相余裕 PM, 立ち上がり時間 Tr の値から，
% 適切なゲイン補償 k を決定してください．
% (問題で設計指針として示した位相余裕は 40 〜 60 [deg] です)
%
% この手順は，講義スライド 13 か，テキスト P.159 の
% 位相遅れ補償の設計手順 [ステップ1] です．
%
% 以上で，ゲイン補償は終了です．
% 次いで，位相遅れ補償の設計に移ります．
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pause % press any key to continue.
