figure
Z = 1:100;  % 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
[mac, CEdge] = PhotonAttenuationQ(Z, E);
colormap(jet(128)) % use hi-res color palette
imagesc(log10(mac)); 
grid on;
axis xy; % put small numbers on y axis on the bottom
title('Photon Mass Attenuation Coefficients (in cm^2/g) and Compton Edges');
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'     ,     -1:4 ); % The image is a log10 of the MAC ...
set(cbar_axes,'YTickLabel',10.^(-1:4)); % ... so add proper labels
hold on
% Add Conpton Edges to the plot
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),'b','LineWidth',2); %Plot 3 L Compton edges
M=plot(ed(:, 5: 9),'g','LineWidth',2); %Plot 5 M Compton edges
N=plot(ed(:,10:14),'Y','LineWidth',2); %Plot first 5 N Compton edges
legend([K(1),L(1),M(1),N(1)], {'K','L','M','N'}, 'Location', 'NorthWest');