This worked for me.
function test1
vn=1
vs=1 ;
M=100 ;
th1=1000;
Dr1=0;
D1=0
y = rand(1,M);
x = rand(1,M);
D=distortion(vn,vs,M,th1,x,y,Dr1,D1)
end
function D=distortion(vn,vs,M,th1,x,y,Dr1,D1)
for i=1:M;
for j=1:M;
if i~=j
d1(i,j)=sqrt(((y(j)-y(i))^2+(x(j)-x(i))^2));
Dr1=exp(-d1(i,j)/th1)+Dr1;
end
end
end
for i=1:M
ds1(i)=sqrt(((y(j)-0)^2+(x(j)-0)^2));% point source at 0 0
cor2=exp(-ds1(i)/th1);
D1=cor2+D1;
end
D=(vs^2)-((vs^4)/(M*((vs^2)+(vn^2))))*(2*(D1-1))+((vs^6)/((M^2)*((vs^2)+(vn^2)))*Dr1);
end
What did you pass in for X and Y?
