% lookmcsubfluor.m
% Plots the output of running gomcsub*
% compiled as cc -o gomcsub mcsubLIB.c mcsub.c
%   which uses mcsub.h.
%
% Step 1: Plot excitation output in mcOUT101.dat
%         Fig. 1: point spread fuction of escape,
%        	J(r) [1/cm^2], for excitation and emission.
%         Fig. 2: point spread function of excitation transport, 
%           Tex [1/cm^2].
%
% Step 2: Plot background fluorescence due to uniform fluorophore
%         Fig. 3: point spread function of emission transport, 
%           Tem [1/cm^2], for homogenous background fluorophore.
%         Fig. 4: Ratio of escape functions Jem(r)/Jex(r).
%         Fig. 5: J
% Step 3: Plot fluorescence from fluorescent heterogeneity
%
% Steven L. Jacques, Oct 24, 2007.
clear all
home

global Nfile J F S A E r z
global mua mus g n1 n2
global mcflag radius waist xs ys zs
global NR NZ dr dz Nphotons

P = 1; % W

p = pwd;
commandwindow

filename = 'mcOUT101.dat';
getmcsub(filename);
Jex = J;
u = (1:NR-1);


filename = 'mcOUT102.dat';
getmcsub(filename);
Jem = J;


figure(1);clf
semilogy(r(u), P*Jex(u),'b-','linewidth',2)
hold on
plot(r(u), P*Jem(u),'r-','linewidth',2)
plot(-r(u), P*Jex(u),'b-','linewidth',2)
plot(-r(u), P*Jem(u),'r-','linewidth',2)
set(gca,'fontsize',18)
xlabel('r [cm]')
ylabel('J(r) [1/cm^-^2]')
legend('excitation','emission')

figure(2);clf
y = Jem./Jex;
plot(r(u), y(u),'r-','linewidth',2)
hold on
semilogy(-r(u), y(u),'r-','linewidth',2)
set(gca,'fontsize',18)
xlabel('r [cm]')
ylabel('Jem/Jex [-]')
axis([-max(r) max(r) 0 max(y(u))*1.2])

figure(3);clf
FF = zeros(NZ,2*NR);
FF(:,NR+1:2*NR) = F;
FF(:,1:NR) = F(:,NR:-1:1);
rr = zeros(2*NR,1);
rr(1:NR) = -r(NR:-1:1);
rr(NR+1:2*NR) = r;
disp(sprintf('total J = %0.4f', sum(J.*r*2*pi*dr)))
imagesc(rr, z, log10(P*FF) )
colormap(makec2f)
colorbar
set(gca,'fontsize',18)
xlabel('r [cm]')
ylabel('z [cm]')
title('log_1_0(Tex [1/cm^2]')
axis equal image


figure(4);clf
FF = zeros(NZ,2*NR);
FF(:,NR+1:2*NR) = F;
FF(:,1:NR) = F(:,NR:-1:1);
rr = zeros(2*NR,1);
rr(1:NR) = -r(NR:-1:1);
rr(NR+1:2*NR) = r;
disp(sprintf('total J = %0.4f', sum(J.*r*2*pi*dr)))
imagesc(rr, z, log10(P*FF) )
colormap(makec2f)
colorbar
set(gca,'fontsize',18)
xlabel('r [cm]')
ylabel('z [cm]')
title('log_1_0(Tem [1/cm^2]')
axis equal image

%%
figure(5);clf
u = (1:NR-1); % omit overflow bin J(NR)
filename = 'mcOUT103.dat'; % positive x range of J
getmcsub(filename);
FF(:,NR+1:2*NR) = F(:,1:NR);
semilogy(r(u), P*J(u),'b-','linewidth',2)
hold on
%
filename = 'mcOUT104.dat'; % negative x range of J
getmcsub(filename);
FF(:,1:NR) = F(:,NR:-1:1);
plot(-r(u), P*J(u),'c--','linewidth',2)
set(gca,'fontsize',18)
xlabel('r [cm]')
ylabel('J(r) [1/cm^-^2]')

%%
figure(6);clf
imagesc(rr,z,log10(FF))
colorbar
colormap(makec2f)
set(gca,'fontsize',18)
xlabel('r [cm]')
ylabel('z [cm]')
title('log_1_0(T_e_xT_e_m [1/cm^4]')
axis equal image

ifigs=1:6;
