%An example MATLAB script for loading the focus fields and plotting them.
%Important point to note in this cript is the way the data are loaded and indexed in MATLAB.
%The same procedure applies for other matrices generated by focus3d and dipole_radiation objects.
%Usage: plot_focus

clear all;

I = sqrt(-1);

%the directory where the data files are
inputdir = ['example_results/successive/focus2_lambda_0.8/data_files/'];

%the directory where the image files will be stored
outputdir = [inputdir '../images/'];
if(exist(outputdir) ~= 7) 			%create output directory if does not exist
        mkdir(outputdir);
end

x = load([inputdir 'x.dat']);
y = load([inputdir 'y.dat']);
z = load([inputdir 'z.dat']);

Nx = length(x);
Ny = length(y);
Nz = length(z);

tempr = load([inputdir 'E_focus_x_real.dat']);
tempi = load([inputdir 'E_focus_x_imag.dat']);
E_focus_x1 = tempr + I*tempi;

tempr = load([inputdir 'E_focus_y_real.dat']);
tempi = load([inputdir 'E_focus_y_imag.dat']);
E_focus_y1 = tempr + I*tempi;

tempr = load([inputdir 'E_focus_z_real.dat']);
tempi = load([inputdir 'E_focus_z_imag.dat']);
E_focus_z1 = tempr + I*tempi;

% Note the way in which the matrices are reshaped. First index corresponds to z and the last index corresponds to x. This is beacuse the way in which the data files are stored by the focus3d object. The same format holds for the files stored by dipole_radiation objects
E_focus_x1 = reshape(E_focus_x1,[Nz,Ny,Nx]);
E_focus_y1 = reshape(E_focus_y1,[Nz,Ny,Nx]);
E_focus_z1 = reshape(E_focus_z1,[Nz,Ny,Nx]);

% the rest of the below code is plotting the figures and storing the results in eps images.

fig_num = 1;

figure(fig_num);
imagesc(x,z,squeeze(abs(E_focus_x1(:,(Ny+1)/2,:)))); axis image; colorbar;
title (['y = ' num2str(y((Ny+1)/2),'%6.2f')]);
xlabel('x (µm)'); ylabel('z (µm)');
saveas(gcf,[outputdir 'xz_amplitude_xcomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num);
imagesc(y,z,squeeze(abs(E_focus_x1(:,:,(Nx+1)/2)))); axis image; colorbar;
title (['x = ' num2str(x((Nx+1)/2),'%6.2f')]);
xlabel('y (µm)'); ylabel('z (µm)');
saveas(gcf,[outputdir 'yz_amplitude_xcomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num);
imagesc(x,z,squeeze(abs(E_focus_y1(:,(Ny+1)/2,:)))); axis image; colorbar;
title (['y = ' num2str(y((Ny+1)/2),'%6.2f')]);
xlabel('x (µm)'); ylabel('z (µm)');
saveas(gcf,[outputdir 'xz_amplitude_ycomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num);
imagesc(y,z,squeeze(abs(E_focus_y1(:,:,(Nx+1)/2)))); axis image; colorbar;
title (['x = ' num2str(x((Nx+1)/2),'%6.2f')]);
xlabel('y (µm)'); ylabel('z (µm)');
saveas(gcf,[outputdir 'yz_amplitude_ycomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num);
imagesc(x,z,squeeze(abs(E_focus_z1(:,(Ny+1)/2,:)))); axis image; colorbar;
title (['y = ' num2str(y((Ny+1)/2),'%6.2f')]);
xlabel('x (µm)'); ylabel('z (µm)');
saveas(gcf,[outputdir 'xz_amplitude_zcomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num);
imagesc(y,z,squeeze(abs(E_focus_z1(:,:,(Nx+1)/2)))); axis image; colorbar;
title (['x = ' num2str(x((Nx+1)/2),'%6.2f')]);
xlabel('y (µm)'); ylabel('z (µm)');
saveas(gcf,[outputdir 'yz_amplitude_zcomp.eps'],'epsc');
fig_num = fig_num+1;

for k = (Nz+1)/2%1:Nz
    figure(fig_num);
    imagesc(y,x,squeeze(abs(E_focus_x1(k,:,:)))); axis image; colorbar;
    title (['z = ' num2str(z(k),'%6.2f')]);
    xlabel('x (µm)'); ylabel('y (µm)');
    saveas(gcf,[outputdir 'xy_amplitude_xcomp.eps'],'epsc');
    fig_num = fig_num+1;

    figure(fig_num);
    imagesc(y,x,squeeze(angle(E_focus_x1(k,:,:)))/pi); axis image; set(gca,'clim',[-1,1]); colorbar;
    title (['z = ' num2str(z(k),'%6.2f')]);
    xlabel('x (µm)'); ylabel('y (µm)');
    saveas(gcf,[outputdir 'xy_phase_xcomp.eps'],'epsc');
    fig_num = fig_num+1;

    figure(fig_num);
    imagesc(y,x,squeeze(abs(E_focus_y1(k,:,:)))); axis image; colorbar;
    title (['z = ' num2str(z(k),'%6.2f')]);
    xlabel('x (µm)'); ylabel('y (µm)');
    saveas(gcf,[outputdir 'xy_amplitude_ycomp.eps'],'epsc');
    fig_num = fig_num+1;

    figure(fig_num);
    imagesc(y,x,squeeze(angle(E_focus_y1(k,:,:)))/pi); axis image; set(gca,'clim',[-1,1]); colorbar;
    title (['z = ' num2str(z(k),'%6.2f')]);
    xlabel('x (µm)'); ylabel('y (µm)');
    saveas(gcf,[outputdir 'xy_phase_ycomp.eps'],'epsc');
    fig_num = fig_num+1;

    figure(fig_num);
    imagesc(y,x,squeeze(abs(E_focus_z1(k,:,:)))); axis image; colorbar;
    title (['z = ' num2str(z(k),'%6.2f')]);
    xlabel('x (µm)'); ylabel('y (µm)');
    saveas(gcf,[outputdir 'xy_amplitude_zcomp.eps'],'epsc');
    fig_num = fig_num+1;

    figure(fig_num);
    imagesc(y,x,squeeze(angle(E_focus_z1(k,:,:)))/pi); axis image; set(gca,'clim',[-1,1]); colorbar;
    title (['z = ' num2str(z(k),'%6.2f')]);
    xlabel('x (µm)'); ylabel('y (µm)');
    saveas(gcf,[outputdir 'xy_phase_zcomp.eps'],'epsc');
    fig_num = fig_num+1;
end

figure(fig_num); plot(x,squeeze(abs(E_focus_x1((Nz+1)/2,(Ny+1)/2,:))));
xlabel('x (µm)'); ylabel('Amplitude');
saveas(gcf,[outputdir 'x_profile_xcomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num); plot(y,squeeze(abs(E_focus_x1((Nz+1)/2,:,(Nx+1)/2))));
xlabel('y (µm)'); ylabel('Amplitude');
saveas(gcf,[outputdir 'y_profile_xcomp.eps'],'epsc');
fig_num = fig_num+1;

figure(fig_num); plot(z,squeeze(abs(E_focus_x1(:,(Ny+1)/2,(Nx+1)/2))));
xlabel('z (µm)'); ylabel('Amplitude');
saveas(gcf,[outputdir 'z_profile_xcomp.eps'],'epsc');
fig_num = fig_num+1;