function fig3

% These are graphic display commands
% The computation engine is in the functions onesidedOpDyn and boundConf

ee=[10 20 50 100];
clf
for i=1:4
	subplot(2,2,i)
	e = ee(i);
	box on
	X = [0;1];X(:,2)=1; % this initial profile will be extended in function onesidedOpDyn
	t = 9;
	Y = onesidedOpDyn(X,e+10^-10,t); % e+10^-10 must be chosen due to numerical vulnerablity 
                                     % of abs(X-Xnew)<=e when abs(X-Xnew)==e
    n = size(Y,1);
	c = colormap(hsv(n));
	set(gca,'Layer','top','nextplot','replacechildren','ColorOrder',c)
	axis([0 t 0 3.2*e])
	plot(0:t,Y')
    title(['\epsilon = ' num2str(e)])
    if i==3 | i==4
        xlabel('t')
    end
    if i==1 | i==3
        ylabel('agents')
    end
end
set(gcf,'Position',[104    82   600   430],'PaperPositionMode','auto')
print -dtiff -r96 fig2



function Y = onesidedOpDyn(X,e,t)

% Y = onesidedOpDyn(X,e)
% X    : Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
% e    : bound of confidence
% t    : time step until we compute the dynamic
%
% Y(:,k) is the opinion profile of time step k
%
% The opinion profile will be extended in the process

Y = X(:,1);

for i=1:t
	% pou = point (index) of uniformity
    n = size(X,1);
	pou = max(find(X(2:n,1)-X(1:n-1,1) ~= 1))+1;
	if isempty(pou)
        pou=0;
	end
	
	X(pou+1:pou+2+2*e,1) = [[pou+1:pou+2+2*e]-1]';	% add enough length
	X(pou+1:pou+2+2*e,2) = 1;
	X = X(1:pou+2+2*e,:);
    X = reduceOp(X); % to spare computation time
    
	X = boundConf(X,e);
    X = allOp(X);
    X(:,2) = 1;
	X(pou+1+e:pou+2+2*e,1) = [[pou+1+e:pou+2+2*e]-1]'; % flat tail to uniformity
    
    Y = [[Y; [size(Y,1):size(X,1)-1]'*ones(1,size(Y,2))] X(:,1)];
end


function X = boundConf(X,e)

% X = boundConf(X,e)
% X    : Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
% e    : bound of confidence
%
% boundConf calculates one step in the HK process of OD
% 
% 02/2003 Jan Lorenz

n = size(X,1); % number of opinions

A=[]; % will get the confidence matrix

% produce rows of the confidence matrix
for i = 1:n
   d = abs(X(:,1)-X(i,1));
   a = [d <= e ]'; % a is confidence set of opinion i
   a = a.*X(:,2)'; % a weighted with number of agents of opinions
   a = a/sum(a); 
   A = [A; a]; 
end
X(:,1)=A*X(:,1);


function X = reduceOp(X)

% X = reduceOp(X)
% X    : Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
%
% reduceOp collects all opinions with the same value in one row
% and adds the numbers of agents with this opinion
%
% 09/2002 Jan Lorenz
% Version 2.0

k=1; % row counter

while size(X,1) > k
   I = ( find(X(:,1) == X(k,1)) ); % row indices equal to X(k,1)
   X(k,2) = sum(X(I,2));
   X = X( sort([k setdiff([1:size(X,1)],I)]) ,:); % delete rows
   k=k+1; 
end


function X = allOp(X)

% X = allOp(X)
%
% X (in): Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
% X (out): column Vector
%
% allOp blows up an opinion profile with collected equal opinions
% to a vector with single opinions, for plotting reasons

M = [];
for i = 1:size(X,1)
    m = zeros(X(i,2),size(X,1));
    m(:,i) = ones(X(i,2),1);
    M = [M;m];
end
X = M*X;
X = X(:,1);
    