From 9e7b6eb4f9ebb5101e6d21347484779958325557 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Tue, 10 Feb 2026 18:45:46 +0900 Subject: [PATCH 1/2] Implement MTCKD4.3 --- .../predefined_absorption_models.cc | 16 + src/core/predefined/CMakeLists.txt | 3 +- src/core/predefined/MT_CKD430.cc | 335 ++++++++++++++++++ src/core/predefined/no_aer.cc | 10 + src/core/predefined/predef.h | 17 + src/core/predefined/predef_data.cc | 62 ++-- src/core/predefined/predef_data.h | 80 ++++- src/core/spec/isotopologues.h | 2 + src/m_predefined_absorption_models.cc | 44 +++ src/python_interface/py_predefined.cpp | 108 ++++++ src/workspace_methods.cpp | 65 +++- tests/aer/ckdmt430.nolgpl.py | 28 ++ 12 files changed, 718 insertions(+), 52 deletions(-) create mode 100644 src/core/predefined/MT_CKD430.cc create mode 100644 tests/aer/ckdmt430.nolgpl.py diff --git a/src/core/absorption/predefined_absorption_models.cc b/src/core/absorption/predefined_absorption_models.cc index ad9bc59eb1..b0dfa5cbc8 100644 --- a/src/core/absorption/predefined_absorption_models.cc +++ b/src/core/absorption/predefined_absorption_models.cc @@ -59,6 +59,22 @@ bool compute_selection( atm_point, std::get(predefined_model_data.data)); return true; + case "H2O-ForeignContCKDMT430"_isot_index: + if constexpr (not check_exist) + MT_CKD430::compute_foreign_h2o( + pm, + f, + atm_point, + std::get(predefined_model_data.data)); + return true; + case "H2O-SelfContCKDMT430"_isot_index: + if constexpr (not check_exist) + MT_CKD430::compute_self_h2o( + pm, + f, + atm_point, + std::get(predefined_model_data.data)); + return true; case "O2-MPM2020"_isot_index: if constexpr (not check_exist) MPM2020::compute(pm, f, atm_point); return true; diff --git a/src/core/predefined/CMakeLists.txt b/src/core/predefined/CMakeLists.txt index 9df721013c..0d6baf8475 100644 --- a/src/core/predefined/CMakeLists.txt +++ b/src/core/predefined/CMakeLists.txt @@ -16,7 +16,8 @@ if (NOT ENABLE_ARTS_LGPL) MT_CKD252.cc CKDMT350.cc CKDMT320.cc - MT_CKD400.cc) + MT_CKD400.cc + MT_CKD430.cc) target_link_libraries(predef_aer PUBLIC matpack rtepack species_tags atm) target_include_directories(predef_aer PRIVATE ${ARTS_SOURCE_DIR}/3rdparty/) target_include_directories(predef_aer PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/predefined/MT_CKD430.cc b/src/core/predefined/MT_CKD430.cc new file mode 100644 index 0000000000..91705724f6 --- /dev/null +++ b/src/core/predefined/MT_CKD430.cc @@ -0,0 +1,335 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "predef.h" + +/*! Ported from Fortran90 code that contains the following statement: + +! -------------------------------------------------------------------------- +! | Copyright ©, Atmospheric and Environmental Research, Inc., 2022 | +! | | +! | All rights reserved. This source code was developed as part of the | +! | LBLRTM software and is designed for scientific and research purposes. | +! | Atmospheric and Environmental Research Inc. (AER) grants USER the right | +! | to download, install, use and copy this software for scientific and | +! | research purposes only. This software may be redistributed as long as | +! | this copyright notice is reproduced on any copy made and appropriate | +! | acknowledgment is given to AER. This software or any modified version | +! | of this software may not be incorporated into proprietary software or | +! | commercial software offered for sale without the express written | +! | consent of AER. | +! | | +! | This software is provided as is without any express or implied | +! | warranties. | +! -------------------------------------------------------------------------- +*/ + +namespace Absorption::PredefinedModel::MT_CKD430 { +namespace { +Numeric RADFN_FUN(const Numeric XVI, const Numeric XKT) noexcept { + // ---------------------------------------------------------------------- B18060 + // LAST MODIFICATION: 12 AUGUST 1991 B17940 + // B17950 + // IMPLEMENTATION: R.D. WORSHAM B17960 + // B17970 + // ALGORITHM REVISIONS: S.A. CLOUGH B17980 + // R.D. WORSHAM B17990 + // J.L. MONCET B18000 + // B18010 + // B18020 + // ATMOSPHERIC AND ENVIRONMENTAL RESEARCH INC. B18030 + // 840 MEMORIAL DRIVE, CAMBRIDGE, MA 02139 B18040 + // B18050 + // B18070 + // WORK SUPPORTED BY: THE ARM PROGRAM B18080 + // OFFICE OF ENERGY RESEARCH B18090 + // DEPARTMENT OF ENERGY B18100 + // B18110 + // B18120 + // SOURCE OF ORIGINAL ROUTINE: AFGL LINE-BY-LINE MODEL B18130 + // B18140 + // FASCOD3 B18150 + // B18160 + // ---------------------------------------------------------------------- B18060 + // B18170 + // IN THE SMALL XVIOKT REGION 0.5 IS REQUIRED + + if (XKT > 0.0) { + const Numeric XVIOKT = XVI / XKT; + + if (XVIOKT <= 0.01) return 0.5 * XVIOKT * XVI; + + if (XVIOKT <= 10) { + Numeric EXPVKT = std::expm1(-XVIOKT); + return -XVI * EXPVKT / (2 + EXPVKT); + } + + return XVI; + } + return XVI; +} + +/*! Single XINT function to compute the interpolation of some data to a point X + + @param[in] P (X1 - X) / (X1 - X0) if A is on X0, X1, X2, X3 and X is the target of the interpolation + @param[in] A The value of the original function at four positions (at X0, X1, X2, X3) +*/ +constexpr Numeric XINT_FUN(const Numeric P, + const std::array& A) noexcept { + const Numeric C = (3 - 2 * P) * P * P; + const Numeric B = 0.5 * P * (1 - P); + const Numeric B1 = B * (1 - P); + const Numeric B2 = B * P; + + return -A[0] * B1 + A[1] * (1 - C + B2) + A[2] * (C + B1) - A[3] * B2; +} + +void check(const WaterData& data) { + bool c = data.self_absco_ref.size() == 0 or data.for_absco_ref.size() == 0 or + data.wavenumbers.size() == 0 or data.self_texp.size() == 0 or + data.for_closure_absco_ref.size() == 0; + ARTS_USER_ERROR_IF(c, "No data") +} +} // namespace + +void compute_foreign_closure_h2o(PropmatVector& propmat_clearsky, + const Vector& f_grid, + const AtmPoint& atm_point, + const WaterData& data) { + using Conversion::freq2kaycm; + + const Numeric P = atm_point.pressure; + const Numeric T = atm_point.temperature; + const Numeric vmrh2o = atm_point["H2O"_spec]; + + // Perform checks to ensure the calculation data is good + check(data); + + const Size n = f_grid.size(); + if (n == 0) return; + if (freq2kaycm(f_grid[0]) > data.wavenumbers.back()) return; + + // Constants + constexpr Numeric RADCN2 = 1.4387752; + + // Data constants + const Numeric last_wavenumber = data.wavenumbers.back(); + const Numeric dvc = data.wavenumbers[1] - data.wavenumbers[0]; + const Numeric recdvc = 1 / dvc; + const auto data_size = Index(data.wavenumbers.size()); + + // Level constants + const Numeric P0 = Conversion::bar2pa(1e-3 * data.ref_press); + const Numeric T0 = data.ref_temp; + const Numeric xkt = T / RADCN2; + const Numeric rho_rat = (P / P0) * (T0 / T); + const Numeric num_den_cm2 = 1e-6 * vmrh2o * P / (Constant::k * T); + + // Scaling logic + auto scl = [vmrh2o, rho_rat, xkt](auto& input_value, auto& input_wavenumber) { + return input_value * (1.0 - vmrh2o) * rho_rat * + RADFN_FUN(input_wavenumber, xkt); + }; + + // Compute data + Index cur = std::distance(data.wavenumbers.begin(), + std::lower_bound(data.wavenumbers.begin(), + data.wavenumbers.end(), + freq2kaycm(f_grid[0]) - 2 * dvc)); + const Numeric* v = data.wavenumbers.data_handle() + cur; + const Numeric* y = data.for_closure_absco_ref.data_handle() + cur; + std::array k{0, 0, 0, 0}; + + //! Follow the Fortran idea (and old MT CKD implementations in ARTS to mirror the zero-frequency values) + for (auto i = -1; i < 3 and cur + i < data_size; i++) { + if (i < 0 and cur == 0) + k[i + 1] = scl(*(y + i + 2), *(v + i + 2)); + else + k[i + 1] = scl(*(y + i), *(v + i)); + } + + // Compute loop + for (Size s = 0; s < n; ++s) { + if (f_grid[s] < 0) continue; + auto x = freq2kaycm(f_grid[s]); + if (x > last_wavenumber) return; + + while (x > *(v + 1)) { + std::shift_left(k.begin(), k.end(), 1); + + k.back() = data_size > cur + 3 ? scl(*(y + 3), *(v + 3)) : 0; + + cur++; + v++; + y++; + } + + auto out = 1e2 * num_den_cm2 * XINT_FUN(recdvc * (x - *v), k); + propmat_clearsky[s].A() += out >= 0 ? out : 0; + } +} + +void compute_foreign_h2o(PropmatVector& propmat_clearsky, + const Vector& f_grid, + const AtmPoint& atm_point, + const WaterData& data) { + using Conversion::freq2kaycm; + + const Numeric P = atm_point.pressure; + const Numeric T = atm_point.temperature; + const Numeric vmrh2o = atm_point["H2O"_spec]; + + // Perform checks to ensure the calculation data is good + check(data); + + const Size n = f_grid.size(); + if (n == 0) return; + if (freq2kaycm(f_grid[0]) > data.wavenumbers.back()) return; + + // Constants + constexpr Numeric RADCN2 = 1.4387752; + + // Data constants + const Numeric last_wavenumber = data.wavenumbers.back(); + const Numeric dvc = data.wavenumbers[1] - data.wavenumbers[0]; + const Numeric recdvc = 1 / dvc; + const auto data_size = Index(data.wavenumbers.size()); + + // Level constants + const Numeric P0 = Conversion::bar2pa(1e-3 * data.ref_press); + const Numeric T0 = data.ref_temp; + const Numeric xkt = T / RADCN2; + const Numeric rho_rat = (P / P0) * (T0 / T); + const Numeric num_den_cm2 = 1e-6 * vmrh2o * P / (Constant::k * T); + + // Scaling logic + auto scl = [vmrh2o, rho_rat, xkt](auto& input_value, auto& input_wavenumber) { + return input_value * (1.0 - vmrh2o) * rho_rat * + RADFN_FUN(input_wavenumber, xkt); + }; + + // Compute data + Index cur = std::distance(data.wavenumbers.begin(), + std::lower_bound(data.wavenumbers.begin(), + data.wavenumbers.end(), + freq2kaycm(f_grid[0]) - 2 * dvc)); + const Numeric* v = data.wavenumbers.data_handle() + cur; + const Numeric* y = data.for_absco_ref.data_handle() + cur; + std::array k{0, 0, 0, 0}; + + //! Follow the Fortran idea (and old MT CKD implementations in ARTS to mirror the zero-frequency values) + for (auto i = -1; i < 3 and cur + i < data_size; i++) { + if (i < 0 and cur == 0) + k[i + 1] = scl(*(y + i + 2), *(v + i + 2)); + else + k[i + 1] = scl(*(y + i), *(v + i)); + } + + // Compute loop + for (Size s = 0; s < n; ++s) { + if (f_grid[s] < 0) continue; + auto x = freq2kaycm(f_grid[s]); + if (x > last_wavenumber) return; + + while (x > *(v + 1)) { + std::shift_left(k.begin(), k.end(), 1); + + k.back() = data_size > cur + 3 ? scl(*(y + 3), *(v + 3)) : 0; + + cur++; + v++; + y++; + } + + auto out = 1e2 * num_den_cm2 * XINT_FUN(recdvc * (x - *v), k); + propmat_clearsky[s].A() += out >= 0 ? out : 0; + } +} + +void compute_self_h2o(PropmatVector& propmat_clearsky, + const Vector& f_grid, + const AtmPoint& atm_point, + const WaterData& data) { + using Conversion::freq2kaycm; + + const Numeric P = atm_point.pressure; + const Numeric T = atm_point.temperature; + const Numeric vmrh2o = atm_point["H2O"_spec]; + + // Perform checks to ensure the calculation data is good + check(data); + + const Size n = f_grid.size(); + if (n == 0) return; + if (freq2kaycm(f_grid[0]) > data.wavenumbers.back()) return; + + // Constants + constexpr Numeric RADCN2 = 1.4387752; + + // Data constants + const Numeric last_wavenumber = data.wavenumbers.back(); + const Numeric dvc = data.wavenumbers[1] - data.wavenumbers[0]; + const Numeric recdvc = 1 / dvc; + const auto data_size = Index(data.wavenumbers.size()); + + // Level constants + const Numeric P0 = Conversion::bar2pa(1e-3 * data.ref_press); + const Numeric T0 = data.ref_temp; + const Numeric xkt = T / RADCN2; + const Numeric rho_rat = (P / P0) * (T0 / T); + const Numeric num_den_cm2 = 1e-6 * vmrh2o * P / (Constant::k * T); + + // Scaling logic + auto scl = [vmrh2o, rho_rat, xkt, r = T0 / T](auto& input_value, + auto& input_exponent, + auto& input_wavenumber) { + return input_value * vmrh2o * rho_rat * std::pow(r, input_exponent) * + RADFN_FUN(input_wavenumber, xkt); + }; + + // Compute data + Index cur = std::distance(data.wavenumbers.begin(), + std::lower_bound(data.wavenumbers.begin(), + data.wavenumbers.end(), + freq2kaycm(f_grid[0]) - 2 * dvc)); + const Numeric* v = data.wavenumbers.data_handle() + cur; + const Numeric* y = data.self_absco_ref.data_handle() + cur; + const Numeric* e = data.self_texp.data_handle() + cur; + std::array k{0, 0, 0, 0}; + + //! Follow the Fortran idea (and old MT CKD implementations in ARTS to mirror the zero-frequency values) + for (auto i = -1; i < 3 and cur + i < data_size; i++) { + if (i < 0 and cur == 0) + k[i + 1] = scl(*(y + i + 2), *(e + i + 2), *(v + i + 2)); + else + k[i + 1] = scl(*(y + i), *(e + i), *(v + i)); + } + + // Compute loop + for (Size s = 0; s < n; ++s) { + if (f_grid[s] < 0) continue; + auto x = freq2kaycm(f_grid[s]); + if (x > last_wavenumber) return; + + while (x > *(v + 1)) { + *std::shift_left(k.begin(), k.end(), 1) = data_size > cur + 3 ? scl(*(y + 3), *(e + 3), *(v + 3)) : 0; + + cur++; + v++; + y++; + e++; + } + + auto out = 1e2 * num_den_cm2 * XINT_FUN(recdvc * (x - *v), k); + propmat_clearsky[s].A() += out >= 0 ? out : 0; + } +} +} // namespace Absorption::PredefinedModel::MT_CKD430 \ No newline at end of file diff --git a/src/core/predefined/no_aer.cc b/src/core/predefined/no_aer.cc index bbaca8dcb8..0cb08ae504 100644 --- a/src/core/predefined/no_aer.cc +++ b/src/core/predefined/no_aer.cc @@ -20,6 +20,16 @@ void compute_self_h2o(PropmatVector &, const AtmPoint &, const WaterData &) NOPE; } // namespace MT_CKD400 +namespace MT_CKD430 { +void compute_foreign_h2o(PropmatVector &, + const Vector &, + const AtmPoint &, + const WaterData &) NOPE; +void compute_self_h2o(PropmatVector &, + const Vector &, + const AtmPoint &, + const WaterData &) NOPE; +} // namespace MT_CKD430 namespace MT_CKD100 { void oxygen_cia(PropmatVector &, const Vector &, const AtmPoint &) NOPE; diff --git a/src/core/predefined/predef.h b/src/core/predefined/predef.h index 0830b2c128..f372dcddf5 100644 --- a/src/core/predefined/predef.h +++ b/src/core/predefined/predef.h @@ -150,4 +150,21 @@ void compute_self_h2o(PropmatVector& propmat_clearsky, const WaterData& data); } // namespace MT_CKD400 +namespace MT_CKD430 { +void compute_foreign_h2o(PropmatVector& propmat_clearsky, + const Vector& f_grid, + const AtmPoint& atm_point, + const WaterData& data); + +void compute_foreign_closure_h2o(PropmatVector& propmat_clearsky, + const Vector& f_grid, + const AtmPoint& atm_point, + const WaterData& data); + +void compute_self_h2o(PropmatVector& propmat_clearsky, + const Vector& f_grid, + const AtmPoint& atm_point, + const WaterData& data); +} // namespace MT_CKD430 + } // namespace Absorption::PredefinedModel diff --git a/src/core/predefined/predef_data.cc b/src/core/predefined/predef_data.cc index 6acee24b44..a3f8283dab 100644 --- a/src/core/predefined/predef_data.cc +++ b/src/core/predefined/predef_data.cc @@ -22,6 +22,21 @@ std::vector WaterData::sizes() const { }; } // namespace MT_CKD400 +namespace MT_CKD430 { +void WaterData::resize(const std::vector &inds) { + ARTS_USER_ERROR_IF(inds.size() not_eq 1, "Expects only one size") + self_absco_ref.resize(inds.front()); + for_absco_ref.resize(inds.front()); + for_closure_absco_ref.resize(inds.front()); + wavenumbers.resize(inds.front()); + self_texp.resize(inds.front()); +} + +std::vector WaterData::sizes() const { + return {static_cast(self_absco_ref.size())}; +}; +} // namespace MT_CKD430 + std::string_view model_name(const ModelVariant &data) { if (std::holds_alternative(data.data)) { return "ModelName"; @@ -31,6 +46,10 @@ std::string_view model_name(const ModelVariant &data) { return "MT_CKD400::WaterData"; } + if (std::holds_alternative(data.data)) { + return "MT_CKD430::WaterData"; + } + throw std::runtime_error("Unspecificed model type, this is a developer bug"); } @@ -43,50 +62,15 @@ ModelVariant model_data(const std::string_view name) { return ModelVariant{.data = MT_CKD400::WaterData{}}; } + if (name == "MT_CKD430::WaterData") { + return ModelVariant{.data = MT_CKD430::WaterData{}}; + } + throw std::runtime_error(std::format( R"(Unknown model name: "{}". Are all models defined?)", name)); } } // namespace Absorption::PredefinedModel -void xml_io_stream::write( - std::ostream &os, - const Absorption::PredefinedModel::MT_CKD400::WaterData &x, - bofstream *pbofs, - std::string_view name) { - XMLTag tag(type_name, "name", name); - tag.write_to_stream(os); - - xml_write_to_stream(os, x.ref_press, pbofs); - xml_write_to_stream(os, x.ref_temp, pbofs); - xml_write_to_stream(os, x.ref_h2o_vmr, pbofs); - xml_write_to_stream(os, x.self_absco_ref, pbofs); - xml_write_to_stream(os, x.for_absco_ref, pbofs); - xml_write_to_stream(os, x.wavenumbers, pbofs); - xml_write_to_stream(os, x.self_texp, pbofs); - - tag.write_to_end_stream(os); -} - -void xml_io_stream::read( - std::istream &is, - Absorption::PredefinedModel::MT_CKD400::WaterData &x, - bifstream *pbifs) { - XMLTag tag; - tag.read_from_stream(is); - tag.check_name(type_name); - - xml_read_from_stream(is, x.ref_press, pbifs); - xml_read_from_stream(is, x.ref_temp, pbifs); - xml_read_from_stream(is, x.ref_h2o_vmr, pbifs); - xml_read_from_stream(is, x.self_absco_ref, pbifs); - xml_read_from_stream(is, x.for_absco_ref, pbifs); - xml_read_from_stream(is, x.wavenumbers, pbifs); - xml_read_from_stream(is, x.self_texp, pbifs); - - tag.read_from_stream(is); - tag.check_end_name(type_name); -} - void xml_io_stream::write( std::ostream &os, const Absorption::PredefinedModel::ModelName &, diff --git a/src/core/predefined/predef_data.h b/src/core/predefined/predef_data.h index 532ba09ee3..bea7b06975 100644 --- a/src/core/predefined/predef_data.h +++ b/src/core/predefined/predef_data.h @@ -26,13 +26,29 @@ struct WaterData { }; } // namespace MT_CKD400 +namespace MT_CKD430 { +struct WaterData { + Numeric ref_press; + Numeric ref_temp; + Vector self_absco_ref; + Vector for_absco_ref; + Vector for_closure_absco_ref; + Vector wavenumbers; + Vector self_texp; + + void resize(const std::vector &); + [[nodiscard]] std::vector sizes() const; +}; +} // namespace MT_CKD430 + struct ModelName { static std::vector sizes() { return {}; }; static constexpr void resize(const std::vector &) {} }; struct ModelVariant { - using var_t = std::variant; + using var_t = + std::variant; var_t data; }; @@ -94,6 +110,41 @@ struct std::formatter { } }; +template <> +struct std::formatter { + format_tags tags; + + [[nodiscard]] constexpr auto &inner_fmt() { return *this; } + [[nodiscard]] constexpr auto &inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context &ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format( + const Absorption::PredefinedModel::MT_CKD430::WaterData &v, + FmtContext &ctx) const { + const std::string_view sep = tags.sep(); + tags.add_if_bracket(ctx, '['); + tags.format(ctx, + v.ref_temp, + sep, + v.ref_press, + sep, + v.self_absco_ref, + sep, + v.for_absco_ref, + sep, + v.for_closure_absco_ref, + sep, + v.self_texp); + tags.add_if_bracket(ctx, ']'); + return ctx.out(); + } +}; + template <> struct std::formatter { format_tags tags; @@ -133,17 +184,25 @@ struct std::formatter { }; template <> -struct xml_io_stream { - static constexpr std::string_view type_name = "PredefWaterData"sv; +struct xml_io_stream_name { + static constexpr std::string_view name = "PredefWaterDataMTCKD400"; +}; - static void write(std::ostream &os, - const Absorption::PredefinedModel::MT_CKD400::WaterData &x, - bofstream *pbofs = nullptr, - std::string_view name = ""sv); +template <> +struct xml_io_stream_aggregate< + Absorption::PredefinedModel::MT_CKD400::WaterData> { + static constexpr bool value = true; +}; - static void read(std::istream &is, - Absorption::PredefinedModel::MT_CKD400::WaterData &x, - bifstream *pbifs = nullptr); +template <> +struct xml_io_stream_name { + static constexpr std::string_view name = "PredefWaterDataMTCKD430"; +}; + +template <> +struct xml_io_stream_aggregate< + Absorption::PredefinedModel::MT_CKD430::WaterData> { + static constexpr bool value = true; }; template <> @@ -159,6 +218,7 @@ struct xml_io_stream { Absorption::PredefinedModel::ModelName &x, bifstream *pbifs = nullptr); }; + template <> struct xml_io_stream_name { static constexpr std::string_view name = "PredefinedModelData"sv; diff --git a/src/core/spec/isotopologues.h b/src/core/spec/isotopologues.h index 965d741efa..5a2d93d811 100644 --- a/src/core/spec/isotopologues.h +++ b/src/core/spec/isotopologues.h @@ -96,6 +96,7 @@ inline constexpr std::array Isotopologues{ Isotope("H2O"_spec, "ForeignContCKDMT320"), Isotope("H2O"_spec, "ForeignContCKDMT350"), Isotope("H2O"_spec, "ForeignContCKDMT400"), + Isotope("H2O"_spec, "ForeignContCKDMT430"), Isotope("H2O"_spec, "ForeignContStandardType"), Isotope("H2O"_spec, "MPM89"), Isotope("H2O"_spec, "PWR2021"), @@ -104,6 +105,7 @@ inline constexpr std::array Isotopologues{ Isotope("H2O"_spec, "SelfContCKDMT320"), Isotope("H2O"_spec, "SelfContCKDMT350"), Isotope("H2O"_spec, "SelfContCKDMT400"), + Isotope("H2O"_spec, "SelfContCKDMT430"), Isotope("H2O"_spec, "SelfContStandardType"), /** Water species **/ diff --git a/src/m_predefined_absorption_models.cc b/src/m_predefined_absorption_models.cc index 13158eca35..1f34b0ea09 100644 --- a/src/m_predefined_absorption_models.cc +++ b/src/m_predefined_absorption_models.cc @@ -108,6 +108,50 @@ void abs_predef_dataAddWaterMTCKD400(PredefinedModelData& abs_predef_data, abs_predef_data["H2O-SelfContCKDMT400"_isot].data = x; } +void abs_predef_dataAddWaterMTCKD430(PredefinedModelData& abs_predef_data, + const Numeric& ref_temp, + const Numeric& ref_press, + const Vector& self_absco_ref, + const Vector& for_absco_ref, + const Vector& for_closure_absco_ref, + const Vector& wavenumbers, + const Vector& self_texp) { + ARTS_TIME_REPORT + + const auto sz = self_absco_ref.size(); + + ARTS_USER_ERROR_IF( + sz not_eq for_absco_ref.size() or sz not_eq wavenumbers.size() or + sz not_eq self_texp.size() or sz not_eq for_closure_absco_ref.size(), + "Mismatching size, all vector inputs must be the same length") + ARTS_USER_ERROR_IF(sz < 4, "It makes no sense to have input shorter than 4") + ARTS_USER_ERROR_IF(not AscendingGrid::is_sorted(wavenumbers), + "The wavenumbers must be increasing in a regular manner") + + using Model = Absorption::PredefinedModel::MT_CKD430::WaterData; + Model x; + x.ref_temp = ref_temp; + x.ref_press = ref_press; + x.self_absco_ref.resize(sz); + x.for_absco_ref.resize(sz); + x.for_closure_absco_ref.resize(sz); + x.wavenumbers.resize(sz); + x.self_texp.resize(sz); + + std::copy( + self_absco_ref.begin(), self_absco_ref.end(), x.self_absco_ref.begin()); + std::copy( + for_absco_ref.begin(), for_absco_ref.end(), x.for_absco_ref.begin()); + std::copy(for_closure_absco_ref.begin(), + for_closure_absco_ref.end(), + x.for_closure_absco_ref.begin()); + std::copy(wavenumbers.begin(), wavenumbers.end(), x.wavenumbers.begin()); + std::copy(self_texp.begin(), self_texp.end(), x.self_texp.begin()); + + abs_predef_data["H2O-ForeignContCKDMT430"_isot].data = x; + abs_predef_data["H2O-SelfContCKDMT430"_isot].data = x; +} + /* Workspace method: Doxygen documentation will be auto-generated */ void spectral_propmatAddPredefined(PropmatVector& spectral_propmat, PropmatMatrix& spectral_propmat_jac, diff --git a/src/python_interface/py_predefined.cpp b/src/python_interface/py_predefined.cpp index bc963e3e09..414d2cfb55 100644 --- a/src/python_interface/py_predefined.cpp +++ b/src/python_interface/py_predefined.cpp @@ -15,6 +15,113 @@ #include "species_tags.h" namespace Python { +void internalCKDMT430(py::module_& m) { + py::class_ mm( + m, "MTCKD430WaterData"); + generic_interface(mm); + mm.def(py::init()) + .def_rw("ref_temp", + &Absorption::PredefinedModel::MT_CKD430::WaterData::ref_temp, + "Reference temperature\n\n.. :class:`Numeric`") + .def_rw("ref_press", + &Absorption::PredefinedModel::MT_CKD430::WaterData::ref_press, + "Reference pressure\n\n.. :class:`Numeric`") + .def_rw( + "self_absco_ref", + &Absorption::PredefinedModel::MT_CKD430::WaterData::self_absco_ref, + "Self absorption\n\n.. :class:`Vector`") + .def_rw("for_absco_ref", + &Absorption::PredefinedModel::MT_CKD430::WaterData::for_absco_ref, + "Foreign absorption\n\n.. :class:`Vector`") + .def_rw("for_closure_absco_ref", + &Absorption::PredefinedModel::MT_CKD430::WaterData::for_closure_absco_ref, + "Foreign absorption closure\n\n.. :class:`Vector`") + .def_rw("wavenumbers", + &Absorption::PredefinedModel::MT_CKD430::WaterData::wavenumbers, + "Wavenumbers\n\n.. :class:`Vector`") + .def_rw("self_texp", + &Absorption::PredefinedModel::MT_CKD430::WaterData::self_texp, + "Self temperature exponent\n\n.. :class:`Vector`") + + .doc() = "Water data representation for the MT CKD 4.3 model"; + + m.def( + "get_foreign_h2o_ckdmt430", + [](const Vector& f, + const AtmPoint& atm, + PredefinedModelData& data) -> Vector { + PropmatVector pm(f.size()); + Absorption::PredefinedModel::MT_CKD430::compute_foreign_h2o( + pm, + f, + atm, + std::get( + data.at("H2O-ForeignContCKDMT430"_isot).data)); + Vector out(pm.size()); + std::transform(pm.begin(), pm.end(), out.begin(), [](auto& prop) { + return prop.A(); + }); + return out; + }, + "f_grid"_a, + "atm"_a, + "predefined_model_data"_a, + R"--(Computes foreign absorption using MT CKD Hitran version 4.3 + +Parameters +---------- +f_grid : ~pyarts3.arts.Vector + Frequency grid [Hz] +atm : AtmPoint + The atmospheric state +predefined_model_data : ~pyarts3.arts.PredefinedModelData + As WSV + +Returns +------- +abs_coef : ~pyarts3.arts.Vector + Absorption coefficients +)--"); + + m.def( + "get_self_h2o_ckdmt430", + [](const Vector& f, + const AtmPoint& atm, + PredefinedModelData& data) -> Vector { + PropmatVector pm(f.size()); + Absorption::PredefinedModel::MT_CKD430::compute_self_h2o( + pm, + f, + atm, + std::get( + data.at("H2O-SelfContCKDMT430"_isot).data)); + Vector out(pm.size()); + std::transform(pm.begin(), pm.end(), out.begin(), [](auto& prop) { + return prop.A(); + }); + return out; + }, + "f_grid"_a, + "atm"_a, + "predefined_model_data"_a, + R"--(Computes self absorption using MT CKD Hitran version 4.3 + +Parameters +---------- +f_grid : ~pyarts3.arts.Vector + Frequency grid [Hz] +atm : AtmPoint + The atmospheric state +predefined_model_data : ~pyarts3.arts.PredefinedModelData + As WSV + +Returns +------- +abs_coef : ~pyarts3.arts.Vector + Absorption coefficients +)--"); +} + void internalCKDMT400(py::module_& m) { py::class_ mm( m, "MTCKD400WaterData"); @@ -1000,6 +1107,7 @@ spectral_propmat : PropmatVector internalCKDMT252(predef); internalCKDMT100(predef); internalCKDMT400(predef); + internalCKDMT430(predef); internalMPM89(predef); internalMPM93(predef); internalPWR98(predef); diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index 5c7126a615..97e5f183c2 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -1387,6 +1387,55 @@ This is based on the works cited here: https://hitran.org/mtckd/ R"--(Self temperature exponent [-])--"}, }; + wsm_data["abs_predef_dataAddWaterMTCKD430"] = { + .desc = R"--(Sets the data for MT CKD 4.3 Water model + +Note that the vectors must have the same length, and that wavenumbers must be growing +at a constant rate. The minimum length is 4. + +Note also that as this is predefined model data, the units of the values of the vectors +must be as described by each vector. + +This is based on the works cited here: https://hitran.org/mtckd/ + +.. note:: + The method itself is implemented from scratch. Using any version of + data after version 4.3 is supported by this method - all that changes + are the values of the vectors. +)--", + .author = {"Richard Larsson"}, + .out = {"abs_predef_data"}, + .in = {"abs_predef_data"}, + .gin = {"ref_temp", + "ref_press", + "self_absco_ref", + "for_absco_ref", + "for_closure_absco_ref", + "wavenumbers", + "self_texp"}, + .gin_type = {"Numeric", + "Numeric", + "Vector", + "Vector", + "Vector", + "Vector", + "Vector"}, + .gin_value = {std::nullopt, + std::nullopt, + std::nullopt, + std::nullopt, + std::nullopt, + std::nullopt, + std::nullopt}, + .gin_desc = {R"--(Reference temperature)--", + R"--(Reference pressure)--", + R"--(Self absorption [1/(cm-1 molecules/cm^2])--", + R"--(Foreign absorption [1/(cm-1 molecules/cm^2)])--", + R"--(Foreign closure absorption [1/(cm-1 molecules/cm^2)])--", + R"--(Wavenumbers [cm-1])--", + R"--(Self temperature exponent [-])--"}, + }; + wsm_data["abs_predef_dataInit"] = { .desc = R"--(Initialize the predefined model data )--", @@ -1524,7 +1573,13 @@ Available models * - ``H2O-ForeignContCKDMT400`` - MT CKD 4 foreign continua. |br| - Use water cutoff of 25 cm-1 and ``H2O-SelfContCKDMT350``. |br| + Use water cutoff of 25 cm-1 and ``H2O-SelfContCKDMT400``. |br| + Requires an external data source. + - :cite:t:`Mlawer2012` and :cite:p:`MTCKD` + + * - ``H2O-ForeignContCKDMT430`` + - MT CKD 4.3 foreign continua. |br| + Use water cutoff of 25 cm-1 and ``H2O-SelfContCKDMT430``. |br| Requires an external data source. - :cite:t:`Mlawer2012` and :cite:p:`MTCKD` @@ -1562,7 +1617,13 @@ Available models * - ``H2O-SelfContCKDMT400`` - MT CKD 4 self continua. |br| - Use water cutoff of 25 cm-1 and ``H2O-SelfContCKDMT350``. |br| + Use water cutoff of 25 cm-1 and ``H2O-ForeignContCKDMT400``. |br| + Requires an external data source. + - :cite:t:`Mlawer2012` and :cite:p:`MTCKD` + + * - ``H2O-SelfContCKDMT430`` + - MT CKD 4.3 self continua. |br| + Use water cutoff of 25 cm-1 and ``H2O-ForeignContCKDMT430``. |br| Requires an external data source. - :cite:t:`Mlawer2012` and :cite:p:`MTCKD` diff --git a/tests/aer/ckdmt430.nolgpl.py b/tests/aer/ckdmt430.nolgpl.py new file mode 100644 index 0000000000..058289d1bf --- /dev/null +++ b/tests/aer/ckdmt430.nolgpl.py @@ -0,0 +1,28 @@ +import pyarts3 as pyarts +import numpy as np + +f = pyarts.arts.convert.kaycm2freq(np.linspace(1, 21000, 101)) + +atm = pyarts.arts.AtmPoint() +atm.pressure = 1e4 +atm.temperature = 250 +atm["H2O"] = 1e-2 +atm["O2"] = 0.21 +atm["N2"] = 0.79 +atm["CO2"] = 400e-6 + +data = pyarts.arts.PredefinedModelData.fromcatalog( + "predef/", ["H2O-ForeignContCKDMT430", "H2O-SelfContCKDMT430"]) +self_abs430 = pyarts.arts.predef.get_self_h2o_ckdmt430(f, atm, data) + +self_ref430 = np.array([9.189769485196123e-08, 0.00015672443964545764, 2.828567174136117e-05, 6.82934630694487e-06, 2.31784464198272e-06, 9.338268580595634e-07, 1.0150564159771669e-06, 1.3946352148593543e-05, 3.35576519056182e-05, 3.155020968821078e-06, 3.3110859725448607e-07, 1.1828277073222903e-07, 7.970201899486387e-08, 7.040290750406673e-08, 6.865340996887612e-08, 2.641005338718162e-07, 3.5717897145869486e-07, 4.7076306017478115e-06, 3.331919198050135e-05, 3.0987459088601393e-06, 3.405220909158323e-07, 1.0449903360012003e-07, 5.669348482088639e-08, 3.976407475026572e-08, 1.258345493843198e-07, 2.1135362269250405e-06, 1.5449048918439632e-06, 1.1304034976872388e-07, 1.381234985540996e-08, 4.684258924732471e-09, 1.7726928396765135e-09, 2.8902352810231995e-09, 2.9352272376466614e-08, 2.041011824740944e-07, 1.11615891644693e-06, 1.3959529838216589e-06, 3.757626351318345e-08, 9.392836675790163e-09, 3.1335903767654105e-09, 1.891029307896446e-09, 6.334118887208509e-09, 1.541906176479611e-08, 1.434519916280068e-07, 1.7949581469070985e-08, 1.8374221406216735e-09, 3.034352571384538e-10, 1.256664878354076e-10, 1.7914391968122892e-10, + 4.695406649008256e-10, 9.644526379116789e-09, 3.991877107055666e-08, 5.577032679652794e-08, 9.776385876462862e-09, 2.681467780647672e-09, 2.4594165638143186e-10, 8.075010034164904e-11, 1.7158714232860604e-10, 8.453038131779653e-10, 5.964159820651672e-09, 7.736888970912498e-10, 4.5668251448705797e-10, 2.272329180160134e-11, 8.54905969452266e-12, 1.2467607396074576e-11, 9.415417117784018e-11, 1.2582833391642096e-09, 5.260224157224619e-09, 5.787001976490031e-10, 9.117492050126759e-10, 7.984672010826418e-11, 1.1636144604285995e-11, 5.73396009633933e-12, 4.1935292651322185e-11, 4.0174862858364136e-10, 7.64266593823464e-11, 5.246484283138208e-11, 1.2275152960302711e-11, 1.834850320760551e-12, 2.633074946638328e-12, 2.500272873850917e-11, 4.326785318794734e-10, 2.843254127172161e-10, 2.8687728193111663e-11, 7.728086311487542e-11, 2.919402391994982e-11, 2.844196687045464e-12, 3.475921372335562e-12, 4.452886673681736e-11, 5.0370099330552684e-11, 3.7923483221681875e-12, 5.051790656967964e-12, 4.302202552617765e-12, 1.2762346648918725e-12, 1.7784717281522425e-11, 2.3544525876029757e-10, 6.980840078950065e-11, 0, 0, 0, 0, 0]) + +assert np.allclose(self_abs430, self_ref430) + +foreign_abs430 = pyarts.arts.predef.get_foreign_h2o_ckdmt430(f, atm, data) + +foreign_ref430 = np.array([7.902266623188642e-08, 0.0006804823452410723, 2.0650827472669016e-05, 1.5494470651488873e-06, 2.5342231495509236e-07, 5.967529722919267e-08, 5.815802888579568e-07, 9.904785169274852e-05, 0.00019247804755044107, 6.549311134655636e-06, 1.675557700066427e-07, 7.84695662567049e-09, 6.105296568147661e-10, 1.7527915546953446e-08, 4.430398440459728e-08, 9.023142674259312e-07, 5.566073153956493e-07, 1.8528804208045185e-05, 0.0001495671313017644, 3.677807698290886e-06, 2.2059324732365718e-07, 4.897731096016647e-08, 2.3098324858751777e-08, 9.421673949704309e-09, 6.65688595136048e-08, 1.2655474801617041e-05, 9.250995620191857e-06, 7.700748348142169e-08, 6.3951274391930446e-09, 1.663946175343862e-09, 6.469483812292643e-10, 1.1793993788786671e-09, 1.5608740386529944e-07, 1.007438289565876e-06, 8.320868435168207e-06, 1.0990592003915888e-05, 8.995168782438871e-08, 4.517692598010068e-09, 8.131593496846035e-10, 5.629520619439137e-09, 2.9053744485723974e-08, 3.6156695835010906e-08, 7.560311192203557e-07, 3.9082335659897955e-08, 1.8146220815061607e-09, 8.108050525642898e-11, 5.1101029727081085e-11, 6.031839093432586e-10, + 5.284200635836903e-10, 3.689465532566e-08, 2.397017019526437e-07, 3.8074241944156685e-07, 4.186473530174785e-08, 3.996397860832967e-08, 3.491776191286313e-10, 4.476140038136932e-11, 8.27111519245472e-10, 2.346524342356196e-09, 3.085759166469117e-08, 1.6341807382578892e-09, 2.3556072044686797e-09, 5.1813725236032335e-11, 4.647791447739354e-12, 1.6596224053830763e-11, 1.0274165828068133e-10, 6.356948828449977e-09, 2.644954523789432e-08, 8.788315971000534e-10, 4.615461603953291e-09, 4.554902147533082e-10, 2.6296699683976694e-11, 5.463385386907692e-12, 1.0915201724324435e-10, 1.9647916317184522e-09, 1.8962542685691598e-10, 2.8742125364309394e-10, 1.630068788521999e-10, 1.8915259645837643e-12, 9.720265339143444e-13, 4.239783567472467e-11, 2.680059704535967e-09, 1.3445087334163242e-09, 4.781183249844617e-11, 4.086265358766565e-10, 7.517071282845343e-11, 9.347318522372578e-12, 3.2018861813102984e-12, 1.632918025252812e-10, 2.16510087407733e-10, 9.357561495983835e-12, 1.8271742567814066e-11, 9.875231659355349e-12, 1.6577468424120544e-12, 1.365083919554613e-11, 4.4075812656506807e-10, 5.577567615725717e-11, 0, 0, 0, 0, 0]) + +assert np.allclose(foreign_abs430, foreign_ref430) From 32dea6646379c851f36f3c45f05d4a8fbdc32797 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 11 Feb 2026 11:45:18 +0900 Subject: [PATCH 2/2] Also pull data --- testdata/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/testdata/CMakeLists.txt b/testdata/CMakeLists.txt index 63c7e1cc68..91682d6b04 100644 --- a/testdata/CMakeLists.txt +++ b/testdata/CMakeLists.txt @@ -12,6 +12,8 @@ else() lines/H2O-161.xml predef/H2O-SelfContCKDMT400.xml predef/H2O-ForeignContCKDMT400.xml + predef/H2O-SelfContCKDMT430.xml + predef/H2O-ForeignContCKDMT430.xml xsec/O3-XFIT.xml xsec/O3-XFIT.xml.bin )