function [y,uc,t]=simu(yr,h,R,S,T)
% simulation of a process
%
% [y,u,t]=simu(u,h)
% simulate the open loop system with 
%	input signal	u 
%	sampling period	h
%	output signal	y
% the vector t is the time vector
%
% [y,uc,t]=simu(yr,h,R,S,T)
% simulate closed loop system with
%	reference signal	yr
%	controller polynomials	R,S,T
% the signal uc is the process input


% (c) Lennart Andersson & Ulf Jönsson 1995
% last update 941228


if nargin==0,
  help simu
  return
end
[N,m]=size(yr);
if m>1,
  error('signals should be column vectors');
end
if nargin==1,
  error('sampling period unspecified');
end
if nargin<5,
  R=1;
  S=0;
  T=1;
end


% flexible servo parameters

J1=22e-6;
J2=60e-6;
kf=4e-3;
d1=3e-5;
d2=3e-5;
km=0.1;
ki=0.25;
kw=0.1;

% statespace representation, sampling and transformation to transfer function

A=[-d1/J1 0 -kf/J1;0 -d2/J2 kf/J2; 1 -1 0];
B=[km*ki/J1;0;0];
C=[0 kw 0];
[Phi,Gamma]=c2d(A,B,h);
[b,a]=ss2tf(Phi,Gamma,C,0);


% noise parameters

omega=25;		% noise filter bandwidth
zeta=1;
sigmae=5*0.01/sqrt(h); % white noise level;
sigmaoutl=1;		% outlier level


% discrete time noise description (approximate !)

[Aw,Bw,Cw,Dw]=tf2ss(omega^2,[1 2*zeta*omega omega^2]);
[Phiw,Gammaw]=c2d(Aw,Bw,h);

[bw,aw]=ss2tf(Phiw,Gammaw,Cw,Dw);
w=filter(bw,aw,sigmae*randn(N,1));

% load disturbance parameter

ld=0.1;

% cloosed loop transfer functions

at=conv(a,R)+[zeros(1,length(a)+length(R)-length(b)-length(S)) conv(b,S)];
bry=conv(b,T);
bly=conv(b,R);
bwy=conv(a,R);
bru=conv(a,T);
blu=conv(a,R);
bwu=conv(a,-S);

% simulate closed loop

uc=filter(bru,at,yr)+filter(bwu,at,w);
y=filter(bry,at,yr)+filter(bly,at,ld*ones(N,1))+filter(bwy,at,w);
t=h.*(0:N-1);