% 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