%%% photobleachDemo.m
% Demonstrate photobleaching of a chromophore.
% Plots the fluence rate (phi) versus depth (z)
% in a non-scattering gel containing a chromophore
% at an initial concentration that yields an
% absorption coefficient (mua0).
%    Problem is 1-D model for chromophore 
% in a non-scattering gel during exposure to
% a photobleaching irradiance (1 W/cm^2).
%
% The absorption within each voxel drops with 
% each incremental time (dt) of exposure:
%   mua(z,t) = mua(z,t)*exp(-phi(z,t)*dt/Hth)
% where
%   mua = absorption coefficient [cm^-1] of choromophore
%   phi = fluence rate [W/cm^2]
%   t   = time of exposure [s]
%   z   = depth [cm]
%   Hth = 1/e photobleaching dose
%
% Steven L. Jacques,
% Student Workshop, Oulu, Finland, June 5, 2010
%%%
clear
home

%% axial position within gel
dz   = 0.1; % [cm]
zmax = 10;
z    = [dz:dz:zmax]';
Nz   = length(z);
%%

%% initial optics
% mua0 = epsilon*Conc*ln(10)
%   epsilon [cm^-1 M^-1], extinction coefficient (literature)
%   Conc [M] = [mole/liter], concentration 
%   ln(10), since 10^(epsilon*C*L) = exp(-mua*L)
mua0 = 1; % [cm^-1], initial absorption coefficient 
mua  = zeros(Nz,1) + mua0;  %<--- set mua for each voxel ----

OD   = cumsum(mua.*dz); % <------ cumulative optical depth --
phi  = exp(-OD);        % <------ fluence rate phi(z) -------
%%

%% Draw initial fluence rate vs depth
figure(2);clf
set(figure(2),'color','w')
plot(z,phi,'r-','linewidth',2)
hold on
set(gca,'fontsize',18)
xlabel('z [cm]')
ylabel('\phi [W/cm^2]')
%%

%% Draw phi(z) during photobleaching.
dt = 1; % incremental timestep [s]
Hth = 50; % 1/e photobleaching dose [J/cm^2]
for it=1:600 % timesteps [s]
    mua = mua.*exp(-phi*dt/Hth);
    OD  = cumsum(mua.*dz);
    phi = exp(-OD);
    
    % plot phi(z) at time t(it)
    if mod(it,60)==0
        plot(z,phi,'b-','linewidth',2)
    elseif mod(it,10)==0
        plot(z,phi,'b-','linewidth',1)
        drawnow
    end
end
%%

disp('done')
