## 3 Computational Results
### Optical Depth Discretizations  
To solve the two challenges mentioned above, which are: (1) Increased error around 90\({}^{\circ}\) with higher OD and (2) increased error from modeling media with complex phase functions, various discretization investigations were performed. Before optimizing the discretization schemes, a theoretical explanation for interpreting the significant error near 90\({}^{\circ}\) was necessary. We considered a single photon that was scattered, then transmitted without being scattered again (the last scattering event). For this photon, if it is assumed that the optical depth from the scattering plane to the transmitted plane is to be \(OD_{t}\), then the probability of this photon being transmitted can be expressed as:  
\[P(OD_{t})=\exp(-\frac{OD_{t}}{\cos\theta_{t}}) \tag{7}\]  
Additionally, if we assumed that a disturbance/error \(\varepsilon\) of \(OD_{t}\) takes place during the simulation, the relative change caused by this disturbance can be written as:  
\[\frac{P(OD_{t}+\varepsilon)}{P(OD_{t})}=\exp(-\frac{OD_{t}+\varepsilon}{\cos \theta_{t}})/\exp(-\frac{OD_{t}}{\cos\theta_{t}})=\exp(-\frac{\varepsilon}{ \cos\theta_{t}}) \tag{8}\]  
In Markov chain discretizations, the disturbance \(\varepsilon\) would have a maximum value equal to the optical depth of the specific layer. As shown in Equation (8), when \(\varepsilon\) << \(\cos\theta_{t}\), \(\exp(-\varepsilon/\cos\theta_{t})\)-1.0, the OD discretization error is less significant. However, when \(\theta_{t}\) is close to 90\({}^{\circ}\), \(\cos\theta_{t}\)-0.0, and \(\exp(-\varepsilon/\cos\theta_{t})\) no longer yields a value close to 1.0, this would cause significant errors in this range. It is also worth noting that when \(\cos\theta_{t}\)-0.0, the value of \(P(OD_{t})\) will also be \(\sim\)0.0 except when \(OD_{t}\) has a sufficiently small value. Therefore, the probability for photons that scatter far away from the exit plane with a \(\theta_{t}\)-90\({}^{\circ}\), then this is too small to introduce significant errors in the angular distribution predictions. Such observations imply transmitted photons with \(\theta_{t}\)-90\({}^{\circ}\) are mostly likely from the last layer (transmitted photons) or the first layer (reflected photons). Therefore, care should be taken when modeling these two cases. In this work, two different non-uniform discretization schemes were examined. In the uniform OD discretization, each layer had 1% of the total OD (100 turbid slabs assumed). Non-uniform Equation (9) emphasized on setting the OD value of the first layer and the last layer, which can be expressed by:  
\[\left\{ \begin{array}{l}OD_{i}=a\,\ i=1,100\\ OD_{i}=(OD_{total}-2\cdot a)/98,\ i\in[2,99]\end{array} \right. \tag{9}\]  
Non-uniform Equation (10) assumed an exponential OD discretization to vary the whole discretization, which can be expressed as:  
\[\left\{\begin{array}{l}OD_{i}=\frac{OD_{total}\cdot(\  
0.3% for OD = 2 cases, and it dropped from 6% to 0.8% for OD = 5 cases. For discretization scheme 2, the error dropped from 6% to 0.3% for OD = 5 cases. These observations verified that proper discretization can help improve computational accuracy. Figure 5: Absolute relative error of Markov chain angular distribution predictions compared to Monte Carlo, with OD discretization scheme 1. Panel (**a**). Phase function 1, OD = 2; Panel (**b**). Phase function 1, OD = 5; Panel (**c**). Phase function 2, OD = 2; Panel (**d**). Phase function 2, OD = 5; Panel (**e**). Phase function 3, OD = 2; Panel (**f**). Phase function 3, OD = 5. Red dash lines indicate an error threshold of 1.0%, and blue dashed lines indicate a threshold of 0.1%. Upon closer examination, it could be readily observed that with the optimal value of either \(a\) or \(b\), the relative error was minimized when the OD in the first layer and last layer was close to 0.022 (\(b\) = 0.03 yields an OD close to 0.022 in the first and last layer).