smallopt3

function xy = solver(kd,bx)

% catch trivial case:
if size(kd,1) < 2
    xy = [mean(bx(1:2)) mean(bx(3:4))];
    return
end

%intialize xy
for iter = 1:10
    best = Inf;
    thisxy = [rand(length(kd),1)*(bx(2)-bx(1))+bx(1) rand(length(kd),1)*(bx(4)-bx(3))+bx(3)];
    s=findstrain(thisxy(:),kd,bx);
    if s < best
        xy = thisxy;
        best=s;
    end
end

[s,M] = findstrain(xy(:),kd,bx);
xyold=Inf*xy;

if s==0
    return
end

iter = 1;
while iter
    fx=0;fy=0; % intialize force vectors
    for indx = 1:length(kd)
        dx=xy(:,1)-xy(indx,1);
        dy=xy(:,2)-xy(indx,2);
        N = sqrt(dx.^2 + dy.^2);
        N(N==0)=1;
        dx=dx./N;
        dy=dy./N;
        % now have direction vectors (dx,dy) from xy(indx,:) to each other pt
        dx = dx .* M(:,indx);
        dy = dy .* M(:,indx);
        dx = sum(dx);
        dy = sum(dy);
        % now have strain vectors
        xy(indx,:) = xy(indx,:) + .05*[dx dy]; % move point
        % clip if outside box
        xy(indx,1) = min(max(xy(indx,1),bx(1)),bx(2));
        xy(indx,2) = min(max(xy(indx,2),bx(3)),bx(4));
        
    end

    iter = iter+1;
    Nx = abs(xy(:,1)-xyold(:,1))/(bx(2)-bx(1));
    Ny = abs(xy(:,2)-xyold(:,2))/(bx(4)-bx(3));
    N = max([Nx;Ny]);
    xyold=xy;
    if iter > 50 | N < .0001
        iter = 0;
    end
    
end

function [s,strainMatrix]=findstrain(inp,kd,bx)
xy = reshape(inp,[length(inp)/2 2]);
% Calculate distance matrix
x=xy(:,1);
y=xy(:,2);

% catch points outside box
if any(x<bx(1) | x>bx(2) | y<bx(3) | y> bx(4))
    s = Inf;
    return
end

[XY1,XY2] = meshgrid(x,y);
dx = XY1 - XY1';
dy = XY2 - XY2';
dist = sqrt(dx.^2 + dy.^2); 

% Calculate strain matrix
strainMatrix = dist - kd;
strainMatrix(kd < 0) = 0;
s = sum(abs(strainMatrix(:)));
return