% -----------------------------------
% Copyright: Andreas Walther, Lund University, 2016
% 
% Note, this code is a research project and is thus permantly under development.
% This means that it may contain bugs and parts of the code may not run optimally, 
% in particular it will only work as intended in versions of matlab later than 2016.
%
% Any user of this code should be aware of this, but should also cite the original
% when using parts of the code for their own work.
%
% -----------------------------------

function Signal_vs_depth
clear all
close all

% options
plotDeltaMu = 1; % if 1, uses manual logplot with semitransparent regions
screen = 1; % screen for plots, 1 for left screen of 2 for right, use 1 if you only have one screen
r_log = 0; %=0 for linear distance, =1 for logarithmic (currently not updated)
r_max_plot = 8; % cm for max r in plot, write 0 for using entire r vector
plot_paper = 1; % removes long title and resizes window

% mu coefficients
mu = struct('a',0); % initilize struct
mu.a = 0.2; % cm^-1, 
mu.s = 5; % cm^-1, 5.2@880 for 'muscle'
mu.a_kon = mu.a;  % contrast, mu.a=4-6 for fully oxygenated blood
mu.range = [0.5 50]; % percent of mu as contrast range
mu.a_sig = logspace(log10(mu.a_kon*(1+mu.range(1)/100)),log10(mu.a_kon*(1+mu.range(2)/100)),100);

% Parameters
params = struct('det_area',1); % initialize
params.r_sig = 0.45; % path length through signal spot (cm), spot size = sqrt(r_sig/mu.s), 0.45 -> size=0.3
params.att = 1e-6; % Filter background attenuation
params.SL_att = 1e-2; % additional attenuation from time gating with slow light effect
params.UOT_pos = 4; % detector position from sender (cm), used for background
params.PAT_pos = 0; % note: does not affect PAT due to time syncing, should this be 0?
params.det_bw = 2e6; % detection bandwidth. depth x freq = 30cm MHz, i.e. 10 cm depth = 3 MHz sound. BW = 0.7*freq = 2 MHz
params.det_area = 1; % detector area in cm^2
params.nbr_vox = 30; % number of voxels (volumes elements)
params.nbr_avg = 1; % set to higher than 1 for number of 250 ms segments to average

% Simulation depths
if(r_log)
    r = logspace(-1.5,1.2,500); % cm
else
    r = linspace(0.2,9,500); % note, some max depths give odd x-axis ticks
end

% Call calculations script and find CNR=1 depths
for i = 1:length(mu.a_sig)
    [r,UOT_out(i,:),UOT_bak_out(i,:),PAT_out(i,:)] = expected_photons_com(mu,params,i,r);
    [~,UOT_ind(i)] = find(UOT_out(i,:)<1,1);
    [~,UOT_bak_ind(i)] = find(UOT_bak_out(i,:)<1,1);
    [~,PAT_ind(i)] = find(PAT_out(i,:)<1,1);
    [~,r_ind] = find(r>2,1);
    kon_UOT(i) = UOT_out(i,r_ind);
    kon_UOT_bak(i) = UOT_bak_out(i,r_ind);
    kon_PAT(i) = PAT_out(i,r_ind);
end
depth_UOT = r(UOT_ind);
depth_UOT_bak = r(UOT_bak_ind);
depth_PAT = r(PAT_ind);

% Plots
% --- SNR vs depth
figure(1)
hold on
set(gcf,'Units','normalized')
if(plot_paper)
    set(gcf,'OuterPosition',[screen-1+0.05 0.1 0.3 0.45]) 
else
    set(gcf,'OuterPosition',[screen-1+0.05 0.1 0.45 0.7]) 
end
set(gca,'FontSize',14)

mid = round(length(mu.a_sig)/2);

