adding control on experiment detector field configuration

Author: MartinPdeS

PyMieSim/_version.py

@@ -21,4 +21,4 @@
 __version__ = version = '5.1.0'
 __version_tuple__ = version_tuple = (5, 1, 0)
 
-__commit_id__ = commit_id = 'g3fac92db6'
+__commit_id__ = commit_id = 'g5b180cc9b'

PyMieSim/cpp/experiment/detector_set/detector_set.cpp

@@ -9,6 +9,7 @@ PhotodiodeSet::PhotodiodeSet(
     const std::vector<double> &gamma_offset,
     const PolarizationSet &polarization_filter_set,
     const MediumSet &medium_set,
+    const std::vector<std::vector<complex128>> &angular_weights,
     const bool is_sequential
 )
 :   sampling(sampling),
@@ -17,8 +18,32 @@ PhotodiodeSet::PhotodiodeSet(
     phi_offset(phi_offset),
     gamma_offset(gamma_offset),
     polarization_filter_set(polarization_filter_set),
-    medium(medium_set)
+    medium(medium_set),
+    angular_weights(angular_weights),
+    has_angular_weights(!angular_weights.empty())
     {
+        if (this->has_angular_weights) {
+            if (!is_sequential) {
+                throw std::invalid_argument(
+                    "angular_weights are currently supported only for sequential PhotodiodeSet builds."
+                );
+            }
+
+            if (this->angular_weights.size() != this->sampling.size()) {
+                throw std::invalid_argument(
+                    "angular_weights must provide one weight vector per sequential detector."
+                );
+            }
+
+            for (size_t index = 0; index < this->angular_weights.size(); ++index) {
+                if (this->angular_weights[index].size() != this->sampling[index]) {
+                    throw std::invalid_argument(
+                        "Each angular weight vector must match the detector sampling at the same sequential index."
+                    );
+                }
+            }
+        }
+
         this->is_sequential = is_sequential;
         this->is_empty = false;
         this->update_shape();
@@ -51,6 +76,12 @@ PhotodiodeSet::get_detector_by_index(long long flat_index) const {
         this->medium[indices[6]]
     );
 
+    if (this->has_angular_weights) {
+        throw std::invalid_argument(
+            "angular_weights are currently supported only in sequential PhotodiodeSet mode."
+        );
+    }
+
     detector->indices = indices;
 
     return detector;
@@ -67,6 +98,11 @@ PhotodiodeSet::get_detector_by_index_sequential(size_t index) const {
         this->polarization_filter_set.polarization_states[index],
         this->medium[index]
     );
+
+    if (this->has_angular_weights) {
+        detector->set_angular_weights(this->angular_weights[index]);
+    }
+
     detector->indices = {index};
     return detector;
 }

PyMieSim/cpp/experiment/detector_set/detector_set.h

@@ -47,6 +47,8 @@ class PhotodiodeSet : public BaseDetectorSet
         std::vector<double> gamma_offset;
         PolarizationSet polarization_filter_set;
         MediumSet medium;
+        std::vector<std::vector<complex128>> angular_weights;
+        bool has_angular_weights = false;
 
         PhotodiodeSet() = default;
 
@@ -58,6 +60,7 @@ class PhotodiodeSet : public BaseDetectorSet
             const std::vector<double> &gamma_offset,
             const PolarizationSet &polarization_filter_set,
             const MediumSet &medium,
+            const std::vector<std::vector<complex128>> &angular_weights,
             const bool is_sequential
         );
 

PyMieSim/cpp/experiment/detector_set/interface.cpp

@@ -14,6 +14,105 @@
 
 namespace py = pybind11;
 
