end
		 end
		end
	end

% Next we will save the number of mfs that are on and the pointer
% e_int as to which rules were on.  This vector of length
% 10 saves the last 10 values of e_count and e_int as time
% progresses (hence, it is a moving window).
% These will be used by the FMRLC knowledge-base modifier.

meme_count=[e_count meme_count(1:9)];
meme_int=[e_int meme_int(1:9)];

% The following if-then structure fills the vector mfc with the
% certainty of each membership fucntion of c for the current
% value of c (to understand this part of the code see the above
% similar code for computing mfe)

	if c<=cc(1)	% Takes care of saturation of left-most mf
         mfc=[1 0 0 0 0 0 0 0 0 0 0];
	 c_count=c_count+1;
	 c_int=1;
	elseif c>=cc(numc)
				% Takes care of saturation of the right-most mf
	 mfc=[0 0 0 0 0 0 0 0 0 0 1];
	 c_count=c_count+1;
	 c_int=numc;
	else
		for i=1:numc
		 if c<=cc(i)
		  mfc(i)=max([0,1+(c-cc(i))/wc]);
			if mfc(i)~=0
			 c_count=c_count+1;
			 c_int=i;  	% This term holds last entry
			 			% with a non-zero term
			end
		 else
		  mfc(i)=max([0,1+(cc(i)-c)/wc]);
			if mfc(i)~=0
			 c_count=c_count+1;
			 c_int=i;	% This term holds last entry
			 			% with a non-zero term
			end
		 end
		end
	end
% Next we will save the number of mfs that are on and the pointer
% c_int as to which rules were on.  This vector of length 10
% saves the last 10 values of c_count and e_int as time progresses
% (hence, it is a moving window). These will be used by the FMRLC
% knowledge-base modifier.

memc_count=[c_count memc_count(1:9)];
memc_int=[c_int memc_int(1:9)];

% These for loops calculate the crisp output using only the non-
% zero premise of error,e, and change in error,c.  This cuts down
% computation time since we will only compute
% the contribution from the rules that are on (i.e., a maximum of
% four rules for our case).  Minimum is used for the premise
% and implication.

