Access kappa and molecular weights from AeroData (#336)

Author: jcurtis2

src/aero_data.F90

@@ -185,6 +185,28 @@ subroutine f_aero_data_get_species_density(ptr_c, idx, val) bind(C)
 
   end subroutine
 
+  subroutine f_aero_data_get_species_kappa(ptr_c, idx, val) bind(C)
+    type(aero_data_t), pointer :: ptr_f => null()
+    type(c_ptr), intent(in) :: ptr_c
+    integer(c_int), intent(in) :: idx
+    real(c_double), intent(out) :: val
+
+    call c_f_pointer(ptr_c, ptr_f)
+    val = ptr_f%kappa(idx+1)
+
+  end subroutine
+
+  subroutine f_aero_data_get_species_molecular_weight(ptr_c, idx, val) bind(C)
+    type(aero_data_t), pointer :: ptr_f => null()
+    type(c_ptr), intent(in) :: ptr_c
+    integer(c_int), intent(in) :: idx
+    real(c_double), intent(out) :: val
+
+    call c_f_pointer(ptr_c, ptr_f)
+    val = ptr_f%molec_weight(idx+1)
+
+  end subroutine
+
   subroutine f_aero_data_n_source(ptr_c, n_source) bind(C)
     type(aero_data_t), pointer :: ptr_f => null()
     type(c_ptr), intent(in) :: ptr_c

src/aero_data.hpp

@@ -28,6 +28,8 @@ extern "C" void f_aero_data_vol2rad(const void *ptr, const double*, double*) noe
 extern "C" void f_aero_data_diam2vol(const void *ptr, const double*, double*) noexcept;
 extern "C" void f_aero_data_vol2diam(const void *ptr, const double*, double*) noexcept;
 extern "C" void f_aero_data_get_species_density(const void *ptr, const int *idx, double *val) noexcept;
+extern "C" void f_aero_data_get_species_kappa(const void *ptr, const int *idx, double *val) noexcept;
+extern "C" void f_aero_data_get_species_molecular_weight(const void *ptr, const int *idx, double *val) noexcept;
 extern "C" void f_aero_data_source_name_by_index(const void *ptr, const int *i_source,
     char *name_data) noexcept;
 extern "C" void f_aero_data_spec_name_by_index(const void *ptr, const int *i_spec,
@@ -252,5 +254,41 @@ struct AeroData {
         }
         return names;
     }
+
+    static auto kappa(const AeroData &self) {
+        int len;
+        f_aero_data_len(
+            self.ptr.f_arg(),
+            &len
+        );
+        std::valarray<double> data(len);
+
+        for (int idx = 0; idx < len; idx++) {
+             f_aero_data_get_species_kappa(
+                 self.ptr.f_arg(),
+                 &idx,
+                 &data[idx]
+            );
+        }
+        return data;
+    }
+
+    static auto molecular_weights(const AeroData &self) {
+        int len;
+        f_aero_data_len(
+            self.ptr.f_arg(),
+            &len
+        );
+        std::valarray<double> data(len);
+
+        for (int idx = 0; idx < len; idx++) {
+             f_aero_data_get_species_molecular_weight(
+                 self.ptr.f_arg(),
+                 &idx,
+                 &data[idx]
+            );
+        }
+        return data;
+    }
 };
 

src/pypartmc.cpp

@@ -97,6 +97,10 @@ PYBIND11_MODULE(_PyPartMC, m) {
             "Radius of primary particles (m)")
         .def_property_readonly("densities", &AeroData::densities,
             "Return array of aerosol species densities")
+        .def_property_readonly("kappa", &AeroData::kappa,
+            "Returns array of aerosol species hygroscopicity parameter kappa")
+        .def_property_readonly("molecular_weights", &AeroData::molecular_weights,
+            "Returns array of aerosol species molecular weights")
         .def("density", &AeroData::density, "Return density of an aerosol species")
         .def("rad2vol", AeroData::rad2vol,
             "Convert geometric radius (m) to mass-equivalent volume (m^3).")

tests/test_aero_data.py

@@ -339,6 +339,100 @@ def test_aero_data_densities():
         # assert
         assert aero_data_densities == densities
 
+    @staticmethod
+    def test_aero_data_kappa():
+        # arrange
+        kappa = [0.65, 0.1, 0.0, 0.0]
+        sut = ppmc.AeroData(
+            (
+                {
+                    "SO4": [
+                        1800.0 * si.kg / si.m**3,
+                        0,
+                        96.0 * si.g / si.mol,
+                        kappa[0],
+                    ]
+                },
+                {
+                    "OC": [
+                        1400.0 * si.kg / si.m**3,
+                        0,
+                        1.0 * si.g / si.mol,
+                        kappa[1],
+                    ]
+                },
+                {
+                    "BC": [
+                        1800.0 * si.kg / si.m**3,
+                        0,
+                        1.0 * si.g / si.mol,
+                        kappa[2],
+                    ]
+                },
+                {
+                    "H2O": [
+                        1000.0 * si.kg / si.m**3,
+                        0,
+                        18.0 * si.g / si.mol,
+                        kappa[3],
+                    ]
+                },
+            )
+        )
+
+        # act
+        aero_data_kappa = sut.kappa
+
+        # assert
+        assert aero_data_kappa == kappa
+
+    @staticmethod
+    def test_aero_data_molecular_weight():
+        # arrange
+        molec_weight = [96.0, 1.0, 1.0, 18.0]
+        sut = ppmc.AeroData(
+            (
+                {
+                    "SO4": [
+                        1800.0 * si.kg / si.m**3,
+                        0,
+                        molec_weight[0] * si.kg / si.mol,
+                        0.65,
+                    ]
+                },
+                {
+                    "OC": [
+                        1400.0 * si.kg / si.m**3,
+                        0,
+                        molec_weight[1] * si.kg / si.mol,
+                        0.1,
+                    ]
+                },
+                {
+                    "BC": [
+                        1800.0 * si.kg / si.m**3,
+                        0,
+                        molec_weight[2] * si.kg / si.mol,
+                        0.0,
+                    ]
+                },
+                {
+                    "H2O": [
+                        1000.0 * si.kg / si.m**3,
+                        0,
+                        molec_weight[3] * si.kg / si.mol,
+                        0.0,
+                    ]
+                },
+            )
+        )
+
+        # act
+        aero_data_molec_weight = sut.molecular_weights
+
+        # assert
+        assert aero_data_molec_weight == molec_weight
+
     @staticmethod
     def test_aero_data_density():
         # arrange