function ClosePoints_mzip
%   MATLAB zip file, a self-extracting MATLAB archive.
%   Usage: Run this file to recreate the original directory.

    fname = mfilename;
    fin = fopen([fname '.m'],'r');
    dname = fname(1:find(fname=='_',1,'last')-1);
    mkdir(dname);
    mkdir([dname filesep 'lib'])
    addpath(dname)
    addpath([dname filesep 'lib'])    
    
    L = fgetl(fin);
    while length(L) < 2 || ~isequal(L(1:2),'%%')
        L = fgetl(fin);
    end
    while ~isequal(L,'%% EOF')
        F = [dname filesep L(4:end)];
        disp(F)
        fout = fopen(F,'w');
        L = fgetl(fin);
        while length(L) < 2 || ~isequal(L(1:2),'%%')
            fprintf(fout,'%s\n',L);
            L = fgetl(fin);
        end
        fclose(fout);
    end
    fclose(fin);
end
%% Center.m
function [d,kout,jout] = Center(z,din)
% Center(z,d) is used by DivCon to examine the
% strip of half-width d about the center point.
    
    n = length(z);
    m = floor(n/2);
    xh = real(z(m));
    d = din;

    [~,p] = sort(imag(z));
    z = z(p);

    s = [];
    ks = [];
    ns = 0;
    for k = 1:n
        if abs(real(z(k)) - xh) <  d
            ns = ns+1;
            s(ns,1) = z(k);
            ks(ns,1) = k;
        end
    end

    kout = NaN;
    jout = NaN;
    for k = 1:ns
        for j = k+1:ns
            if (imag(s(j)) - imag(s(k))) < d && abs(s(k) - s(j)) < d
                d = abs(s(k) - s(j));
                kout = p(ks(k));
                jout = p(ks(j));
            end
        end
    end
end
               
%% DivCon.m
function [d,kout,jout] = DivCon(z,~)
% [d,k,j] = DivCon(z) finds indices k and j so that
% z(k) and z(j) minimizes d = abs(z(k)-z(j)).
%
% DivCon employs a recursive divide and conquer algorithm with
% complexity O(n*log(n)) where n = length(z).
% [d,k,j] = DivCon(z,true) is a recursive call with ascending real(z).
%
% Center(z,d) examines the central strip of half-width d.
% Pairs(z) computes the same thing as DivCon, but has complexity O(n^2).

    n = length(z);
    if n <= 3
        [d,kout,jout] = Pairs(z);
    else
        if nargin < 2
            [~,p] = sort(real(z));
            z = z(p);
        else
            p = (1:n)';
        end

        m = floor(n/2);

        % Left half
        [dl,kl,jl] = DivCon(z(1:m),true);

        % Right half
        [dr,kr,jr] = DivCon(z(m+1:end),true);

        if dl < dr
            d = dl; k = kl; j = jl;
        else
            d = dr; k = kr+m; j = jr+m;
        end

        % Center strip
        [ds,ks,js] = Center(z,d);
        if ds < d
            d = ds; k = ks; j = js;
        end

        kout = min(p(k),p(j));
        jout = max(p(k),p(j));
      
    end
end
%% Pairs.m
function [d,kout,jout] = Pairs(z)
% [d,k,j] = Pairs(z) finds indices k and j so that
% z(k) and z(j) minimizes d = abs(z(k)-z(j)).
%
% Pairs employs a brute force algorithm that examines all possible 
% pairs.  This has complexity O(n^2) where n = length(z).
% DivCon(z) computes the same thing, but has complexity O(n*log(n))
% and consequently is much faster.  

    n = length(z);
    d = Inf; 
    for k = 1:n
        for j = k+1:n
            if abs(z(k) - z(j)) < d
                d = abs(z(k) - z(j));
                kout = k;
                jout = j;
            end    
        end
    end
end
%% xDivCon.m
function xDivCon(z,seed)
    if nargin == 0   
        disp('------')
        for n = 4:20
            disp(n)
            for seed = 0:20
                rng(seed)
                z = rand(n,1) + rand(n,1)*1i;
                xDivCon(z,seed)
            end
        end
        return
    end

    [db,kb,jb] = Pairs(z);

    [d,k,j] = DivCon(z);
 
    if ~(d == db && k == kb && j == jb)
        n = length(z);
        disp([n seed kb jb k j])
        fprintf('Pairs:  %8.4f  %6.4f+%6.4fi  %6.4f %6.4fi\n', ...
            db,real(z(kb)),imag(z(kb)),real(z(jb)),imag(z(jb)))
        fprintf('DivCon: %8.4f  %6.4f+%6.4fi  %6.4f %6.4fi\n', ...
            d,real(z(k)),imag(z(k)),real(z(j)),imag(z(j)))
    end
    
end
%% EOF