% [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.