% draw W map
% Create W from C for +/- r vs z
% Plot in Figure 2

maxC             = max(max(C));
temp             = C(2:Nz-1, 2:Nr-1)*exp(-mua*c*time)/Co;
Nw = 2*(Nr-2);
W(:,Nw/2+1:Nw)   = temp;
W(:,Nw/2:-1:1)   = temp;
rw = (-rmax:dr:rmax);
zw = (0:dz:zmax)';
[Nwz Nwr] = size(W);

% select scan lines
h = 1;       % @z
v = Nwr/2+1; % @r

%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
figure(1);clf
set(figure(1), 'position', [460 100 500 600])

imagesc(rw,zw,log10(W),[-12 0])
hold on
plot([-rmax rmax],[zw(h) zw(h)], 'w-') % horizontal line
plot([rw(v) rw(v)], [0 zmax],'w-') % vertical line
xlabel('r [cm]')
xlabel('z [cm]')
axis([-rmax rmax 0 zmax])
axis('equal')

summ = 0;
Vol  = 0;
for ir=v:Nwr
	rr = (ir - v + 0.5)*dr;
	for iz=1:Nwz
		summ = summ + W(iz,ir)*2*pi*rr*dr*dz;
		Vol  = Vol  + 2*pi*rr*dr*dz;
	end
end
title(sprintf('total U = %1.3f = %1.3f. Vol = %1.3f = %1.3f',summ,Uo,Vol, zmax*pi*rmax^2),'fontsize',18)

%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
figure(2);clf
set(figure(2), 'position', [30 100 400 700])

% C vs r
subplot(2,1,1)
semilogy(rw(1:Nwr)+dr/2,W(h,:), 'r-')
hold on
set(gca,'fontsize',24)
xlabel('r [cm]')
ylabel('C [J/cm^3]')
axis('square')
axis([-rmax rmax 1e-9 1])

if time<1000e-12
	title(sprintf('time = %5.0f [ps]',time*1e12),'fontsize',36)
else
	title(sprintf('time = %5.3f [ns]',time*1e9),'fontsize',36)
end

% C vs z
subplot(2,1,2)
semilogy(zw(1:38)+dz/2,W(:,v), 'r-')
hold on
set(gca,'fontsize',24)
xlabel('z [cm]')
ylabel('C [J/cm^3]')
axis('square')
axis([0 rmax 1e-9 1])

drawTRDiffusion

