293 implement dust attempt2

docs/src/4_extra_options/including_dust_continuum.rst

@@ -0,0 +1,41 @@
+Including dust or continuum sources in Magritte
+===============================================
+
+Starting from Magritte version 0.9.11, you can include additional continuum sources in your model.
+These are treated in LTE; therefore you need to provide a dust/continuum opacity table (per point, per frequency) and temperature (per point).
+This can be done using:
+
+.. code-block:: python
+
+    model.dust.dust_opacities.set(continuum_opacity_table)# in 1/m dims: [N_points, N_continuum_frequencies]
+    model.dust.dust_frequencies.set(frequency_table)# in Hz dims: [N_continuum_frequencies]
+    model.dust.dust_temperatures.set(temperature_table)# in K dims: [N_points]
+
+Evidently, we have also provides some tools to help calculate the dust opacities from tabulated data.
+In particular, if you have a table from the Jena dust database (https://www2.astro.uni-jena.de/Laboratory/OCDB/), you can use the following commands:
+
+.. code-block:: python
+
+    import magritte.tools as tools
+    continuum_reference_density = astropy.Quantity # in mass density [e.g. g/cm^3]; given by the continuum data; used for rescaling the values
+    frequency_grid, complex_refractive_index = tools.read_dust_opacity_table(continuum_file_location, ASTROPY_UNIT_OF_FREQ_WAVELENGTH_...)
+    # reads the file and returns the frequency grid (first column), multiplied whatever astropy unit you specified, and the complex refractive index (unitless, third column) as a numpy array
+    frequency_grid_Hz, absorption_coeff_per_density = tools.convert_dust_opacity_table_to_SI(frequency_grid, complex_refractive_index, continuum_ref_density)
+    # returns the frequency grid in Hz and absorption coeff in m^2/kg (needs to be multiplied with the mass density to get the absorption coeff per m)
+    #for this example, assume we scale the continuum/dust number density with the H2 number density
+    continuum_opacity_table = (absorption_coeff_per_density/1000) * fraction_density_dust * nH2 * 2.016 / 6.022e23 #H2: 2.016 g/mol, 6.022e23 particles/mol -> particles/m3 to g/m3
+
+In which the last step converts the absorption coefficient per mass density (m^2/kg) to an absorption coefficient per number density (1/m).
+In this release, we also included options to image the model using a custom frequency grid (given that equidistant grids might not be optimal for imaging the continuum).
+
+.. code-block:: python
+
+    # Instead of model.compute_spectral_discretization(...)
+    model.set_custom_spectral_discretization(frequency_grid_Hz)# np.array, in Hz, dims: [N_frequencies_to_image]
+    model.compute_image_new(...) #as usual
+
+
+.. Warning::
+
+    Even though you might have a model for which you only have continuum sources (no line sources), you still need to include a dummy species producing lines.
+    You can even set its abundance to zero, but Magritte needs at least one species to properly initialize a model. See tests/benchmarks/analytic/continuum_implementation.py for an example.

docs/src/4_extra_options/index.rst

@@ -11,3 +11,4 @@ we will explain a few interesting options for our users.
    other_solvers
    adaptive_ng_acceleration
    NLTE_in_low_density
+   including_dust_continuum

magritte/tools.py

@@ -3,6 +3,7 @@
 
 from datetime          import datetime
 from time              import perf_counter
+import pandas as pd
 from astropy.io        import fits
 from astropy           import units, constants
 from scipy.interpolate import griddata, interp1d
@@ -775,3 +776,97 @@ def check_one_line_approximation(model):
     contribution = np.exp(-diff_max**2)
 
     return contribution
+
+
+def convert_dust_opacity_table_to_SI(wavelengths_wavenumbers_wavefrequencies: units.Quantity, complex_refractory_index: np.ndarray, reference_density: units.Quantity, new_frequency_grid: units.Quantity=None, left_fill_value: float=0.0, right_fill_value: float=0.0):
+    """Convert a dust opacity table (wavelengths, wavenumbers or frequencies) with complex refractive index to absorption coefficient in SI units (m^-1) and rescaled to unit density (kg/m^3).
+
+    Args:
+        wavelengths_wavenumbers_wavefrequencies (Astropy.units.Quantity): 1D array containing the wavelengths, wavenumbers or frequencies in compatible units (m, m^-1 or Hz).
+        complex_refractory_index (np.ndarray): 1D array containing the complex refractive index values (unitless). Same length as `wavelengths_wavenumbers_wavefrequencies`.
+        reference_density (astropy.Quantity): Refernce density in units compatible with kg/m^3, used for rescaling the absorption coefficient.
+        new_frequency_grid (Astropy.units.Quantity, optional): New frequency grid to interpolate values upon. If None, use `wavelengths_wavenumbers_wavefrequencies' instead. Defaults to None.
+        left_fill_value (float, optional): Fill value for dust opacity (in m^2/kg) for frequencies lower than the input frequency range. Defaults to 0.0.
+        right_fill_value (float, optional): Fill value for dust opacity (in m^2/kg) for frequencies higher than the input frequency range. Defaults to 0.0.
+
+    Raises:
+        ValueError: Incompatible units for `wavelengths_wavenumbers_wavefrequencies` or `reference_density` or `new_frequency_grid'.
+        TypeError: `wavelengths_wavenumbers_wavefrequencies` or `reference_density` or `new_frequency_grid' is not an Astropy Quantity.
+
+    Returns:
+        tuple(Astropy.units.Quantity, Astropy.units.Quantity): frequencies in Hz, rescaled absorption coefficient in m^2/kg (needs to be multiplied by density in kg/m^3 to get absorption coefficient in m^-1).
+    """
+
+    frequencies_in_hz = None
+    wavelengths_in_m = None
+    reference_density_in_kg_per_m3 = None
+    #first check whether the wavelengths_or_wavenumbers are have astropy units compatible with m (wavelength), m^-1 (wavenumber) of Hz (frequency)
+    #if so, convert to Hz
+    if isinstance(wavelengths_wavenumbers_wavefrequencies, units.Quantity):
+        if wavelengths_wavenumbers_wavefrequencies.unit.is_equivalent(units.m) or wavelengths_wavenumbers_wavefrequencies.unit.is_equivalent(1.0/units.m) or wavelengths_wavenumbers_wavefrequencies.unit.is_equivalent(units.Hz):
+            # Convert to Hz
+            frequencies_in_hz = wavelengths_wavenumbers_wavefrequencies.to(units.Hz, equivalencies=units.spectral())
+            wavelengths_in_m = wavelengths_wavenumbers_wavefrequencies.to(units.m, equivalencies=units.spectral())
+            pass
+        else:
+            raise ValueError("Input must be in units compatible to m or Hz.")
+    else:
+        raise TypeError("Input must be an astropy Quantity. This is for ease of unit conversion.")
+
+    if isinstance(reference_density, units.Quantity):
+        if not reference_density.unit.is_equivalent(units.kg/units.m**3):
+            raise ValueError("Reference density must be in units compatible to kg/m^3.")
+        # Convert to kg/m^3
+        reference_density_in_kg_per_m3 = reference_density.to(units.kg/units.m**3)
+    else:
+        raise TypeError("Reference density must be an astropy Quantity.")
+
+
+    # Convert complex refractory index to absorption coefficient
+    absorption_coefficient = 4*np.pi * complex_refractory_index / wavelengths_in_m #[m^-1]
+    #rescale to unit density (kg/m^3)
+    rescaled_absorption_coefficient = absorption_coefficient / reference_density_in_kg_per_m3 #[m^2/kg]
+
+    #co-sort frequencies and rescaled absorption coefficient (because of possible unit conversion)
+    sorted_indices = np.argsort(frequencies_in_hz)
+    frequencies_in_hz = frequencies_in_hz[sorted_indices]
+    rescaled_absorption_coefficient = rescaled_absorption_coefficient[sorted_indices]
+
+    if new_frequency_grid == None:
+        return frequencies_in_hz, rescaled_absorption_coefficient
+    else: 
+        new_frequency_grid_in_hz = None
+        # Convert new frequency grid to Hz if it is not already
+        if isinstance(new_frequency_grid, units.Quantity):
+            if new_frequency_grid.unit.is_equivalent(units.m) or new_frequency_grid.unit.is_equivalent(1.0/units.m) or new_frequency_grid.unit.is_equivalent(units.Hz):
+                # Convert to Hz
+                new_frequency_grid_in_hz = new_frequency_grid.to(units.Hz, equivalencies=units.spectral())
+            else:
+                raise ValueError("New frequency grid must be in units compatible to m, m^-1 or Hz.")
+        else:
+            raise TypeError("New frequency grid must be an astropy Quantity in Hz.")
+        # Interpolate the data onto the new frequency grid
+        # interpolation = interp1d(frequencies_in_hz, rescaled_absorption_coefficient, bounds_error=False, fill_value=0.0)
+        interpolated_rescaled_absorption_coefficient = np.interp(new_frequency_grid_in_hz, frequencies_in_hz, rescaled_absorption_coefficient, left=left_fill_value, right=right_fill_value)
+        return new_frequency_grid, interpolated_rescaled_absorption_coefficient
+
+
+def read_dust_opacity_table(filename: str, wavelength_frequency_unit: units.Unit):
+    """Simple script for reading headerless dust opacity data files (with file format like the Jena dust database).
+    Returns the first column (wavelengths, wavenumbers or frequencies) and the third column (complex refractive index).
+
+    Args:
+        filename (str): Location of the dust opacity table file.
+        wavelength_frequency_unit (Astropy.units.Unit): Unit of the first column (wavelengths, wavenumbers or frequencies).
+
+    Returns:
+        tuple(Astropy.units.Quantity, np.ndarray): wavelengths or wavenumbers or frequencies, complex refractive index
+    """
+    reader = pd.read_csv(filename, sep = '\s+', header=None, comment='#')
+    wavelength_frequency = np.array(reader[0]) * wavelength_frequency_unit # wavelengths or wavenumbers or frequencies
+    complex_refractive_index = np.array(reader[2]) # complex refractive index; unitless
+
+
+    return wavelength_frequency, complex_refractive_index
+
+

src/CMakeLists.txt

@@ -26,6 +26,7 @@ set (SOURCE_FILES
     model/radiation/radiation.cpp
     model/radiation/frequencies/frequencies.cpp
     model/image/image.cpp
+    model/dust/dust.cpp
     solver/solver.cpp
 )
 

src/bindings/CMakeLists.txt

@@ -28,6 +28,7 @@ set (SOURCE_FILES
     ../model/radiation/radiation.cpp
     ../model/radiation/frequencies/frequencies.cpp
     ../model/image/image.cpp
+    ../model/dust/dust.cpp
     ../solver/solver.cpp
 )
 

src/bindings/pybindings.cpp

@@ -129,6 +129,7 @@ PYBIND11_MODULE(core, module) {
         .def_readonly("thermodynamics", &Model::thermodynamics)
         .def_readonly("radiation", &Model::radiation)
         .def_readonly("images", &Model::images)
+        .def_readonly("dust", &Model::dust)
         .def_readwrite("eta", &Model::eta)
         .def_readwrite("chi", &Model::chi)
         .def_readonly("S_ray", &Model::S_ray)
@@ -179,6 +180,12 @@ PYBIND11_MODULE(core, module) {
             "images with the given min and max frequency. Can also specify the "
             "amount of frequency bins to use (instead of defaulting to "
             "parameters.nfreqs).")
+        .def("set_custom_spectral_discretization",
+            (int(Model::*)(
+                const py::array_t<Real, py::array::c_style | py::array::forcecast> frequency_grid))
+                & Model::set_custom_spectral_discretization,
+            "Set a custom spectral discretization for the model, using the "
+            "given frequency grid.")
         .def("compute_LTE_level_populations", &Model::compute_LTE_level_populations,
             "Compute the level populations for the model assuming local "
             "thermodynamic equilibrium (LTE).")
@@ -680,6 +687,21 @@ PYBIND11_MODULE(core, module) {
         .def("read", &Frequencies::read, "Read object from file.")
         .def("write", &Frequencies::write, "Write object to file.");
 
+    // Dust
+    py::class_<Dust>(module, "Dust", "Class containing the dust/continuum properties.")
+        // attributes
+        .def_readonly("n_dust_frequencies", &Dust::n_dust_frequencies,
+            "Number of dust/continuum frequency bins.")
+        .def_readwrite("dust_frequencies", &Dust::dust_frequencies,
+            "Array with the dust/continuum frequency bins.")
+        .def_readwrite(
+            "dust_opacities", &Dust::dust_opacities, "Dust/continuum opacity per frequency bin.")
+        .def_readwrite("dust_temperature", &Dust::dust_temperature,
+            "Dust/continuum temperature per model point.")
+        // functions
+        .def("read", &Dust::read, "Read object from file.")
+        .def("write", &Dust::write, "Write object to file.");
+
     // Vector <Size>
     py::class_<Vector<Size>>(module, "VSize", py::buffer_protocol())
         // buffer

src/model/dust/dust.cpp

@@ -0,0 +1,43 @@
+#include "dust.hpp"
+
+/// Reader for dust data
+///    @param[in] io: io data object
+void Dust::read(const Io& io) {
+    cout << "Reading dust..." << endl;
+
+    // by default initialized to zero
+    n_dust_frequencies = 0;
+    // Note: dust is optional to include in the model, so first check if data is present
+    int err = io.read_number("dust/.n_dust_frequencies", n_dust_frequencies);
+    if (err == 0 && n_dust_frequencies > 0) {
+        dust_temperature.resize(parameters->npoints());
+        dust_frequencies.resize(n_dust_frequencies);
+        dust_opacities.resize(parameters->npoints(), n_dust_frequencies);
+
+        io.read_list("dust/dust_temperature", dust_temperature);
+        io.read_list("dust/dust_frequencies", dust_frequencies);
+        io.read_list("dust/dust_opacities", dust_opacities);
+
+        cout << "  Read dust/continuum opacities at " << n_dust_frequencies << " frequencies."
+             << endl;
+    } else {
+        // in this case, n_dust_frequencies remains zero, so when calculating dust properties, this
+        // must be checked to skip that step in case of no dust present in the model.
+        cout << "No dust/continuum data found." << endl;
+    }
+}
+
+/// Writer for dust data
+///    @param[in] io: io data object
+void Dust::write(const Io& io) const {
+    cout << "Writing dust..." << endl;
+
+    // Note: data will be written either way, even if it all zeros
+    io.write_list("dust/dust_temperature", dust_temperature);
+    io.write_list("dust/dust_frequencies", dust_frequencies);
+    io.write_list("dust/dust_opacities", dust_opacities);
+
+    Size temp_n_dust_frequencies = dust_frequencies.size();
+
+    io.write_number("dust/.n_dust_frequencies", temp_n_dust_frequencies);
+}

src/model/dust/dust.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include "io/io.hpp"
+#include "model/parameters/parameters.hpp"
+#include "tools/types.hpp"
+
+/// Data structure for dust or continuum radiation in general; Note: these opacity sources will be
+/// treated in LTE, using a single temperature for all at once
+struct Dust {
+    Size n_dust_frequencies; ///< number of frequencies at which dust opacities are defined
+    std::shared_ptr<Parameters> parameters; ///< data structure containing model parameters
+    Vector<Real> dust_temperature;          ///< dust temperature at each point
+    Vector<Real> dust_frequencies;          ///< frequency grid at which dust opacities are defined
+    Matrix<Real> dust_opacities;            ///< dust opacities (point, dust frequency index)
+    // Matrix<Real> dust_emissivities; ///< = dust opacities * B_nu(T_dust)
+    //  Note: we cannot precalculate dust emissivities, as it depends on the exact frequency
+
+    Dust(std::shared_ptr<Parameters> params) : parameters(params){};
+
+    void read(const Io& io);
+    void write(const Io& io) const;
+
+    void compute_dust_opacity_emissivity(
+        Size p, Real freq, Real& dust_opacity, Real& dust_emissivity) const;
+
+    inline Real planck(Real temp, Real freq) const;
+};
+
+#include "dust.tpp"

src/model/dust/dust.tpp

@@ -0,0 +1,59 @@
+#include "tools/constants.hpp"
+#include "tools/types.hpp"
+
+///  Planck function: copied from solver.cpp
+///    @param[in] temp : temperature of the corresponding
+///    black body
+///    @param[in] freq : frequency at which to evaluate the
+///    function
+///    @return Planck function evaluated at this frequency
+///////////////////////////////////////////////////////////////////////////
+inline Real Dust::planck(Real temp, Real freq) const {
+    return TWO_HH_OVER_CC_SQUARED * (freq * freq * freq)
+         / expm1(HH_OVER_KB * freq / temp); // best not to use float version here
+}
+
+///  Computes the dust opacity by interpolating the precalculated dust opacities
+///    @param[in] p : index of the cell
+///    @param[in] freq : Comoving frame frequency at which to evaluate the dust opacity
+///    @param[out] dust_opacity : dust opacity at this point, at this frequency
+///    @param[out] dust_emissivity : dust emissivity at this point, at this frequency
+inline void Dust::compute_dust_opacity_emissivity(
+    Size p, Real freq, Real& dust_opacity, Real& dust_emissivity) const {
+    // Return 0 if no dust present
+    if (this->n_dust_frequencies == 0) {
+        dust_opacity    = 0.0;
+        dust_emissivity = 0.0;
+        return;
+    }
+    // If frequency is outside the range of the dust opacities, return the boundary value
+    if (freq < dust_frequencies[0]) {
+        dust_opacity    = dust_opacities(p, 0);
+        dust_emissivity = dust_opacity * planck(dust_temperature[p], freq);
+        return;
+    }
+    if (freq > dust_frequencies[dust_frequencies.size() - 1]) {
+        dust_opacity    = dust_opacities(p, dust_frequencies.size() - 1);
+        dust_emissivity = dust_opacity * planck(dust_temperature[p], freq);
+        return;
+    }
+
+    // find the right bin using binary search
+    Size f_low  = 0;
+    Size f_high = dust_frequencies.size() - 1;
+    while (f_high - f_low > 1) {
+        Size f_mid = (f_high + f_low) / 2;
+        if (dust_frequencies[f_mid] > freq) {
+            f_high = f_mid;
+        } else {
+            f_low = f_mid;
+        }
+    }
+
+    // and do linear interpolation
+    Real slope = (dust_opacities(p, f_high) - dust_opacities(p, f_low))
+               / (dust_frequencies[f_high] - dust_frequencies[f_low]);
+    dust_opacity    = dust_opacities(p, f_low) + slope * (freq - dust_frequencies[f_low]);
+    dust_emissivity = dust_opacity * planck(dust_temperature[p], freq);
+    return;
+}