+namespace {
+std::vector<std::vector<complex128>> cast_py_to_sequential_angular_weights(
+    const py::object& angular_weights,
+    const std::vector<unsigned>& sampling,
+    const size_t target_size
+) {
+    std::vector<std::vector<complex128>> weights_by_detector;
+
+    if (angular_weights.is_none()) {
+        return weights_by_detector;
+    }
+
+    if (py::isinstance<py::array>(angular_weights)) {
+        py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> array =
+            angular_weights.cast<py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast>>();
+
+        if (array.ndim() == 1) {
+            const size_t sampling_size = static_cast<size_t>(array.shape(0));
+
+            for (const unsigned detector_sampling : sampling) {
+                if (sampling_size != detector_sampling) {
+                    throw py::value_error(
+                        "One-dimensional angular_weights arrays require a uniform detector sampling."
+                    );
+                }
+            }
+
+            std::vector<complex128> shared_weight_vector;
+            vector_assign_from_numpy(shared_weight_vector, array);
+            weights_by_detector.assign(target_size, shared_weight_vector);
+            return weights_by_detector;
+        }
+
+        if (array.ndim() == 2) {
+            if (static_cast<size_t>(array.shape(0)) != target_size) {
+                throw py::value_error(
+                    "Two-dimensional angular_weights arrays must have shape (target_size, sampling)."
+                );
+            }
+
+            const size_t sampling_size = static_cast<size_t>(array.shape(1));
+            for (const unsigned detector_sampling : sampling) {
+                if (sampling_size != detector_sampling) {
+                    throw py::value_error(
+                        "Two-dimensional angular_weights arrays require a uniform sampling shared by all detectors."
+                    );
+                }
+            }
+
+            auto unchecked = array.unchecked<2>();
+            weights_by_detector.resize(target_size);
+
+            for (size_t detector_index = 0; detector_index < target_size; ++detector_index) {
+                weights_by_detector[detector_index].resize(sampling_size);
+
+                for (size_t sample_index = 0; sample_index < sampling_size; ++sample_index) {
+                    weights_by_detector[detector_index][sample_index] = unchecked(
+                        static_cast<py::ssize_t>(detector_index),
+                        static_cast<py::ssize_t>(sample_index)
+                    );
+                }
+            }
+
+            return weights_by_detector;
+        }
+
+        throw py::value_error("angular_weights must be one- or two-dimensional.");
+    }
+
+    py::sequence sequence = py::reinterpret_borrow<py::sequence>(angular_weights);
+
+    if (static_cast<size_t>(py::len(sequence)) != target_size) {
+        throw py::value_error(
+            "Sequence angular_weights inputs must provide one weight vector per sequential detector."
+        );
+    }
+
+    weights_by_detector.resize(target_size);
+
+    for (size_t detector_index = 0; detector_index < target_size; ++detector_index) {
+        py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast> vector =
+            py::cast<py::array_t<std::complex<double>, py::array::c_style | py::array::forcecast>>(sequence[detector_index]);
+
+        std::vector<complex128> cast_weights;
+        vector_assign_from_numpy(cast_weights, vector);
+
+        if (cast_weights.size() != sampling[detector_index]) {
+            throw py::value_error(
+                "Each angular weight vector must match the detector sampling at the same sequential index."
+            );
+        }
+
+        weights_by_detector[detector_index] = std::move(cast_weights);
+    }
+
+    return weights_by_detector;
+}
+}
+
 
 PYBIND11_MODULE(detector_set, module) {
     py::object ureg = get_shared_ureg();
@@ -109,6 +208,7 @@ PYBIND11_MODULE(detector_set, module) {
                         Casting::cast_py_to_vector<double>(gamma_offset, "radian"),
                         Casting::Polarization::cast_py_to_polarization_set(polarization_filter),
                         Casting::Material::create_material_set_from_pyobject<MediumSet, double, BaseMedium, ConstantMedium>(medium, "medium"),
+                        std::vector<std::vector<complex128>>{},
                         false
                     );
                 }
