
function [D]= dist(C, im, r1, c1, S, m)

[rows, cols, chan] = size(im);
[x,y] = meshgrid(c1:(c1+cols-1), r1:(r1+rows-1));
x = x-C(4);  % x and y dist from cluster centre
y = y-C(5);
ds2 = x.^2 + y.^2; 

deltaLbarg = C(1) - im(:,:,1);
lbar = (im(:,:,1)+ C(1))/2;
clba = sqrt(im(:,:,2).^2+im(:,:,2).^2);
c2ba =sqrt(C(2).^2+C(3).^2);

cinversebarg= (clba  +c2ba)/2;
delCbar = clba - clba ;

G = 0.5*(1- sqrt((cinversebarg.^(7))./((cinversebarg.^(7))+25^7)));
a1barg = (1+G).* im(:,:,2);
a2barg = (1+G).* C(2);
c1barg = sqrt((a1barg.^2) + (im(:,:,3).^2));
c2barg = sqrt(a2barg.^2) + (C(3).^2);

Cpprod = (c2barg.*c1barg);
zcidx = find(Cpprod == 0);


h1bar = atan2d(im(:,:,3),a1barg);
h1bar= h1bar +2*pi*(h1bar < 0);  % rollover ones that come -ve
h1bar(find( (abs(a1barg)+abs(im(:,:,3)))== 0) ) = 0;

h2bar= atan2d(C(3),a2barg);
h2bar = h2bar+2*pi*(h2bar < 0);
h2bar(find( (abs(a2barg)+abs(C(3)))==0) ) = 0;


deltahbarg = (h2bar-h1bar);
deltahbarg = deltahbarg - 2*pi* (deltahbarg > pi );
deltahbarg = deltahbarg + 2*pi* (deltahbarg < (-pi) );
deltahbarg(zcidx ) = 0;
delHbarg = 2*sqrt(Cpprod).*sin(deltahbarg/2);


Hinversebarg = (h1bar + h2bar)/2;
Hinversebarg = Hinversebarg - ( abs(h1bar-h2bar)  > pi ) *pi;
Hinversebarg = Hinversebarg + (Hinversebarg < 0)*2*pi;
Hinversebarg(zcidx) = h2bar(zcidx)+h1bar(zcidx);


T = 1- 0.17.*cosd(Hinversebarg - 30) + 0.24 .*cosd(2.*Hinversebarg) + 0.32.*cosd(3.*Hinversebarg +6) - 0.20.* cosd(4.*Hinversebarg - 63);

sl = 1+ (0.015.* (lbar-50).^2./ sqrt(20+ (lbar - 50).^2));
sc = 1+ 0.045.*cinversebarg;
sh = 1+ 0.015.*cinversebarg.*T;
rt = -2 .* sqrt(cinversebarg.^(7)./(cinversebarg.^(7)+ 25.^7)) .* sind(60 .* exp(-1.* ((Hinversebarg - 275)./25).^2));


d1 = (deltaLbarg./sl).^2 + (delCbar./sc).^2 + (delHbarg ./ sh).^2 + rt.*(delCbar./sc).*(delHbarg ./ sh); %d1


D = sqrt(d1 + ds2/S^2*m^2);