// sfb 30 Aug 2009; updated 14 Sep 2009, 27 Jan 2010
lines(0);

function circle(x, y, r, colorIndex)
	xfarc(x-r,y+r,2*r,2*r,0,64*360);
	harc = gce();
	harc.background=colorIndex;
endfunction

function circles(x, y, r, c)
// x y r columns c row.
	a1 = zeros(x);
	a2 = 64*360*ones(x);
	arcs = [(x-r)';(y+r)';(2*r)';(2*r)';a1';a2'];
	xfarcs(arcs,c);
endfunction;

// theta and colors are column vectors
function track(theta, colors)
	x = cos(theta);
	y = sin(theta);
	r = 0.05*ones(theta);
	circles(x,y,r,colors');
endfunction;

function blink(theta, colors)
	r = max(0,sin(theta))/9;
	circles(positions(:,1), positions(:,2), r, colors')
endfunction

function y = initPositions(n) // blinking fire flys
	y = zeros(n,2);
	col = 1;
	row = 1;
	for i = 1:n,
		y(i,1) = -1 + col/4;
		y(i,2) = 3 - row/4;
		col = col + 1;
		if col == 8,
			col = 1;
			row = row + 1;
		end
	end;
endfunction;

function [theta, omega, colors] = init(n)
	rand('uniform')
	theta = 2*%pi*rand(n,1);
	rand('normal')
	omega = sort(1.5 + 0.5*rand(n,1));
	colors = ones(theta)*color('red');
	colors((1+floor(n/4)):(1+floor(n/2))) =color('orange');
	colors((1+floor(n/2)):(1+floor(3*n/4))) =color('green');
	colors((1+floor(3*n/4)):n) =color('blue');
endfunction;

function y = step(theta, deltat, K)
	n = size(theta,1);
	thetadot = omega + (K/n)*sum(sin(ones(n,1)*theta' - theta*ones(1,n)),2);
	y = theta + thetadot*deltat;
endfunction;

function blackball(theta)
	n = size(theta,1);
	x = sum(cos(theta))/n;
	y = sum(sin(theta))/n;
	b = color('black');
	circle(x,y,0.05*sqrt(2),b);
endfunction;

function phi = centerofmass(theta)
	n = size(theta,1);
	x = sum(cos(theta))/n;
	y = sum(sin(theta))/n;
	phi = atan(y,x);
endfunction;

// the whole thing
function bigshow(n, K)
deltat = 0.05;
[theta, omega, colors] = init(n);
positions = initPositions(n);
scf(0);
clf;
f=gcf();
f.pixmap='on';
isoview(-1.1,-1.1,1.1,3.1);
square(-1.1,-1.1,1.1,3.1);
track(theta, colors);
blink(theta, colors);
show_pixmap();
for i = 1:1000,
	clf();
	isoview(-1.1,-1.1,1.1,3.1);
	square(-1.1,-1.1,1.1,3.1);
	theta  = step(theta, deltat, K);
	// ORACLE CHANGES HERE 
	track(theta, colors);
	blink(theta, colors);
	show_pixmap();
end;
endfunction;

function phi = jogger(omega1, omega2, K, style)
deltat = 0.05;
omega = [max(omega1, omega2); min(omega1, omega2)]
theta = [0; %pi]
colors = [color('red'); color('blue')];
scf(1);
clf;
f=gcf();
f.pixmap='on';
t=0:%pi/40:2*%pi;
plot(cos(t),sin(t),'k-');
isoview(-1.1,-1.1,1.1,1.1);
square(-1.1,-1.1,1.1,1.1);
track(theta, colors);
show_pixmap();
cc = 500;
zip = zeros(10, cc/10);
for i = 1:cc,
	clf();
	plot(cos(t),sin(t),'k-');
	isoview(-1.1,-1.1,1.1,1.1);
	square(-1.1,-1.1,1.1,1.1);
	theta  = step(theta, deltat, K);
	if style == 0,
		track(theta, colors);
	elseif style == 1;
		track(theta, colors);
		blackball(theta);
	elseif style == 2;
		phi = centerofmass(theta); // center of mass angle
		track(theta - phi + %pi/2, colors);
	end;
	show_pixmap();
	zip(i) = theta(1) - theta(2);
end;
zip = zip'
phi = zip(:,10); // prints a hundred lines.
endfunction;
