function [ConInd, ConRob, ConDes] = conicity(m,range,controller,params)
%CONICITY compute the conicity stabity criterium
%    
%    This method is	applicable only if the plant system is
%    has single takagi sugeno  mapping for all the outputs or
%    it is a linear one.  In addition the controller must be a
%    static one. The outputs of  the function are: CONICITY
%    INDEX = CONIC DESVIATION / CONIC ROBUSTNESS  CONIC
%    DEVIATION = GAIN(CONTROLLER - CENTRE) CONIC ROBUSTNESS =
%    1 / GAIN(FEEDB((A,B,C,D),CENTRE) where CENTRE is the
%    center of the cone which  is calculated with the method
%    specified in the field	SEARCH_METHOD of the PARAMS
%    structure. Possible values for this field are:
%    
%      LINEARIZATION: which uses a linarization technique
%          for finding  the center of the cone
%    
%    if no SEARCH_METHOD is specified the linearization
%    techique is used.
%    
%    CONICITY(PLANT,RANGE,PARAMS,CONTROLLER) coputes the
%    conicity stability criterium for the plant PLANT in the
%    range specified  by the array RANGE with the static
%    controller  CONTROLLER.
%    
%    CONICITY(PLANT,RANGE,CONTROLLER)
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose  compute the conicity stabity criterium. This method is 
%   applicable only if the plant system is has single takagi sugeno 
%   mapping for all the outputs or it is a linear one. 
%   In addition the controller must be a static one. The outputs of 
%   the function are:
%
%   CONICITY INDEX = CONIC DESVIATION / CONIC ROBUSTNESS 
%
%	CONIC DEVIATION = GAIN(CONTROLLER - CENTRE)
%
%   CONIC ROBUSTNESS = 1 / GAIN(FEEDB((A,B,C,D),CENTRE)
% 
%   where <code>centre</code> is the center of the cone which 
%   is calculated with the method specified in the field 
%   <code>search_method</code> of the <code>params</code> structure.
%   Possible values for this field are:
%<dl>
%<dt> linearization <dd> which uses a linarization technique for finding 
%                        the center of the cone
% </dl>  
%   if no <code>search_method</code> is specified the linearization
%   techique is used.
%
%   @synopsis   conicity(plant,range,params,controller)
%   @description coputes the conicity stability criterium
%   for the plant <code>plant</code> in the range specified 
%   by the array <code>range</code> with the static controller 
%   <code>controller</code>.
%
%   @synopsis  conicity(plant,range,controller)
%   @description
%
%MAN_PAGE_END

error(nargchk(2,4,nargin));

%check if the system has a single  taksug mapping or it is a linear one

links=get(m,'links');

if any(links(:,1)~=links(1,1)),
	error('The conicity criterium is applicable only to systems with sinlge mapping');
end

statemappings=get(m,'statemapping');
plant = statemappings{links(1,1)};

if not(isa(plant,'taksug') | isa(plant,'linear')),
	error('The conicity criterium is applicable only to taksug or linear systems');
end   

%if exist('controller') & ~(is_static(controller)),
%	error('The conicity criterium is applicable only with static controllers');
%end

%separate the linear part of the system

plant_states = get(m,'n_state');
lin=get(plant,'linears');

if isa(plant,'taksug'),
	A=sum(lin(:,1:plant_states,:),3);
else 
    A=lin(:,1:plant_states);
end 
B=eye(plant_states);
C=eye(plant_states);
D=zeros(plant_states);

%prepare the grid with all the combination of the states

totalPoints=10;		% Number of state vector generated
num_variables=size(range,2);
interval= ((range(2,:)-range(1,:)).*ones(1,num_variables))/...
				((totalPoints/num_variables)^(1/num_variables));
states=[];

for i=1:size(range,2) 
	if isempty(states),
		states=[range(1,i):interval(i):range(2,i)]';
	else
  	    xx=states;
	    yy=[range(1,i):interval(i):range(2,i)]';
	    nx=size(xx,1);
	    ny=size(yy,1);
	    xx = repmat(xx,ny,1);
	    yy=reshape(repmat(yy,1,nx)',nx*ny,1);
	    states=[xx yy];
  	end
end

if exist('controller'),
	feedback_values=eval(controller,states);
end

if isa(plant,'taksug'),
	w = membership(plant,states)';
	
	outputs = zeros(size(states,1),plant_states);
	for j=1:plant_states,
	
		%calculate the outputs given the grid of the states
	    
	    lin_j=squeeze(lin(j,:,:))';
	    
	    lin_a=lin_j(:,1:plant_states);
	    lin_b=lin_j(:,plant_states+1: end-1);
	    lin_c=lin_j(:,end);
	    
	    sf=0;
	    sc=0;
	    sb=0;
	    
	    for i=1:get(plant,'n_rules'),
	    	sf= sf + sum(repmat(lin_a(i,:),size(states,1),1).* ...
	    	             states,2).* (1-w(:,i));
	    	
	    	if exist('controller')
	    		sb= sb + sum(repmat(lin_b(i,:),size(feedback_values,1),1).* ...
	    		             feedback_values,2).* w(:,i);
	    	end
	    	
	    	sc = sc + repmat(lin_c(i,:),size(states,1),1).* w(:,i);
	    end
	    
	    outputs(:,j)=-sf+sb+sc;
	
	end
else
	outputs = zeros(size(states,1),plant_states);
	for j=1:plant_states,
	    sc=0;
	    sb=0;
	    
	    if exist('controller')
	    	sb= sum(repmat(lin(j,plant_states+1:end-1),size(feedback_values,1),1).* ...
	    		             feedback_values,2);
	    end
	    
	    sc = repmat(lin(j,end),size(states,1),1);
	
		outputs(:,j)=sb+sc;
	end
end

keyboard;

states=states';
outputs=outputs';


%calculate the center of the cone
if exist('params') and isfield(params, 'search_method'),
	switch params.search_method
		case 'linearization',
			center=conse_cc(states,outputs);
    %add here the other methods
    end
else
	center=conse_cc(states,outputs);
end		

%calcuate the inverse of conic robustness
ltiGain=x_conlti(A,B,C,D,center);	

%and the conic deviation
nlsGain=x_connls(states,outputs,center);		

%calculates the indeces
ConInd=nlsGain*ltiGain		% Conicity Index
ConRob=1/ltiGain
ConDes=nlsGain

