%% This code was used to generate Fig. 4.16-a.
%  It implements the delay-and-sum (DAS) method described in Sections 4.5.2,
%  4.5.3, and 4.5.4
%  Version Date: Aug. 28, 2016
%  Author: Natalia K. Nikolova 
%% get excitation pulse and generate time-shifted versions
clear all;
close all;
c = 2.99792458e+8; % speed of light in m/s
degtorad = pi/180; % conversion coefficient from degrees to radians
radtodeg = 180/pi; % conversion coefficient from  radians to degress
%% READ DATA
Nskipline = 7; % skip first 7 lines in data files (no data there)
DataFormat = '%d %g %g';
sizeA = [3 Inf];
ext = '.txt';
mi(1)=7; % mapping indeces, simulated data idx -> circ idx 
mi(2)=3;
mi(3)=5;
mi(4)=1;
mi(5)=2;
mi(6)=6;
mi(7)=8;
mi(8)=4;
%% coordinates of point sources in 8-element circular array
ca=zeros(8,2);
ca(1,:)=[100 0]; % in mm
ca(2,:)=[70.5 70.5]; % in mm
ca(3,:)=[0 100]; % in mm
ca(4,:)=[-70.5 70.5]; % in mm
ca(5,:)=[-100 0]; % in mm
ca(6,:)=[-70.5 -70.5]; % in mm
ca(7,:)=[0 -100]; % in mm
ca(8,:)=[70.5 -70.5]; % in mm
ca=ca*1e-3; % convert to m
%% read data
disp('loading data...');
NtMax=15001; % number of temporal samples
DeltaT=3.3356409519815207e-012; % temporal sampling step
PB=zeros(NtMax,8,8); % baseline (background, RO) signals
P=zeros(NtMax,8,8); % signals obtained with the OUT (obect under test)
for i=1:8 % Rx index
    idx1 = num2str(i); % Rx index sweeps from 1 to 8
    idx2 = num2str(1); % Tx index is 1
    fname = strcat('p',idx1,idx2,'b',ext); % full file name to read baseline data from
    fid = fopen(fname);
    for k=1:Nskipline
        tline = fgets(fid); % first lines do not contain data
    end
% read baseline data when source 1 excited
    A = fscanf(fid,DataFormat,sizeA);
    A = A';
% extract the vector of the number of time samples and the vector of actual
% temporal values
    if (i==1) TS=A(:,1); T=A(:,2);
    end
% populate baseline response when source 1 is excited
    newi = mi(i); % Rx index
    newj = mi(1); % Tx index is 1
    PB(:,newi,newj) = A(:,3); % baseline response at receiver i (newi)
    fclose(fid);
    for j=1:8 % Tx index
        idx2 = num2str(j);
        fname = strcat('p',idx1,idx2,ext); % full file name to read OUT data from
        fid = fopen(fname);
        for k=1:Nskipline
            tline = fgets(fid); % first lines do not contain data
        end
% read OUT data
        A = fscanf(fid,DataFormat,sizeA);
        A = A';
        newj = mi(j);
% populate OUT data
        P(:,newi,newj) = A(:,3); % OUT probe data at i-th Rx with j-th excitation
        fclose(fid);
    end
end
%% populate baseline data for all sources using the circular symmetry
j=mi(1);
for n=1:7
    jj=j+1; % Tx index
    if (jj>8) jj=1;
    end
    for i=1:8
        ii=i+1; % Rx index
        if (ii>8) ii=1;
        end
        PB(:,ii,jj)=PB(:,i,j);
    end
    j=j+1;
    if (j>8) j=1;
    end
end
%% subtract baseline data from OUT data
P=P-PB;
%% generate DAS image
mindelay=3.3356e-10; % min expected delay in sec
maxdelay=1.1341e-09; % max expected delay in sec
ds = 2e-3; % discretization step in m
extenthalf = 70e-3; % half spatial extent in m
nshalf = round(extenthalf/ds); % number of discrete steps along half-axis
I=zeros(2*nshalf+1); % image matrix
xindex=0;
for ix=-nshalf:nshalf % x-coordinate focus-point sweep
    xp=ix*ds; % x-coordinate of focus point
    xindex=xindex+1; % x-coordinate focus-point index
    yindex=0;
    for iy=-nshalf:nshalf % y-coordinate focus-point sweep
        PFS=zeros(NtMax,1); % initialize sum of filtered responses vector
        yp=iy*ds; % y-coordinate focus point
        yindex=yindex+1; % y-coordinate focus-point index
        for j=1:8 % Tx index
            xt=ca(j,1);yt=ca(j,2); % Tx coordinates
            txd=sqrt((xp-xt)^2+(yp-yt)^2); % Tx-to-focus-point distance
            for i=1:8 % Rx index
                xr=ca(i,1);yr=ca(i,2); % Rx coordinates
                rxd=sqrt((xp-xr)^2+(yp-yr)^2); % Rx-to-focus-point distance
                delay=(txd+rxd)/c; % signal travel time
                if (delay > mindelay & delay < maxdelay)
                    delaysteps=round(delay/DeltaT);
                    gain=sqrt(txd*rxd); % gain factor due to 2D spread loss
                    for n = 1:NtMax % filter response
                        np = n + delaysteps;
                        if (np > NtMax) 
                            PF(n,i,j)=0;
                        else
                            PF(n,i,j) = P(np,i,j);
                            %PF(n,i,j) = gain*P(np,i,j); % optional use of
                            %gain function
                        end
                    end
                    PFS=PFS+PF(:,i,j);
                end
            end
        end
        I(xindex,yindex)=norm(PFS)^2;
    end
    jobdone = xindex / (2*nshalf) * 100;
    clc;
    disp('processing... ');
    disp_string = strcat(num2str(jobdone),'%');
    disp(disp_string)
end
% normalize I
[M,i] = max(I(:));
I = I/M;
I=I';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% plots %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
myfont = 'Times New Roman';
myfontsz = 14;
step = ds*1000; % step in mm
xmin = -nshalf*step;
xmax =  nshalf*step;
%[xmesh,ymesh] = meshgrid(xmin:step:xmax,xmin:step:xmax);
x=[xmin xmax];
y=[xmin xmax];
figure;
myfont = 'Times New Roman';
myfontsz = 14;
imagesc(x,y,I);
colormap bone;
%shading interp;
xlabel('{\it x} (mm)','FontName',myfont,'FontSize',myfontsz,'Interpreter','latex');% label x
ylabel('{\it y} (mm)','FontName',myfont,'FontSize',myfontsz,'Interpreter','latex');% label y
set(gca,'FontName',myfont,'FontSize',myfontsz);
% save image
save('image.mat','x','y','I');
