function particles( length, time, win, np, ntr)
%PARTICLES
% Implementation of a "particle"-like on an initially
% Reversible Cellular Automaton with a nbits-period
% clock corresponding to a "Rotating" Lattice
% Sequential application of several arbitrary permutation
% operators guarantees step-by-step reversibility of nbit-substrings.
% Periodic Boundary Conditions are naturally employed
% Long term reversibility is currently broken due to the
% movement & collision rules employed
% Stochasticity of the "particle" trajectories is the
% result of randomness in the in. condition and
% the mixing of states during nbits-step time evolution
% Individually tracked particle trajectories exhibit a strongly
% non-stationary statistics for small times
% To avoid overcrowding, only half of the lattice positions can
% be 'excited' with a p-state at anyone time.
%
% INPUT :
% length = total string array length in multiples of 2*nbits
% time = total simulation time
% win = time window for graphics
% np = no. of particle states (<= max. saturation level = length/6)
% ntrack = ptr of part. of which the trajectory is to be tracked (in [1, np])
%
%TEST :
% single part. : try particles(1000,2400,600,1,1);
% many part. : try particles(1000,2400,600,5,1);
%
% T. E. Raptis, DAT-NCSRD (c) 2010
% Reference : http://cag.dat.demokritos.gr/publications/RCAposter.pdf
% rtheo@dat.demokritos.gr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc, close all
nbits = 3;
nrules = 4;
dim = 2^nbits;
disp('"PARTICLES" ');
disp('========================================');
% Rules for substring operators defined as permutations over
% the integers : Cardinality of the rule set increases as (2^nbits)!
% we can apply as many rules we want in any sequence
% we like.
rules(1,:) = [0,1,2,7,4,5,6,3]; % Reversible AND
rules(2,:) = [0,5,6,3,4,1,2,7]; % Reversible XOR
rules(3,:) = [4,5,6,7,0,1,2,3];
rules(4,:) = [7,1,6,2,5,3,4,0];
% Expand rules in the extended alphabet to hold "particle" states
rules = [rules,dim+rules];
% locate non-invariant substrings
index = 0:2*dim-1;
disp('Langton Parameters :');
for i=1:nrules
mask(i,:) = 1-(rules(i,:) == index);
lp = sum(mask(i,:))/(2*dim);
disp(['Rule ',num2str(i),' : l = ',num2str(lp)]);
end
% store all substring binary patterns
bin = generator(nbits);
%Expand patterns in the extended alphabet of p-states
tmp = bin; tmp(3,:) = 2;
bin = [bin,tmp];
% create initial condition
l0 = mod(length,2*nbits);
if l0~=0
length = length - l0 + 2*nbits;
disp('WARNING!');
disp('Total length adapted to lattice constrains');
disp(['New lattice length : ',num2str(length),' cells']);
if length < 30 disp('Lattice too small! Initialization may not end'); end
end
L = length/3;
if np<round(length/2)
% create full length random initial condition for v-states
s = zeros(1,length);
for i=1:length
if rand>0.5 s(i) = 1; end
end
% fill in with particle states
aux = zeros(1,L/2);
while sum(aux)<np
index = round(rand*(L/2-1))+1;
aux(index) = 1;
end
p = [];
for i=1:L/2
if aux(i)==0
p = [p,[0,0]];
else
p = [p,[1,0]];
end
end
ptr0 = track(p,ntr);
if ntr>1
f = find(p==1);old = diff(f);
end
% setup random spin configuration
spin = zeros(1,L);
for i=1:L
if p(i)==1
if rand>0.5 spin(i) = 1;
else spin(i) = -1;
end
end
end
if sum(p)~=np disp('Incorrect Particle no.! Exiting'), return, end
s(3*find(p==1)) = 2;
disp(['Particle Density = ',num2str(100*np/L),' %']);
else
disp('Saturation level reached!');
disp('Increase Lattice length or decrease no. of particles')
disp('Exiting');
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAIN LOOP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k = 1; ridx = 1;
CA = zeros(length,win); % hold binary image of ca evolution
TV = zeros(length/3,win); % hold nbits-colored image
PO = zeros(length/3,win); % hold particle trajectories
Inv = ones(length,win); % auxilliary array for bit inversion
CInv = ones(length/3,win);
b = 2.^(0:nbits-1); % alphabet base for representation change
scrsz = get(0,'ScreenSize');
for t = 1:time
s = circshift(s',-1)'; % Left cyclic shift of lattice sites
[s,spin,ptr] = move(s, spin, nbits, ntr, t); % swap "excited" states for particle motion
LR(t) = ptr0-ptr; ptr0 = ptr; % record tracked particle motion
v = transform(s,b,L); % change alphabet - decrease length to L
v = 1+rules(ridx,v); % apply rule mapping
s = decode(v,bin,L); % change back to binary - increase length to 3L
p = (v>8); % extract part. positions
if ntr>1
[dp, old] = rdist(p, old, L, t); % find part. relative distances
end
CA(:,k) = s;
TV(:,k) = v;
PO(:,k) = p;
if sum(p)~=np
% plotter(Inv-CA,TV,CInv-PO,scrsz);
error(['Particle conservation violation at step : ',num2str(t)]);
end
if mod(t,win)==0 k=plotter(Inv-CA,TV,CInv-PO,scrsz);pause(.1);end
ridx = mod(ridx,4)+1;
k = k + 1;
end
% reconstruct tracked particle signal
x0 = 0;
for i=1:time
tr(i) = x0 + LR(i);x0 = tr(i);
end
disp(['Neg. steps = ',num2str(sum(LR<0)),' Pos. steps = ',num2str(sum(LR>0))])
figure(4)
set(gcf, 'Position',[scrsz(3)/2,1,(scrsz(3)-10)/2,(scrsz(4)-10)/2]);
hist(tr,30);
title(['"Statistics of part. ',num2str(ntr),' : Mean = ',num2str(mean(tr)),' St. Dev. = ',num2str(std(tr))]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% END MAIN
%%%%%%%%%%%%%%%%%%%%%%%%%
function v = transform(s,b,L)
j = 1;
for i=1:L
v(i) = 1+dot(b,s(j:j+2));
j = j + 3;
end
end
function s = decode(v,b,L)
s = [];
for i=1:L
s = [s,b(:, v(i))'];
end
end
function [s,p,ptr] = move(s, sp, n, ntr, t)
len = length(s);
s3 = [s(1:3:len);s(2:3:len)>1;s(3:3:len)];
% - hold prev. particle pos.
p = s3(2,:); n0 = sum(p);
L = length(p);
% implement periodic b. c. for collision mask
s3 = [s3(:,L),s3,s3(:,1)];
% clean up previous 'excited' p-states
s = ceil(s/2);
% track all possible movements
for i=1:L
if s3(2+sp(i),i)
idx = i+1; if i==L idx=1; end
if p(i) && p(idx) == 0
p(i)=0; p(idx)=1;
sp(idx) = sp(i); sp(i) = 0;
end
else
idx = i-1; if i==1 idx=L; end
if p(i) && p(idx) == 0
p(i)=0;p(idx)=1;
sp(idx) = sp(i); sp(i) = 0;
end
end
end
for i=1:L
if sp(i)==0 p(i)=0;end
end
%if sum(p)~=n0 disp(['conservation violation in Move : step=',num2str(t)]), end
ptr = track(p,ntr);
% recreate new 'excited' p-states
for i=1:L
if p(i) s(3*i) = 2; end
end
end
function ptr = track(p, ntr)
f = find(p==1); ptr = f(ntr);
end
function [dp, op] = rdist(v, op, L, t)
p = (v>8);
dp = zeros(1,L-1);
f = find(p==1);
dp = diff(f);
dd = abs(dp - op);
op = dp;
if any(dd>2)
disp(['Anomalous transition : step = ',num2str(t)]);
end
end
function ss=generator(n)
% generates all binary patterns in the range [0, (2^n)-1]
p2=1;
len = 2^n;
ss = zeros(n,len);
for i=1:n
ss(i,:) = floor(mod((0:len-1)/p2,2));
p2 = 2*p2;
end
end
function k=plotter(CA,TV,PO,scr)
figure(1)
set(gcf, 'Position',[1,(scr(4)-10)/2,(scr(3)-10)/2,(scr(4)-10)/2]);
imagesc(CA);
colormap gray;
figure(2)
set(gcf, 'Position',[scr(3)/2,(scr(4)-10)/2,(scr(3)-10)/2,(scr(4)-10)/2]);
imagesc(TV);
figure(3)
set(gcf, 'Position',[1,1,(scr(3)-10)/2,(scr(4)-10)/2]);
imagesc(PO)
colormap gray
k = 0;
end
