%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Redbird-m Demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath(fullfile(pwd, '../matlab'));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Domain Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all

[cfg.node,cfg.face, cfg.elem ]= meshabox([0 0 0], [70, 60, 50], 15);    % Iso2Mesh command for meshing a box. meshabox(min, max, max elemnt volume)
cfg.seg = ones(size(cfg.elem,1),1);     % Single segmentation/"tissue" type

xi = [20:10:50];
cfg.srcpos = [xi(:),ones(length(xi(:)),1).*20,ones(length(xi(:)),1)];   % Declaring source location 
[xi,yi] = meshgrid(20:10:50,20:10:40);
cfg.detpos = [xi(:),yi(:),ones(length(xi(:)),1).*50];   % Detector location
cfg.srcdir = [0 0 1];   % Declare source direction

cfg.prop = [0 0 1 1;0.005 1 0 1.37];    % Optical properties: mua, mus, g, n. mua and mus are given in mm-1
cfg.omega = 0;  % Modulation frequnecy for frequency domain measurements. 0 = CW.

%=====Exercise : Try Changing Source/Detector Layouts or Properties=====%

cfg = rbmeshprep(cfg);

%% Forward Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[detphi,phi] = rbrunforward(cfg);

%=====This Should Provide Same Results as Following Section=====%

%% Forward Simulation v2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Build Left Hand Side (FEM Stiffness Matrix)
%%%%%
[Amat,deldotdel] = rbfemlhs(cfg); 

% Build Right Hand Side
%%%%%
[rhs,loc,bary ] = rbfemrhs(cfg);


% Solve for Solutions at All Free Nodes
%%%%%%
tic
fprintf(1,'Solving ... \n')
phi = rbfemsolve(Amat,rhs);
toc

% Extract Detector Readings
%%%%%%
detphi = rbfemgetdet(phi, cfg, loc, bary);

%% Figure 1 - Source Fluence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure,plotmesh([cfg.node log10(abs(phi(:,1)))],cfg.elem,'y=20','facecolor','interp','linestyle','none');
hold on
plotmesh([cfg.node log10(abs(phi(:,1)))],cfg.elem,'x=20','facecolor','interp','linestyle','none');
view(3);
title('Log(abs(fluence)) for Source 1');

%% Figure 2 - Detector Fluence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure,plotmesh([cfg.node log10(abs(phi(:,14)))],cfg.elem,'y=20','facecolor','interp','linestyle','none');
hold on
plotmesh([cfg.node log10(abs(phi(:,14)))],cfg.elem,'x=50','facecolor','interp','linestyle','none');
view(3);
title('Log(abs(fluence)) for Detector 10');

%% Figure 3 - Detector Measurements from Source 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure,imagesc(reshape(detphi(:,2),3,4));
title('Detector fluence measurements for Source 3');


%=====Check Other Measurements by Changing wWhich Column of phi/detphi You
%Look At=====%

%====Alternative, Plot Rows to See Source Inten by Detector=====%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Multiple Options for Building Jacobians
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sd,dist] = rbsdmap(cfg);   %Returns channels and channel distances for input cfg

%=====Run One of the Below Options for Building Jacobian=====%

%% Native MATLAB code
tic;fprintf('Native code timing... \n')
[Jmua,Jmua_elem,JD,JD_elem] = rbjac(sd,phi,cfg.deldotdel,cfg.elem,cfg.evol);
toc


%% .mex to build Jacobian
tic; fprintf('Mex timing... \n')
Jmua = rbfemmatrix(cfg,sd,phi);
toc


%% Nodal-adjoint method
tic; fprintf('Nodal-adjoint method... \n')
Jmua = rbjacmuafast(sd,phi,cfg.nvol);
toc

%=====Try Different Options to Check Speed and Results (next)=====%


%% Figure 4 - Sensitivity Profile
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure,plotmesh([cfg.node log10(abs(Jmua(3,:)'))],cfg.elem,'x=20','facecolor','interp','linestyle','none');
view([1 0 0])
title('Log(abs(Jmua)) : Source 1, Detector 3')

