Diffing "Stuning02" and "Stuning03"

function xy = solver(kd,bx)
npts0 = length(kd);
dx=bx(2);
dy=bx(4);
D=dx+1i*dy;
pm=D/2;
xy=zeros(1,npts0)+pm;

d=0.1*max(dx,dy);
x=0:d:dx;
y=0:d:dy;
[xx,yy] = ndgrid(x,1i*y);
XX=xx(:)+yy(:);

dr=abs(D);
kd(find(kd>dr))=dr;
mx=max(kd-eye(npts0));

ismart=find(mx>=0);
if isempty(ismart)|dx+dy<1e-3
	xy=zeros(npts0,2);
	return;
elseif length(ismart)==2
	xy=zeros(npts0,2);
	f=kd(ismart(1),ismart(2))/dr;
	xy0=[0;D*f];
	xy(ismart,:)=[real(xy0) imag(xy0)];
	return
end
npts=length(ismart);
kd=kd(ismart,ismart);

S=(kd-eye(npts))>=0;
nX=length(XX);
xOnes=ones(nX,1);
xZeros=zeros(nX,1);
sCut = 1;

olds = 1.0e30;
[dum,p]=sort(-max(kd));
S=S(p,p);
kd=kd(p,p);
[v,w]=sort(p);
if npts<50
	xxx=11;
	xxx=12;
else
	xxx=10;
end;    
for iter=1:xxx
    for index = 1:npts
        I = find(S(:,index)); 
        if ~isempty(I)
            st = sum(abs(abs(XX(:,ones(size(I)))-xy(xOnes,I))-kd(xZeros+index,I)),2);
            [null,minloc] = min(st);
            xy(index) = XX(minloc);
        end
    end
    [s,M]=findstrain(S,xy,kd);
    if((abs(s-olds)/(s+1)) < 0.002) | s < sCut
        break;
    end
    olds = s;
end
xy=xy.';
xy=xy(w,:);
kd=kd(w,w);
S=S(w,w);

olds = 1.0e50;
oldxy = xy;
alpha = 1;

for iter=1:40
    [s,M]=findstrain(S,xy,kd);
    if( iter > 5 & abs(s-olds) / (s+1) < 0.0001 ) | s < sCut
        xy = oldxy;
        break;
    end
    
    if( s > olds )
        xy = oldxy;
        break;
    end
    oldxy = xy;
    olds = s;
    for indx = 1:npts
        I = find(S(:,index)>0);
        dX=xy(I)-xy(indx);
        Dx=mean(dX.*(1-kd(I,index)./(abs(dX)+eps)));
%            Dx=mean(dX./(abs(dX)+eps) .* M(I,indx));
        
        xy(indx) = xy(indx) + alpha*Dx;
        xy(indx) = min(dx,max(0,real(xy(indx))))+1i*min(dy,max(0,imag(xy(indx))));
    end
end

alpha = 0.1;
ngrad = 40;

obj = 1.0e20*ones(ngrad,1);
xybig = zeros(npts,ngrad);

[objnew,strainMatrix]=findstrain(S,xy,kd);
obj(1) = objnew;
xybig(:,1) = xy;
gradient = zeros(npts,1);

for ij=2:ngrad
   [s,sM,dM]=findstrain(S,xy,kd);
   gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM));
   
   xy = xy - alpha  * gradient;
   xy = min(max(real(xy),0),dx) + j * min(max(imag(xy),0),dy);
   [objnew,strainMatrix]=findstrain(S,xy,kd);
   
   obj(ij) = objnew;
   xybig(:,ij) = xy;
   if( obj(ij) > obj(ij-1) )
       xy = xybig(:,ij-1);
       alpha = alpha / 3;
   end
   [bestobj,index] = min(obj);
   if( ij > 3 & index < 2 ) | bestobj < sCut
      break
  end
end
[bestobj,index] = min(obj);
xy = xybig(:,index);
xy0 = zeros(npts0,1);
xy0(ismart)=xy;
xy=[real(xy0) imag(xy0)];

function [s,strainMatrix,dX]=findstrain(S,xy,kd)
sI=find(S);
strainMatrix=S;
dX=S;
[X1,X2]=meshgrid(xy); 
dX(sI)=X1(sI)-X2(sI);
strainMatrix(sI) = abs(dX(sI))-kd(sI);
s = sum(abs(strainMatrix(sI)));