space-domain back-projection algorithm implementation

hello,

I haved designed a millimeter wave imaging system. I have data about received electrical field at each receiver point resulted from each transmitter at each particular frequency. my system has 16 transmitters and 88 receivers and measurements has been done at 11 uniformly-spaced frequencies between 10 to 20 GHz. I have attached Efield file (needed to implement following code). I just want to know that what is problem with my code? its formula is really straightforward but I do not know why I can not implement it.

also here is the link to another file containing nearfield textfiles extracted from feko and saveData.m file needed to create Efield data:

https://drive.google.com/file/d/1uhwtq0UoonoiMw6fwzRpyRhelRSmbB1k/view?usp=sharing

Here it is my code related to back-projection algorithm:

load Efield

% Efield is a 3-dimentional 1408*6*11 matrix which columns 1,2,3 show x,y,z position of
% our receiver and column 4,5,6 show electrical field in x,y,z direction.
% third dimention represents frequency.
% (pay attention that for each receiver would be 16 row in order to have
% all information from each transmitter at each frequency)

Tx(1)=29.250000e+00;
Ty(1)=5.250000e+00;

Tx(2)=29.250000e+00;
Ty(2)=-5.2500000e+00;

Tx(3)=-29.250000e+00;
Ty(3)=5.250000e+00;

Tx(4)=-29.250000e+00;
Ty(4)=-5.250000e+00;

Tx(5)=24.750000e+00;
Ty(5)=15.750000e+00;

Tx(6)=24.750000e+00;
Ty(6)=-15.750000e+00;

Tx(7)=-24.750000e+00;
Ty(7)=15.750000e+00;

Tx(8)=-24.7500000e+00;
Ty(8)=-15.750000e+00;

Tx(9)=15.7500000e+00;
Ty(9)=24.750000e+00;

Tx(10)=15.750000e+00;
Ty(10)=-24.750000e+00;

Tx(11)=-15.750000e+00;
Ty(11)=24.750000e+00;

Tx(12)=-15.750000e+00;
Ty(12)=-24.750000e+00;

Tx(13)=5.250000e+00;
Ty(13)=29.250000e+00;

Tx(14)=5.250000e+00;
Ty(14)=-29.250000e+00;

Tx(15)=-5.250000e+00;
Ty(15)=29.250000e+00;

Tx(16)=-5.250000e+00;
Ty(16)=-29.250000e+00;


numF = 11;
x = (-15:0.2:15);
y = (-15:0.2:15);
coor = zeros(numel(x)*numel(y),2);
count = 1;

for x1 = x,
    for y1 = y,
       coor(count,1) = x1;
       coor(count,2) = y1;
       count = count + 1;
    end
end

freq = 10e9:(10e9)/(numF-1):20e9;
or = zeros(numel(x)*numel(y),3);
c = 3e8;

for i = 1:numel(x)*numel(y),    
    for tx = 1:16  ,
        rt = sqrt((Tx(tx)-coor(i,1))^2 + (Ty(tx)-coor(i,2))^2 + 15^2);
        for rx = 1:88
            for f = 1:numF,
                rr = sqrt((real(Efield((tx-1)*88 + rx,1,f))-coor(i,1))^2 + (real(Efield((tx-1)*88 + rx,2,f))-coor(i,2))^2 + 15^2);
                r = (rr + rt)/100;
                
                or(i,1) = or(i,1) + Efield((tx-1)*88 + rx,4,f)*exp(1j*(2*pi*freq(f)/c)*r);
                or(i,2) = or(i,2) + Efield((tx-1)*88 + rx,5,f)*exp(1j*(2*pi*freq(f)/c)*r);
                or(i,3) = or(i,3) + Efield((tx-1)*88 + rx,6,f)*exp(1j*(2*pi*freq(f)/c)*r);
            end
        end
    end    
end

image = zeros(numel(x)*numel(y),1);

image(:) = nthroot(abs(or(:,1)).^2 + abs(or(:,2)).^2 + abs(or(:,3)).^2 , 2);

scatter(coor(:,1),coor(:,2),[],image(:),'filled');
colormap(gray);
colorbar;

xlabel('X position');
ylabel('Y position');
title('2D Intensity plot');