num=0;
den=0;
	for k=(e_int-e_count+1):e_int
						% Scan over e indices ofmfs that are on
		for l=(c_int-c_count+1):c_int
						% Scan over c indices ofmfs that are on
		  prem=min([mfe(k) mfc(l)]);
		  				% Value of premisemembership function
 % This next calculation of num adds up the numerator for the
 % defuzzification formula. fuzzyrules(k,l) is the output center
 % for the rule.base*(prem-(prem)^2/2 is the area of a symmetric
 % triangle with base width "base" and that is chopped off at
 % a height of prem (since we use minimum to represent the
 % implication). Computation of den is similar but without
 % fuzzyrules(k,l).
		  num=num+fuzzyrules(k,l)*base*(prem-(prem)^2/2);
		  den=den+base*(prem-(prem)^2/2);
		end
	end

   u(index)=num/den;
   			% Crisp output of fuzzy controller that is the input
   			% to the plant

% Next, we perform computations for the fuzzy inverse model.

ye=ym(index)-y(index);		% Calculates ye
yc=(ye-yeold)/step;			% Calculates yc
yeold=ye;					% Saves the value of ye for use the
							% next time

ye_count=0;,yc_count=0;    	% Counts the non-zero certainties
							% of ye and ycsimilar to how we did
							% for the fuzzycontroller

% The following if-then structure fills the vector mfye with the
% certainty of each membership fucntion of ye (similar to the
% fuzzy controller).  Notice that we use the same number of
% input membership functions as in the fuzzy controller -
% just for convenience - it does not have to be this way

	if ye<=cye(1)			% Takes care of saturation
         mfye=[1 0 0 0 0 0 0 0 0 0 0];
	 ye_count=ye_count+1;,ye_int=1;
	elseif ye>=cye(numye)	% Takes care of saturation
	 mfye=[0 0 0 0 0 0 0 0 0 0 1];
	 ye_count=ye_count+1;,ye_int=numye;
	else
		for i=1:numye
		 if ye<=cye(i)
		  mfye(i)=max([0 1+(ye-cye(i))/wye]);
			if mfye(i)~=0
			 ye_count=ye_count+1;
			 ye_int=i;	% This term holds last entry with
			 			% a non-zero term
			end
		 else
		  mfye(i)=max([0,1+(cye(i)-ye)/wye]);
			if mfye(i)~=0
			 ye_count=ye_count+1;
			 ye_int=i;	% This term holds last entry with
			 			% a non-zero term
			end
		 end
		end
	end

% The following if-then structure fills the vector mfyc with the
% certainty of each membership fucntion of yc

	if yc<=cyc(1)
         mfyc=[1 0 0 0 0 0 0 0 0 0 0];
	 yc_count=yc_count+1;,yc_int=1;
	elseif yc>=cyc(numyc)
	 mfyc=[0 0 0 0 0 0 0 0 0 0 1];
	 yc_count=yc_count+1;,yc_int=numyc;
	else
		for i=1:numyc
		 if yc<=cyc(i)
		  mfyc(i)=max([0 1+(yc-cyc(i))/wyc]);
			if mfyc(i)~=0
			 yc_count=yc_count+1;
			 yc_int=i;
			end
		 else
		  mfyc(i)=max([0,1+(cyc(i)-yc)/wyc]);
			if mfyc(i)~=0
			 yc_count=yc_count+1;
			 yc_int=i;
			end
		 end
		end
	end

% These for loops calculate the crisp output of the inverse model
% using only the non-zero premise of error,ye, and change in
% error,yc.  This cuts down computation time (similar to the
% fuzzy controller). Minimum is used for the premise and the
% implication.  To understand this part of the code see the
% similar code for the fuzzy controller.

invnum=0;
invden=0;
	for k=(ye_int-ye_count+1):ye_int
		for l=(yc_int-yc_count+1):yc_int
		  prem=min([mfye(k) mfyc(l)]);
		  invnum=invnum+inverrules(k,l)*invbase*(prem-(prem)^2/2);
		  invden=invden+invbase*(prem-(prem)^2/2);
		end
	end

% Next we compute the output of the fuzzy inverse model.
% If you want to let gp=0 to test the fuzzy controller by itself
% then you will have invden=0 and you will not be able to compute p.
% To make it possible to let gp=0 we put in the following
% if-then rule.

if gp==0
p=0;
else
p=invnum/invden;   			% Crisp output of inverse model
if abs(p)<.05, p=0; end     % robustification term

end

% Next we modify the centers of ouput membership functions
% associated with the rules that were on d steps ago
% The indexing sheme is similar to the one used in the
% computation of the outputs of the fuzzy controller and
% fuzzy inverse model.

	for k=(meme_int(d)-meme_count(d)+1):meme_int(d)
		for l=(memc_int(d)-memc_count(d)+1):memc_int(d)

		  fuzzyrules(k,l)=fuzzyrules(k,l)+p;
		end
	end

% Next, the Runge-Kutta equations to find the next state.

	time(index)=t;
	F=A*x+B*u(index);
	k1=step*F;
	xnew=x+k1/2;
	F=A*xnew+B*u(index);
	k2=step*F;
	xnew=x+k2/2;
	F=A*xnew+B*u(index);
	k3=step*F;
	xnew=x+k3;
	F=A*xnew+B*u(index);
	k4=step*F;
	x=x+(1/6)*(k1+2*k2+2*k3+k4); % Calculated next state

t=t+step;  			% Increments time
index=index+1;	 	% Increments the indexing term so that
					% index=1 corresponds to time t=0.

end % This end statement goes with the first "while" statement
    % in the program

%
% Next, we provide plots of the input and output of the plant
% and the output of the reference model
%

subplot(211)
plot(time,y,'-',time,ym,'--')
xlabel('Time (sec)')
title('Output of plant (solid) and reference model (dashed)')
subplot(212)
plot(time,u)
xlabel('Time (sec)')
title('Output of fuzzy controller (input to the plant)')