if(plotDeltaMu==1)
    if(r_log)
        rplot = log10(r);
    else
        rplot = r;
    end
    plot(rplot,log10(UOT_out(1,:)),'b'); hold on %plot UOT
    plot(rplot,log10(UOT_out(end,:)),'b')
    fhb = fill([rplot  fliplr(rplot)],[log10(UOT_out(1,:)) fliplr(log10(UOT_out(end,:)))],'b','faceAlpha',0.4);
%     plot(rplot,log10(UOT_bak_out(1,:)),'g')
%     plot(rplot,log10(UOT_bak_out(end,:)),'g')
%     fhg = fill([rplot  fliplr(rplot)],[log10(UOT_bak_out(1,:)) fliplr(log10(UOT_bak_out(end,:)))],'g','faceAlpha',0.4);
    plot(rplot,log10(PAT_out(1,:)),'r') % plot PAT
    plot(rplot,log10(PAT_out(end,:)),'r')
    fhr = fill([rplot  fliplr(rplot)],[log10(PAT_out(1,:)) fliplr(log10(PAT_out(end,:)))],'r','faceAlpha',0.4);
    if(plot_paper==0) % title shows current values
        title({['Coefficients: \mu_s=' num2str(mu.s) ' cm^-^1, \mu_a=' num2str(mu.a) ' cm^-^1,']...
               ['              \Delta\mu_a_,_s_i_g=' num2str(mu.range(1),2) '%...' num2str(mu.range(2),2)...
               '% , r_d_e_t=' num2str(params.r_sig*10) ' mm, det. area: ' num2str(params.det_area) ' cm^2']... 
               ['UOT: Separation: ' num2str(params.UOT_pos) ' cm, Abs. filter: 10^' num2str(log10(1/params.att))...
               ', Slow Light filter: 10^' num2str(log10(1/params.SL_att))]...
               ['PAT: Signal to thermal NEP, Separation: ' num2str(params.PAT_pos) ' cm, line of ' num2str(params.nbr_vox) ' pts']})
    end
    if(r_max_plot>0) % cap x-axis for plotting
        xlim([min(rplot),r_max_plot])
    else
        xlim([min(rplot),max(rplot)])
    end
    ylim([log10(1e-1), log10(5e4)])
    if(r_log) %the same for now
        xscale = rplot(round(linspace(1,length(r),10))); % note, the fixed number of ticks gives some issues
    else
        xscale = rplot(round(linspace(1,length(r),10)));
    end
    yscale = log10(logspace(-1,7,9)); %values for y-scale, may need different number of ticks for some r vectors
    set(gca,'xtick',xscale)
    set(gca,'ytick',yscale)
    if(r_log) % convert whole vector of x ticks to strings at once
        set(gca,'xticklabel',cellfun(@num2str,num2cell(10.^xscale),num2cell(1*ones(1,length(xscale))),'uniformoutput',false))
    else
        set(gca,'xticklabel',cellfun(@num2str,num2cell(xscale),num2cell(1*ones(1,length(xscale))),'uniformoutput',false))
    end
    yscale_c = cell(1,length(yscale));
    for c = 1:length(yscale) % convert y ticks to strings, but no 10^ notation for 1 and 10
        if(yscale(c)>=-1 && yscale(c)<=1)
            yscale_c(c) = cellstr(num2str(10^yscale(c),2));
        else
            yscale_c(c) = cellstr([num2str(yscale(c),'10^%d')]); %'10^'
        end
    end
    set(gca,'yticklabel',yscale_c) %cellfun(@num2str,num2cell(10.^yscale),num2cell(2*ones(1,length(yscale)),'uniformoutput',false))
    plot([min(rplot) max(rplot)],[0 0], 'k--') % line when CNR is 1
    leg = legend([fhb fhr],'UOT','PAT','location','northeast');
    % set(leg,'fontsize',12)
