%Propagation.m sample program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                               %
%                     PULSE PROPAGATION                         %
%                                                               %
%        by E. Amador (emmanuel.amador@insa-rennes.fr)          %
%                         IETR/DGA                              %
%                                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

global Lt c R N POS

tic
c = 299792458;%
Lt = 1e-6; %Time-window length in seconds
nbre_elements = 1; %number of radiating elements

%loading the Position matrix from the image generator
filename = sprintf('%delem_%dns1s8.mat',nbre_elements,round(Lt/(1e-9)));
load(filename)

%Pulsed signal
tau = .1e-6; %length of the pulse in seconds
f0 = .5e9; %monochromatic pulse frequency

N = round(5*Lt*f0) %number of points for the chosen time-window (Lt)
t = 0:Lt/(N-1):Lt; %time scale
x = 0:1/((N-1)/Lt):tau;
s = sin(2*pi*f0*x); %pulsed signal

%Loss coefficient
R = 0.98; %lossy cavity 0.98 instead of 0.998

%Dimensions of the chamber
l = 8.7;
p = 3.7;
h = 2.9;


%Reception point rectangular coordinates
delta=.1; %reduce the value to increase the resolution and the computing time



X = 0:delta:l;
Y = 0:delta:p;
Z = 2; %the plane chosen is at z=2 m

for i=1:length(X)
    
    for j=1:length(Y)
        
        
        [Sx,Sy,Sz] = CIR(X(i),Y(j),Z);
        
        %Convolution of the CIRs with the chosen pulsed signal
        Six = conv(Sx,s);
        Siy = conv(Sy,s);
        Siz = conv(Sz,s);
        Signalfinalx(i,j,:) = Six(1:N);
        Signalfinaly(i,j,:) = Siy(1:N);
        Signalfinalz(i,j,:) = Siz(1:N);
    end
    disp([num2str(round(i/length(X)*100)),' %'])
end

%total power vs. time
PdB=10*log10(Signalfinalx.^2+Signalfinaly.^2+Signalfinalz.^2);
Pmax=max(max(max(PdB(:,:,:))));

figure(1)
for k=1:5:length(t)
    imagesc(X,Y,PdB(:,:,k)')
    colorbar
    caxis([-30 Pmax])
    axis equal
    axis tight
    xlabel('x [m]')
    ylabel('y [m]')
    title(['Total power in dBm, t=',num2str(round(t(k)/1e-9)),'ns'])
    getframe;
    
    %     %Uncomment if you want to save the animation as an image sequence
    %     file=sprintf('Propag_%d.png',k)
    %     saveas(gcf,file)
    
end

