This Notebook was written under Mathematica 3.01 to calculate the total transmission
and reflection of a stack of an arbitrary number of layers of different thickness, absorption
and scattering coefficient as well as scattering anisotropy by Martin Hammer (hammer@bach.med.uni-jena.de).
The basic algorithm is the Adding-Doubling method first described by H.C.van de Hulst (1963, A New Look at
Multiple Scattering. Unnumbered mimeographed report, NASA Institute for Space SciencesNew York). 
This Program follows the description by Scott A. Prahl (in: Optical-
Thermal Response of Laser-Irradiated Tissue. ed.:A.J. Welch and M.J.C. van Gemert,
Plenum Press New York, 1995). The Henyey-Greenstein function is used as phase function.
In this version of the program neither internal sources nor 
reflecting layers (i.e. differences of the real part of the refractive index of the layers) are
included. In the example given here, the optical properties of the ocular fundus at 560 nm 
(Hammer et al., Optical properties of ocular fundus tissue..., Phys. Med. Biol. 40 (1995)
pp. 963-978) are used.

n = 8     Quadrature = 1     layers = 4

<< NumericalMath`GaussianQuadrature`

\(\( (*calculation\ of\ the\ quadrature\ angles\ and\ weights*) 

    If[Quadrature == 0, 
      G = GaussianQuadratureWeights[n, \ 0, \ 1], {
        X = N[Solve[
              LegendreP[n - 1, x] + 
                  \(x - 1\)\/n*Dt[LegendreP[n - 1, x], x] == 0, \ x]], 
        G = Table[{If[i < n, \ 1/2 - 1/2\ Re[x] /. X[\([n - i]\)], 1], 
              1\/\(2\ \((1 - Re[x])\)\ Dt[LegendreP[n - 1, Re[x]], x]^2\) /. 
                If[i < n, X[\([n - i]\)], x -> \(-1. \)]}, \ {i, n}]}]\)\)

<< DiscreteMath`KroneckerDelta`

II = Table[KroneckerDelta[i - j] ,  {i, n},  {j, n}]EE = 
      Table[KroneckerDelta[i - j]/(2 G[[i, 1]] G[[i, 2]]) , 
         {i, n},  {j, n}]F = 
      Table[KroneckerDelta[i - j] (2 G[[i, 1]] G[[i, 2]]) , 
         {i, n},  {j, n}]c = 
      Table[KroneckerDelta[i - j] G[[i, 2]] ,  {i, n},  {j, n}]

\(
 (*
      Optical\ properties\ of\ the\ single\ layers\ of\ the\ 
        stack . 
			For\ each\ layer\ one\ entry\ has\ to\ be\ given\ in\ 
        the\ vectors
			d, \ μs, \ μa, \ 
      and\ g . \ The\ first\ elements\ of\ these\ vectors
			belongs\ to\ 
        the\ first\ layer\ of\ the\ stack\ which\ is\ 
			exposed\ to\ the
        \ \(irradiance . \)*) 
d = {0.2, \ 0.01, \ 0.25, \ 0.7}\)μs = {2.4,  144. ,  46. ,  70. }μa = {0.11,  17. ,  3.5,  0.18}g = {0.97,  0.84,  0.945,  0.9}a = μs/(μs + μa)tau = (μs + μa)*dttau = (1 - a (g^n)) tauaa = ((a (1 - (g^n)))/(1 - a (g^n)))deltatau = ttau

m = Table[0,  {i, layers}]

