A Sample Adding-Doubling Calculation Here are numerical results from my adding-doubling chapter S. A. Prahl, ``The Adding-Doubling Method,'' in Optical Thermal Response of Laser Irradiated Tissue, edited by A. J. Welch and M. J. C. van Gemert, Plenum Press, New York, pp. 101--129, 1995. I have supplemented the example in section 5.3.5 with a few extra intermediate results --- particularly for helping to generate the initial R and T matrices for a thin layer. There are are few discrepancies between these values and those in the book (5.64 and 5.72). These were pointed out to me by Martin Hammer (thanks!). The mistake in 5.64 was caused by too quick translation of the submatrices of the full redistribution matrix. The mistake in 5.72 was the result of insufficient checking of the manuscript before publication. If you find any others, let me know. Hope this helps, Scott Prahl https://omlc.org/~prahl/ The quadrature angles are (cf 5.63) 0.15751 0.58784 0.83024 1.00000 The quadrature weights are (cf 5.63) 0.15751 0.58784 0.83024 1.00000 angles 0.37268 0.37268 0.19098 0.06366 The h-- redistribution matrix is (cf 5.64) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.60813 1.30369 0.66305 -0.03647 0.96564 0.58784 1.30369 1.94346 1.95748 1.73691 1.84650 0.83024 0.66305 1.95748 3.15547 4.23620 2.47555 1.00000 -0.03647 1.73691 4.23620 6.84908 2.97218 flux 0.96564 1.84650 2.47555 2.97218 2.08589 The h+- redistribution matrix is (cf 5.64) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.35031 0.65833 0.23115 -0.03455 0.51588 0.58784 0.65833 0.05804 0.08635 0.34517 0.17405 0.83024 0.23115 0.08635 0.12037 0.15326 0.12265 1.00000 -0.03455 0.34517 0.15326 -0.37395 0.14817 flux 0.51588 0.17405 0.12265 0.14817 0.19459 The h-+ redistribution matrix is (cf 5.64) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.35031 0.65833 0.23115 -0.03455 0.51588 0.58784 0.65833 0.05804 0.08635 0.34517 0.17405 0.83024 0.23115 0.08635 0.12037 0.15326 0.12265 1.00000 -0.03455 0.34517 0.15326 -0.37395 0.14817 flux 0.51588 0.17405 0.12265 0.14817 0.19459 The h++ redistribution matrix is (cf 5.64) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.60813 1.30369 0.66305 -0.03647 0.96564 0.58784 1.30369 1.94346 1.95748 1.73691 1.84650 0.83024 0.66305 1.95748 3.15547 4.23620 2.47555 1.00000 -0.03647 1.73691 4.23620 6.84908 2.97218 flux 0.96564 1.84650 2.47555 2.97218 2.08589 A from equation 5.55 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.25138 -0.05967 -0.01555 0.00029 -0.00153 0.58784 -0.01599 0.06324 -0.01230 -0.00364 0.02147 0.83024 -0.00576 -0.01700 0.04761 -0.00628 0.00618 1.00000 0.00026 -0.01252 -0.01565 0.04275 -0.00498 flux 0.02072 0.01372 0.00589 0.00189 0.01055 B from equation 5.55 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.06180 0.03013 0.00542 -0.00027 0.02214 0.58784 0.00807 0.00071 0.00054 0.00072 0.00152 0.83024 0.00201 0.00075 0.00054 0.00023 0.00076 1.00000 -0.00025 0.00249 0.00057 -0.00046 0.00118 flux 0.01140 0.00440 0.00112 0.00030 0.00366 Inverse of G from equation 5.56 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.24804 -0.06120 -0.01584 0.00027 0.11472 0.58784 0.01314 1.06224 -0.01255 -0.00363 0.46253 0.83024 0.00470 0.01632 1.04733 -0.00634 0.33903 1.00000 -0.00020 0.01177 0.01508 1.04261 0.14266 flux 0.15374 0.46491 0.32669 0.12918 0.34180 G from equation 5.56 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.80193 0.04637 0.01267 0.00003 0.11849 0.58784 0.01242 0.94163 0.01133 0.00335 0.41805 0.83024 0.00469 0.01565 0.95490 0.00581 0.31097 1.00000 0.00002 0.01153 0.01447 0.95913 0.13176 flux 0.10108 0.42445 0.31111 0.12543 0.31247 R for dtau=0.25 before dividing by 2mu[i]w[i] (eq 5.54) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.08054 0.05021 0.01014 -0.00013 0.03466 0.58784 0.01345 0.00279 0.00137 0.00132 0.00340 0.83024 0.00376 0.00189 0.00113 0.00044 0.00169 1.00000 -0.00012 0.00453 0.00110 -0.00080 0.00222 flux 0.01653 0.00829 0.00229 0.00060 0.00638 R for dtau=0.25 after incorporating 2mu[i]w[i] (cf 5.65) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.68600 0.11460 0.03199 -0.00103 0.14076 0.58784 0.11460 0.00637 0.00432 0.01033 0.01893 0.83024 0.03199 0.00432 0.00357 0.00348 0.00723 1.00000 -0.00103 0.01033 0.00348 -0.00630 0.00471 flux 0.14076 0.01893 0.00723 0.00471 0.02771 T for dtau=0.25 before dividing by 2mu[i]w[i] (eq 5.54) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.60386 0.09274 0.02535 0.00005 0.11957 0.58784 0.02485 0.88326 0.02266 0.00670 0.39796 0.83024 0.00938 0.03130 0.90980 0.01162 0.30481 1.00000 0.00005 0.02305 0.02893 0.91826 0.13620 flux 0.08476 0.41075 0.30510 0.12354 0.30241 T for dtau=0.25 after incorporating 2mu[i]w[i] (cf 5.66) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 5.14345 0.21165 0.07993 0.00042 0.72199 0.58784 0.21165 2.01587 0.07144 0.05262 0.93746 0.83024 0.07993 0.07144 2.86893 0.09123 0.96210 1.00000 0.00042 0.05262 0.09123 7.21208 0.97029 flux 0.72199 0.93746 0.96210 0.97029 0.92416 R for dtau=0.5 after incorporating 2mu[i]w[i] (cf 5.67) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.95214 0.18793 0.05639 0.00160 0.21221 0.58784 0.18793 0.01727 0.01016 0.01908 0.03528 0.83024 0.05639 0.01016 0.00755 0.00684 0.01434 1.00000 0.00160 0.01908 0.00684 -0.01098 0.00932 flux 0.21221 0.03528 0.01434 0.00932 0.04611 T for dtau=0.5 after incorporating 2mu[i]w[i] (cf 5.68) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 3.15128 0.32286 0.12977 0.00829 0.55364 0.58784 0.32286 1.78949 0.13128 0.09694 0.87595 0.83024 0.12977 0.13128 2.61440 0.16846 0.92329 1.00000 0.00829 0.09694 0.16846 6.62644 0.94056 flux 0.55364 0.87595 0.92329 0.94056 0.86135 R for dtau=1.0 after incorporating 2mu[i]w[i] (cf 5.69) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.10688 0.26132 0.08567 0.00942 0.27282 0.58784 0.26132 0.04220 0.02348 0.03360 0.06089 0.83024 0.08567 0.02348 0.01612 0.01343 0.02717 1.00000 0.00942 0.03360 0.01343 -0.01668 0.01796 flux 0.27282 0.06089 0.02717 0.01796 0.06961 T for dtau=1.0 after incorporating 2mu[i]w[i] (cf 5.70) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.23674 0.38732 0.17810 0.03193 0.37545 0.58784 0.38732 1.42637 0.22054 0.16560 0.76146 0.83024 0.17810 0.22054 2.18142 0.28772 0.84595 1.00000 0.03193 0.16560 0.28772 5.60397 0.88106 flux 0.37545 0.76146 0.84595 0.88106 0.75816 R01 for the boundary is (cf 5.71) 0.15751 0.58784 0.83024 1.00000 angles 0.00000 0.00000 0.02393 0.00509 0.00000 (natural) 0.00000 0.00000 0.07545 0.04000 0.00000 *2aw R10 should be the same 0.15751 0.58784 0.83024 1.00000 angles 0.11740 0.43815 0.02393 0.00509 0.00000 (natural) 1.00000 1.00000 0.07545 0.04000 0.00000 *2aw T01 for the boundary is (cf 5.72) 0.15751 0.58784 0.83024 1.00000 angles 0.00000 0.00000 0.29320 0.12223 0.00000 (natural) 0.00000 0.00000 0.92455 0.96000 0.00000 *2aw T01 for the boundary is (cf 5.72) 0.15751 0.58784 0.83024 1.00000 angles 0.00000 0.00000 2.91544 7.53994 0.00000 (natural) 0.00000 0.00000 9.19343 59.21944 0.00000 *2aw My program's T01 is divided by 2 nu w 0.15751 0.58784 0.83024 1.00000 angles 0.00000 0.00000 0.92455 0.96000 0.00000 (natural) 0.00000 0.00000 2.91544 7.53994 0.00000 *2aw R for dtau=1.0 after adding boundaries (cf 5.73) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.00000 0.00000 0.00000 0.00000 0.00000 0.58784 0.00000 0.00000 0.00000 0.00000 0.00000 0.83024 0.00000 0.00000 0.40698 0.07282 0.13833 1.00000 0.00000 0.00000 0.07282 0.47810 0.08396 flux 0.00000 0.00000 0.13833 0.08396 0.05456 T for dtau=1.0 after adding boundaries (cf 5.74) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.00000 0.00000 0.00000 0.00000 0.00000 0.58784 0.00000 0.00000 0.00000 0.00000 0.00000 0.83024 0.00000 0.00000 1.91968 0.29193 0.64594 1.00000 0.00000 0.00000 0.29193 5.19418 0.75391 flux 0.00000 0.00000 0.64594 0.75391 0.30083 R for dtau=0.25 before dividing by 2mu[i]w[i] (eq 5.54) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.08054 0.05021 0.01014 -0.00013 0.03466 0.58784 0.01345 0.00279 0.00137 0.00132 0.00340 0.83024 0.00376 0.00189 0.00113 0.00044 0.00169 1.00000 -0.00012 0.00453 0.00110 -0.00080 0.00222 flux 0.01653 0.00829 0.00229 0.00060 0.00638 R for dtau=0.25 after incorporating 2mu[i]w[i] (cf 5.65) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.68600 0.11460 0.03199 -0.00103 0.14076 0.58784 0.11460 0.00637 0.00432 0.01033 0.01893 0.83024 0.03199 0.00432 0.00357 0.00348 0.00723 1.00000 -0.00103 0.01033 0.00348 -0.00630 0.00471 flux 0.14076 0.01893 0.00723 0.00471 0.02771 T for dtau=0.25 before dividing by 2mu[i]w[i] (eq 5.54) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.60386 0.09274 0.02535 0.00005 0.11957 0.58784 0.02485 0.88326 0.02266 0.00670 0.39796 0.83024 0.00938 0.03130 0.90980 0.01162 0.30481 1.00000 0.00005 0.02305 0.02893 0.91826 0.13620 flux 0.08476 0.41075 0.30510 0.12354 0.30241 T for dtau=0.25 after incorporating 2mu[i]w[i] (cf 5.66) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 5.14345 0.21165 0.07993 0.00042 0.72199 0.58784 0.21165 2.01587 0.07144 0.05262 0.93746 0.83024 0.07993 0.07144 2.86893 0.09123 0.96210 1.00000 0.00042 0.05262 0.09123 7.21208 0.97029 flux 0.72199 0.93746 0.96210 0.97029 0.92416 R for dtau=0.5 after incorporating 2mu[i]w[i] (cf 5.67) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.95214 0.18793 0.05639 0.00160 0.21221 0.58784 0.18793 0.01727 0.01016 0.01908 0.03528 0.83024 0.05639 0.01016 0.00755 0.00684 0.01434 1.00000 0.00160 0.01908 0.00684 -0.01098 0.00932 flux 0.21221 0.03528 0.01434 0.00932 0.04611 T for dtau=0.5 after incorporating 2mu[i]w[i] (cf 5.68) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 3.15128 0.32286 0.12977 0.00829 0.55364 0.58784 0.32286 1.78949 0.13128 0.09694 0.87595 0.83024 0.12977 0.13128 2.61440 0.16846 0.92329 1.00000 0.00829 0.09694 0.16846 6.62644 0.94056 flux 0.55364 0.87595 0.92329 0.94056 0.86135 R for dtau=1.0 after incorporating 2mu[i]w[i] (cf 5.69) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.10688 0.26132 0.08567 0.00942 0.27282 0.58784 0.26132 0.04220 0.02348 0.03360 0.06089 0.83024 0.08567 0.02348 0.01612 0.01343 0.02717 1.00000 0.00942 0.03360 0.01343 -0.01668 0.01796 flux 0.27282 0.06089 0.02717 0.01796 0.06961 T for dtau=1.0 after incorporating 2mu[i]w[i] (cf 5.70) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 1.23674 0.38732 0.17810 0.03193 0.37545 0.58784 0.38732 1.42637 0.22054 0.16560 0.76146 0.83024 0.17810 0.22054 2.18142 0.28772 0.84595 1.00000 0.03193 0.16560 0.28772 5.60397 0.88106 flux 0.37545 0.76146 0.84595 0.88106 0.75816 T01 for the boundary is (BIG DISCREPANCY!) (cf 5.72) 0.15751 0.58784 0.83024 1.00000 angles 0.00000 0.00000 0.29320 0.12223 0.00000 (natural) 0.00000 0.00000 0.92455 0.96000 0.00000 *2aw R01 for the boundary is (cf 5.71) 0.15751 0.58784 0.83024 1.00000 angles 0.00000 0.00000 0.02393 0.00509 0.00000 (natural) 0.00000 0.00000 0.07545 0.04000 0.00000 *2aw ** Note that 0.07545+0.92455=1.0000 and 0.04000+0.96000=1.0000 R for dtau=1.0 after adding boundaries (cf 5.73) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.00000 0.00000 0.00000 0.00000 0.00000 0.58784 0.00000 0.00000 0.00000 0.00000 0.00000 0.83024 0.00000 0.00000 0.40698 0.07282 0.13833 1.00000 0.00000 0.00000 0.07282 0.47810 0.08396 <=== eq 5.75 flux 0.00000 0.00000 0.13833 0.08396 0.05456 The zeros on the diagonal in R are not important because those are for virtual angles anyhow. I must have changed the default behavior of my program to set them to zero instead of letting them be unity. (I may very well have eliminated calculating them to improve efficiency.) T for dtau=1.0 after adding boundaries (cf 5.74) 0 0.15751 0.58784 0.83024 1.00000 flux 0.15751 0.00000 0.00000 0.00000 0.00000 0.00000 0.58784 0.00000 0.00000 0.00000 0.00000 0.00000 0.83024 0.00000 0.00000 1.91968 0.29193 0.64594 1.00000 0.00000 0.00000 0.29193 5.19418 0.75391 <===eq 5.76 flux 0.00000 0.00000 0.64594 0.75391 0.30083