%CavityModes.m sample program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                               %
%                       Cavity Modes                            %
%                                                               %
%        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)

f0 = .5e9; %monochromatic pulse frequency

N = round(5*Lt*f0) %number of points for the chosen time-window (Lt)

%Loss coefficient
R = 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);
        
        %Frequency response
        Fs = N/Lt; %sampling frequency
        T = 1/Fs;                     % sample length
        L = N;                  %Number of points
        tt = (0:L-1)*T;                % time...
        
        NFFT = 2^nextpow2(N); % Next power of 2 from length of y
        %Yx = fft((Sx),NFFT);
        %Yy = fft((Sy),NFFT);
        Yz = fft((Sz),NFFT);
        f = Fs/2*linspace(0,1,NFFT/2);
        
        % Three axis frequency response
        %FFTx(i,j,:) = Yx(1:NFFT/2);
        %FFTy(i,j,:) = Yy(1:NFFT/2);
        FFTz(i,j,:) = Yz(1:NFFT/2);
        freq = f;
        
    end
    disp([num2str(round(i/length(X)*100)),' %'])
end



figure(1)
for k=1:1:length(freq)
    imagesc(X,Y,abs(FFTz(:,:,k))')
    colorbar
    axis equal
    axis tight
    xlabel('x [m]')
    ylabel('y [m]')
    title(['Vertical component of the E-field, f=',num2str(round(freq(k)/1e6)),' MHz'])
    getframe;
    
    %     %Uncomment if you want to save the animation as an image sequence
    %     file=sprintf('CavityModes_%d.png',k)
    %     saveas(gcf,file)
end