else % older code
    loglog(r,UOT_out(mid,:),'b'); hold on
    loglog(r,UOT_bak_out(mid,:),'g')
    loglog(r,PAT_out(mid,:),'r')
    
    text(5e-2,3e5,'UOT SNR','color','b')
    text(5e-2,6e6,'UOT SBR','color','g')
    text(5e-2,3e4,'PAT SNR','color','r','rotation',-30)
    plot([min(r) max(r)],[1 1], 'k--')
    xlim([min(r),max(r)])
    ylim([1e-1, 1e7])
    title({['Coefficients: \mu_s=10 cm^-^1, \mu_a=0.1 cm^-^1,']... % note, still hardcoded values...
           ['              \Delta\mu_a_,_s_i_g=' num2str(mu.a_sig(mid),3) ' cm^-^1 , US_f_o_c_u_s=3 mm']... 
           ['UOT: Detectors 4 cm apart, Abs. filter: 10^6, Slow Light filter: 10^2']...
           ['PAT: Signal to 0.1 mbar background pressure, 1 cm^2 det. area']})
end
xlabel('Imaging depth (cm)')
ylabel('Contrast-to-noise ratio (CNR)')


% fig 2 --- mu.a vs depth
% O2Sat = 100 - 270*dmu
% O2res = 2.7*dmu   % @880 nm 
% O2res = 0.19*dmu  % @700 nm
% Hbres =  1.0*dmu @ 880 and blood vs no blood
konv = 2.7; % conversion factor delta_mu -> blood O2%
dmu = (mu.a_sig-mu.a_kon)/mu.a_kon; % relative mu difference
figure(2) 
axl = semilogy(depth_UOT,dmu,'b','linewidth',1.5); % plot UOT depth
axl.Parent.YTick = [0.01 0.03 0.1 0.3];
axl.Parent.YTickLabel = {'1%' '3%' '10%' '30%'};
set(gcf,'Units','normalized')
set(gcf,'OuterPosition',[screen-1+0.5 0.5 0.3 0.45]) % [x y dx dy] (from bottom)
set(gca,'FontSize',14)
yyaxis right 
ax = semilogy(depth_UOT,konv*dmu,'b'); % establish right axis as O2%
ax.Parent.YTick = [0.01 0.03 0.1 0.3];
ax.Parent.YTickLabel = {'1%' '3%' '10%' '30%'};
ax.Parent.YColor = [0 0 0];
ylabel('Resolution of oxygen saturation')
yyaxis left
hold on
semilogy(depth_PAT,dmu,'r','linewidth',1.5) % plot PAT depth, since same scale as UOT, right axis is auto valid
% semilogy([min([depth_UOT depth_PAT]) max([depth_UOT depth_PAT])],[0.1/konv 0.1/konv],'k--') % min([depth_UOT depth_PAT]) max([depth_UOT depth_PAT]), [1,7]
% text(max([depth_UOT depth_PAT])/3,0.034,'Required to detect heart ischemia') % 0.034
% text(1.5,0.042,'10% reduction in O_2 saturation')

ylabel('\Delta\mu_a/\mu_a ')
xlabel('Imaging depth (cm)')
xax = gca;
xax.XTick = [1:8]; % only use integer cm's
xax.XTickLabel = [cellfun(@num2str,num2cell(1:8),'uniformoutput',false)];
title('Depth where CNR=1')
axis tight; % required to match the left and right axes
xlim([min(rplot),r_max_plot-1]) % if max is too low, the right/left axis becomes mismatched
legend('UOT','PAT','location','northwest')

end




%% ------------------------------------------

function [r,UOT_out,UOT_bak_out,PAT_out] = expected_photons_com(mu,params,loop,r)

% clear all
% close all

% Options
PATint = 3; % at the pulsed energy limit, 10 Hz can be used, giving 3 shots in 250 ms.
UOTint = 2.5*10^4; % a 2W laser gives 2uJ/1us_pulse. Limit of 200mW gives 100kHz rate and 2.5*10^4 pulses in 250 ms.
nbr_vox = params.nbr_vox; % number of voxels. Reduces the nbr photons for UOT but not PAT, since optical here is all at once.

