Feature/197 powerlawscatterer option to define lambda0

src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs

@@ -7,6 +7,9 @@ namespace Vts.Test.Modeling.Spectroscopy
     [TestFixture]
     public class PowerLawScattererTests
     {
+        /// <summary>
+        /// Test constructor that specifies power law coefficients A, B, C, D
+        /// </summary>
         [Test]
         public void Test_power_law_scatterer_constructor()
         {
@@ -18,13 +21,34 @@ public void Test_power_law_scatterer_constructor()
             Assert.That(scatterer.D, Is.EqualTo(0.0));
         }
 
+        /// <summary>
+        /// Test constructor that specifies power law coefficients A, B, C, D,
+        /// and lambda0 the wavelength normalization factor
+        /// </summary>
+        [Test]
+        public void Test_power_law_scatterer_constructor_with_lambda0_specification()
+        {
+            var scatterer = new PowerLawScatterer(1, 0.1, 0.0, 0.0, 1000);
+            Assert.That(scatterer, Is.InstanceOf<PowerLawScatterer>());
+            Assert.That(scatterer.A, Is.EqualTo(1));
+            Assert.That(scatterer.B, Is.EqualTo(0.1));
+            Assert.That(scatterer.C, Is.EqualTo(0.0));
+            Assert.That(scatterer.D, Is.EqualTo(0.0));
+        }
+
+        /// <summary>
+        /// Test that wrong specification of Tissue Type throws exception
+        /// </summary>
         [Test]
         public void Test_set_tissue_type_undefined()
         {
             var scatterer = new PowerLawScatterer();
             Assert.Throws<ArgumentOutOfRangeException>(() => scatterer.SetTissueType((TissueType) 100));
         }
 
+        /// <summary>
+        /// Test ability to call GetMusp method with and without lambda0
+        /// </summary>
         [Test]
         public void Verify_user_ability_to_specify_lambda0()
         {
@@ -38,7 +62,6 @@ public void Verify_user_ability_to_specify_lambda0()
             // set up call to GetMusp with another lambda0
             musp = scatterer.GetMusp(1000, 800);
             Assert.That(musp, Is.Not.EqualTo(4.0));
-
         }
     }
 }

src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs

@@ -1,5 +1,4 @@
 using NUnit.Framework;
-using Vts.MonteCarlo;
 using Vts.SpectralMapping;
 
 namespace Vts.Test.Modeling.Spectroscopy
@@ -9,6 +8,9 @@ public class TissueTests
     {
         private Tissue _tissue;
 
+        /// <summary>
+        /// Instantiate Tissue with predefined chromophores and scatterer
+        /// </summary>
         [OneTimeSetUp]
         public void One_time_setup()
         {
@@ -28,6 +30,9 @@ public void One_time_setup()
                 n);
         }
 
+        /// <summary>
+        /// Test setup specification of chromophores and optical properties
+        /// </summary>
         [Test]
         public void Test_tissue_constructor()
         {
@@ -38,6 +43,9 @@ public void Test_tissue_constructor()
             Assert.That(_tissue.Absorbers[3].Concentration, Is.EqualTo(0.87));
         }
 
+        /// <summary>
+        /// Test instantiation of Tissue with TissueTypes
+        /// </summary>
         [Test]
         public void Test_tissue_constructor_tissue_type()
         {
@@ -52,40 +60,58 @@ public void Test_tissue_constructor_tissue_type()
             Assert.That(tissue.ScattererType, Is.EqualTo(ScatteringType.PowerLaw));
         }
 
+        /// <summary>
+        /// Test GetMua method
+        /// </summary>
         [Test]
         public void Test_get_mua()
         {
             var mua = _tissue.GetMua(1000);
             Assert.That(mua, Is.EqualTo(0.067854).Within(0.000001));
         }
 
+        /// <summary>
+        /// Test GetMusp method
+        /// </summary>
         [Test]
         public void Test_get_musp()
         {
             var musp = _tissue.GetMusp(1000);
             Assert.That(musp, Is.EqualTo(0.839999).Within(0.000001));
         }
 