\(\( (*\ Calculation\ of\ the\ optical\ depth\ of\ the\ initial\ layers, 
      
	the\ doubling\ steps\ per\ layer, \ the\ redistribution\ matrices, 
      \ and\ the\ reflection\ and\ transmission\ matrices\ of\ the\ single\ 
        layers*) 
	
    For[l = 1, \ l <= layers, \ \(l++\), 
      
{While[deltatau[\([l]\)] >= G[\([1, 1]\)], {deltatau[\([l]\)] *= 1/2, 
            \ \(m[\([l]\)]++\)}], 

        hpp = Table[
            ∑\+\(k = 0\)\%\(n - 1\)\((2\ k + 1)\)*
                \((g[\([l]\)]^k - g[\([l]\)]^n)\)/\((1 - g[\([l]\)]^n)\)*
                LegendreP[k, G[\([i, 1]\)]]*LegendreP[k, G[\([j, 1]\)]], 
            \ {i, n}, \ {j, n}], 

        hpm = Table[
            ∑\+\(k = 0\)\%\(n - 1\)\((2\ k + 1)\)*
                \((g[\([l]\)]^k - g[\([l]\)]^n)\)/\((1 - g[\([l]\)]^n)\)*
                LegendreP[k, G[\([i, 1]\)]]*LegendreP[k, \(-G[\([j, 1]\)]\)], 
            \ {i, n}, \ {j, n}], 

        A = Table[KroneckerDelta[i - j]/G[\([i, 1]\)], \ {i, n}, \ {j, n}] . 
              \((II - aa[\([l]\)]\/2\ hpp . c)\)*deltatau[\([l]\)]/2, 

        B = aa[\([l]\)]\/2*
            Table[KroneckerDelta[i - j]/G[\([i, 1]\)], \ {i, n}, \ {j, n}] . 
              hpm . c*deltatau[\([l]\)]/2, 

        GG = Inverse[II + A - B . Inverse[II + A] . B], 

        R1 = 2*GG . B . Inverse[II + A], 

        R = Table[R1[\([i, j]\)]/\((2\ G[\([j, 1]\)]\ G[\([j, 2]\)])\), 
            \ {i, n}, {j, n}], 
T1 = 2\ GG - II, 

        T = Table[T1[\([i, j]\)]/\((2\ G[\([j, 1]\)]\ G[\([j, 2]\)])\), 
            \ {i, n}, {j, n}], 

        For[i = 1, \ i <= m[\([l]\)], \ \(i++\), 

          	{R1 = R + T . Inverse[EE - R . F . R] . R . F . T, 
\ \ 
            T1 = T . Inverse[EE - R . F . R] . T, 
	R = R1, 
	T = T1}], 

        TT[l - 1, l] = \(TT[l, l - 1] = T1\), 

        RR[l - 1, l] = \(RR[l, l - 1] = R1\), 
}]\)\)

\(\( (*
      Calculating\ the\ reflection\ and\ transmission\ matrices
			of\ the
        \ whole\ stack\ by\ adding\ the\ single\ layers*) 

    For[i = 2, \ i <= layers, \ \(i++\), 

      	{RR[i, 0] = 
          RR[i, i - 1] + 
            TT[i - 1, i] . Inverse[EE - RR[i - 1, 0] . F . RR[i - 1, i]] . 
              RR[i - 1, 0] . F . TT[i, i - 1], 
\ \ 
        RR[0, i] = 
          RR[0, i - 1] + 
            TT[i - 1, 0] . Inverse[EE - RR[i - 1, i] . F . RR[i - 1, 0]] . 
              RR[i - 1, i] . F . TT[0, i - 1], 
	\ 
        TT[i, 0] = 
          TT[i - 1, 0] . Inverse[EE - RR[i - 1, i] . F . RR[i - 1, 0]] . 
            TT[i, i - 1], 
	\ 
        TT[0, i] = 
          TT[i - 1, i] . Inverse[EE - RR[i - 1, 0] . F . RR[i - 1, i]] . 
            TT[0, i - 1]}]\)\)

\( (*Calculation\ of\ the\ total\ reflection\ and\ transmission\ of

        			the\ stack\ for\ irradiance\ normal\ to\ the\ surface\ 
        \((Radau - 
				quadrature)\)\ or\ under\ the\ smallest\ 
        quadrature\ angle\ 
			\((Gauss - quadrature)\)*) 
Rc = 0\)For[i = 1,  i <= n,  i++,  
      Rc += 2 G[[i, 1]] G[[i, 2]] RR[0, layers][[i, n]]]

Rc

Tc = 0For[i = 1,  i <= n,  i++,  
      Tc += 2 G[[i, 1]] G[[i, 2]] TT[layers, 0][[i, n]]]

Tc