function  [yout,x] = flsim(a, b, c, d, u, t, x0, rules)
%
%  FLSIM
%         ... = FLSIM(..., u,t,x0,RULES)
%
%         The function FLSIM is a modification of the Matlab front-end
%       function LSIM.  FLSIM works in an analgous manner, except that
%       closed-loop control is facilitated through the specification of
%       the RULES string.  This string contains a function name for use
%       by FEVAL which specifies u as a function of the original input,
%       u and of the current state, x.  
%         This modification allows FLSIM to be used with a rule-base file
%       in order to simulate fuzzy control where the input u is adaptively
%       changed in real time.  
%          For help in construction of the RULE-BASE, see the demo function
%       URUL1.M   See also, FLTITR, where RULES is used
%         The main difference between LSIM and FLSIM is the addition of the
%       rules argument and the use of the external function FLTITR.M rather
%       than the internal LTITR.  This slows down the simulation, so LSIM
%       should be used when simulating deterministic systems (such as the 
%       results of CLOOP).
%         Since it is expected that the nature of a fuzzy feedback signal
%       will be somewhat discrete, the zero-order-hold function is always
%       used.  If you thind a first-order might be more appropriate, you
%       should explore the rest of the program.
%
%   NOTE:  Be careful with model formulation so that needed
%          variables are accessible and meaningful for the
%          rules function!
%                                                 Roger McNichols 2/28/95


%	LSIM normally linearly interpolates the input (using a first order hold)
%	which is more accurate for continuous inputs. For discrete inputs such 
%	as square waves LSIM tries to detect these and uses a more accurate 
%	zero-order hold method. LSIM can be confused and for accurate results
%	a small time interval should be used.

%  NOTE:  FLSIM always uses a first-order hold!

%	J.N. Little 4-21-85
%       Revised A.C.W.Grace 8-27-89 (added first order hold)
%	                    1-21-91 (test to see whether to use foh or zoh)
%	Copyright (c) 1986-93 by the MathWorks, Inc.

error(nargchk(5,8,nargin));

if (nargin==5),			% transfer function description 
	[num,den] = tfchk(a,b);
	u = c;
	t = d;
	% Convert to state space
	[a,b,c,d] = tf2ss(num,den);

elseif (nargin==6),
	 error('Wrong number of input arguments.');
else
	error(abcdchk(a,b,c,d));
end

if (nargin==7),
    rules = x0;
end

[ns,n] = size(a);
if (nargin==7)|(nargin==5),
	x0 = zeros(ns,1);
end

[p,m]=size(d);
if p*m==0, x=[]; if nargout~=0, yout=[]; end; return; end

if m==1, u=u(:); end, [nu,mu]=size(u);
t=t(:); nt = length(t);

% Make sure u has the right number of columns and rows.
if mu~=m, error('U must have the same number of columns as inputs.'); end
if nu~=nt, error('U must have the same number of rows as the length of T.'); end
	
TS=t(2)-t(1);
% First Order Hold Approximation 

% Try to find out whether the input is step-like or continuous

% Differentiate u
udiff = u(2:nu,:)-u(1:nu-1,:);

% Find all transition points of the input i.e., all changing inputs
uf = find(udiff ~= 0);

% Find out if changes in the input occur after a number of constant inputs.
% If they do then assume that it is a step-like input and use a zero 
% order hold method. Otherwise assume that the input is continuous and 
% linearly interpolate using first order hold method.
usteps = find(diff(uf) == 1);
%  foh = (length(usteps) > 0);

foh = 0;  %  **** This mandates that a zero-order hold be used!!!


if foh 
	% Use first order hold approximation 

% For first order hold approximation first add m integrators in series
	[a,b,c,d]=series(zeros(m),eye(m),eye(m),zeros(m),a,b,c,d);

% For first order hold add (z-1)/TS in series
% This is equivalent to differentiating u.
% Transfer first sample to initial conditions.

	x0=[zeros(m,1);x0(:)]+(b*u(1,:).');
	u(1:nu-1,:) = udiff/TS;
	u(nu,:)=u(nu-1,:);
end

% Get equivalent zero order hold discrete system
[A,B] = c2d(a,b,t(2)-t(1));


x = fltitr(A,B,u,x0,rules);
%  y = x * c.' + u * d.';

y = x*c.';   %  NO FEED FORWARD because of we assume feedback.


if foh
	% Remove the integrator state
	x=x(:,1+m:n+m);
end

if nargout==0		% If no output arguments, plot graph
	plot(t,y), xlabel('Time (secs)'), ylabel('Amplitude')
	return % Suppress output
end
yout = y;