%% Plot mean free path of photons in different media
% Plot as function of energy and atomic number of media/element
figure
Z = 1:99;  % elements with Z in 1-100 range - number of columns
nr = 500;   % number of rows to use in the plot
E = logspace(log10(0.001), log10(20), 500);  % define energy grid
mfp = PhotonAttenuation(Z, E, 'mean free path');
colormap(jet(128)) % use hi-res color palette
imagesc(log10(mfp)); 
grid on;
axis xy; % put small numbers on y axis on the bottom
title('Mean free path of photons in different media (in cm)');
xlabel('Atomic number of elements');
ylabel('Energy in MeV');
% Add X-Axis 
EPos = [6 16 26 35 47 55 65 74 82 94];   % define array to store label location
ELab = { '6-C','16-S','26-Fe','35-Br','47-Ag','55-Cs','65-Tb','74-W','82-Pb','94-Pu'}; %Define Energy labels for y-axis
set(gca,'XTick'     ,EPos);
set(gca,'XTickLabel',ELab);
% Add Y-Axis
ELab = [0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]; %Define Energy labels for y-axis
EPos = size(ELab);   % define array to store label location
for i=1:length(ELab), [tmp EPos(i)]=min(abs(E-ELab(i))); end
set(gca,'YTick'     ,EPos);
set(gca,'YTickLabel',ELab);
% add Colorbar
cbar_axes = colorbar;
set(cbar_axes,'YTick'     ,     -5:5 ); % The image is a log10 of the MAC ...
set(cbar_axes,'YTickLabel',10.^(-5:5)); % ... so add proper labels
hold on
% Add Compton Edges to the plot
[tmp, CEdge] = PhotonAttenuationQ(Z, E);
ed = accumarray([CEdge(:,1),CEdge(:,2)],CEdge(:,3)); % get per element energies of 14 compton edges 
ed = 500*(log(ed')-log(0.001))/(log(20)-log(0.001)); % convert energy to row numbers of the image
K=plot(ed(:,1)    ,'m','LineWidth',3); %Plot K Compton edge
L=plot(ed(:, 2: 4),'r','LineWidth',2); %Plot 3 L Compton edges
M=plot(ed(:, 5: 9),'g','LineWidth',2); %Plot 5 M Compton edges
N=plot(ed(:,10:end),'y','LineWidth',2); %Plot first 5 N Compton edges
legend([K(1),L(1),M(1),N(1)], {'K','L','M','N'}, 'Location', 'NorthEast');