function y=dmgrade(x,mfs)
%MAN_PAGE_BEGIN
%
%   @purpose	Compute the derivatives of a fuzzy set.
%
%   @synopsis	y=dmgrade(x,mfs)
%   @description  
%
%   @see mgrade
%
%MAN_PAGE_END

% By Antoine Duchateau. All Rights Reserved.
% Written: 18/05/99

% Copyright (c) Robert Babuska, 1992-98

if ~iscell(mfs),		% standard situation
	[k,l]=size(mfs);
	x = x(:)';
	lx = length(x);
	
	mx = zeros(1,lx);	% <<<<<<<
	
	if ~all(mfs(:,2)<=mfs(:,3) & mfs(:,3)<=mfs(:,4) & mfs(:,4)<=mfs(:,5)),
		error('Illegal membership function parameters');
	end;
	
	for j=1:k
		y(j,:) = mx;
		
		if mfs(j,1) == 0,                    % ANY
		elseif mfs(j,1) == 1,                % trapezoidal
			if mfs(j,3)>mfs(j,2),
				y(j,:) = ((x > mfs(j,2)) & (x < mfs(j,3)))/(mfs(j,3) - mfs(j,2));
			end
         	if mfs(j,5) > mfs(j,4),
               y(j,:) = - ((x > mfs(j,4)) & (x < mfs(j,5)))/(mfs(j,5) - mfs(j,4));
         	end;
      	elseif mfs(j,1) == 2,                % exponential
         	wl = mfs(j,3) - mfs(j,2);
         	wr = mfs(j,5) - mfs(j,4);
         	if wl > 0
            	y(j,:) = (x<mfs(j,3)) .* exp(-7*(x-mfs(j,3)).^2/wl^2) .* (x-mfs(j,3)) * (-14/wl^2);
         	end;
         	if wr > 0
            	y(j,:) = (x>mfs(j,4)) .* exp(-7*(x-mfs(j,4)).^2/wr^2) .* (x-mfs(j,4)) * (-14/wr^2);
         	end;
      	elseif mfs(j,1)==3,                  % sigmoidal
			wl = (mfs(j,3) - mfs(j,2))/2;
			wr = (mfs(j,5) - mfs(j,4))/2;
			if wl > 0
				y(j,:) = ((x>mfs(j,2)) & (x<=mfs(j,2)+wl)) .* ((x - mfs(j,2))/wl^2) ...
					   + ((x>mfs(j,2)+wl) & (x<mfs(j,3))) .* ((mfs(j,3) - x)/wl^2);
			end;
			if wr > 0
				y(j,:) = ((x>mfs(j,4)) & (x<=mfs(j,4)+wl)) .* ((mfs(j,4) - x)/wr^2) ...
					   + ((x>mfs(j,4)+wl) & (x<mfs(j,5))) .* ((x - mfs(j,5))/wr^2);
			end;
		end;
	 end;
	y = y';   
else			% cell array supplied
   for k = 1 : size(x,1),
      a = mgrade(x(k,1),mfs{1});
      sb = length(a);
      for i = 2 : length(mfs),
         b = dmgrade(x(k,i),mfs{i});
         sb(i) = length(b);
         a = a(:);
         %	   a = min(a*ones(1,sb(i)),ones(size(a))*b(:)');
         a = a*b(:)';
      end;   
      %   y = zeros(sb);       % ND array
      %   y(:) = a(:);
      y(k,:) = a(:)';
	end;   
end;   
