
% CompMethods() 
% 
% This script compares various methods for the forward problem 
% using an infinite medium since the exact solution for a sphere is 
% only accurate for an infinite medium. 
% 
% Methods include: Exact, 1st Born, Nth Born, Full Born, Ext Born, 
% and Nth Ext Born 
% 
% Feel free to change the parameters to explore different optical 
% contrast in the object. The code is set for an absorption 
% perturbation. If you include scattering then don't forget to 
% switch the ds.Fwd.Method.ObjVec_musp flag to 1. 
% 
% Note that at present the Full Born and Ext Born methods do not 
% work for a scattering perturbation, but Nth Born and Nth Ext Born 
% do. 
% 
% Note that small differences between the exact solution and the 
% discretized solutions are expected because of small differences 
% in the volume of the object. 
% 

clear ds 
ds.Debug = 0; %Show Debug Information 

% 
% OBJECT POSITION AND RADIUS 
% 
xo = 0; 
yo = 0; 
zo = 3; 
ro = 0.4; 
Thickness = 6; 

% 
% SET THE FORWARD MODEL STRUCTURE 
% 
ds.Fwd.Lambda = [780]; 
ds.Fwd.idxRefr = [1.37]; 
ds.Fwd.Mu_s = [100]; 
ds.Fwd.g = [0.9]; 

ds.Fwd.Mu_sp = ds.Fwd.Mu_s .* (1 - ds.Fwd.g); 
ds.Fwd.v = 2.99e10 ./ ds.Fwd.idxRefr; 
ds.Fwd.Mu_a = [0.05]; 

x = [0:0.4:6]; 
ds.Fwd.Src = SetOptode( 0, 0, 0, 1e-3 ); 
ds.Fwd.Det = SetOptode( x, 0, Thickness, 1e-3 ); 

ds.Fwd.ModFreq = 100; 

ds.Fwd.Boundary.Geometry = 'infinite'; 

% ONLY INCLUDE THE COMPUTATIONAL VOLUME AROUND THE SPHERE 
step = ro/5; 
ds.Fwd.CompVol.Type = 'uniform'; 
ds.Fwd.CompVol.ZStep = step; 
ds.Fwd.CompVol.Z=[zo-ro+step/2:step:zo+ro-step/2]; 
ds.Fwd.CompVol.XStep = step; 
ds.Fwd.CompVol.X=[xo-ro+step/2:step:xo+ro-step/2]; 
ds.Fwd.CompVol.YStep = step; 
ds.Fwd.CompVol.Y=[yo-ro+step/2:step:yo+ro-step/2]; 


ds.Fwd.MeasList = genMeasList('all', ds.Fwd); 
ds.Inv = ds.Fwd; 

% 
% SET THE OBJECT PARAMETERS 
% 
ds.Object = cell(1); 

ds.Object{1}.Type = 'Sphere'; 
ds.Object{1}.Pos = [xo yo zo]; 
ds.Object{1}.Radius = ro; 
ds.Object{1}.Mu_a = [0.3]; 
ds.Object{1}.Mu_s = [100]; 
ds.Object{1}.g = [0.9]; 
ds.Object{1}.Mu_sp = ds.Object{1}.Mu_s .* (1 - ds.Object{1}.g); 


% 
% Perform the Forward Calculations 
% 
ds.Fwd.Method.ObjVec_mua = 1; 
ds.Fwd.Method.ObjVec_musp = 0; 
ds.Fwd.Method.ObjVec_sd = 0; 

ds.Fwd.Method.Type = 'BornN'; 
ds.Fwd.Method.Born_Order = 10; 
ds.Fwd = genFwdMat(ds.Fwd); 
ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); 
rbn = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... 
ds.Fwd.P.PhiInc, 1 ); 

ds.Fwd.Method.Type = 'FullBorn'; 
ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); 
rbf = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... 
ds.Fwd.P.PhiInc, 1 ); 

ds.Fwd.Method.Type = 'Born'; 
ds.Fwd = genFwdMat(ds.Fwd); 
ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); 
rb1 = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... 
ds.Fwd.P.PhiInc, 1 ); 

ds.Fwd.Method.Type = 'ExtBorn'; 
ds.Fwd.Method.ExtBorn_Radius = 0.9; 
ds.Fwd = genFwdMat(ds.Fwd); 
ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); 
reba = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... 
ds.Fwd.P.PhiInc, 1 ); 

ds.Fwd.Method.Type = 'ExtBornN'; 
ds.Fwd.Method.ExtBorn_Radius = 0.9; 
ds.Fwd = genFwdMat(ds.Fwd); 
ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); 
rebna = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... 
ds.Fwd.P.PhiInc, 1 ); 

ds.Fwd.Method.Type = 'Spherical'; 
ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); 
re = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... 
ds.Fwd.P.PhiInc, 1 ); 

figure 
plot(x,abs(re),x,abs(rb1),x,abs(rbn),x,abs(reba),x,abs(rbf),x,abs(rebna)) 
legend('Exact','Born1','BornN','ExtBorn','FullBorn','ExtBornN') 
title('Amplitude Perturbation') 
xlabel('X Position of Detector (cm)') 
ylabel('Relative change') 

if ds.Fwd.ModFreq > 0 
figure 
plot(x,angle(re)*180/pi, ... 
x,angle(rb1)*180/pi, ... 
x,angle(rbn)*180/pi, ... 
x,angle(reba)*180/pi, ... 
x,angle(rbf)*180/pi, ... 
x,angle(rebna)*180/pi ) 
legend('Exact','Born1','BornN','ExtBorn','FullBorn','ExtBornN') 
title('Phase Perturbation') 
xlabel('X Position of Detector (cm)') 
ylabel('Phase change (degrees)') 

end