@@ -243,7 +343,8 @@ PYBIND11_MODULE(detector_set, module) {
                 const py::object& sampling,
                 const py::object& cache_numerical_aperture,
                 const py::object& polarization_filter,
-                const py::object& medium
+                const py::object& medium,
+                const py::object& angular_weights
             ) {
                 std::vector<unsigned> sampling_value =
                     Casting::cast_py_to_broadcasted_vector<unsigned>(
@@ -295,6 +396,13 @@ PYBIND11_MODULE(detector_set, module) {
                         target_size
                 );
 
+                std::vector<std::vector<complex128>> angular_weights_value =
+                    cast_py_to_sequential_angular_weights(
+                        angular_weights,
+                        sampling_value,
+                        target_size
+                    );
+
                 return std::make_shared<PhotodiodeSet>(
                     sampling_value,
                     numerical_aperture_value,
@@ -303,6 +411,7 @@ PYBIND11_MODULE(detector_set, module) {
                     gamma_offset_value,
                     polarization_filter_set,
                     medium_set,
+                    angular_weights_value,
                     true
                 );
             },
@@ -314,6 +423,7 @@ PYBIND11_MODULE(detector_set, module) {
             py::arg("cache_numerical_aperture") = py::float_(0.0),
             py::arg("polarization_filter") = PolarizationState(),
             py::arg("medium") = ConstantMedium(1.0),
+            py::arg("angular_weights") = py::none(),
             R"pdoc(
                 Build a sequential photodiode detector set.
 
@@ -344,6 +454,12 @@ PYBIND11_MODULE(detector_set, module) {
                     Surrounding medium. This can be provided either as medium
                     objects or as refractive index values. Default is
                     ``ConstantMedium(1.0)``.
+                angular_weights : array-like, optional
+                    Explicit complex angular weight vectors for each sequential
+                    detector. This can be a single one-dimensional array shared
+                    by all detectors, a two-dimensional array of shape
+                    ``(target_size, sampling)``, or a list of one-dimensional
+                    arrays when detector sampling varies by sequential index.
 
                 Returns
                 -------

PyMieSim/cpp/single/detector/base.h

@@ -98,8 +98,8 @@ class BaseDetector {
         gamma_offset(_gamma_offset),
         polarization_filter(polarization_filter),
         medium(std::move(_medium)),
-        angular_weights(_sampling, complex128(1.0, 0.0)),
-        mean_coupling(_mean_coupling)
+        mean_coupling(_mean_coupling),
+        angular_weights(_sampling, complex128(1.0, 0.0))
     {}
 
     void set_angular_weights(const std::vector<complex128>& weights) {

docs/examples/extras/README.rst

@@ -14,3 +14,4 @@ Contents
 - ``plot_Qsca_vs_permittivity_vs_size_parameter.py`` - inspect scattering as material and size vary.
 - ``array_scattering.py`` - compute fields and Stokes parameters for arbitrary angles.
 - ``photodiode_angular_weights.py`` - apply a custom angular mask on a single photodiode detector.
+- ``experiment_photodiode_angular_weights.py`` - use explicit angular weight vectors in a sequential experiment detector set.

tests/validation/test_sequential.py

@@ -88,6 +88,50 @@ def test_sequential_vs_standard_no_detector():
         data_sequential == data_standard
     ), "Mismatch between sequential and standard computed data."
 
+def test_sequential_photodiode_angular_weights_reduce_coupling():
+    size = 2
+    sampling = 128
+
+    source_sequential = GaussianSet.build_sequential(
+        target_size=size,
+        wavelength=[1000, 1000] * ureg.nanometer,
+        polarization=PolarizationSet(angles=[0, 0] * ureg.degree),
+        optical_power=[1, 1] * ureg.watt,
+        numerical_aperture=[0.2, 0.2],
+    )
+
+    scatterer_sequential = SphereSet.build_sequential(
+        target_size=size,
+        diameter=[1000, 1000] * ureg.nanometer,
+        material=[1.5 + 0.5j, 1.5 + 0.5j],
+        medium=[1.0, 1.0],
+    )
+
+    angular_weights = numpy.ones((size, sampling), dtype=numpy.complex128)
+    angular_weights[1, :] = 0.0
+
+    detector_sequential = PhotodiodeSet.build_sequential(
+        target_size=size,
+        phi_offset=[90, 90] * ureg.degree,
+        gamma_offset=[0, 0] * ureg.degree,
+        numerical_aperture=[0.1, 0.1],
+        cache_numerical_aperture=[0.0, 0.0],
+        sampling=[sampling, sampling],
+        polarization_filter=[0, 0] * ureg.degree,
+        medium=[1.0, 1.0],
+        angular_weights=angular_weights,
+    )
+
+    setup_sequential = Setup(
+        scatterer_set=scatterer_sequential,
+        source_set=source_sequential,
+        detector_set=detector_sequential,
+    )
+
+    coupling = setup_sequential.get_sequential("coupling")
+
+    assert coupling[0] > 0
+    assert numpy.isclose(coupling[1], 0.0, atol=1e-20)
 
 def test_sequential_vs_standard_detector():
     source_standard = GaussianSet(