+        /// <summary>
+        /// Test GetMus method
+        /// </summary>
         [Test]
         public void Test_get_mus()
         {
             var mus = _tissue.GetMus(1000);
             Assert.That(mus, Is.EqualTo(4.2).Within(0.000001));
         }
 
+        /// <summary>
+        /// Test GetG method
+        /// </summary>
         [Test]
         public void Test_get_g()
         {
             var g = _tissue.GetG(1000);
             Assert.That(g, Is.EqualTo(0.8));
         }
 
+        /// <summary>
+        /// Test ToString method
+        /// </summary>
         [Test]
         public void Test_to_string()
         {
             Assert.That(_tissue.ToString(), Is.EqualTo("test_tissue"));
         }
 
+        /// <summary>
+        /// Test GetOpticalProperties method with single wavelength
+        /// </summary>
         [Test]
         public void Test_get_optical_properties()
         {
@@ -97,6 +123,25 @@ public void Test_get_optical_properties()
             Assert.That(opticalProperties.G, Is.EqualTo(0.8));
         }
 
+        /// <summary>
+        /// Test GetOpticalProperties method with single wavelength and lambda0 specification.
+        /// Lambda0 is the wavelength normalization specified in PowerLawScatterer
+        /// </summary>
+        [Test]
+        public void Test_get_optical_properties_with_lambda0_specification()
+        {
+            const double lambda0 = 1000;
+            var opticalProperties = _tissue.GetOpticalProperties(1000, lambda0);
+            Assert.That(opticalProperties.N, Is.EqualTo(1.4));
+            Assert.That(opticalProperties.Mua, Is.EqualTo(0.067854).Within(0.000001));
+            Assert.That(opticalProperties.Mus, Is.EqualTo(4.2).Within(0.000001));
+            Assert.That(opticalProperties.Musp, Is.EqualTo(0.84).Within(0.000001));
+            Assert.That(opticalProperties.G, Is.EqualTo(0.8));
+        }
+
+        /// <summary>
+        /// Test GetOpticalProperties with array of wavelengths
+        /// </summary>
         [Test]
         public void Test_get_optical_properties_wavelength_array()
         {
@@ -119,5 +164,33 @@ public void Test_get_optical_properties_wavelength_array()
             Assert.That(opticalPropertyArray[2].Musp, Is.EqualTo(0.84).Within(0.000001));
             Assert.That(opticalPropertyArray[2].G, Is.EqualTo(0.8));
         }
+
+
+        /// <summary>
+        /// Test GetOpticalProperties with array of wavelengths and lambda0 specification
+        /// </summary>
+        [Test]
+        public void Test_get_optical_properties_wavelength_array_and_lambda0_specification()
+        {
+            const double lambda0 = 1000;
+            var wavelengths = new double[] { 600, 700, 1000 };
+            var opticalPropertyArray = _tissue.GetOpticalProperties( 
+                wavelengths, lambda0);
+            Assert.That(opticalPropertyArray[0].N, Is.EqualTo(1.4));
+            Assert.That(opticalPropertyArray[0].Mua, Is.EqualTo(0.314619).Within(0.000001));
+            Assert.That(opticalPropertyArray[0].Mus, Is.EqualTo(5.562449).Within(0.000001));
+            Assert.That(opticalPropertyArray[0].Musp, Is.EqualTo(1.112489).Within(0.000001));
+            Assert.That(opticalPropertyArray[0].G, Is.EqualTo(0.8));
+            Assert.That(opticalPropertyArray[1].N, Is.EqualTo(1.4));
+            Assert.That(opticalPropertyArray[1].Mua, Is.EqualTo(0.036097).Within(0.000001));
+            Assert.That(opticalPropertyArray[1].Mus, Is.EqualTo(5.110287).Within(0.000001));
+            Assert.That(opticalPropertyArray[1].Musp, Is.EqualTo(1.022057).Within(0.000001));
+            Assert.That(opticalPropertyArray[1].G, Is.EqualTo(0.8));
+            Assert.That(opticalPropertyArray[2].N, Is.EqualTo(1.4));
+            Assert.That(opticalPropertyArray[2].Mua, Is.EqualTo(0.067854).Within(0.000001));
+            Assert.That(opticalPropertyArray[2].Mus, Is.EqualTo(4.2).Within(0.000001));
+            Assert.That(opticalPropertyArray[2].Musp, Is.EqualTo(0.84).Within(0.000001));
+            Assert.That(opticalPropertyArray[2].G, Is.EqualTo(0.8));
+        }
     }
 }

