%%4D CF from 'Data_3GHz_1cm.txt'

clear 
clc

%% 1cm 
A = dlmread('Data_3GHz_Hy_probe10_1cm.txt'); %read correlation data from file 
% A(theta,Hx)

L=36; % number of paddle rotation

Hx=(A(:,2));

nx=60;
ny=60;

n=1
for k=1:L    
    for x = 1:nx
        for y = 1:ny
        
        data(x,y,k)=Hx(n);
       
        n=n+1
           
        end
    end
end

for t=1:L
    t
    for x1=1:nx
        for y1=1:ny
            for x2=1:nx
                for y2=1:ny
                     
                    CF(x1,y1,x2,y2,t)=data(x1,y1,t)*conj(data(x2,y2,t));
        
                end
            end
        end
    end
end

for x1=1:nx
    x1
    for y1=1:ny
        for x2=1:nx
            for y2=1:ny
                     
                    R1(x1,y1,x2,y2)=mean(squeeze(CF(x1,y1,x2,y2,:))); % CF4D = R(x1,y1,x2,y2)
            
            end
        end
    end
end

figure
imagesc(squeeze(abs((R1(30,30,:,:))))'); %CF2(x1,y1,x2,y2)
title('3GHz - 4D CF(y1,x1,30,30) probe10 at 1cm')
set(gca,'YDir','normal')
colormap('jet');

%% 10cm 
A = dlmread('Data_3GHz_Hy_probe10_10cm.txt'); %read correlation data from file 
% A(theta,Hx)

%L=1; % number of paddle rotation

Hx=(A(:,2));

nx=60;
ny=60;

n=1
for k=1:L    
    for x = 1:nx
        for y = 1:ny
        
        data(x,y,k)=Hx(n);
       
        n=n+1
           
        end
    end
end

for t=1:L
    t
    for x1=1:nx
        for y1=1:ny
            for x2=1:nx
                for y2=1:ny
                     
                    CF(x1,y1,x2,y2,t)=data(x1,y1,t)*conj(data(x2,y2,t));
        
                end
            end
        end
    end
end

for x1=1:nx
    x1
    for y1=1:ny
        for x2=1:nx
            for y2=1:ny
                     
                    R3(x1,y1,x2,y2)=mean(squeeze(CF(x1,y1,x2,y2,:))); % CF4D = R(x1,y1,x2,y2)
            
            end
        end
    end
end

figure
imagesc(squeeze(abs((R3(30,30,:,:))))'); %CF2(x1,y1,x2,y2)
title('3GHz - 4D CF(y1,x1,30,30) probe10 at 10cm')
set(gca,'YDir','normal')
colormap('jet');