% Parameters
% medically allowed pulse: 1.1e4*dur^(1/4)*1e-4/(6e-34*3e8/800e-9); 20 mJ/cm2 for a pulse
% continuous: 200 mW/cm2 giving a rep rate of only 10 Hz if max peak energy is used (Yao2014)
P = 10^17; % number of photons in a pulse /cm^2
mu_a = mu.a; % cm^-1 %
mu_s = mu.s; % cm^-1 %
mu_a_sig = mu.a_sig(loop);
mu_a_kon = mu.a_kon;
r_sig = params.r_sig; % path length through signal spot (cm) (note this is actual path length, larger than size b/c of scattering)
att = params.att; % Filter background attenuation
SL_att = params.SL_att; %1e-2; % additional attenuation from time gating with slow light effect
UOT_pos = params.UOT_pos; % detector position from sender (cm), used for background
PAT_pos = params.PAT_pos; % 
det_bw = params.det_bw; % detection bandwidth. depth x freq = 30cm MHz, i.e. 10 cm depth = 3 MHz sound. BW = 0.7*freq = 2 MHz
det_area = params.det_area; % detector area in cm^2
det_eff = 0.5; % detector efficiency
mu_eff = sqrt(3*mu_a.*(mu_a+mu_s)); % note mu_s here is actually the reduced coefficient mu_s' = (1-g)*mu_s


% Calculations

% UOT Signal
r1 = r - 1/mu_s; % boundary effects, distance to positive source
r2 = r + 5/mu_s; % distance to negative source, using n_tissue = 1.4 -> z_b = 2/mu_s, see Haskell1994
r_UOT1 = sqrt(r1.^2 + (UOT_pos/2)^2); % actual travelled light distance with fixed detectors
r_UOT2 = sqrt(r2.^2 + (UOT_pos/2)^2); % actual travelled light distance with fixed detectors
phi_uot = P/UOTint/nbr_vox*mu_eff.^2./(4*pi*mu_a).*(1./r_UOT1.*exp(-mu_eff*r_UOT1)-1./r_UOT2.*exp(-mu_eff*r_UOT2)); % photons /cm^2 /pulse reaching the sound wave
phi_uot = phi_uot*0.1; % sound diffraction efficiency
phi_uot_sig = phi_uot*exp(-mu_a_sig*r_sig); % fluence loss due to Lambert-Beer absorption in the small signal region
phi_uot_dif = phi_uot*exp(-mu_a_kon*r_sig); % contrast difference
phi_uot_sig = phi_uot_sig*mu_eff.^2./(4*pi*mu_a).*(1./r_UOT1.*exp(-mu_eff*r_UOT1)-1./r_UOT2.*exp(-mu_eff*r_UOT2)); % light from sound to other side, using same boundary conditions as at entry
phi_uot_dif = phi_uot_dif*mu_eff.^2./(4*pi*mu_a).*(1./r_UOT1.*exp(-mu_eff*r_UOT1)-1./r_UOT2.*exp(-mu_eff*r_UOT2)); % light from sound to other side
P_end_sig = phi_uot_sig*det_area*det_eff; % total photons collected in 1 cm^2 in a 1 us pulse * detector efficiency;
P_end_dif = phi_uot_dif*det_area*det_eff; % total photons collected in 1 cm^2 in a 1 us pulse * detector efficiency;

% UOT Background
r_bak = UOT_pos*ones(1,length(r));
r_bak1 = sqrt(r_bak.^2 - (1/mu_s)^2); % boundary effects
r_bak2 = sqrt(r_bak.^2 + (5/mu_s)^2);
phi_bak = P/UOTint/nbr_vox*mu_eff.^2./(4*pi*mu_a).*(1./r_bak1.*exp(-mu_eff*r_bak1)-1./r_bak2.*exp(-mu_eff*r_bak2)); % photons /cm^2 going directly
phi_bak = phi_bak*att; % attenuation from filter
P_bak = phi_bak*det_area; % total background in 1 cm^2 in one pulse;
P_bak_SL = P_bak*SL_att; % attenuation from slow light time gating

