# doc-cache created by Octave 10.2.0
# name: cache
# type: cell
# rows: 3
# columns: 47
# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
besselhprime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 436

 hp=besselhprime(n,k,z)

 Hankel function first order derivative 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Hankel function
     k: kind of the Hankel function
     z: input variable

 output:
     hn: Hankel function first order derivative 

 example:
     hn=besselhprime(0,1,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 24

 hp=besselhprime(n,k,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
besseljprime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 435

 jp=besseljprime(n,z)

 Bessel function (Bessel first kind) first order derivative 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Hankel function
     z: input variable

 output:
     jp: Bessel function (Bessel first kind) first order derivative

 example:
     jp=besseljprime(0,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 22

 jp=besseljprime(n,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
besselyprime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 440

 yp=besselyprime(n,z)

 Neumann function (Bessel second kind) first order derivative 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Hankel function
     z: input variable

 output:
     yp: Neumann function (Bessel second kind) first order derivative 

 example:
     yp=besselyprime(0,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 22

 yp=besselyprime(n,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 14
cart2sphorigin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 589

 [T,P,R]=cart2sphorigin(xi,yi,zi,x0,y0,z0)

 Converting Cartesian coordinates to spherical with a reset of origin

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     xi,yi,zi: list of Cartesian coordinates
     x0,y0,z0: new origin in Cartesian coordinates

 output:
     T: theta of the output spherical coordinates
     P: phi of the output spherical coordinates
     R: R of the output spherical coordinates

 example:
     [T,P,R]=cart2sphorigin(1:5,1:5,1:5,2,2,2);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 43

 [T,P,R]=cart2sphorigin(xi,yi,zi,x0,y0,z0)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
genT5mesh


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 467
 tess_lat: simplicial tessellation of a rectangular lattice
 usage: [tessellation,vertices]=genT5mesh(v1,v2,v3,...)

 arguments: input
 v1,v2,v3,... - numeric vectors defining the lattice in
 each dimension.
 Each vector must be of length >= 1

 arguments: (output)
 vertices - factorial lattice created from (v1,v2,v3,...)
 Each row of this array is one node in the lattice
 tess - integer array defining simplexes as references to
 rows of "vertices".



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 tess_lat: simplicial tessellation of a rectangular lattice
 usage: [tessell...



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
genT6mesh


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 544
 tess_lat: simplicial tessellation of a rectangular lattice
 usage: [tessellation,vertices]=genT6mesh(v1,v2,v3,...)

 URL: http://www.mathkb.com/Uwe/Forum.aspx/matlab/50484/Constant-surfaces

 arguments: input
 v1,v2,v3,... - numeric vectors defining the lattice in
 each dimension.
 Each vector must be of length >= 1

 arguments: (output)
 vertices - factorial lattice created from (v1,v2,v3,...)
 Each row of this array is one node in the lattice
 tess - integer array defining simplexes as references to
 rows of "vertices".



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 tess_lat: simplicial tessellation of a rectangular lattice
 usage: [tessell...



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
generate_g1


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 2043

   [tau,g1]=generate_g1(fhist,tau, disp_model, DV, lambda, format)

   Compute simulated electric-field auto-correlation function using
   simulated photon pathlengths and scattering momentum transfer

   author: Stefan Carp (carp <at> nmr.mgh.harvard.edu)

   input:
       fhist:      the file name of the output .mch file 
       tau:        correlation times at which to compute g1 
                   (default: 1e-7 to 1e-1 seconds, log equidistant)
       disp_model: displacement model ('brownian', 'random_flow', <custom>)
                   (default: brownian, see further explanation below)
       disp_var:   value of displacement variable using mm as unit of
                   length and s as unit of time
                   (default: 1e-7 mm^2/s, see further explanation below)
       lambda:     wavelenght of light used in nm
                   (default: 785)
       format:     the format used to save the .mch file 
                   (default: 'float')

   output:

       tau:        correlation times at which g1 was computed provided for
                   convenience (copied from input if set, otherwise 
                   outputs default)
       g1:         field auto-correlation curves, one for each detector

   The displacement model indicates the formula used to compute the root
   mean square displacement of scattering particles during a given delay
   
   brownian:       RMS= 6 * DV * tau; 
                   DV(displacement variable)=Db (brownian diffusion coeff)
   random_flow:    RMS= DV^2 * tau^2; 
                   DV = V (first moment of velocity distribution)
   <custom>:       any string other than 'brownian' or 'random_flow' will
                   be evaluate as is using Matlab evalf, make sure it uses
                   'DV' as the flow related independent variable, tau is
                   indexed as tau(J). Any additional parameters can be 
                   sent via "varargin"

   This file is part of Mesh-Based Monte Carlo
   License: GPLv3, see http://mcx.sf.net for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 67

   [tau,g1]=generate_g1(fhist,tau, disp_model, DV, lambda, format)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
load_mc_prop


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 273

   [mua,mus,g,n]=load_mc_prop(fname)

   Loads the absorption and scattering coefficients, as well as the 
   scattering anisotropy and index of refraction from the MMC optical
   property file given as an argument
  
   author: Stefan Carp (carp <at> nmr.mgh.harvard.edu



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37

   [mua,mus,g,n]=load_mc_prop(fname)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
loadmch


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 2202

    [data, header]=loadmch(fname,format,endian)

    author: Qianqian Fang (q.fang <at> neu.edu)

    input:
        fname: the file name to the output .mch file
        format:a string to indicate the format used to save
               the .mch file; if omitted, it is set to 'float'
        endian: optional, specifying the endianness of the binary file
               can be either 'ieee-be' or 'ieee-le' (default)

    output:
        data:   the output detected photon data array
                data has at least M*2+2 columns (M=header.medium), the first column is the 
                ID of the detector; columns 2 to M+1 store the number of 
                scattering events for every tissue region; the following M
                columns are the partial path lengths (in mm) for each medium type;
                the last column is the initial weight at launch time of each detecetd
                photon; when the momentum transfer is recorded, M columns of
                momentum tranfer for each medium is inserted after the partial path;
                when the exit photon position/dir are recorded, 6 additional columns
                are inserted before the last column, first 3 columns represent the
                exiting position (x/y/z); the next 3 columns are the dir vector (vx/vy/vz).
                in other words, data is stored in the follow format
                    [detid(1) nscat(M) ppath(M) mom(M) p(3) v(3) w0(1)]
        header: file header info, a structure has the following fields
                [version,medianum,detnum,recordnum,totalphoton,detectedphoton,
                 savedphoton,lengthunit,seedbyte,normalizer,respin,srcnum,savedetflag]
        photonseed: (optional) if the mch file contains a seed section, this
                returns the seed data for each detected photon. Each row of 
                photonseed is a byte array, which can be used to initialize a  
                seeded simulation. Note that the seed is RNG specific. You must use
                the an identical RNG to utilize these seeds for a new simulation.

    this file is part of Monte Carlo eXtreme (MCX)
    License: GPLv3, see http://mcx.sf.net for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 48

    [data, header]=loadmch(fname,format,endian)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
mmc2json


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1066

 Format:
    mmc2json(cfg,filestub)

 Save MCXLAB simulation configuration to a JSON file for MCX binary

 Author: Qianqian Fang <q.fang at neu.edu>

 Input:
    cfg: a struct defining the parameters associated with a simulation. 
         Please run 'help mcxlab' or 'help mmclab' to see the details.
         mcxpreview supports the cfg input for both mcxlab and mmclab.
    filestub: the filestub is the name stub for all output files,including
         filestub.json: the JSON input file
         filestub_vol.bin: the volume file if cfg.vol is defined
         filestub_shapes.json: the domain shape file if cfg.shapes is defined
         filestub_pattern.bin: the domain shape file if cfg.pattern is defined

 Dependency:
    this function depends on the savejson/saveubjson functions from the 
    Iso2Mesh toolbox (http://iso2mesh.sf.net) or JSONlab toolbox 
    (http://iso2mesh.sf.net/jsonlab)

 This function is part of Monte Carlo eXtreme (MCX) URL: http://mcx.space

 License: GNU General Public License version 3, please read LICENSE.txt for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 36

 Format:
    mmc2json(cfg,filestub)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
mmcadddet


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1125

 [newnode,newelem]=mmcadddet(node,elem,det,opt)
   or
 [newnode,newelem]=mmcadddet(node,elem,det,'Param1',value1, 'Param2',value2, ...)

 Adding an external wide-field (polyhedral) detector domain to an
 existing tetrahedral mesh

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     see mmcaddsrc for details.

 output:
     see mmcaddsrc for details.

 example:
     [node,face,elem]=meshasphere([0 0 0],24,1,2);
     elem(:,5)=1;

     cfg=struct('srctype','cone','srcpos',[0 0 28],'srcdir',[0 0 -1]);
     [nodesrc,elemsrc]=mmcaddsrc(node,elem,cfg);

     % example with additional options
     cfg=struct('srctype','planar','srcpos',[-30 -5 5],'srcdir',[1 0 0],...
                'srcparam1',[0 10 0 0],'srcparam2',[0 0 8 0]);
     [nodedet,elemdet]=mmcadddet(nodesrc,elemsrc,cfg);
     plotmesh(nodedet,elemdet);

     cfg=struct('srctype','disk','srcpos',[-30 0 0],'srcdir',[1 0 0],'srcparam1',[4 0 0 0]);
     [nodedet,elemdet]=mmcadddet(nodesrc,elemsrc,cfg);
     figure;plotmesh(nodedet,elemdet);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80

 [newnode,newelem]=mmcadddet(node,elem,det,opt)
   or
 [newnode,newelem]=mmc...



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
mmcaddsrc


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 2370

 [newnode,newelem]=mmcaddsrc(node,elem,src,opt)
   or
 [newnode,newelem]=mmcaddsrc(node,elem,src,'Param1',value1, 'Param2',value2, ...)

 Adding an external wide-field (polyhedral) source domain to an
 existing tetrahedral mesh

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     node: the node list of the original mesh
     elem: the elem list of the original mesh
     src:  the coordinates of the vertices that define the source domain
           if src is a struct, it will be treated as an MMCLAB input
           structure and the source vertices will be computed
           automatically
     opt: optional, a struct defining additional options, the possible fields
              are identical to those supported by meshrefine.m in iso2mesh.
              the opt struct can be replaced by a list of 'param','value' pairs

 output:
     newnode: the node list of the modified mesh, by default, 
              newnode is a concatination of src at the end of node
     newelem: the elem list of the modified mesh, by default, the elements
              between the convex hull of the combined node/src and
              the original mesh surface are labeled as 0; the elements
              sharing the source domain are labeled as -1.

 example:
     [node,face,elem]=meshasphere([0 0 0],24,1,2);
     elem(:,5)=1;

     source=double([-5 -5;-5 5;5 5;5 -5]);
     source(:,3)=30;
     [newnode,newelem]=mmcaddsrc(node,elem,source);
     plotmesh(newnode,newelem,'x>-5');

     % example with additional options
     [newnode,newelem]=mmcaddsrc(node,elem,source,'extcmdopt','-Y -a1','extlabel',-2);

     % example with mmclab cfg input
     cfg=struct('srctype','cone','srcpos',[0 0 28],'srcdir',[0 0 -1]);
     [newnode,newelem]=mmcaddsrc(node,elem,cfg);
     figure; plotmesh(newnode,newelem);

     % for box-shaped objects, you have to use meshabox with matlab 7.11
     % or newer, older matlab versions will fail to get the convex hull
     [node,face,elem]=meshabox([-10 -10 -10],[10 10 0],1,2);
     elem(:,5)=1;

     cfg=struct('srctype','fourier','srcpos',[-5 -5 10],'srcdir',...
        [0 0 -1],'srcparam1',[10 0 0 3],'srcparam2',[0 8 0 1]);
     [newnode,newelem]=mmcaddsrc(node,elem,cfg);
     figure; plotmesh(newnode,newelem);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80

 [newnode,newelem]=mmcaddsrc(node,elem,src,opt)
   or
 [newnode,newelem]=mmc...



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
mmcdettime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 768

 dett=mmcdettime(detp,prop,unitinmm)

 Recalculate the detected photon time using partial path data and 
 optical properties (for perturbation Monte Carlo or detector readings)

 author: Qianqian Fang (q.fang <at> neu.edu)
	  Ruoyang Yao (yaor <at> rpi.edu) 

 input:
     detp: the 2nd output from mmclab. detp must be a struct
     prop: optical property list, as defined in the cfg.prop field of mmclab's input
     unitinmm: voxel edge-length in mm, should use cfg.unitinmm used to generate detp; 
           if ignored, assume to be 1 (mm)

 output:
     dett: re-caculated detected photon time based on the partial path data and optical property table

 this file is copied from Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mmc.space/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37

 dett=mmcdettime(detp,prop,unitinmm)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
mmcdettpsf


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 819

 tpsf=mmcdettpsf(detp,detnum,prop,time)

 Calculate the temporal point spread function curve of a specified detector
 given the partial path data, optical properties, and distribution of time bins

 author: Qianqian Fang (q.fang <at> neu.edu)
	  Ruoyang Yao (yaor <at> rpi.edu) 

 input:
     detp: the 2nd output from mmclab. detp must be a struct
     detnum: specified detector number
     prop: optical property list, as defined in the cfg.prop field of mmclab's input
     time: distribution of time bins, a 1*3 vector [tstart tend tstep]
           could be defined different from in the cfg of mmclab's input
 output:
     tpsf: caculated temporal point spread function curve of the specified detector

 this file is copied from Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.space/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 40

 tpsf=mmcdettpsf(detp,detnum,prop,time)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
mmcdetweight


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 736

 detw=mmcdetweight(detp,prop,unitinmm)

 Recalculate the detected photon weight using partial path data and 
 optical properties (for perturbation Monte Carlo or detector readings)

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     detp: the 2nd output from mmclab. detp must a struct 
     prop: optical property list, as defined in the cfg.prop field of mmclab's input
     unitinmm: voxel edge-length in mm, should use cfg.unitinmm used to generate detp; 
           if ignored, assume to be 1 (mm)

 output:
     detw: re-caculated detected photon weight based on the partial path data and optical property table

 this file is copied from Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mmc.space/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 39

 detw=mmcdetweight(detp,prop,unitinmm)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
mmcjacobian


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1336

 [Jmua, Jmus]=mmcjacobian(cfg,detp,seeds,detnum) (element based)

 Generate time-domain Jacobians (sensitivity matrix) for absorption and scattering (mus) 
 perturbation of a specified source-detector pair with configurations, detected photon 
 seeds and detector readings from an initial MC simulation

 author: Ruoyang Yao (yaor <at> rpi.edu)
         Qianqian Fang (q.fang <at> neu.edu)

 input:
     cfg: the simulation configuration structure used for the initial MC simulation by mmclab
     detp: detector readings from the initial MC simulation
     seeds: detected photon seeds from the initial MC simulation
     detnum: the detector number whose detected photons will be replayed

 output:
     Jmua: the Jacobian for absorption coefficient for a specified source-detector pair
           also output Jmua as a by-product.
     Jmus: (optional) the Jacobian for scattering perturbation of a specified source detector pair
		  number of rows is the number of the mesh elements
		  number of columns is the number of time gates

 example:
	  [cube,detp,ncfg,seeds]=mmclab(cfg);            % initial MC simulation
	  [Jmua,Jmus] = mmcjacobian(ncfg,detp,seeds,1);  % generate scattering Jacobian of the first detector

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 65

 [Jmua, Jmus]=mmcjacobian(cfg,detp,seeds,detnum) (element based)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
mmcjmua


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1184

 Ja=mmcjmua(cfg,detp,seeds,detnum) (element based)

 Generate a time-domain Jacobian (sensitivity matrix) for absorption (mua) perturbation of a specified detector
 with configurations, detected photon seeds and detector readings from an initial MC simulation

 author: Ruoyang Yao (yaor <at> rpi.edu)
         Qianqian Fang (q.fang <at> neu.edu)

 input:
     cfg: the simulation configuration structure used for the initial MC simulation by mmclab
     detp: detector readings from the initial MC simulation, must be a
	        structure (supported after MMC v2016.4)
     seeds: detected photon seeds from the initial MC simulation
     detnum: the detector number whose detected photons will be replayed

 output:
     Ja: a Jacobian for absorption perturbation of the specified detector
		  number of rows is the number of the mesh nodes or elements
		  number of columns is the number of time gates

 example:
	  [cube,detp,ncfg,seeds]=mmclab(cfg);	% initial MC simulation
	  Jmua = mmcjmua(ncfg,detp,seeds,1);    % generate absorption Jacobian of the first detector

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 51

 Ja=mmcjmua(cfg,detp,seeds,detnum) (element based)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
mmcjmus


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1275

 [Jmus, Jmua]=mmcjmus(cfg,detp,seeds,detnum) (element based)

 Generate a time-domain Jacobian (sensitivity matrix) for scattering (mus) perturbation of a specified detector
 with configurations, detected photon seeds and detector readings from an initial MC simulation

 author: Ruoyang Yao (yaor <at> rpi.edu)
         Qianqian Fang (q.fang <at> neu.edu)

 input:
     cfg: the simulation configuration structure used for the initial MC simulation by mmclab
     detp: detector readings from the initial MC simulation
     seeds: detected photon seeds from the initial MC simulation
     detnum: the detector number whose detected photons will be replayed

 output:
     Jmus: the Jacobian for scattering perturbation of a specified source detector pair
		  number of rows is the number of the mesh elements
		  number of columns is the number of time gates
     Jmua: (optional) because calculating Jmus requires the Jacobian for mua, so one can
           also output Jmua as a by-product.

 example:
	  [cube,detp,ncfg,seeds]=mmclab(cfg);   % initial MC simulation
	  Jmus = mmcjmus(ncfg,detp,seeds,1);    % generate scattering Jacobian of the first detector

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 61

 [Jmus, Jmua]=mmcjmus(cfg,detp,seeds,detnum) (element based)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
mmclab


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 15276


         MMCLAB - Mesh-based Monte Carlo (MMC) for MATLAB/GNU Octave         %
          Copyright (c) 2010-2019 Qianqian Fang <q.fang at neu.edu>          %
                            http://mcx.space/#mmc                            %
                                                                             %
 Computational Optics & Translational Imaging (COTI) Lab- http://fanglab.org %
            Department of Bioengineering, Northeastern University            %
                                                                             %
               Research funded by NIH/NIGMS grant R01-GM114365               %

$Rev::      $2019.4$Date::                       $ by $Author::             $%


 Format:
    [fluence,detphoton,ncfg,seeds]=mmclab(cfg);
          or
    fluence=mmclab(cfg);
    newcfg=mmclab(cfg,'prep');
    [fluence,detphoton,ncfg,seeds]=mmclab(cfg, options);

 Input:
    cfg: a struct, or struct array. Each element in cfg defines 
         a set of parameters for a simulation. 

    option: (optional), options is a string, specifying additional options
         option='preview': this plots the domain configuration using mcxpreview(cfg)
         option='opencl':  force using OpenCL (set cfg.gpuid=1 if not set)
                           instead of SSE on CPUs/GPUs that support OpenCL


    cfg may contain the following fields:

== Required ==
     *cfg.nphoton:     the total number of photons to be simulated (integer)
     *cfg.prop:        an N by 4 array, each row specifies [mua, mus, g, n] in order.
                       the first row corresponds to medium type 0 which is 
                       typically [0 0 1 1]. The second row is type 1, and so on.
     *cfg.node:        node array for the input tetrahedral mesh, 3 columns: (x,y,z)
     *cfg.elem:        element array for the input tetrahedral mesh, 4 columns
     *cfg.elemprop:    element property index for input tetrahedral mesh
     *cfg.tstart:      starting time of the simulation (in seconds)
     *cfg.tstep:       time-gate width of the simulation (in seconds)
     *cfg.tend:        ending time of the simulation (in second)
     *cfg.srcpos:      a 1 by 3 vector, the position of the source in mesh node length unit
     *cfg.srcdir:      if defined as [vx, vy, vy], it specifies the incident vector
                       if defined as [vx, vy, vy, focus], the first 3 elements define
                       the incident vector; focus controls the convergence or 
                       divergence of the beam:
                       focus=0: collimated beam
                       focus<0: diverging beam from an imaginary src at c0-|focus|*[vx vy vz]
                       focus>0: converging beam, focusing to a point at c0+|focus|*[vx vy vz]
                       where c0 is the centroid of the source domain. Setting focus does 
                       not impact pencil/isotropic/cone sources.

== MC simulation settings ==
      cfg.seed:        seed for the random number generator (integer)
                       if set to a uint8 array, the binary data in each column is used 
                       to seed a photon (i.e. the "replay" mode), default value: 1648335518
      cfg.isreflect:   [1]-consider refractive index mismatch, 0-matched index
                       2 - total absorption on exterior surface
                       3 - prefect reflection (mirror) on exterior surface
      cfg.isnormalized:[1]-normalize the output fluence to unitary source, 0-no reflection
      cfg.isspecular:  [1]-calculate specular reflection if source is outside
      cfg.ismomentum:  [0]-save momentum transfer for each detected photon
      cfg.method:      ray-tracing method, ["plucker"]:Plucker, "havel": Havel (SSE4),
                       "badouel": partial Badouel, "elem": branchless Badouel (SSE), 
                       "grid": dual-grid MMC
      cfg.mcmethod:    0 use MCX-styled MC method, 1 use MCML style MC
      cfg.nout:        [1.0] refractive index for medium type 0 (background)
      cfg.minenergy:   terminate photon when weight less than this level (float) [0.0]
      cfg.roulettesize:[10] size of Russian roulette
      cfg.unitinmm:    defines the default length unit (to interpret mesh nodes, src/det positions 
                       the default value is 1.0 (mm). For example, if the mesh node length unit is 
                       in cm, one should set unitinmm to 10.
      cfg.basisorder:  [1]-linear basis, 0-piece-wise constant basis

== Source-detector parameters ==
      cfg.detpos:      an N by 4 array, each row specifying a detector: [x,y,z,radius]
      cfg.srctype:     source type, the parameters of the src are specified by cfg.srcparam{1,2}
                      'pencil' - default, pencil beam, no param needed
                      'isotropic' - isotropic source, no param needed
                      'cone' - uniform cone beam, srcparam1(1) is the half-angle in radian
                      'gaussian' - a gaussian beam, srcparam1(1) specifies the waist radius 
                                (in default length unit); if one specifies a non-zero focal length
                                using cfg.srcdir, the gaussian beam can be converging to or 
                                diverging from the waist center, which is located at srcpos+focus*srcdir;
                                optionally, one can specify the wavelength lambda (in cfg.unitinmm mm), 
                                using srcparam1(2). This will rescale the Gaussian profile according 
                                to w(z)=w0*sqrt(1-(z/z0)^2), where w0 is the waist radius, z is the 
                                distance (in mm) to the waist center (focus), and z0 is the Rayleigh 
                                range (in mm), and z0 is related to w0 by z0=w0^2*pi/lambda
                      'planar' - a 3D quadrilateral uniform planar source, with three corners specified 
                                by srcpos, srcpos+srcparam1(1:3) and srcpos+srcparam2(1:3)
                      'pattern' - a 3D quadrilateral pattern illumination, same as above, except
                                srcparam1(4) and srcparam2(4) specify the pattern array x/y dimensions,
                                and srcpattern is a floating-point pattern array, with values between [0-1]. 
                                if cfg.srcnum>1, srcpattern must be a floating-point array with 
                                a dimension of [srcnum srcparam1(4) srcparam2(4)]
                                Example: <demo_photon_sharing.m>
                      'fourier' - spatial frequency domain source, similar to 'planar', except
                                the integer parts of srcparam1(4) and srcparam2(4) represent
                                the x/y frequencies; the fraction part of srcparam1(4) multiplies
                                2*pi represents the phase shift (phi0); 1.0 minus the fraction part of
                                srcparam2(4) is the modulation depth (M). Put in equations:
                                    S=0.5*[1+M*cos(2*pi*(fx*x+fy*y)+phi0)], (0<=x,y,M<=1)
                      'arcsine' - similar to isotropic, except the zenith angle is uniform
                                distribution, rather than a sine distribution.
                      'disk' - a uniform disk source pointing along srcdir; the radius is 
                               set by srcparam1(1) (in default length unit)
                      'fourierx' - a general Fourier source, the parameters are 
                               srcparam1: [v1x,v1y,v1z,|v2|], srcparam2: [kx,ky,phi0,M]
                               normalized vectors satisfy: srcdir cross v1=v2
                               the phase shift is phi0*2*pi
                      'fourierx2d' - a general 2D Fourier basis, parameters
                               srcparam1: [v1x,v1y,v1z,|v2|], srcparam2: [kx,ky,phix,phiy]
                               the phase shift is phi{x,y}*2*pi
                      'zgaussian' - an angular gaussian beam, srcparam1(1) specifies the variance in  
                               the zenith angle
      cfg.{srcparam1,srcparam2}: 1x4 vectors, see cfg.srctype for details
      cfg.srcpattern: see cfg.srctype for details
      cfg.srcnum:     the number of source patterns that are
                      simultaneously simulated; only works for 'pattern'
                      source, see cfg.srctype='pattern' for details
                      Example <demo_photon_sharing.m>
      cfg.replaydet:  only works when cfg.outputtype is 'jacobian', 'wl', 'nscat', or 'wp' and cfg.seed is an array
                      -1 replay all detectors and save in separate volumes (output has 5 dimensions)
                       0 replay all detectors and sum all Jacobians into one volume
                       a positive number: the index of the detector to replay and obtain Jacobians
      cfg.voidtime:   for wide-field sources, [1]-start timer at launch, 0-when entering
                      the first non-zero voxel

      by default, mmc assumes the mesh and source position settings are all in mm unit.
      if the mesh coordinates/source positions are not in mm unit, one needs to define
      cfg.unitinmm  (in mm) to specify the actual length unit.

== Optional mesh data ==
     -cfg.facenb:      element face neighbohood list (calculated by faceneighbors())
     -cfg.evol:        element volume (calculated by elemvolume() with iso2mesh)
     -cfg.e0:          the element ID enclosing the source, if not defined,
                       it will be calculated by tsearchn(node,elem,srcpos);
                       if cfg.e0 is set as one of the following characters,
                       mmclab will do an initial ray-tracing and move
                       srcpos to the first intersection to the surface:
                       '>': search along the forward (srcdir) direction
                       '<': search along the backward direction
                       '-': search both directions

== Output control ==
      cfg.issaveexit: [0]-save the position (x,y,z) and (vx,vy,vz) for a detected photon
      cfg.issaveref:  [0]-save diffuse reflectance/transmittance on the exterior surfaces.
                      The output is stored as flux.dref in a 2D array of size [#Nf,  #time_gate]
                      where #Nf is the number of triangles on the surface; #time_gate is the
                      number of total time gates. To plot the surface diffuse reflectance, the output
                      triangle surface mesh can be extracted by faces=faceneighbors(cfg.elem,'rowmajor');
                      where 'faceneighbors' can be found in the iso2mesh toolbox.
                      Example: see <demo_mmclab_basic.m>
      cfg.issaveseed:  [0]-save the RNG seed for a detected photon so one can replay
      cfg.isatomic:    [1]-use atomic operations for saving fluence, 0-no atomic operations
      cfg.outputtype:  'flux' - output fluence-rate
                       'fluence' - fluence, 
                       'energy' - energy deposit, 
                       'jacobian' - mua Jacobian (replay mode)
                       'wl'- weighted path lengths to build mua Jacobian (replay mode)
                       'wp'- weighted scattering counts to build mus Jacobian (replay mode)
      cfg.debuglevel:  debug flag string, a subset of [MCBWDIOXATRPE], no space
      cfg.debugphoton: print the photon movement debug info only for a specified photon ID

      fields marked with * are required; options in [] are the default values
      fields marked with - are calculated if not given (can be faster if precomputed)


    type: omit or 'omp' for multi-threading version; 'sse' for the SSE4 MMC,
          the SSE4 version is about 25% faster, but requires newer CPUs; 
          if type='prep' with a single output, mmclab returns ncfg only.

 Output:
      fluence: a struct array, with a length equals to that of cfg.
            For each element of fluence, fluence(i).data is a 2D array with
            dimensions [size(cfg.node,1), total-time-gates] if cfg.basisorder=1,
            or [size(cfg.elem,1), total-time-gates] if cfg.basisorder=0. 
            The content of the array is the normalized fluence-rate (or others 
            depending on cfg.outputtype) at each mesh node and time-gate.
            In the "replay" mode, if cfg.replaydet is set to -1 and multiple 
            detectors exist, fluence.data will add a 5th dimension for the detector number.

            If cfg.issaveref is set to 1, fluence(i).dref is not empty, and stores
            the surface diffuse reflectance (normalized by default). The surface mesh
            that the dref output is attached can be obtained by faces=faceneighbors(cfg.elem,'rowmajor');
      detphoton: (optional) a struct array, with a length equals to that of cfg.
            Starting from v2016.5, the detphoton contains the below subfields:
              detphoton.detid: the ID(>0) of the detector that captures the photon
              detphoton.nscat: cummulative scattering event counts in each medium
              detphoton.ppath: cummulative path lengths in each medium (partial pathlength)
                   one need to multiply cfg.unitinmm with ppath to convert it to mm.
              detphoton.mom: cummulative cos_theta for momentum transfer in each medium
              detphoton.p or .v: exit position and direction, when cfg.issaveexit=1
              detphoton.w0: photon initial weight at launch time
              detphoton.prop: optical properties, a copy of cfg.prop
              detphoton.data: a concatenated and transposed array in the order of
                    [detid nscat ppath mom p v w0]'
              "data" is the is the only subfield in all MMCLAB before 2016.5
      ncfg: (optional), if given, mmclab returns the preprocessed cfg structure,
            including the calculated subfields (marked by "-"). This can be
            used in the subsequent simulations to avoid repetitive preprocessing.
      seeds: (optional), if give, mmclab returns the seeds, in the form of
            a byte array (uint8) for each detected photon. The column number
            of seed equals that of detphoton.

 Example:
      cfg.nphoton=1e5;
      [cfg.node face cfg.elem]=meshabox([0 0 0],[60 60 30],6);
      cfg.elemprop=ones(size(cfg.elem,1),1);
      cfg.srcpos=[30 30 0];
      cfg.srcdir=[0 0 1];
      cfg.prop=[0 0 1 1;0.005 1 0 1.37];
      cfg.tstart=0;
      cfg.tend=5e-9;
      cfg.tstep=5e-10;
      cfg.debuglevel='TP';
      % preprocessing to populate the missing fields to save computation
      ncfg=mmclab(cfg,'prep');

      cfgs(1)=ncfg;   % when using struct array input, all fields must be defined
      cfgs(2)=ncfg;
      cfgs(1).isreflect=0;
      cfgs(2).isreflect=1;
      cfgs(2).detpos=[30 20 0 1;30 40 0 1;20 30 1 1;40 30 0 1];
      % calculate the fluence and partial path lengths for the two configurations
      [fluxs,detps]=mmclab(cfgs);


 This function is part of Mesh-based Monte Carlo (MMC) URL: http://mcx.space/#mmc

 License: GNU General Public License version 3, please read LICENSE.txt for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 0




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
mmcmeanpath


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 533

 avgpath=mmcmeanpath(detp,prop)

 Calculate the average pathlengths for each tissue type for a given source-detector pair

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     detp: the 2nd output from mmclab. detp can be either a struct or an array (detp.data)
     prop: optical property list, as defined in the cfg.prop field of mmclab's input

 output:
     avepath: the average pathlength for each tissue type 

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 32

 avgpath=mmcmeanpath(detp,prop)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
mmcmeanscat


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 559

 avgnscat=mmcmeanscat(detp,prop)

 Calculate the average scattering event counts for each tissue type for a given source-detector pair

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     detp: the 2nd output from mmclab. detp can be either a struct or an array (detp.data)
     prop: optical property list, as defined in the cfg.prop field of mmclab's input

 output:
     avgnscat: the average scattering event count for each tissue type 

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 33

 avgnscat=mmcmeanscat(detp,prop)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
mmcraytrace


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 625

 [p,e0]=mmcraytrace(node,elem,p0,v0,e0)

 Determine the entry position and element for a ray to intersect a mesh

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     node: the mesh coordinate list
     elem: the tetrahedral mesh element list, 4 columns
     p0: origin of the ray
     v0: direction vector of the ray
     e0: search direction: '>' forward search, '<' backward search, '-' bidirectional

 output:
     p: the intersection position
     e0: if found, the index of the intersecting element ID

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 40

 [p,e0]=mmcraytrace(node,elem,p0,v0,e0)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
mmcsrcdomain


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1225

 [srcnode,srcface]=mmcsrcdomain(cfg)

 Defining a source domain (for launching new photons) in the form of 
 polyhedra based on an MMCLAB simulation configuration structure

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     cfg: the simulation configuration structure used by mmclab, the
          following fields are required
               cfg.srcpos, cfg.srcdir, cfg.srctype, cfg.srcparam1,
               cfg.srcparam2
          run "help mmclab" for details.
     meshbbx: the bounding box of the original mesh in the form of
          [min(node); max(node)], i.e. a 2x3 matrix

 output:
     srcnode: the vertices of the source domain
     srcface (optional): the polyhedral patches of the source domain

 example:
     [node,face,elem]=meshasphere([0 0 0],20,1);
     elem(:,5)=1;
     cfg=struct('srctype','cone','srcpos',[0 0 25],'srcdir',[0 0 -1]);
     srccone=mmcsrcdomain(cfg,[min(node);max(node)]);

     cfg=struct('srctype','fourier','srcpos',[-5 -5 25],'srcdir',...
        [0 0 -1],'srcparam1',[10 0 0 3],'srcparam2',[0 8 0 1]);
     srcfourier=mmcsrcdomain(cfg,[min(node);max(node)]);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37

 [srcnode,srcface]=mmcsrcdomain(cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
readmmcelem


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 387

 elem=readmmcelem(filename)

 Loading MMC mesh element file

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     filename: the file name to the element data file

 output:
     elem: the tetrahedral mesh element list 

 example:
     elem=readmmcelem('elem_sph1.dat');

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 28

 elem=readmmcelem(filename)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
readmmcface


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 399

 face=readmmcface(filename)

 Loading MMC surface triangle file

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     filename: the file name to the surface element data file

 output:
     face: the surface triangle element list 

 example:
     face=readmmcface('face_sph1.dat');

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 28

 face=readmmcface(filename)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
readmmcmesh


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 510

 [node elem]=readmmcmesh(key)

 Loading MMC node and element data files

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     key: the file name stub to the node coordinate file. The full 
          file names are {node,elem}_key.dat

 output:
     node: the node coordinate list
     elem: the tetrahedra node index list

 example:
     [node elem]=readmmcmesh('sph1');

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80

 [node elem]=readmmcmesh(key)

 Loading MMC node and element data files
...



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
readmmcnode


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 410

 node=readmmcnode(filename)

 Loading MMC node coordinates data file

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     filename: the file name to the node coordinate file

 output:
     node: the node coordinate list 

 example:
     node=readmmcnode('node_sph1.dat');

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80

 node=readmmcnode(filename)

 Loading MMC node coordinates data file

 ...



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
savemmcmesh


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 802

 savemmcmesh(key,node,elem,face,evol,facenb)

 Export a tetrahedral mesh to the MMC mesh format

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     key: a string included in all exported mesh file names, the 
          format of the files are {node,elem,face,facenb,evol}_key.dat
      node: a node coordinate list, 3 columns for x/y/z
      elem: a tetrahedral element list
      face: a triangular surface face list, if missing, it will be calculated
      evol: the volumes of all tetrahedra, if missing, it will be calculated
      facenb: the 4 neighboring elements for each element, if missing, it will
             be calculated

 example:
     savemmcmesh('sph1',node,elem);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 45

 savemmcmesh(key,node,elem,face,evol,facenb)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
spbesselh


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 409

 hn=spbesselhprime(n,k,z)

 spherical Hankel function 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Hankel function
     k: kind of the Hankel function
     z: input variable

 output:
     hn: spherical Hankel function 

 example:
     hn=spbesselh(0,1,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26

 hn=spbesselhprime(n,k,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 14
spbesselhprime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 460

 hn=spbesselhprime(n,k,z)

 spherical Hankel function first order derivative 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Hankel function
     k: kind of the Hankel function
     z: input variable

 output:
     hn: spherical Hankel function first order derivative 

 example:
     hn=spbesselhprime(0,1,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26

 hn=spbesselhprime(n,k,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
spbesselj


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 385

 jn=spbesselj(n,z)

 spherical Bessel function

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Bessel function
     z: input variable

 output:
     jn: spherical Bessel function first order derivative

 example:
     jn=spbesselj(0,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 19

 jn=spbesselj(n,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 14
spbesseljprime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 418

 jp=spbesseljprime(n,z)

 spherical Bessel function first order derivative

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Bessel function
     z: input variable

 output:
     jp: spherical Bessel function first order derivative

 example:
     jp=spbesseljprime(0,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 24

 jp=spbesseljprime(n,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
spbessely


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 391

 yn=spbessely(n,z)

 spherical Neumann function 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Neumann function
     z: input variable

 output:
     yn:  spherical Neumann function first order derivative 

 example:
     yn=spbessely(0,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 19

 yn=spbessely(n,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 14
spbesselyprime


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 424

 yp=spbesselyprime(n,z)

 spherical Neumann function first order derivative 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     n: order of the spherical Neumann function
     z: input variable

 output:
     yp:  spherical Neumann function first order derivative 

 example:
     yp=spbesselyprime(0,1)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 24

 yp=spbesselyprime(n,z)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
spharmonic


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 439

 Y=spharmonic(l,m,theta,phi)

 spherical harmonic function Y_{m,l}(theta,phi)

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     m: angular index
     l: order
     theta,phi: spherical angular coordinates, can be vectors

 output:
     Y:  the values for Y_{m,l}(theta,phi)

 example:
     Y=spharmonic(l,m,theta,phi)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 29

 Y=spharmonic(l,m,theta,phi)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 13
sphdiffAcoeff


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 444

 A=sphdiffAcoeff(m,l,cfg)

 sphere diffusion exterior solution A coefficients

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     m: angular index
     l: order
     cfg: the problem domain setup, see sphdiffusioninfinite.m

 output:
     res:  the coefficient at the specified order

 example:
     A=sphdiffAcoeff(m,l,cfg)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26

 A=sphdiffAcoeff(m,l,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 13
sphdiffBcoeff


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 444

 B=sphdiffBcoeff(m,l,cfg)

 sphere diffusion exterior solution B coefficients

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     m: angular index
     l: order
     cfg: the problem domain setup, see sphdiffusioninfinite.m

 output:
     res:  the coefficient at the specified order

 example:
     B=sphdiffBcoeff(0,1,cfg)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26

 B=sphdiffBcoeff(m,l,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 13
sphdiffCcoeff


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 444

 C=sphdiffCcoeff(m,l,cfg)

 sphere diffusion interior solution C coefficients

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     m: angular index
     l: order
     cfg: the problem domain setup, see sphdiffusioninfinite.m

 output:
     res:  the coefficient at the specified order

 example:
     C=sphdiffCcoeff(0,1,cfg)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26

 C=sphdiffCcoeff(m,l,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 15
sphdiffexterior


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 955

 res=sphdiffexterior(r,theta,phi,cfg)

 sphere exterior solution (incident+scatter) of the diffusion model

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     r,theta,phi: source position in spherical coordinates.
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both interior and exterior regions

 example:
     phi_ext=sphdiffexterior(30,pi,0,cfg);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38

 res=sphdiffexterior(r,theta,phi,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 15
sphdiffincident


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 938

 phi=sphdiffincident(r,theta,phi,cfg)

 insident field of a sphere with a diffusion model

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     r,theta,phi: source position in spherical coordinates.
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both interior and exterior regions

 example:
     phi_inc=sphdiffincident(30,pi,0,cfg);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38

 phi=sphdiffincident(r,theta,phi,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 15
sphdiffinterior


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 936

 res=sphdiffinterior(r,theta,phi,cfg)

 sphere interior solution of the diffusion model

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     r,theta,phi: source position in spherical coordinates.
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both interior and exterior regions

 example:
     phi_int=sphdiffinterior(30,pi,0,cfg);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38

 res=sphdiffinterior(r,theta,phi,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 14
sphdiffscatter


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 920

 res=sphdiffscatter(r,theta,phi,cfg)

 sphere exterior scattered field 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     r,theta,phi: source position in spherical coordinates.
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both interior and exterior regions

 example:
     phi_scat=sphdiffscatter(30,pi,0,cfg);

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37

 res=sphdiffscatter(r,theta,phi,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
sphdiffusion


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1098

 [res,xi,yi,zi]=sphdiffusion(xrange,yrange,zrange,cfg)

 diffusion solution for a sphere inside the infinite homogeneous medium 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     xrange,yrange,zrange: a vector from where a grid will be created
       		    and the phi values will be calculated
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both interior and exterior regions

 example:
   [phi_ana,xa,ya,za]=sphdiffusion(-30:0.8:30,0,-30:0.8:30);
   contourf(xa,za,log10(abs(phi_ana)),40)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 55

 [res,xi,yi,zi]=sphdiffusion(xrange,yrange,zrange,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 20
sphdiffusioninfinite


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1114

 [res,xi,yi,zi]=sphdiffusioninfinite(xrange,yrange,zrange,cfg)

 diffusion solution for a sphere inside the infinite homogeneous medium 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     xrange,yrange,zrange: a vector from where a grid will be created
       		    and the phi values will be calculated
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both interior and exterior regions

 example:
   [phi_ana,xa,ya,za]=sphdiffusioninfinite(-30:0.8:30,0,-30:0.8:30);
   contourf(xa,za,log10(abs(phi_ana)),40)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 63

 [res,xi,yi,zi]=sphdiffusioninfinite(xrange,yrange,zrange,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 23
sphdiffusionscatteronly


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1100

 [res,xi,yi,zi]=sphdiffusionscatteronly(xrange,yrange,zrange,cfg)

 analytical solution to the diffusion model (scattered field only)

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     xrange,yrange,zrange: a vector from where a grid will be created
       		    and the phi values will be calculated
     cfg: the problem domain setup: 
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both exterior region

 example:
   [phi_ana,xa,ya,za]=sphdiffusionscatteronly(-30:0.8:30,0,-30:0.8:30);
   contourf(xa,za,log10(abs(phi_ana)),40)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 66

 [res,xi,yi,zi]=sphdiffusionscatteronly(xrange,yrange,zrange,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 16
sphdiffusionsemi


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1185

 [res,xi,yi,zi]= sphdiffusionsemi(Reff,xrange,yrange,zrange,cfg)

 diffusion solution for a sphere inside a semi-infinite homogeneous medium 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     Reff:  the effective reflection coeff.
     xrange,yrange,zrange: a vector from where a grid will be created
       		    and the phi values will be calculated
     cfg: domain structure for internal/external parameters
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both the interior and exterior
     regions

 example:
   [phi,xi,yi,zi]=sphdiffusionsemi(0,-30:0.8:30,0,-30:0.8:30);
   contourf(xi,zi,log10(abs(phi)),40)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 65

 [res,xi,yi,zi]= sphdiffusionsemi(Reff,xrange,yrange,zrange,cfg)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 16
sphdiffusionslab


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1225

 [res,xi,yi,zi]= sphdiffusionslab(Reff,h,xrange,yrange,zrange,cfg)

 diffusion solution for a sphere inside an infinite homogeneous slab 

 author: Qianqian Fang (q.fang <at> neu.edu)

 input:
     Reff:  the effective reflection coeff.
     xrange,yrange,zrange: a vector from where a grid will be created
       		    and the phi values will be calculated
     h: the height of the slab
     cfg: domain structure for internal/external parameters
          cfg.v: speed of light in vacuum (mm/s)
          cfg.a: sphere radius (mm)
          cfg.omua: background (outside) mua (1/mm)
          cfg.omusp: background (outside) mus' (1/mm)
          cfg.imua: sphere (inside) mua (1/mm)
          cfg.imusp: sphere (inside) mus' (1/mm)
          cfg.src: spherical source position (R,theta,phi) R in mm
          cfg.maxl: maximum serial expansion terms
          cfg.omega: DPDW modulation frequency

 output:
     res:  the output fluence for both the interior and exterior
     regions

 example:
   [phi_ana,xa,ya,za]=sphdiffusionslab(0,0,60,-30:0.8:30,0,-30:0.8:30);
   contourf(xa,za,log10(abs(phi_ana)),40)

 this file is part of Mesh-based Monte Carlo (MMC)

 License: GPLv3, see http://mcx.sf.net/mmc/ for details




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 67

 [res,xi,yi,zi]= sphdiffusionslab(Reff,h,xrange,yrange,zrange,cfg)