src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs

@@ -4,30 +4,46 @@ namespace Vts.SpectralMapping
 {
     /// <summary>
     /// Returns scattering values based on Steve Jacques' Skin Optics Summary:
-    /// https://omlc.ogi.edu/news/jan98/skinoptics.html
+    /// https://omlc.org/news/jan98/skinoptics.html
     /// This returned reduced scattering follows the approximate formula:
-    /// mus' = A1*lamda(-b1) + A2*lambda(-b2)
+    /// mus' = A1*(lamda/lambda0)(-b1) + A2*(lambda/lambda0)(-b2)
+    /// with default value of lambda0=1000nm
     /// </summary>
     public class PowerLawScatterer : BindableObject, IScatterer
     {
         private double _a;
         private double _b;
         private double _c;
         private double _d;
+        private double _lambda0;
 
         /// <summary>
-        /// Constructs a power law scatterer; i.e. mus' = a*lamda^-b + c*lambda^-d
+        /// Constructs a power law scatterer; i.e. mus' = a*(lamda/lambda0)^-b + c*(lambda/lambda0)^-d
         /// </summary>
         /// <param name="a">The first prefactor</param>
         /// <param name="b">The first exponent</param>
         /// <param name="c">The second prefactor</param>
         /// <param name="d">The second exponent</param>
-        public PowerLawScatterer(double a, double b, double c, double d)
+        /// <param name="lambda0">Wavelength normalization factor</param>
+        public PowerLawScatterer(double a, double b, double c, double d, double lambda0)
         {
             A = a;
             B = b;
             C = c;
             D = d;
+            Lambda0 = lambda0;
+        }
+
+        /// <summary>
+        /// Constructs a power law scatterer; i.e. mus' = a*lamda^-b + c*lambda^-d
+        /// </summary>
+        /// <param name="a">The first prefactor</param>
+        /// <param name="b">The first exponent</param>
+        /// <param name="c">The second prefactor</param>
+        /// <param name="d">The second exponent</param>
+        public PowerLawScatterer(double a, double b, double c, double d)
+            : this(a, b, c, d, 1000.0)
+        {
         }
 
         /// <summary>
@@ -36,7 +52,7 @@ public PowerLawScatterer(double a, double b, double c, double d)
         /// <param name="a">The first prefactor</param>
         /// <param name="b">The first exponent</param>
         public PowerLawScatterer(double a, double b)
-            : this(a,b,0.0,0.0)
+            : this(a,b,0.0,0.0, 1000.0)
         {
         }
 
@@ -70,6 +86,7 @@ public PowerLawScatterer()
         /// <param name="tissueType">Tissue type</param>
         public void SetTissueType(TissueType tissueType)
         {
+            Lambda0 = 1000;
             switch (tissueType)
             {
                 case TissueType.Skin:
@@ -176,6 +193,19 @@ public double D
             }
         }
 
+        /// <summary>
+        /// The wavelength normalization factor
+        /// </summary>
+        public double Lambda0
+        {
+            get => _lambda0;
+            set
+            {
+                _lambda0 = value;
+                OnPropertyChanged(nameof(Lambda0));
+            }
+        }
+
         /// <summary>
         /// Returns mus' based on Steve Jacques' Skin Optics Summary:
         /// https://omlc.ogi.edu/news/jan98/skinoptics.html
@@ -184,7 +214,7 @@ public double D
         /// <returns>The reduced scattering coefficient Mus'</returns>
         public double GetMusp(double wavelength)
         {
-            return A * Math.Pow(wavelength/1000, - B) + C * Math.Pow(wavelength/1000, - D);
+            return GetMusp(wavelength, 1000.0);
         }
 
         /// <summary>