% PAT
r_PAT_det = sqrt(r.^2 + (PAT_pos/2)^2); % distance to PAT detector
r_PAT = sqrt(r.^2 + (0.1)^2); %to avoid an artifact, and only depths longer than 1 mm are used anyway
r_PAT1 = r_PAT - 1/mu_s; % boundary effects
r_PAT2 = r_PAT + 5/mu_s;
phi_pat = P*mu_eff.^2./(4*pi*mu_a).*(1./r_PAT1.*exp(-mu_eff*r_PAT1)-1./r_PAT2.*exp(-mu_eff*r_PAT2)); % photons /cm^2 /pulse reaching the voxel
E_pat_sig = phi_pat *2.5e-19 * (1-exp(-mu_a_sig*r_sig)); % absorbed energy per unit area with mu_a_sig abs coeff
E_pat_dif = phi_pat *2.5e-19 * (1-exp(-mu_a_kon*r_sig)); % contrast difference
dT_sig = E_pat_sig/(1*r_sig*4.18); % dT = E/A * 1/(rho*r*c_v); note: V=A*r, J/cm^2 -> K
dT_dif = E_pat_dif/(1*r_sig*4.18); % dT = E/A * 1/(rho*r*c_v), [c_v] = J/(g K)
th_exp = 2e9*2e-4; % thermal expansion, isothermal bulk modulus K_T [N/m^2] * volume exp alpha [1/K] = Pa/K
P_pat_sig = dT_sig*th_exp*1/4; % dT*(Pa/K)*efficiency = Pressure
P_pat_dif = dT_dif*th_exp*1/4; % dT*(Pa/K)*efficiency = Pressure
dP_pat_sig = P_pat_sig.*r_sig./r_PAT_det; % 1/r linear scaling ; old: detector area at distance r_PAT: sqrt(det_area./(4*pi*r_PAT.^2))
dP_pat_dif = P_pat_dif.*r_sig./r_PAT_det; 
% Gr = 2e-4*2e9/(4180*1000); %  alpha[1/K] * K_T[N/m^2] / (C_v[J/(kg*K)] * rho[kg/m^3)], dimensionless
% P_grun = phi_pat *2.5e-19 *300e-6*2e9/(4.18*1e6) *1/4*1e2 *1e4*1e-2; %phi * J/photon *Gr *eta *mu(m^-1) *m^2/cm^2*mbar/Pa
% dP_grun(j,:) = P_grun*0.1; % 10% detection efficiency; this form is taken from Wang2015, Optica (OSA), same as other form

% PAT averaging
if(PATint>0)
    dP_pat_sig = sqrt(PATint)*dP_pat_sig;
    dP_pat_dif = sqrt(PATint)*dP_pat_dif;
end

% PAT background
kB = 1.38e-23; % J/K = Nm/K
T = 300; % K
F = 2; % preamp noise factor
eta = 0.5; % efficiency, for resonant detectors
Z = 1.5e6; % acoustic impendence for water [Pa s /m = Ns/m^3] (note: wrong unit in winkler2013!, should not be additional /m^2)
A = det_area*(1/100)^2; % m^2
NEP = sqrt(kB*T*(1+F/eta)*Z/A); % see Winkler2013
dP_bak = NEP*sqrt(det_bw)*ones(1,length(r)); % Pa background, ~1e-4

% Output
UOT_bak_out = abs(P_end_sig-P_end_dif)./P_bak_SL;
% UOT signal compared with shot noise from total light, increased if any averaging
UOT_out = sqrt(params.nbr_avg) * sqrt(UOTint) * abs(P_end_sig-P_end_dif)./sqrt(P_end_dif+P_bak_SL);
PAT_out = sqrt(params.nbr_avg) * abs(dP_pat_sig-dP_pat_dif)./dP_bak;

end