function qc = dynamicqc(maxpts, qc, plot)
%
%    Copyright (C) 2005 by  
%    Sasha Gutfraind <ag362@cornell.edu> 
%    Distributed under the terms of the GNU Lesser General Public License 
%    http://www.gnu.org/copyleft/lesser.html 
% 

%todo: introduce layer counter and also an index for the last pt in the previous layer
%both will be modified by the findNextUntapped, but used elsewhere.  each pt will have a field for it later.
%alternatively, perhaps the same can be achieved by going over the qc just before drawing.

if nargin < 1
    maxpts = 2500;
end

if (nargin == 3)
   plotQC(qc);
   return;
end

[qc, metaQC] = setupInitialQC1(maxpts);

while  metaQC.qcSize < maxpts
    %debug
    metaQC.qcSize
   
    metaQC = updateMetaInfo(qc, metaQC);
    
    if metaQC.lastUntapped == -1
        break;
    end
    [qc, metaQC] = addNewILC(qc, metaQC);
end

qc = qc(1:metaQC.qcSize);

plotQC(qc);


function [qc, metaQC] = setupInitialQC1(maxpts, qcSize)
%this function sets up the initial QC with just 1 pt at the origin
%it might not be suitable for all cases, so then use setupInitialQC11 instead
newPt.x = 0;
newPt.y = 0;
newPt.tapped = 0; %don't forget this field!
newPt.parent{1}.x = 0;
newPt.parent{1}.y = 0;
qc(1) = newPt;

metaQC.qcSize             = 1;
metaQC.lastUntapped       = 1; %used to speedup search for new sites
metaQC.endOfLast          = 1; %used to mark end of last symmetric configuration

qc(maxpts).x = 0;  %klugdy initialization

function [qc, metaQC] = setupInitialQC11(maxpts, qcSize)
%this function sets up the initial QC with 11 pts
%it might not be suitable for all cases, so then use setupInitialQC1 instead
newPt.x = 0;
newPt.y = 0;
newPt.tapped = 10; %don't forget this field!
newPt.parent{1}.x = 0;
newPt.parent{1}.y = 0;
qc(1) = newPt;

for i=1:10
    if i <= 5
        sign = 1;
    else
        sign = -1;
    end
    newPt.x = 0 + sign*cos(2*pi/5*i);
    newPt.y = 0 + sign*sin(2*pi/5*i);
    newPt.parent = cell(0,0); %we will get this pt later.
    newPt.tapped = 0;
    qc(1+i) = newPt;
end

metaQC.qcSize             = 11;
metaQC.lastUntapped       = 2; %used to speedup search for new sites
metaQC.endOfLast          = 11; %used to mark end of last symmetric configuration

qc(maxpts).x = 0;  %klugdy initialization


function metaQC = updateMetaInfo(qc, metaQC)

found   = 0;
qcSize  = metaQC.qcSize;

for untapped=metaQC.lastUntapped:qcSize
    if qc(untapped).tapped < 10
        found = 1;
        break;
    end
end

if ~found
    metaQC.lastUntapped = -1
    warning('The QC is fully tapped!');
    return;
end

if untapped == metaQC.endOfLast + 1
    metaQC.endOfLast    = qcSize;
end
    
metaQC.lastUntapped       = untapped;


function [qc, metaQC] = addNewILC(qc, metaQC, qcSize)

startPtIndex = metaQC.lastUntapped;
startPt      = qc(startPtIndex);

%for now, we will just add one pt, and then quit
for i=startPt.tapped+1:10
    if i <= 5
        sign = 1;
    else
        sign = -1;
    end
    newPt.x = startPt.x + sign*cos(2*pi/5*i);
    newPt.y = startPt.y + sign*sin(2*pi/5*i);
    if ~isExistingPt(qc, metaQC.qcSize, newPt)
        newPt.tapped = 0;
        %we will get the parent when we try tapping this point.
        %newPt.parent{1}.x = startPt.x;
        %newPt.parent{1}.y = startPt.y;
        newPt.parent = cell(0,0);
        if acceptable(qc, metaQC, newPt)
            qc(metaQC.qcSize + 1) = newPt;
            metaQC.qcSize  = metaQC.qcSize + 1;
            %break;
        end
    else
        %unfortunately, this creates some duplication as step-son and
        %a pt from another later find each other
        lenPar = length(startPt.parent);
        startPt.parent{lenPar + 1}.x = newPt.x;
        startPt.parent{lenPar + 1}.y = newPt.y;
    end
    startPt.tapped = startPt.tapped + 1;
end
qc(startPtIndex) = startPt;
    

function answer = isExistingPt(qc, qcSize, newPt)
%we assume that newPt is not yet in the QC
%this algorithm is pretty bad - it scales as n^2, unfortunately,
%what we need is more sophisticated data structures, and we don't have them
answer = 0;

for i=1:qcSize
    delta = abs(qc(i).x - newPt.x) + abs(qc(i).y - newPt.y);
    if delta < 0.0001
        answer = 1;
        break;
    end
end

function acceptable = acceptable(qc, metaQC, newPt)

multiplier  = 1;%2.1222;
e0          = 0.3823*multiplier;
theta4      = ( 0.5*(1+sqrt(5)) )^4;

%the symmetrized version
%lastPtIndex = metaQC.endOfLast;

%the fully-dynamic version
lastPtIndex = metaQC.qcSize;

%the fixed-window version
%lastPtIndex = 11;

%dot products of x, y with each one of the mu's
dotProducts = newPt.x * [qc(1:lastPtIndex).x] + newPt.y * [qc(1:lastPtIndex).y];

tmp = exp(2*pi*i*2*theta4 * dotProducts);

%we need to forget about the imaginary part.  there should not be any imaginary part
%in a trully symmetrical acceptance function
f = sum(tmp)/lastPtIndex;

%asympetrical construction
%f = sum(tmp)/(metaQC.qcSize);

acceptable = real(f) >= e0;


function plotQC(qc)
%the iterator business is just for optimization of drawing

plot([qc.x], [qc.y], '.');

lineIterator = 0;
lenQC = length(qc);
xPts  = zeros(2,100);
yPts  = zeros(2,100);
for i = 1:lenQC
    numParents = length(qc(i).parent);
    for j = 1:numParents        
        lineIterator = lineIterator + 1;  
        xPts(:,lineIterator) = [qc(i).x ; qc(i).parent{j}.x];
        yPts(:,lineIterator) = [qc(i).y ; qc(i).parent{j}.y];
        if lineIterator == 100
            line(xPts, yPts, 'Color', 'b');
            xPts  = zeros(2,100);
            yPts  = zeros(2,100);
            lineIterator = 0;
        end 
    end
end
%radius = sqrt((qc(lenQC).x)^2 + (qc(lenQC).y)^2) + 20; %radius plus margin
radius = 25;
set(gca, 'XLim', [-radius radius], 'YLim', [-radius radius])


line([xPts(:,1:lineIterator)], [yPts(:,1:lineIterator,:)], 'Color', 'b');
