Skip to content

Commit

Permalink
Prototype of the scattering species interface. (#852)
Browse files Browse the repository at this point in the history 
Adds a HenyeyGreensteinScatterer that can be used as a scattering
species and provides functionality to calculate bulk scattering
properties.
  • Loading branch information
@riclarsson
riclarsson authored Oct 30, 2024
2 parents f5fd6ed + 4bbfeb0 commit 9402278
Show file tree
Hide file tree
Showing 25 changed files with 1,832 additions and 488 deletions.
222 changes: 170 additions & 52 deletions examples/getting-started/4-scattering/scattering_species.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
"source": [
"# 4. Scattering calculations\n",
"\n",
"This example demonstrates the basic principles of performing scattering calcualtions with ARTS. Populations of particles that scatter radiation are referred to as *scattering species*. In essence, a scattering species defines a mapping from one or several fields of *scattering species properties* to *scattering properties* which form the input to the actual scattering calculation."
"This example demonstrates the basic principles of performing scattering calcualtions with ARTS. Populations of particles that scatter radiation are referred to as *scattering species*. Each scattering species defines a mapping from one or several fields of *scattering species properties* to *scattering properties* which form the input to the actual scattering calculation."
]
},
{
Expand All @@ -25,20 +25,16 @@
"metadata": {},
"outputs": [],
"source": [
"# import numpy as np\n",
"# import pyarts\n",
"import numpy as np\n",
"import pyarts\n",
"\n",
"# ws = pyarts.workspace.Workspace()\n",
"# ws.PlanetSet(option = \"Earth\")\n",
"# ws.atm_fieldInit(toa=100e3)\n",
"# ws.atm_fieldAddGriddedData(\n",
"# key=pyarts.arts.String(\"t\"),\n",
"# data=pyarts.arts.GriddedField3.fromxml(\"planets/Earth/afgl/tropical/t.xml\")\n",
"# )\n",
"# ws.atm_fieldAddGriddedData(\n",
"# key=pyarts.arts.String(\"p\"),\n",
"# data=pyarts.arts.GriddedField3.fromxml(\"planets/Earth/afgl/tropical/p.xml\")\n",
"# )"
"ws = pyarts.Workspace()\n",
"ws.surface_fieldSetPlanetEllipsoid(option=\"Earth\")\n",
"ws.surface_field[pyarts.arts.SurfaceKey(\"t\")] = 295.0\n",
"ws.atmospheric_fieldRead(\n",
" toa=100e3, basename=\"planets/Earth/afgl/tropical/\", missing_is_zero=1\n",
")\n",
"ws.atmospheric_fieldIGRF(time=\"2000-03-11 14:39:37\")"
]
},
{
Expand All @@ -48,7 +44,15 @@
"source": [
"## Add a field of scattering species properties\n",
"\n",
"In order to perform a simulation with ice particles, it is necessary to first define the *scattering species property* that represents the properties of the scattering species in the atmosphere. Below we define the **mass density** of **ice** as the scattering species property we use to describe the distribution of ice hydrometeors in the atmosphere."
"For this example, we will use scatterers with a Henyey-Greenstein phase function to represent ice particles. We will use two *scattering species properties* to represent the scattering species *ice*: The extinction and the single-scattering albed (SSA). The ice extinction and SSA thus become atmospheric fields. Below we define scattering-species-property object that identify these fields."
]
},
{
"cell_type": "markdown",
"id": "fde7da44-6123-44b6-a545-66e6271a892d",
"metadata": {},
"source": [
"### Ice"
]
},
{
Expand All @@ -58,15 +62,16 @@
"metadata": {},
"outputs": [],
"source": [
"# ice_mass_density = pyarts.arts.ScatteringSpeciesProperty(\"ice\", pyarts.arts.ParticulateProperty(\"MassDensity\"))"
"ice_extinction = pyarts.arts.ScatteringSpeciesProperty(\"ice\", pyarts.arts.ParticulateProperty(\"Extinction\"))\n",
"ice_ssa = pyarts.arts.ScatteringSpeciesProperty(\"ice\", pyarts.arts.ParticulateProperty(\"SingleScatteringAlbedo\"))"
]
},
{
"cell_type": "markdown",
"id": "7dfa6a96-2456-42b2-a455-4b03ab254068",
"id": "59bfd385-6125-4751-94da-c06f85a0f599",
"metadata": {},
"source": [
"We then define a GriddedField3 representing the ice mass density profile in the atmosphere and add it to ``atm_field`` of the workspace."
"We then define a GriddedField3 representing the ice extinction and single-scattering albedo and add it to ``atm_field`` of the workspace."
]
},
{
Expand All @@ -76,12 +81,70 @@
"metadata": {},
"outputs": [],
"source": [
"# grids = ws.atm_field[\"t\"].data.grids\n",
"# z = grids[0]\n",
"# ice_mass_density_profile = np.zeros_like(z)\n",
"# ice_mass_density_profile[(z > 10e3) * (z < 15e3)] = 1e-4\n",
"# ice_mass_density_profile = pyarts.arts.GriddedField3(data=ice_mass_density_profile[..., None, None], grids=grids)\n",
"# ws.atm_field[ice_mass_density] = ice_mass_density_profile"
"grids = ws.atmospheric_field[\"t\"].data.grids\n",
"z = grids[0]\n",
"ice_extinction_profile = np.zeros_like(z)\n",
"ice_extinction_profile[(z > 5e3) * (z < 15e3)] = 1.0\n",
"ice_extinction_profile = pyarts.arts.GriddedField3(data=ice_extinction_profile[..., None, None], grids=grids)\n",
"ws.atmospheric_field[ice_extinction] = ice_extinction_profile\n",
"\n",
"ice_ssa_profile = np.zeros_like(z)\n",
"ice_ssa_profile[(z > 5e3) * (z < 15e3)] = 0.5\n",
"ice_ssa_profile = pyarts.arts.GriddedField3(data=ice_ssa_profile[..., None, None], grids=grids)\n",
"ws.atmospheric_field[ice_ssa] = ice_ssa_profile"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9fa4cdf8-3e41-4120-a9db-c95f3973e4e5",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "1523f263-1c70-4c30-b6d4-8c5f8a2ae68b",
"metadata": {},
"source": [
"### Rain"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a94b899d-f78e-4772-b46f-53378b17f276",
"metadata": {},
"outputs": [],
"source": [
"rain_extinction = pyarts.arts.ScatteringSpeciesProperty(\"rain\", pyarts.arts.ParticulateProperty(\"Extinction\"))\n",
"rain_ssa = pyarts.arts.ScatteringSpeciesProperty(\"rain\", pyarts.arts.ParticulateProperty(\"SingleScatteringAlbedo\"))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bee6b305-76f3-4ecc-9a3f-bc7642384c8d",
"metadata": {},
"outputs": [],
"source": [
"rain_extinction_profile = np.zeros_like(z)\n",
"rain_extinction_profile[z < 5e3] = 1.0\n",
"rain_extinction_profile = pyarts.arts.GriddedField3(data=rain_extinction_profile[..., None, None], grids=grids)\n",
"ws.atmospheric_field[rain_extinction] = rain_extinction_profile\n",
"\n",
"rain_ssa_profile = np.zeros_like(z)\n",
"rain_ssa_profile[z < 5e3] = 0.5\n",
"rain_ssa_profile = pyarts.arts.GriddedField3(data=rain_ssa_profile[..., None, None], grids=grids)\n",
"ws.atmospheric_field[rain_ssa] = rain_ssa_profile"
]
},
{
"cell_type": "markdown",
"id": "056c6b3c-7cab-4d61-bb24-0764ed1e86cc",
"metadata": {},
"source": [
"## Create the scattering species"
]
},
{
Expand All @@ -91,65 +154,120 @@
"metadata": {},
"outputs": [],
"source": [
"# import matplotlib.pyplot as plt\n",
"# plt.plot(ice_mass_density_profile[:, 0, 0] * 1e3, z/ 1e3)\n",
"# plt.ylim([0, 40])\n",
"# plt.ylabel(\"Altitude [km]\")\n",
"# plt.xlabel(\"Ice water content [g m$^{-3}$]\")"
"from pyarts.arts import HenyeyGreensteinScatterer, ArrayOfScatteringSpecies\n",
"\n",
"hg_ice = HenyeyGreensteinScatterer(ice_extinction, ice_ssa, 0.5)\n",
"hg_rain = HenyeyGreensteinScatterer(rain_extinction, rain_ssa, 0.0)\n",
"scattering_species = ArrayOfScatteringSpecies()\n",
"scattering_species.add(hg_ice)\n",
"scattering_species.add(hg_rain)"
]
},
{
"cell_type": "markdown",
"id": "d7342431-699a-40da-a183-c5ca3e03a6a1",
"metadata": {},
"source": [
"### Particle size distributions\n",
"\n"
"## Extracting scattering properties\n",
"\n",
"We can now extract bulk scattering properties from the scattering species."
]
},
{
"cell_type": "markdown",
"id": "3093e2be-366b-4cd4-9a35-77392a285bcb",
"metadata": {},
"source": [
"### Ice\n",
"\n",
"We can extract the ice phase function from the combind scattering species by calculating the bulk-scattering properties at an altitude above 5 km. To check the consistency of the Henyey-Greenstein scatterer, we extract the data in gridded and spectral representation and ensure that they are the same when both are converted to gridded representation."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "43d853f8-62f7-4a2f-a740-9a52cbf39bb3",
"metadata": {},
"outputs": [],
"source": [
"atm_pt = ws.atmospheric_field.at(11e3, 0.0, 0.0)\n",
"f_grid = np.array([89e9])\n",
"bsp = scattering_species.get_bulk_scattering_properties_tro_spectral(atm_pt, f_grid, 32, 1)\n",
"pm_spectral = bsp.phase_matrix"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "034fc429-f63c-4586-86e6-b46f5d67c9ad",
"id": "6710e1fc-28e1-43e9-b074-8cbbef6b2b79",
"metadata": {},
"outputs": [],
"source": [
"# psd_abel12 = pyarts.arts.MGDSingleMoment(ice_mass_density, 0.22, 2.2, 0.0, 1.0, 180.0, 270.0, False)\n",
"# psd_wang16 = pyarts.arts.MGDSingleMoment(ice_mass_density, \"Wang16\", 180.0, 270.0, False)\n",
"# psd_field19 = pyarts.arts.MGDSingleMoment(ice_mass_density, \"Field19\", 180.0, 270.0, False)"
"za_scat_grid = pm_spectral.to_gridded().get_za_scat_grid()\n",
"bsp = scattering_species.get_bulk_scattering_properties_tro_gridded(atm_pt, f_grid, za_scat_grid, 1)\n",
"pm_gridded = bsp.phase_matrix"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b9271a49-a0a4-4a13-8567-e7ab42ca61b3",
"id": "e3da481e-1ce8-43ae-96f6-f20192748afb",
"metadata": {},
"outputs": [],
"source": [
"# atm_point = ws.atm_field.at([12e3], [0.0], [0.0])\n",
"# sizes = np.logspace(-6, -2, 41)\n",
"# scat_species_a = 479.9830 \n",
"# scat_species_b = 3.0\n",
"import matplotlib.pyplot as plt\n",
"plt.plot(pm_spectral.to_gridded().flatten(), label=\"Converted from Spectral\")\n",
"plt.plot(pm_gridded.flatten(), ls=\"--\", label=\"Gridded\")\n",
"plt.title(\"Henyey-Greenstein phase function\")\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"id": "e389b4b1-c997-4e50-9d46-d695753e2d1f",
"metadata": {},
"source": [
"### Rain\n",
"\n",
"# y_abel12 = psd_abel12.evaluate(atm_point[0], sizes, scat_species_a, scat_species_b)\n",
"# y_wang16 = psd_wang16.evaluate(atm_point[0], sizes, scat_species_a, scat_species_b)\n",
"# y_field19 = psd_field19.evaluate(atm_point[0], sizes, scat_species_a, scat_species_b)"
"Similarly, we can extract the phase function for rain by calculating the bulk-scattering properties at a lower point in the atmosphere."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0ac9d751-328e-4b00-b611-e56cef020978",
"metadata": {},
"outputs": [],
"source": [
"atm_pt = ws.atmospheric_field.at(4e3, 0.0, 0.0)\n",
"bsp = scattering_species.get_bulk_scattering_properties_tro_spectral(atm_pt, f_grid, 32, 1)\n",
"pm_spectral = bsp.phase_matrix"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b345c993-a997-4925-a3d7-2ad5c1636997",
"metadata": {},
"outputs": [],
"source": [
"za_scat_grid = pm_spectral.to_gridded().get_za_scat_grid()\n",
"bsp = scattering_species.get_bulk_scattering_properties_tro_gridded(atm_pt, f_grid, za_scat_grid, 1)\n",
"pm_gridded = bsp.phase_matrix"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8ba10eb6-dc0a-4bea-9070-3d2a6d4be1b5",
"id": "04a485a5-cd07-4b46-abf0-7db349073c8f",
"metadata": {},
"outputs": [],
"source": [
"# plt.plot(sizes, y_abel12, label=\"Abel 12\")\n",
"# plt.plot(sizes, y_wang16, label=\"Wang 16\")\n",
"# plt.plot(sizes, y_field19, label=\"Field 19\")\n",
"# plt.xscale(\"log\")\n",
"# plt.yscale(\"log\")\n",
"# plt.xlabel(\"Particle size [m]\")\n",
"# plt.ylabel(\"Particle number density [m$^{-4}$]\")"
"import matplotlib.pyplot as plt\n",
"plt.plot(pm_spectral.to_gridded().flatten(), label=\"Converted from Spectral\")\n",
"plt.plot(pm_gridded.flatten(), ls=\"--\", label=\"Gridded\")\n",
"plt.title(\"Henyey-Greenstein phase function\")\n",
"plt.legend()"
]
}
],
Expand All @@ -169,7 +287,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.6"
"version": "3.11.9"
}
},
"nbformat": 4,
Expand Down
2 changes: 2 additions & 0 deletions src/core/options/arts_options.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1070,6 +1070,8 @@ radiation).
Value{"NumberDensity", "n", "Number density in m^{-3}"},
Value{"DMax", "dmax", "Maximum particle diameter in m."},
Value{"DVeq", "dveq", "Volume-equivalent diameter in m."},
Value{"Extinction", "ext", "Extinction in m^{-1}"},
Value{"SingleScatteringAlbedo", "ssa", "Single scattering albedo"},
Value{"ShapeParameter",
"ShapeParameter",
"PSD shape parmeter in arbitary units."},
Expand Down
9 changes: 8 additions & 1 deletion src/core/scattering/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ add_library(scattering STATIC
psd.cc
sht.cc
phase_matrix.cc
henyey_greenstein.cc
)

target_include_directories(scattering PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
Expand All @@ -15,4 +16,10 @@ target_compile_definitions(scattering PUBLIC ARTS_NO_SHTNS)
else()
target_include_directories(scattering PRIVATE "${FFTW_INCLUDE_DIR}" "${CMAKE_SOURCE_DIR}/3rdparty/shtns")
target_link_libraries(scattering PRIVATE fftw3 shtns)
endif()
endif()

add_library(scattering_io STATIC
xml_io_scattering.cc
)
target_include_directories(scattering_io PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
target_link_libraries(scattering_io PUBLIC artsworkspace)
30 changes: 29 additions & 1 deletion src/core/scattering/absorption_vector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ class AbsorptionVectorData<Scalar, Format::TRO, repr, stokes_dim>
absorption::get_n_mat_elems(Format::TRO, stokes_dim);
using CoeffVector = Eigen::Matrix<Scalar, 1, n_stokes_coeffs>;

AbsorptionVectorData() {};

/** Create a new AbsorptionVectorData container.
*
* Creates a container to hold extinction matrix data for the
Expand Down Expand Up @@ -106,6 +108,25 @@ class AbsorptionVectorData<Scalar, Format::TRO, repr, stokes_dim>
{this->extent(0), this->extent(1)});
}

AbsorptionVectorData<Scalar, Format::ARO, repr, stokes_dim>
to_lab_frame(std::shared_ptr<const Vector> za_inc_grid) {
AbsorptionVectorData<Scalar, Format::ARO, repr, stokes_dim> av_new{t_grid_, f_grid_, za_inc_grid};
for (Index t_ind = 0; t_ind < t_grid_->size(); ++t_ind) {
for (Index f_ind = 0; f_ind < f_grid_->size(); ++f_ind) {
for (Index za_inc_ind = 0; za_inc_ind < za_inc_grid->size(); ++za_inc_ind) {
av_new(t_ind, f_ind, za_inc_ind, 0) = this->operator()(t_ind, f_ind, 0);
}
}
}
return av_new;
}

AbsorptionVectorData<Scalar, Format::TRO, Representation::Spectral, stokes_dim> to_spectral() {
AbsorptionVectorData<Scalar, Format::TRO, Representation::Spectral, stokes_dim> avd_new(t_grid_, f_grid_);
reinterpret_cast<matpack::matpack_data<Scalar, 3>&>(avd_new) = *this;
return avd_new;
}

AbsorptionVectorData regrid(const ScatteringDataGrids grids,
const RegridWeights weights) {
AbsorptionVectorData result(grids.t_grid, grids.f_grid);
Expand Down Expand Up @@ -158,6 +179,7 @@ class AbsorptionVectorData<Scalar, Format::ARO, repr, stokes_dim>
absorption::get_n_mat_elems(Format::ARO, stokes_dim);
using CoeffVector = Eigen::Matrix<Scalar, 1, n_stokes_coeffs>;

AbsorptionVectorData() {};
/** Create a new AbsorptionVectorData container.
*
* Creates a container to hold extinction matrix data for the
Expand Down Expand Up @@ -210,9 +232,15 @@ class AbsorptionVectorData<Scalar, Format::ARO, repr, stokes_dim>
{this->extent(0), this->extent(1), this->extent(2)});
}

AbsorptionVectorData<Scalar, Format::ARO, Representation::Spectral, stokes_dim> to_spectral() {
AbsorptionVectorData<Scalar, Format::ARO, Representation::Spectral, stokes_dim> avd_new(t_grid_, f_grid_, za_inc_grid_);
reinterpret_cast<matpack::matpack_data<Scalar, 4>&>(avd_new) = *this;
return avd_new;
}

AbsorptionVectorData regrid(const ScatteringDataGrids& grids,
const RegridWeights& weights) {
AbsorptionVectorData result(grids.t_grid, grids.f_grid, grids.za_inc_grid);
AbsorptionVectorData result(grids.t_grid, grids.f_grid, std::make_shared<Vector>(grid_vector(*grids.za_inc_grid)));
auto coeffs_this = get_coeff_vector_view();
auto coeffs_res = result.get_coeff_vector_view();
for (Index i_t = 0; i_t < static_cast<Index>(weights.t_grid_weights.size()); ++i_t) {
Expand Down
Loading

0 comments on commit 9402278

Please sign in to comment.