% [sig] = slab(tmin, tmax, dtgate, R, mua, musp, dz {, n {, err}})
%
% Computes the time-resolved detected surface flux in multi-layered slab
% due to injection of a point source, in the diffusion approximation.
% The algorithm is capable of arbitrary absorption and scattering
% coefficient functions of depth, but the interface only handles
% multi-layered slabs.
%
% tmin, tmax, dtgate - give signal time-range and sample (gate) period
% R - S-D distance on surface (or row vector of such)
% mua, musp - row vectors, length defines number of layers N_l
% dz - row vector gives layer thicknesses (>0) in mm
% if length N_l then all layers given finite thickness
% if length N_l-1 then last layer is infinite thickness
% n - refractive index of whole system (default 1.4)
%
% err - array of error tolerances, model choices, flags.
% (see slab.c for defaults, beware C array offset of 1)
% err(1) = max est rel err due to box depth z_box
% err(2) = max est rel err due to
% z-lattice (h) at t_min
% (if <0, size sets h directly).
% err(3) = max est rel err due to timestepping
% (if <0, size sets dt [or dt/t for CN]).
% (only CN implemented).
% err(4) = max est rel err, induced by radial box.
% err(5) = max est rel err, trunc (k_max) induced.
% err(6) = calculation mode.
% -2: 2d, no gradient, top layer props.
% -1: full 3d (default)
% 0: 1d, u(z,t) (output test.1d.dat)
% (all R's are duplicates of data)
% 1,2...: only sum contrib due to that
% radial mode (output test.1d.dat)
%
% err(7) = verbosity, integer with single-bit flags:
% 0x1 : stdout basics
% 0x2 : stdout info about modes
% 0x4 : stdout signal data
% 0x8 : stdout mode contribs at each R
% 0x10 : stdout add_as_samples info
% 0x20 : stdout 1d lattice cell props (init_FD)
% 0x40 : stdout 1d lat info per mode (onedim_FD)
% 0x80 : stdout info about lowest mode only
% 0x100 : write 1d lattice datafile test.1dlat.dat
%
% err(8) = timestepping scheme.
% 0: const-dt expl-implicit (Backwards Euler)
% 1: variable-dt CN (w/ initial const-dt FE).
% 2: Forward Euler: fixed dt = dt_courant/2.
% 3: const-dt CN (no kill).
% 4: const-then-var CN (w/ initial CN kill).
% (4 is the default).
%
% err(9) = boundary conditions and detection model
% 0: dirichlet (ZBC), det gradient.
% 1: extrap BC, det value (Haskell94)(R3).
% 2: extrap BC, det gradient (Martelli99,R1).
% 3: extrap BC, det lin comb (Kienle97)(R2).
% (3 not yet implemented).
%
% If err not given as full length vector, remaining elements are set to
% defaults. (this retains backwards-compatibility). Defaults
% are set in slab.h
%
% sig(i,j) - array of signal samples (i=time, j=radius)
%
% All units are mm (length) and ps (time).
%
% Normalization of signal corresponds to launching 1 unit of mass,
% that is, a fluence value of 1 over a (1mm)^3 volume. This is the
% same as the analytic Greens function when expressed in mm and ps units.
% Also same as Kienle97, without the overall factor of c.
%
% The value of asking for multiple R measurement distances in a
% single call is that the computation time is no more than for a
% single R.
%
% Error sizes are discussed in the companion publication,
%
% ``A fast numerical method for time-resolved photon diffusion
% in general stratified turbid media'', Alex Barnett, preprint.
%
% Available at http://www.cims.nyu.edu/~barnett/pubs.html
%
% In brief, a reasonable choice is the values given in 2nd example below.
% err(2)=h and err(3)=beta are most crucial and have large effects on runtime.
% Calling with default settings (ie no err vector) gives of order
% 1% error.
%
% Returned signals are C1-differentiable in material input
% parameters for semi-infinite media. But, for finite thickness, h is tweaked
% to next smallest h such that the thickness is integer multiple of
% h (beware: this non-smoothness may not always be small).
%
% Example usage: a 3-material semi-infinite slab, with top two
% layers 7 mm and 8 mm, refractive index 1.4 everywhere, measuring
% at distance 0 mm and 30 mm, at 0.1 ns intervals from 0.1 to 5 ns:
% s = slab(100,5000,100,[0,30],[0.01,0.005,0.002],[1,2,1.5],[7,8],1.4);
%
% Notice that returned signals are instantaneous at each time; if
% you want signals integrated across timegates (eg for photon
% counting), you will have to approximate the integral yourself.
%
% Controlling error parameters explicitly, requesting diagnostic output:
% s = slab(100,5000,100,[0,30],[0.01,0.005,0.002],[1,2,1.5],[7,8],1.4,...
% [1e-6,-0.4,-0.05,1e-6,1e-12,-1,129,4,1]);
%
% Testing can be explored with the (undocumented) test.m, or the fig_*.m
% files used to generate figures in the paper.
%
% Author: Alex Barnett (Courant Institute, New York University).
% barnett at cims.nyu.edu
%
% Version history:
% v 0.0 (11/18/02)
% v 0.7 (12/3/02) multiple-Rs, 1d: forward-Euler finite-difference (slow!)
% v 1.0 (12/9/02) Variable-dt Crank-Nicholson, problems with bessel cancel
% v 1.1 (1/5/03) fixed-dt expl-implicit for now, typ 1 sec run time.
% v 1.2 (4/19/03) Bessel cancel fix->fast! Extrap BCs and det type options.
% v 1.3 (4/26/03) exact LAPACK solve for lowest 3d eigval, sets mu_shift.
% v 1.4 (5/19/03) new timestepping schemes. CN var as default.
% v 1.5 (8/1/03) Switch to CN const-var as default timestepping.
% v 1.6 (5/22/05) Debug segfault from static alloc to stack, changed to malloc
%
% Future plans: (also see publication)
% * nonuniform (fixed-in-time) z-grid based on error estimates. Faster!
% * output of full 2d slice in time: phi(R,z,t), for Matlab import.
% * adjoint differentiation to get derivative wrt input params.
% * different refractive index for each layer (low priority).
% * z_b (for EBC) a continuous function of top layer refractive
% index. (low priority).
%
% Also see: HOMOG.