ARTS 2.5.9 (git: 825fa5f2)
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#include <predef.h>
32
33#include "debug.h"
34#include "jacobian.h"
35#include "lin_alg.h"
36#include "linescaling.h"
37#include "matpack.h"
38#include "matpackI.h"
39#include "predefined/predef_data.h"
40#include "propagationmatrix.h"
41#include "quantum_numbers.h"
42#include "species.h"
43
61template <bool check_exist>
62bool compute_selection(PropagationMatrix& pm [[maybe_unused]],
63 const SpeciesIsotopeRecord& model,
64 const Vector& f [[maybe_unused]],
65 const Numeric& p [[maybe_unused]],
66 const Numeric& t [[maybe_unused]],
67 const VMRS& vmr [[maybe_unused]],
68 const PredefinedModelData& predefined_model_data
69 [[maybe_unused]]) {
70 switch (Species::find_species_index(model)) {
71 case find_species_index(Species::Species::Water, "ForeignContCKDMT400"):
72 if constexpr (not check_exist) MT_CKD400::compute_foreign_h2o(pm, f, p, t, vmr.H2O, predefined_model_data.get<MT_CKD400::WaterData>());
73 return true;
74 case find_species_index(Species::Species::Water, "SelfContCKDMT400"):
75 if constexpr (not check_exist) MT_CKD400::compute_self_h2o(pm, f, p, t, vmr.H2O, predefined_model_data.get<MT_CKD400::WaterData>());
76 return true;
77 case find_species_index(Species::Species::Oxygen, "MPM2020"):
78 if constexpr (not check_exist) MPM2020::compute(pm, f, p, t, vmr.O2);
79 return true;
80 case find_species_index(Species::Species::Oxygen, "PWR98"):
81 if constexpr (not check_exist) PWR98::oxygen(pm, f, p, t, vmr.O2, vmr.H2O);
82 return true;
83 case find_species_index(Species::Species::Oxygen, "TRE05"):
84 if constexpr (not check_exist) TRE05::oxygen(pm, f, p, t, vmr.O2, vmr.H2O);
85 return true;
86 case find_species_index(Species::Species::Water, "PWR98"):
87 if constexpr (not check_exist) PWR98::water(pm, f, p, t, vmr.H2O);
88 return true;
89 case find_species_index(Species::Species::Oxygen, "MPM89"):
90 if constexpr (not check_exist) MPM89::oxygen(pm, f, p, t, vmr.O2, vmr.H2O);
91 return true;
92 case find_species_index(Species::Species::Water, "MPM89"):
93 if constexpr (not check_exist) MPM89::water(pm, f, p, t, vmr.H2O);
94 return true;
95 case find_species_index(Species::Species::Water, "ForeignContCKDMT350"):
96 if constexpr (not check_exist) CKDMT350::compute_foreign_h2o(pm, f, p, t, vmr.H2O);
97 return true;
98 case find_species_index(Species::Species::Water, "SelfContCKDMT350"):
99 if constexpr (not check_exist) CKDMT350::compute_self_h2o(pm, f, p, t, vmr.H2O);
100 return true;
101 case find_species_index(Species::Species::Water, "ForeignContStandardType"):
102 if constexpr (not check_exist) Standard::water_foreign(pm, f, p, t, vmr.H2O);
103 return true;
104 case find_species_index(Species::Species::Water, "SelfContStandardType"):
105 if constexpr (not check_exist) Standard::water_self(pm, f, p, t, vmr.H2O);
106 return true;
107 case find_species_index(Species::Species::Oxygen, "SelfContStandardType"):
108 if constexpr (not check_exist) Standard::oxygen(pm, f, p, t, vmr.O2, vmr.H2O);
109 return true;
110 case find_species_index(Species::Species::Nitrogen, "SelfContStandardType"):
111 if constexpr (not check_exist) Standard::nitrogen(pm, f, p, t, vmr.N2);
112 return true;
113 case find_species_index(Species::Species::CarbonDioxide, "CKDMT252"):
114 if constexpr (not check_exist) MT_CKD252::carbon_dioxide(pm, f, p, t, vmr.CO2);
115 return true;
116 case find_species_index(Species::Species::Oxygen, "visCKDMT252"):
117 if constexpr (not check_exist) MT_CKD252::oxygen_vis(pm, f, p, t, vmr.O2);
118 return true;
119 case find_species_index(Species::Species::Nitrogen, "CIAfunCKDMT252"):
120 if constexpr (not check_exist) MT_CKD252::nitrogen_fun(pm, f, p, t, vmr.N2, vmr.H2O, vmr.O2);
121 return true;
122 case find_species_index(Species::Species::Nitrogen, "CIArotCKDMT252"):
123 if constexpr (not check_exist) MT_CKD252::nitrogen_rot(pm, f, p, t, vmr.N2, vmr.H2O, vmr.O2);
124 return true;
125 case find_species_index(Species::Species::Oxygen, "CIAfunCKDMT100"):
126 if constexpr (not check_exist) MT_CKD100::oxygen_cia(pm, f, p, t, vmr.O2);
127 return true;
128 case find_species_index(Species::Species::Oxygen, "v0v0CKDMT100"):
129 if constexpr (not check_exist) MT_CKD100::oxygen_v0v0(pm, f, p, t, vmr.O2, vmr.N2);
130 return true;
131 case find_species_index(Species::Species::Oxygen, "v1v0CKDMT100"):
132 if constexpr (not check_exist) MT_CKD100::oxygen_v0v1(pm, f, p, t, vmr.O2);
133 return true;
134 case find_species_index(Species::Species::liquidcloud, "ELL07"):
135 if constexpr (not check_exist) ELL07::compute(pm, f, t, vmr.LWC);
136 return true;
137 case find_species_index(Species::Species::FINAL, "Not A Model"): break;
138 }
139 return false;
140}
141
144 return compute_selection<true>(pm, model, {}, {}, {}, {}, {});
145}
146
153template <bool special>
154constexpr Numeric dvmr_calc(Numeric x) noexcept {
155 if constexpr (special) {
156 return x;
157 } else {
158 constexpr Numeric d = 1e-6;
159 constexpr Numeric l = d * 1e-4;
160 return x < l ? d : x * d;
161 }
162}
163
181template <bool special>
183 const PropagationMatrix& pm,
184 const SpeciesIsotopeRecord& model,
185 const Vector& f,
186 const Numeric& p,
187 const Numeric& t,
188 VMRS vmr,
189 const Species::Species spec,
190 const PredefinedModelData& predefined_model_data
191 [[maybe_unused]]) {
192 Numeric dvmr = 0;
193
194 switch (spec) {
195 case Species::Species::Oxygen:
196 dvmr = dvmr_calc<special>(vmr.O2);
197 if constexpr (not special) vmr.O2 += dvmr;
198 break;
199 case Species::Species::Water:
200 dvmr = dvmr_calc<special>(vmr.H2O);
201 if constexpr (not special) vmr.H2O += dvmr;
202 break;
203 case Species::Species::Nitrogen:
204 dvmr = dvmr_calc<special>(vmr.N2);
205 if constexpr (not special) vmr.N2 += dvmr;
206 break;
207 case Species::Species::CarbonDioxide:
208 dvmr = dvmr_calc<special>(vmr.CO2);
209 if constexpr (not special) vmr.CO2 += dvmr;
210 break;
211 case Species::Species::liquidcloud:
212 dvmr = dvmr_calc<special>(vmr.LWC);
213 if constexpr (not special) vmr.LWC += dvmr;
214 break;
215 default:
216 return false; // Escape mechanism when nothing should be done
217 }
218 static_assert(
219 sizeof(VMRS) / sizeof(Numeric) == 5,
220 "It seems you have changed VMRS. Please check that the derivatives are up-to-date above this assert");
221
222 if constexpr (not special) {
223 dpm.SetZero();
224 compute_selection<false>(dpm, model, f, p, t, vmr, predefined_model_data);
225 dpm -= pm;
226 dpm /= dvmr;
227 } else {
228 if (dvmr not_eq 0) {
229 dpm = pm;
230 dpm /= dvmr;
231 } else {
232 compute_vmr_deriv<false>(
233 dpm, pm, model, f, p, t, vmr, spec, predefined_model_data);
234 }
235 }
236 return true;
237}
238
239void compute(PropagationMatrix& propmat_clearsky,
240 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
241 const SpeciesIsotopeRecord& model,
242 const Vector& f_grid,
243 const Numeric& rtp_pressure,
244 const Numeric& rtp_temperature,
245 const VMRS& vmr,
246 const ArrayOfRetrievalQuantity& jacobian_quantities,
247 const PredefinedModelData& predefined_model_data) {
248 if (not compute_selection<true>(propmat_clearsky,
249 model,
250 f_grid,
251 rtp_pressure,
252 rtp_temperature,
253 vmr,
254 predefined_model_data))
255 return;
256
257 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
258 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
259 const bool do_vmrs_jac =
260 std::any_of(jacobian_quantities.begin(),
261 jacobian_quantities.end(),
262 [](auto& deriv) { return deriv == Jacobian::Line::VMR; }) or
263 std::any_of(jacobian_quantities.begin(),
264 jacobian_quantities.end(),
265 [model](auto& deriv) {
266 return deriv == Jacobian::Special::ArrayOfSpeciesTagVMR and
267 std::any_of(deriv.Target().species_array_id.begin(),
268 deriv.Target().species_array_id.end(),
269 [model](auto& tag) {
270 return tag.Isotopologue() == model;
271 });
272 });
273
274 if (do_freq_jac or do_temp_jac or do_vmrs_jac) {
276 PropagationMatrix pm(f_grid.nelem());
277 PropagationMatrix dpm(f_grid.nelem());
278 compute_selection<false>(pm,
279 model,
280 f_grid,
281 rtp_pressure,
282 rtp_temperature,
283 vmr,
284 predefined_model_data);
285
286 // Add absorption to the forward parameter
287 propmat_clearsky.Kjj() += pm.Kjj();
288
289 if (do_temp_jac) {
290 const Numeric d = temperature_perturbation(jacobian_quantities);
291 ARTS_ASSERT(d not_eq 0)
292
293 dpm.SetZero();
294 compute_selection<false>(dpm,
295 model,
296 f_grid,
297 rtp_pressure,
298 rtp_temperature + d,
299 vmr,
300 predefined_model_data);
301 dpm -= pm;
302 dpm /= d;
303 for (Index iq = 0; iq < dpropmat_clearsky_dx.nelem(); iq++) {
304 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
305 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
306 }
307 }
308 }
309
310 if (do_freq_jac) {
311 const Numeric d = frequency_perturbation(jacobian_quantities);
312 ARTS_ASSERT(d not_eq 0)
313
314 Vector f_grid_d{f_grid};
315 f_grid_d += d;
316
317 dpm.SetZero();
318 compute_selection<false>(dpm,
319 model,
320 f_grid_d,
321 rtp_pressure,
322 rtp_temperature,
323 vmr,
324 predefined_model_data);
325 dpm -= pm;
326 dpm /= d;
327 for (Index iq = 0; iq < dpropmat_clearsky_dx.nelem(); iq++) {
328 if (is_frequency_parameter(jacobian_quantities[iq])) {
329 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
330 }
331 }
332 }
333
334 for (Index iq = 0; iq < dpropmat_clearsky_dx.nelem(); iq++) {
335 auto& deriv = jacobian_quantities[iq];
336 if (deriv == Jacobian::Line::VMR) {
337 if (compute_vmr_deriv<false>(dpm,
338 pm,
339 model,
340 f_grid,
341 rtp_pressure,
342 rtp_temperature,
343 vmr,
344 deriv.QuantumIdentity().Species(),
345 predefined_model_data))
346 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
347 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR and
348 std::any_of(deriv.Target().species_array_id.begin(),
349 deriv.Target().species_array_id.end(),
350 [model](auto& tag) {
351 return tag.Isotopologue() == model;
352 })) {
353 if (compute_vmr_deriv<true>(
354 dpm,
355 pm,
356 model,
357 f_grid,
358 rtp_pressure,
359 rtp_temperature,
360 vmr,
361 deriv.Target().species_array_id.front().Spec(),
362 predefined_model_data))
363 dpropmat_clearsky_dx[iq].Kjj() += dpm.Kjj();
364 }
365 }
366 } else {
367 compute_selection<false>(propmat_clearsky,
368 model,
369 f_grid,
370 rtp_pressure,
371 rtp_temperature,
372 vmr,
373 predefined_model_data);
374 }
375}
376} // namespace Absorption::PredefinedModel
This can be used to make arrays out of anything.
Definition: array.h:48
void SetZero()
Sets all data to zero.
The Vector class.
Definition: matpackI.h:910
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1175
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1150
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1006
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1183
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1191
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
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, const PredefinedModelData &predefined_model_data)
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, const PredefinedModelData &predefined_model_data)
Compute the selected model and returns if it can be computed.
constexpr Numeric dvmr_calc(Numeric x) noexcept
Sets a VMR perturbation.
bool can_compute(const SpeciesIsotopeRecord &model)
Returns true if the model 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, const PredefinedModelData &predefined_model_data)
Compute the predefined model.
constexpr Index find_species_index(const Species spec, const std::string_view isot) noexcept
Stuff related to the propagation matrix.
Contains known required VMR values.
Struct containing all information needed about one isotope.
Definition: isotopologues.h:16
#define d