ARTS 2.5.4 (git: 4c0d3b4d)
predefined_absorption_models.cc
Go to the documentation of this file.
1/* Copyright (C) 2020
2 * Richard Larsson <ric.larsson@gmail.com>
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2, or (at your option) any
7 * later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 * USA. */
18
27
28#include <Faddeeva/Faddeeva.hh>
29#include <algorithm>
30#include <iomanip>
31
32#include "debug.h"
33#include "jacobian.h"
34#include "lin_alg.h"
35#include "linescaling.h"
36#include "matpack.h"
37#include "matpackI.h"
38#include "predefined/predef.h"
39#include "propagationmatrix.h"
40#include "quantum_numbers.h"
41#include "species.h"
42
60template <bool check_exist>
61bool compute_selection(PropagationMatrix& pm [[maybe_unused]],
62 const SpeciesIsotopeRecord& model,
63 const Vector& f [[maybe_unused]],
64 const Numeric& p [[maybe_unused]],
65 const Numeric& t [[maybe_unused]],
66 const VMRS& vmr [[maybe_unused]]) {
67 switch (Species::find_species_index(model)) {
68 case find_species_index(Species::Species::Oxygen, "MPM2020"):
69 if constexpr (not check_exist) MPM2020::compute(pm, f, p, t, vmr.O2);
70 return true;
71 case find_species_index(Species::Species::Water, "ForeignContCKDMT350"):
72 if constexpr (not check_exist)
73 CKDMT350::compute_foreign_h2o(pm, f, p, t, vmr.H2O);
74 return true;
75 case find_species_index(Species::Species::Water, "SelfContCKDMT350"):
76 if constexpr (not check_exist)
77 CKDMT350::compute_self_h2o(pm, f, p, t, vmr.H2O);
78 return true;
79 }
80 return false;
81}
82
89template <bool special>
90constexpr Numeric dvmr_calc(Numeric x) noexcept {
91 if constexpr (special) {
92 return x;
93 } else {
94 constexpr Numeric d = 1e-6;
95 constexpr Numeric l = d * 1e-4;
96 return x < l ? d : x * d;
97 }
98}
99
117template <bool special>
119 const PropagationMatrix& pm,
120 const SpeciesIsotopeRecord& model,
121 const Vector& f,
122 const Numeric& p,
123 const Numeric& t,
124 VMRS vmr,
125 const Species::Species spec) {
126 Numeric dvmr = 0;
127
128 switch (spec) {
129 case Species::Species::Oxygen:
130 dvmr = dvmr_calc<special>(vmr.O2);
131 if constexpr (not special) vmr.O2 += dvmr;
132 break;
133 case Species::Species::Water:
134 dvmr = dvmr_calc<special>(vmr.H2O);
135 if constexpr (not special) vmr.H2O += dvmr;
136 break;
137 default:
138 return false; // Escape mechanism when nothing should be done
139 }
140
141 if constexpr (not special) {
142 dpm.SetZero();
143 compute_selection<false>(dpm, model, f, p, t, vmr);
144 dpm -= pm;
145 dpm /= dvmr;
146 } else {
147 if (dvmr not_eq 0) {
148 dpm = pm;
149 dpm /= dvmr;
150 } else {
151 compute_vmr_deriv<false>(dpm, pm, model, f, p, t, vmr, spec);
152 }
153 }
154 return true;
155}
156
157void compute(PropagationMatrix& propmat_clearsky,
158 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
159 const SpeciesIsotopeRecord& model,
160 const Vector& f_grid,
161 const Numeric& rtp_pressure,
162 const Numeric& rtp_temperature,
163 const VMRS& vmr,
164 const ArrayOfRetrievalQuantity& jacobian_quantities) {
165 if (not compute_selection<true>(
166 propmat_clearsky, model, f_grid, rtp_pressure, rtp_temperature, vmr))
167 return;
168
169 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
170 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
171 const bool do_vmrs_jac =
172 std::any_of(jacobian_quantities.begin(),
173 jacobian_quantities.end(),
174 [](auto& deriv) { return deriv == Jacobian::Line::VMR; }) or
175 std::any_of(jacobian_quantities.begin(),
176 jacobian_quantities.end(),
177 [model](auto& deriv) {
178 return deriv == Jacobian::Special::ArrayOfSpeciesTagVMR and
179 std::any_of(deriv.Target().species_array_id.begin(),
180 deriv.Target().species_array_id.end(),
181 [model](auto& tag) {
182 return tag.Isotopologue() == model;
183 });
184 });
185
186 if (do_freq_jac or do_temp_jac or do_vmrs_jac) {
189 f_grid.nelem()); // NOTE: Change if ever stokes_dim not_eq 1
191 f_grid.nelem()); // NOTE: Change if ever stokes_dim not_eq 1
192 compute_selection<false>(
193 pm, model, f_grid, rtp_pressure, rtp_temperature, vmr);
194
195 // Add absorption to the forward parameter
196 propmat_clearsky.Kjj() += pm.Kjj();
197
198 if (do_temp_jac) {
199 const Numeric d = temperature_perturbation(jacobian_quantities);
200 ARTS_ASSERT(d not_eq 0)
201
202 dpm.SetZero();
203 compute_selection<false>(
204 dpm, model, f_grid, rtp_pressure, rtp_temperature + d, vmr);
205 dpm -= pm;
206 dpm /= d;
207 for (Index iq = 0; iq < dpropmat_clearsky_dx.nelem(); iq++) {
208 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
209 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
210 }
211 }
212 }
213
214 if (do_freq_jac) {
215 const Numeric d = temperature_perturbation(jacobian_quantities);
216 ARTS_ASSERT(d not_eq 0)
217
218 Vector f_grid_d{f_grid};
219 f_grid_d += d;
220
221 dpm.SetZero();
222 compute_selection<false>(
223 dpm, model, f_grid_d, rtp_pressure, rtp_temperature, vmr);
224 dpm -= pm;
225 dpm /= d;
226 for (Index iq = 0; iq < dpropmat_clearsky_dx.nelem(); iq++) {
227 if (is_frequency_parameter(jacobian_quantities[iq])) {
228 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
229 }
230 }
231 }
232
233 for (Index iq = 0; iq < dpropmat_clearsky_dx.nelem(); iq++) {
234 auto& deriv = jacobian_quantities[iq];
235 if (deriv == Jacobian::Line::VMR) {
236 if (compute_vmr_deriv<false>(dpm,
237 pm,
238 model,
239 f_grid,
240 rtp_pressure,
241 rtp_temperature,
242 vmr,
243 deriv.QuantumIdentity().Species()))
244 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
245 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR and
246 std::any_of(deriv.Target().species_array_id.begin(),
247 deriv.Target().species_array_id.end(),
248 [model](auto& tag) {
249 return tag.Isotopologue() == model;
250 })) {
251 if (compute_vmr_deriv<true>(
252 dpm,
253 pm,
254 model,
255 f_grid,
256 rtp_pressure,
257 rtp_temperature,
258 vmr,
259 deriv.Target().species_array_id.front().Spec()))
260 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
261 }
262 }
263 } else {
264 compute_selection<false>(
265 propmat_clearsky, model, f_grid, rtp_pressure, rtp_temperature, vmr);
266 }
267}
268} // namespace Absorption::PredefinedModel
This can be used to make arrays out of anything.
Definition: array.h:108
void SetZero()
Sets all data to zero.
The Vector class.
Definition: matpackI.h:908
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1166
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1145
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1001
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1174
Routines for setting up the jacobian.
Linear algebra functions.
Constains various line scaling functions.
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
constexpr Numeric dvmr_calc(Numeric x) noexcept
Sets a VMR perturbation.
bool compute_vmr_deriv(PropagationMatrix &dpm, const PropagationMatrix &pm, const SpeciesIsotopeRecord &model, const Vector &f, const Numeric &p, const Numeric &t, VMRS vmr, const Species::Species spec)
Compute the partial VMR derivative.
bool compute_selection(PropagationMatrix &pm, const SpeciesIsotopeRecord &model, const Vector &f, const Numeric &p, const Numeric &t, const VMRS &vmr)
Compute the selected model and returns if it can be computed.
void compute(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const SpeciesIsotopeRecord &model, const Vector &f_grid, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const VMRS &vmr, const ArrayOfRetrievalQuantity &jacobian_quantities)
Compute the predefined model.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ArrayOfSpeciesTagVMR
Definition: jacobian.h:101
Temperature
Definition: jacobian.h:58
constexpr Index find_species_index(const Species spec, const std::string_view isot) noexcept
Stuff related to the propagation matrix.
#define d
Contains known required VMR values.
Struct containing all information needed about one isotope.
Definition: isotopologues.h:16