ARTS 2.5.10 (git: 2f1c442c)
m_abs.cc
Go to the documentation of this file.
1
2/* Copyright (C) 2000-2012
3 Stefan Buehler <sbuehler@ltu.se>
4 Patrick Eriksson <patrick.eriksson@chalmers.se>
5 Axel von Engeln <engeln@uni-bremen.de>
6 Thomas Kuhn <tkuhn@uni-bremen.de>
7
8 This program is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 2, or (at your option) any
11 later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
21 USA. */
22
23//
24
33#include <algorithm>
34#include <cmath>
35#include <cstddef>
36#include <memory>
37#include <utility>
38
39#include "absorption.h"
40#include "absorptionlines.h"
41#include "agenda_class.h"
42#include "agenda_set.h"
43#include "array.h"
44#include "arts.h"
45#include "arts_constants.h"
46#include "arts_omp.h"
47#include "artstime.h"
48#include "auto_md.h"
49#include "check_input.h"
50#include "debug.h"
51#include "depr.h"
52#include "file.h"
53#include "global_data.h"
54#include "hitran_species.h"
55#include "jacobian.h"
56#include "lineshape.h"
57#include "m_xml.h"
58#include "math_funcs.h"
59#include "matpackI.h"
60#include "messages.h"
61#include "methods.h"
62#include "montecarlo.h"
63#include "optproperties.h"
64#include "parameters.h"
65#include "physics_funcs.h"
66#include "rte.h"
67#include "species_tags.h"
68#include "xml_io.h"
69
70#ifdef ENABLE_NETCDF
71#include <netcdf.h>
72
73#include "nc_io.h"
74#endif
75
78inline constexpr Numeric PI=Constant::pi;
81
82/* Workspace method: Doxygen documentation will be auto-generated */
83void AbsInputFromRteScalars( // WS Output:
84 Vector& abs_p,
85 Vector& abs_t,
86 Matrix& abs_vmrs,
87 // WS Input:
88 const Numeric& rtp_pressure,
89 const Numeric& rtp_temperature,
90 const Vector& rtp_vmr,
91 const Verbosity&) {
92 // Prepare abs_p:
93 abs_p.resize(1);
94 abs_p = rtp_pressure;
95
96 // Prepare abs_t:
97 abs_t.resize(1);
98 abs_t = rtp_temperature;
99
100 // Prepare abs_vmrs:
101 abs_vmrs.resize(rtp_vmr.nelem(), 1);
102 abs_vmrs = rtp_vmr;
103}
104
105/* Workspace method: Doxygen documentation will be auto-generated */
107 ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
108 // WS Input:
109 const ArrayOfAbsorptionLines& abs_lines,
110 const ArrayOfArrayOfSpeciesTag& abs_species,
111 const Verbosity&) {
112 // Size is set but inner size will now change from the original definition of species tags...
113 abs_lines_per_species.resize(abs_species.nelem());
114
115 // The inner arrays need to be emptied, because they may contain lines
116 // from a previous calculation
117 for (auto& lines : abs_lines_per_species) lines.resize(0);
118
119#pragma omp parallel for schedule(dynamic) if (!arts_omp_in_parallel())
120 for (Index ilines = 0; ilines < abs_lines.nelem(); ilines++) {
121 AbsorptionLines lines = abs_lines[ilines];
122
123 // Skip empty lines
124 if (lines.NumLines() == 0) continue;
125
126 // Loop all the tags
127 for (Index i = 0; i < abs_species.nelem() and lines.NumLines(); i++) {
128 for (auto& this_tag : abs_species[i]) {
129 // Test isotopologue, we have to hit the end of the list for no isotopologue or the exact value
130 if (not same_or_joker(this_tag.Isotopologue(), lines.Isotopologue()))
131 continue;
132
133 // If there is a frequency range, we have to check so that only selected lines are included
134 if (this_tag.lower_freq >= 0 or this_tag.upper_freq >= 0) {
135 const Numeric low = (this_tag.lower_freq >= 0)
136 ? this_tag.lower_freq
137 : std::numeric_limits<Numeric>::lowest();
138 const Numeric upp = (this_tag.upper_freq >= 0)
139 ? this_tag.upper_freq
140 : std::numeric_limits<Numeric>::max();
141
142 // Fill up a copy of the line record to match with the wished frequency criteria
143 AbsorptionLines these_lines = lines;
144 these_lines.lines.resize(0);
145 for (Index k = lines.NumLines() - 1; k >= 0; k--)
146 if (low <= lines.lines[k].F0 and upp >= lines.lines[k].F0)
147 these_lines.AppendSingleLine(lines.PopLine(k));
148
149 // Append these lines after sorting them if there are any of them
150 if (these_lines.NumLines()) {
151 these_lines.ReverseLines();
152#pragma omp critical
153 abs_lines_per_species[i].push_back(these_lines);
154 }
155
156 // If this means we have deleted all lines, then we leave
157 if (lines.NumLines() == 0) goto leave_inner_loop;
158 } else {
159#pragma omp critical
160 abs_lines_per_species[i].push_back(lines);
161 goto leave_inner_loop;
162 }
163 }
164 }
165 leave_inner_loop : {}
166 }
167
168 abs_lines_per_species.shrink_to_fit();
169 for (auto& spec_band : abs_lines_per_species)
170 std::sort(spec_band.begin(), spec_band.end(), [](auto& a, auto& b) {
171 return a.lines.size() and b.lines.size() and
172 a.lines.front().F0 < b.lines.front().F0;
173 });
174}
175
176/* Workspace method: Doxygen documentation will be auto-generated */
179 Index& propmat_clearsky_agenda_checked,
180 // Control Parameters:
181 const String& basename,
182 const Verbosity& verbosity) {
184
185 // Invalidate agenda check flags
186 propmat_clearsky_agenda_checked = false;
187
188 // We want to make lists of included and excluded species:
189 ArrayOfString included(0), excluded(0);
190
191 tgs.resize(0);
192
193 for (Index i = 0; i < Index(Species::Species::FINAL); ++i) {
194 const String specname = Species::toShortName(Species::Species(i));
195
196 String filename = basename;
197 if (basename.length() && basename[basename.length() - 1] != '/')
198 filename += ".";
199 filename += specname;
200
201 try {
202 find_xml_file(filename, verbosity);
203 // Add to included list:
204 included.push_back(specname);
205
206 // Add this tag group to tgs:
207 tgs.emplace_back(ArrayOfSpeciesTag(specname));
208 } catch (const std::runtime_error& e) {
209 // The file for the species could not be found.
210 excluded.push_back(specname);
211 }
212 }
213
214 // Some nice output:
215 out2 << " Included Species (" << included.nelem() << "):\n";
216 for (Index i = 0; i < included.nelem(); ++i)
217 out2 << " " << included[i] << "\n";
218
219 out2 << " Excluded Species (" << excluded.nelem() << "):\n";
220 for (Index i = 0; i < excluded.nelem(); ++i)
221 out2 << " " << excluded[i] << "\n";
222}
223
224/* Workspace method: Doxygen documentation will be auto-generated */
225void abs_speciesDefineAll( // WS Output:
226 ArrayOfArrayOfSpeciesTag& abs_species,
227 Index& propmat_clearsky_agenda_checked,
228 // Control Parameters:
229 const Verbosity& verbosity) {
230 // Species lookup data:
231
232 // We want to make lists of all species
233 ArrayOfString specs(0);
234 for (Index i = 0; i < Index(Species::Species::FINAL); ++i) {
235 if (Species::Species(i) not_eq Species::Species::Bath) {
236 specs.emplace_back(Species::toShortName(Species::Species(i)));
237 }
238 }
239
240 // Set the values
241 abs_speciesSet(abs_species,
242 propmat_clearsky_agenda_checked,
243 specs,
244 verbosity);
245}
246
247/* Workspace method: Doxygen documentation will be auto-generated */
248void AbsInputFromAtmFields( // WS Output:
249 Vector& abs_p,
250 Vector& abs_t,
251 Matrix& abs_vmrs,
252 // WS Input:
253 const Index& atmosphere_dim,
254 const Vector& p_grid,
255 const Tensor3& t_field,
256 const Tensor4& vmr_field,
257 const Verbosity&) {
258 // First, make sure that we really have a 1D atmosphere:
259 ARTS_USER_ERROR_IF(1 != atmosphere_dim,
260 "Atmospheric dimension must be 1D, but atmosphere_dim is ",
261 atmosphere_dim,
262 ".")
263
264 abs_p = p_grid;
265 abs_t = t_field(joker, 0, 0);
266 abs_vmrs = vmr_field(joker, joker, 0, 0);
267}
268
269//======================================================================
270// Methods related to continua
271//======================================================================
272
273/* Workspace method: Doxygen documentation will be auto-generated */
275 StokesVector& nlte_source,
276 ArrayOfStokesVector& dnlte_source_dx,
277 // WS Input:
278 const ArrayOfMatrix& src_coef_per_species,
279 const ArrayOfMatrix& dsrc_coef_dx,
280 const ArrayOfRetrievalQuantity& jacobian_quantities,
281 const Vector& f_grid,
282 const Numeric& rtp_temperature,
283 const Verbosity&) {
284 // nlte_source has format
285 // [ abs_species, f_grid, stokes_dim ].
286 // src_coef_per_species has format ArrayOfMatrix (over species),
287 // where for each species the matrix has format [f_grid, abs_p].
288
289 Index n_species = src_coef_per_species.nelem(); // # species
290
291 ARTS_USER_ERROR_IF(not n_species, "Must have at least one species.")
292
293 Index n_f = src_coef_per_species[0].nrows(); // # frequencies
294
295 // # pressures must be 1:
296 ARTS_USER_ERROR_IF(1 not_eq src_coef_per_species[0].ncols(),
297 "Must have exactly one pressure.")
298
299 // Check frequency dimension of propmat_clearsky
300 ARTS_USER_ERROR_IF(nlte_source.NumberOfFrequencies() not_eq n_f,
301 "Frequency dimension of nlte_source does not\n"
302 "match abs_coef_per_species.")
303
304 const Vector B = planck(f_grid, rtp_temperature);
305
306 StokesVector sv(n_f, nlte_source.StokesDimensions());
307 for (Index si = 0; si < n_species; ++si) {
308 sv.Kjj() = src_coef_per_species[si](joker, 0);
309 sv *= B;
310 nlte_source.Kjj() += sv.Kjj();
311 }
312
313 // Jacobian
314 for (Index ii = 0; ii < jacobian_quantities.nelem(); ii++) {
315 const auto& deriv = jacobian_quantities[ii];
316
317 if (not deriv.propmattype()) continue;
318
319 if (deriv == Jacobian::Atm::Temperature) {
320 const Vector dB = dplanck_dt(f_grid, rtp_temperature);
321
322 for (Index si = 0; si < n_species; ++si) {
323 sv.Kjj() = src_coef_per_species[si](joker, 0);
324 sv *= dB;
325 dnlte_source_dx[ii].Kjj() += sv.Kjj();
326 }
327
328 sv.Kjj() = dsrc_coef_dx[ii](joker, 0);
329 sv *= B;
330 dnlte_source_dx[ii].Kjj() += sv.Kjj();
331 } else if (is_frequency_parameter(deriv)) {
332 const Vector dB = dplanck_df(f_grid, rtp_temperature);
333
334 for (Index si = 0; si < n_species; ++si) {
335 sv.Kjj() = src_coef_per_species[si](joker, 0);
336 sv *= dB;
337 dnlte_source_dx[ii].Kjj() += sv.Kjj();
338 }
339
340 sv.Kjj() = dsrc_coef_dx[ii](joker, 0);
341 sv *= B;
342 dnlte_source_dx[ii].Kjj() += sv.Kjj();
343 } else {
344 sv.Kjj() = dsrc_coef_dx[ii](joker, 0);
345 sv *= B;
346 dnlte_source_dx[ii].Kjj() += sv.Kjj();
347 }
348 }
349}
350
351/* Workspace method: Doxygen documentation will be auto-generated */
352void propmat_clearskyInit( //WS Output
353 PropagationMatrix& propmat_clearsky,
354 StokesVector& nlte_source,
355 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
356 ArrayOfStokesVector& dnlte_source_dx,
357 //WS Input
358 const ArrayOfRetrievalQuantity& jacobian_quantities,
359 const Vector& f_grid,
360 const Index& stokes_dim,
361 const Index& propmat_clearsky_agenda_checked,
362 const Verbosity&) {
363 const Index nf = f_grid.nelem();
364 const Index nq = jacobian_quantities.nelem();
365
367 !propmat_clearsky_agenda_checked,
368 "You must call *propmat_clearsky_agenda_checkedCalc* before calling this method.")
369
370 ARTS_USER_ERROR_IF(not nf, "No frequencies");
371
372 ARTS_USER_ERROR_IF(stokes_dim < 1 or stokes_dim > 4,
373 "stokes_dim not in [1, 2, 3, 4]");
374
375 // Set size of propmat_clearsky or reset it's values
376 if (propmat_clearsky.StokesDimensions() == stokes_dim and
377 propmat_clearsky.NumberOfFrequencies() == nf) {
378 propmat_clearsky.SetZero();
379 } else {
380 propmat_clearsky = PropagationMatrix(nf, stokes_dim);
381 }
382
383 // Set size of dpropmat_clearsky_dx or reset it's values
384 if (dpropmat_clearsky_dx.nelem() not_eq nq) {
385 dpropmat_clearsky_dx =
386 ArrayOfPropagationMatrix(nq, PropagationMatrix(nf, stokes_dim));
387 } else {
388 for (auto& pm : dpropmat_clearsky_dx) {
389 if (pm.StokesDimensions() == stokes_dim and
390 pm.NumberOfFrequencies() == nf) {
391 pm.SetZero();
392 } else {
393 pm = PropagationMatrix(nf, stokes_dim);
394 }
395 }
396 }
397
398 // Set size of nlte_source or reset it's values
399 if (nlte_source.StokesDimensions() == stokes_dim and
400 nlte_source.NumberOfFrequencies() == nf) {
401 nlte_source.SetZero();
402 } else {
403 nlte_source = StokesVector(nf, stokes_dim);
404 }
405
406 // Set size of dnlte_source_dx or reset it's values
407 if (dnlte_source_dx.nelem() not_eq nq) {
408 dnlte_source_dx = ArrayOfStokesVector(nq, StokesVector(nf, stokes_dim));
409 } else {
410 for (auto& pm : dnlte_source_dx) {
411 if (pm.StokesDimensions() == stokes_dim and
412 pm.NumberOfFrequencies() == nf) {
413 pm.SetZero();
414 } else {
415 pm = StokesVector(nf, stokes_dim);
416 }
417 }
418 }
419}
420
421/* Workspace method: Doxygen documentation will be auto-generated */
423 PropagationMatrix& propmat_clearsky,
424 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
425 const Index& stokes_dim,
426 const Index& atmosphere_dim,
427 const Vector& f_grid,
428 const ArrayOfArrayOfSpeciesTag& abs_species,
429 const ArrayOfSpeciesTag& select_abs_species,
430 const ArrayOfRetrievalQuantity& jacobian_quantities,
431 const Vector& rtp_vmr,
432 const Vector& rtp_los,
433 const Vector& rtp_mag,
434 const Verbosity&) {
435 Index ife = -1;
436 for (Index sp = 0; sp < abs_species.nelem() && ife < 0; sp++) {
437 if (abs_species[sp].FreeElectrons()) {
438 ife = sp;
439 }
440 }
441
442 ARTS_USER_ERROR_IF(ife < 0,
443 "Free electrons not found in *abs_species* and "
444 "Faraday rotation can not be calculated.");
445
446 // Allow early exit for lookup table calculations
447 if (select_abs_species.nelem() and select_abs_species not_eq abs_species[ife]) return;
448
449 // All the physical constants joined into one static constant:
450 // (abs as e defined as negative)
451 static const Numeric FRconst =
455
456 const bool do_magn_jac = do_magnetic_jacobian(jacobian_quantities);
457 const Numeric dmag = magnetic_field_perturbation(jacobian_quantities);
458
460 stokes_dim < 3,
461 "To include Faraday rotation, stokes_dim >= 3 is required.")
463 atmosphere_dim == 1 && rtp_los.nelem() < 1,
464 "For applying propmat_clearskyAddFaraday, los needs to be specified\n"
465 "(at least zenith angle component for atmosphere_dim==1),\n"
466 "but it is not.\n")
468 atmosphere_dim > 1 && rtp_los.nelem() < 2,
469 "For applying propmat_clearskyAddFaraday, los needs to be specified\n"
470 "(both zenith and azimuth angle components for atmosphere_dim>1),\n"
471 "but it is not.\n")
472
473 const Numeric ne = rtp_vmr[ife];
474
475 if (ne != 0 && (rtp_mag[0] != 0 || rtp_mag[1] != 0 || rtp_mag[2] != 0)) {
476 // Include remaining terms, beside /f^2
477 const Numeric c1 =
478 2 * FRconst *
480 rtp_los, rtp_mag[0], rtp_mag[1], rtp_mag[2], atmosphere_dim);
481
482 Numeric dc1_u = 0.0, dc1_v = 0.0, dc1_w = 0.0;
483 if (do_magn_jac) {
484 dc1_u = (2 * FRconst *
485 dotprod_with_los(rtp_los,
486 rtp_mag[0] + dmag,
487 rtp_mag[1],
488 rtp_mag[2],
489 atmosphere_dim) -
490 c1) /
491 dmag;
492 dc1_v = (2 * FRconst *
493 dotprod_with_los(rtp_los,
494 rtp_mag[0],
495 rtp_mag[1] + dmag,
496 rtp_mag[2],
497 atmosphere_dim) -
498 c1) /
499 dmag;
500 dc1_w = (2 * FRconst *
501 dotprod_with_los(rtp_los,
502 rtp_mag[0],
503 rtp_mag[1],
504 rtp_mag[2] + dmag,
505 atmosphere_dim) -
506 c1) /
507 dmag;
508 }
509
510 for (Index iv = 0; iv < f_grid.nelem(); iv++) {
511 const Numeric f2 = f_grid[iv] * f_grid[iv];
512 const Numeric r = ne * c1 / f2;
513 propmat_clearsky.AddFaraday(r, iv);
514
515 // The Jacobian loop
516 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
517 if (is_frequency_parameter(jacobian_quantities[iq]))
518 dpropmat_clearsky_dx[iq].AddFaraday(-2.0 * ne * r / f_grid[iv], iv);
519 else if (jacobian_quantities[iq] == Jacobian::Atm::MagneticU)
520 dpropmat_clearsky_dx[iq].AddFaraday(ne * dc1_u / f2, iv);
521 else if (jacobian_quantities[iq] == Jacobian::Atm::MagneticV)
522 dpropmat_clearsky_dx[iq].AddFaraday(ne * dc1_v / f2, iv);
523 else if (jacobian_quantities[iq] == Jacobian::Atm::MagneticW)
524 dpropmat_clearsky_dx[iq].AddFaraday(ne * dc1_w / f2, iv);
525 else if (jacobian_quantities[iq] == Jacobian::Atm::Electrons)
526 dpropmat_clearsky_dx[iq].AddFaraday(r, iv);
527 else if (jacobian_quantities[iq] == abs_species[ife])
528 dpropmat_clearsky_dx[iq].AddFaraday(r, iv);
529 }
530 }
531 }
532}
533
534/* Workspace method: Doxygen documentation will be auto-generated */
536 // WS Output:
537 PropagationMatrix& propmat_clearsky,
538 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
539 // WS Input:
540 const Index& stokes_dim,
541 const Index& atmosphere_dim,
542 const Vector& f_grid,
543 const ArrayOfArrayOfSpeciesTag& abs_species,
544 const ArrayOfSpeciesTag& select_abs_species,
545 const ArrayOfRetrievalQuantity& jacobian_quantities,
546 const Vector& rtp_vmr,
547 const Vector& rtp_los,
548 const Numeric& rtp_temperature,
549 const ArrayOfArrayOfSingleScatteringData& scat_data,
550 const Index& scat_data_checked,
551 const Index& use_abs_as_ext,
552 // Verbosity object:
553 const Verbosity& verbosity) {
555
556 ARTS_USER_ERROR_IF(select_abs_species.nelem(), R"--(
557 We do not yet support select_abs_species for lookup table calculations
558 )--")
559
560 // (i)yCalc only checks scat_data_checked if cloudbox is on. It is off here,
561 // though, i.e. we need to check it here explicitly. (Also, cloudboxOff sets
562 // scat_data_checked=0 as it does not check it and as we ususally don't need
563 // scat_data for clearsky cases, hence don't want to check them by
564 // scat_data_checkedCalc in that case. This approach seems to be the more
565 // handy compared to cloudboxOff setting scat_data_checked=1 without checking
566 // it assuming we won't use it anyways.)
567 ARTS_USER_ERROR_IF(scat_data_checked != 1,
568 "The scat_data must be flagged to have "
569 "passed a consistency check (scat_data_checked=1).")
570
571 const Index ns = TotalNumberOfElements(scat_data);
572 Index np = 0;
573 for (Index sp = 0; sp < abs_species.nelem(); sp++) {
574 if (abs_species[sp].Particles()) {
575 np++;
576 }
577 }
578
580 np == 0,
581 "For applying propmat_clearskyAddParticles, *abs_species* needs to"
582 "contain species 'particles', but it does not.\n")
583
585 ns != np,
586 "Number of 'particles' entries in abs_species and of elements in\n"
587 "*scat_data* needs to be identical. But you have ",
588 np,
589 " 'particles' entries\n"
590 "and ",
591 ns,
592 " *scat_data* elements.\n")
593
595 atmosphere_dim == 1 && rtp_los.nelem() < 1,
596 "For applying *propmat_clearskyAddParticles*, *rtp_los* needs to be specified\n"
597 "(at least zenith angle component for atmosphere_dim==1),\n"
598 "but it is not.\n")
600 atmosphere_dim > 1 && rtp_los.nelem() < 2,
601 "For applying *propmat_clearskyAddParticles*, *rtp_los* needs to be specified\n"
602 "(both zenith and azimuth angle components for atmosphere_dim>1),\n"
603 "but it is not.\n")
604
605 // Use for rescaling vmr of particulates
606 Numeric rtp_vmr_sum = 0.0;
607
608 // Tests and setup partial derivatives
609 const bool do_jac_temperature = do_temperature_jacobian(jacobian_quantities);
610 const bool do_jac_frequencies = do_frequency_jacobian(jacobian_quantities);
611 const Numeric dT = temperature_perturbation(jacobian_quantities);
612
613 const Index na = abs_species.nelem();
614 Vector rtp_los_back;
615 mirror_los(rtp_los_back, rtp_los, atmosphere_dim);
616
617 // 170918 JM: along with transition to use of new-type (aka
618 // pre-f_grid-interpolated) scat_data, freq perturbation switched off. Typical
619 // clear-sky freq perturbations yield insignificant effects in particle
620 // properties. Hence, this feature is neglected here.
621 if (do_jac_frequencies) {
622 out1 << "WARNING:\n"
623 << "Frequency perturbation not available for absorbing particles.\n";
624 }
625
626 // creating temporary output containers
627 ArrayOfArrayOfTensor5 ext_mat_Nse;
628 ArrayOfArrayOfTensor4 abs_vec_Nse;
629 ArrayOfArrayOfIndex ptypes_Nse;
630 Matrix t_ok;
631
632 // preparing input in format needed
633 Vector T_array;
634 if (do_jac_temperature) {
635 T_array.resize(2);
636 T_array = rtp_temperature;
637 T_array[1] += dT;
638 } else {
639 T_array.resize(1);
640 T_array = rtp_temperature;
641 }
642 Matrix dir_array(1, 2);
643 dir_array(0, joker) = rtp_los_back;
644
645 // ext/abs per scat element for all freqs at once
646 opt_prop_NScatElems(ext_mat_Nse,
647 abs_vec_Nse,
648 ptypes_Nse,
649 t_ok,
650 scat_data,
651 stokes_dim,
652 T_array,
653 dir_array,
654 -1);
655
656 const Index nf = abs_vec_Nse[0][0].nbooks();
657 Tensor3 tmp(nf, stokes_dim, stokes_dim);
658
659 // Internal computations necessary since it relies on zero start
660 PropagationMatrix internal_propmat(propmat_clearsky.NumberOfFrequencies(),
661 propmat_clearsky.StokesDimensions());
662
663 // loop over the scat_data and link them with correct vmr_field entry according
664 // to the position of the particle type entries in abs_species.
665 Index sp = 0;
666 Index i_se_flat = 0;
667 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++) {
668 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
669 // forward to next particle entry in abs_species
670 while (sp < na && not abs_species[sp].Particles()) sp++;
671 internal_propmat.SetZero();
672
673 // running beyond number of abs_species entries when looking for next
674 // particle entry. shouldn't happen, though.
675 ARTS_ASSERT(sp < na);
677 rtp_vmr[sp] < 0.,
678 "Negative absorbing particle 'vmr' (aka number density)"
679 " encountered:\n"
680 "scat species #",
681 i_ss,
682 ", scat elem #",
683 i_se,
684 " (vmr_field entry #",
685 sp,
686 ")\n")
687
688 if (rtp_vmr[sp] > 0.) {
689 ARTS_USER_ERROR_IF(t_ok(i_se_flat, 0) < 0.,
690 "Temperature interpolation error:\n"
691 "scat species #",
692 i_ss,
693 ", scat elem #",
694 i_se,
695 "\n")
696 if (use_abs_as_ext) {
697 if (nf > 1)
698 for (Index iv = 0; iv < f_grid.nelem(); iv++)
699 internal_propmat.AddAbsorptionVectorAtPosition(
700 abs_vec_Nse[i_ss][i_se](iv, 0, 0, joker), iv);
701 else
702 for (Index iv = 0; iv < f_grid.nelem(); iv++)
703 internal_propmat.AddAbsorptionVectorAtPosition(
704 abs_vec_Nse[i_ss][i_se](0, 0, 0, joker), iv);
705 } else {
706 if (nf > 1)
707 for (Index iv = 0; iv < f_grid.nelem(); iv++)
708 internal_propmat.SetAtPosition(
709 ext_mat_Nse[i_ss][i_se](iv, 0, 0, joker, joker), iv);
710 else
711 for (Index iv = 0; iv < f_grid.nelem(); iv++)
712 internal_propmat.SetAtPosition(
713 ext_mat_Nse[i_ss][i_se](0, 0, 0, joker, joker), iv);
714 }
715 propmat_clearsky += rtp_vmr[sp] * internal_propmat;
716 }
717
718 // For temperature derivatives (so we don't need to check it in jac loop)
719 if (do_jac_temperature) {
721 t_ok(i_se_flat, 1) < 0.,
722 "Temperature interpolation error (in perturbation):\n"
723 "scat species #",
724 i_ss,
725 ", scat elem #",
726 i_se,
727 "\n")
728 }
729
730 // For number density derivatives
731 if (jacobian_quantities.nelem()) rtp_vmr_sum += rtp_vmr[sp];
732
733 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
734 const auto& deriv = jacobian_quantities[iq];
735
736 if (not deriv.propmattype()) continue;
737
738 if (deriv == Jacobian::Atm::Temperature) {
739 if (use_abs_as_ext) {
740 tmp(joker, joker, 0) = abs_vec_Nse[i_ss][i_se](joker, 1, 0, joker);
741 tmp(joker, joker, 0) -= abs_vec_Nse[i_ss][i_se](joker, 0, 0, joker);
742 } else {
743 tmp = ext_mat_Nse[i_ss][i_se](joker, 1, 0, joker, joker);
744 tmp -= ext_mat_Nse[i_ss][i_se](joker, 0, 0, joker, joker);
745 }
746
747 tmp *= rtp_vmr[sp];
748 tmp /= dT;
749
750 if (nf > 1)
751 for (Index iv = 0; iv < f_grid.nelem(); iv++)
752 if (use_abs_as_ext)
753 dpropmat_clearsky_dx[iq].AddAbsorptionVectorAtPosition(
754 tmp(iv, joker, 0), iv);
755 else
756 dpropmat_clearsky_dx[iq].AddAtPosition(tmp(iv, joker, joker),
757 iv);
758 else
759 for (Index iv = 0; iv < f_grid.nelem(); iv++)
760 if (use_abs_as_ext)
761 dpropmat_clearsky_dx[iq].AddAbsorptionVectorAtPosition(
762 tmp(0, joker, 0), iv);
763 else
764 dpropmat_clearsky_dx[iq].AddAtPosition(tmp(0, joker, joker),
765 iv);
766 }
767
768 else if (deriv == Jacobian::Atm::Particulates) {
769 for (Index iv = 0; iv < f_grid.nelem(); iv++)
770 dpropmat_clearsky_dx[iq].AddAtPosition(internal_propmat, iv);
771 }
772
773 else if (deriv == abs_species[sp]) {
774 dpropmat_clearsky_dx[iq] += internal_propmat;
775 }
776 }
777 sp++;
778 i_se_flat++;
779 }
780 }
781
782 //checking that no further 'particle' entry left after all scat_data entries
783 //are processes. this is basically not necessary. but checking it anyway to
784 //really be safe. remove later, when more extensively tested.
785 while (sp < na) {
786 ARTS_ASSERT(abs_species[sp][0].Type() != Species::TagType::Particles);
787 sp++;
788 }
789
790 if (rtp_vmr_sum != 0.0) {
791 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
792 const auto& deriv = jacobian_quantities[iq];
793
794 if (not deriv.propmattype()) continue;
795
796 if (deriv == Jacobian::Atm::Particulates) {
797 dpropmat_clearsky_dx[iq] /= rtp_vmr_sum;
798 }
799 }
800 }
801}
802
804 const Vector& f_grid,
805 const Numeric& sparse_df,
806 const String& speedup_option,
807 // Verbosity object:
808 const Verbosity&) {
809 // Return empty for nothing
810 if (not f_grid.nelem()) {
811 sparse_f_grid.resize(0);
812 return;
813 };
814
815 switch (Options::toLblSpeedupOrThrow(speedup_option)) {
816 case Options::LblSpeedup::LinearIndependent:
817 sparse_f_grid = LineShape::linear_sparse_f_grid(f_grid, sparse_df);
819 break;
820 case Options::LblSpeedup::QuadraticIndependent:
821 sparse_f_grid = LineShape::triple_sparse_f_grid(f_grid, sparse_df);
822 break;
823 case Options::LblSpeedup::None:
824 sparse_f_grid.resize(0);
825 break;
826 case Options::LblSpeedup::FINAL: { /* Leave last */
827 }
828 }
829}
830
832 const Numeric& sparse_df,
833 const String& speedup_option,
834 // Verbosity object:
835 const Verbosity& verbosity) {
836 Vector sparse_f_grid;
838 sparse_f_grid, f_grid, sparse_df, speedup_option, verbosity);
839 return sparse_f_grid;
840}
841
842/* Workspace method: Doxygen documentation will be auto-generated */
843void propmat_clearskyAddLines( // Workspace reference:
844 // WS Output:
845 PropagationMatrix& propmat_clearsky,
846 StokesVector& nlte_source,
847 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
848 ArrayOfStokesVector& dnlte_source_dx,
849 // WS Input:
850 const Vector& f_grid,
851 const ArrayOfArrayOfSpeciesTag& abs_species,
852 const ArrayOfSpeciesTag& select_abs_species,
853 const ArrayOfRetrievalQuantity& jacobian_quantities,
854 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
855 const SpeciesIsotopologueRatios& isotopologue_ratios,
856 const Numeric& rtp_pressure,
857 const Numeric& rtp_temperature,
858 const EnergyLevelMap& rtp_nlte,
859 const Vector& rtp_vmr,
860 const Index& nlte_do,
861 const Index& lbl_checked,
862 // WS User Generic inputs
863 const Numeric& sparse_df,
864 const Numeric& sparse_lim,
865 const String& speedup_option,
866 const Index& robust,
867 // Verbosity object:
868 const Verbosity& verbosity) {
869 // Size of problem
870 const Index nf = f_grid.nelem();
871 const Index nq = jacobian_quantities.nelem();
872 const Index ns = abs_species.nelem();
873
874 // Possible things that can go wrong in this code (excluding line parameters)
875 ARTS_USER_ERROR_IF(not lbl_checked, "Must check LBL calculations")
876 check_abs_species(abs_species);
877 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
878 "*rtp_vmr* must match *abs_species*")
879 ARTS_USER_ERROR_IF(propmat_clearsky.NumberOfFrequencies() not_eq nf,
880 "*f_grid* must match *propmat_clearsky*")
881 ARTS_USER_ERROR_IF(nlte_source.NumberOfFrequencies() not_eq nf,
882 "*f_grid* must match *nlte_source*")
884 not nq and (nq not_eq dpropmat_clearsky_dx.nelem()),
885 "*dpropmat_clearsky_dx* must match derived form of *jacobian_quantities*")
887 not nq and bad_propmat(dpropmat_clearsky_dx, f_grid),
888 "*dpropmat_clearsky_dx* must have frequency dim same as *f_grid*")
890 nlte_do and (nq not_eq dnlte_source_dx.nelem()),
891 "*dnlte_source_dx* must match derived form of *jacobian_quantities* when non-LTE is on")
893 nlte_do and bad_propmat(dnlte_source_dx, f_grid),
894 "*dnlte_source_dx* must have frequency dim same as *f_grid* when non-LTE is on")
896 "Negative frequency (at least one value).")
897 ARTS_USER_ERROR_IF((any_cutoff(abs_lines_per_species) or speedup_option not_eq "None") and not is_increasing(f_grid),
898 "Must be sorted and increasing if any cutoff or speedup is used.")
900 "Negative VMR (at least one value).")
902 "Negative NLTE (at least one value).")
903 ARTS_USER_ERROR_IF(rtp_temperature <= 0, "Non-positive temperature")
904 ARTS_USER_ERROR_IF(rtp_pressure <= 0, "Non-positive pressure")
906 sparse_lim > 0 and sparse_df > sparse_lim,
907 "If sparse grids are to be used, the limit must be larger than the grid-spacing.\n"
908 "The limit is ",
909 sparse_lim,
910 " Hz and the grid_spacing is ",
911 sparse_df,
912 " Hz")
913
914 if (not nf) return;
915
916 // Deal with sparse computational grid
917 const Vector f_grid_sparse = create_sparse_f_grid_internal(
918 f_grid, sparse_df, speedup_option, verbosity);
919 const Options::LblSpeedup speedup_type =
920 f_grid_sparse.nelem() ? Options::toLblSpeedupOrThrow(speedup_option)
921 : Options::LblSpeedup::None;
923 sparse_lim <= 0 and speedup_type not_eq Options::LblSpeedup::None,
924 "Must have a sparse limit if you set speedup_option")
925
926 // Calculations data
927 LineShape::ComputeData com(f_grid, jacobian_quantities, nlte_do);
928 LineShape::ComputeData sparse_com(
929 f_grid_sparse, jacobian_quantities, nlte_do);
930
931 if (arts_omp_in_parallel()) {
932 for (Index ispecies = 0; ispecies < ns; ispecies++) {
933 if (select_abs_species.nelem() and
934 select_abs_species not_eq abs_species[ispecies])
935 continue;
936
937 // Skip it if there are no species or there is Zeeman requested
938 if (not abs_species[ispecies].nelem() or abs_species[ispecies].Zeeman() or
939 not abs_lines_per_species[ispecies].nelem())
940 continue;
941
942 for (auto& band : abs_lines_per_species[ispecies]) {
944 sparse_com,
945 band,
946 jacobian_quantities,
947 rtp_nlte,
948 band.BroadeningSpeciesVMR(rtp_vmr, abs_species),
949 abs_species[ispecies],
950 rtp_vmr[ispecies],
951 isotopologue_ratios[band.Isotopologue()],
952 rtp_pressure,
953 rtp_temperature,
954 0,
955 sparse_lim,
957 speedup_type,
958 robust not_eq 0);
959 }
960 }
961 } else { // In parallel
962 const Index nbands = [](auto& lines) {
963 Index n = 0;
964 for (auto& abs_lines : lines) n += abs_lines.nelem();
965 return n;
966 }(abs_lines_per_species);
967
968 std::vector<LineShape::ComputeData> vcom(
971 f_grid, jacobian_quantities, static_cast<bool>(nlte_do)});
972 std::vector<LineShape::ComputeData> vsparse_com(
975 f_grid_sparse, jacobian_quantities, static_cast<bool>(nlte_do)});
976
977#pragma omp parallel for schedule(dynamic)
978 for (Index i = 0; i < nbands; i++) {
979 const auto [ispecies, iband] =
980 flat_index(i, abs_species, abs_lines_per_species);
981
982 if (select_abs_species.nelem() and
983 select_abs_species not_eq abs_species[ispecies])
984 continue;
985
986 // Skip it if there are no species or there is Zeeman requested
987 if (not abs_species[ispecies].nelem() or abs_species[ispecies].Zeeman() or
988 not abs_lines_per_species[ispecies].nelem())
989 continue;
990
991 auto& band = abs_lines_per_species[ispecies][iband];
993 vsparse_com[arts_omp_get_thread_num()],
994 band,
995 jacobian_quantities,
996 rtp_nlte,
997 band.BroadeningSpeciesVMR(rtp_vmr, abs_species),
998 abs_species[ispecies],
999 rtp_vmr[ispecies],
1000 isotopologue_ratios[band.Isotopologue()],
1001 rtp_pressure,
1002 rtp_temperature,
1003 0,
1004 sparse_lim,
1006 speedup_type,
1007 robust not_eq 0);
1008 }
1009
1010 for (auto& pcom: vcom) com += pcom;
1011 for (auto& pcom: vsparse_com) sparse_com += pcom;
1012 }
1013
1014 switch (speedup_type) {
1015 case Options::LblSpeedup::LinearIndependent:
1016 com.interp_add_even(sparse_com);
1017 break;
1018 case Options::LblSpeedup::QuadraticIndependent:
1019 com.interp_add_triplequad(sparse_com);
1020 break;
1021 case Options::LblSpeedup::None: /* Do nothing */
1022 break;
1023 case Options::LblSpeedup::FINAL: { /* Leave last */
1024 }
1025 }
1026
1027 // Sum up the propagation matrix
1028 propmat_clearsky.Kjj() += com.F.real();
1029
1030 // Sum up the Jacobian
1031 for (Index j = 0; j < nq; j++) {
1032 if (not jacobian_quantities[j].propmattype()) continue;
1033 dpropmat_clearsky_dx[j].Kjj() += com.dF.real()(joker, j);
1034 }
1035
1036 if (nlte_do) {
1037 // Sum up the source vector
1038 nlte_source.Kjj() += com.N.real();
1039
1040 // Sum up the Jacobian
1041 for (Index j = 0; j < nq; j++) {
1042 if (not jacobian_quantities[j].propmattype()) continue;
1043 dnlte_source_dx[j].Kjj() += com.dN.real()(joker, j);
1044 }
1045 }
1046}
1047
1048/* Workspace method: Doxygen documentation will be auto-generated */
1050 const Vector& f_grid,
1051 const Index& stokes_dim,
1052 const Verbosity&) {
1053 propmat_clearsky = PropagationMatrix(f_grid.nelem(), stokes_dim);
1054}
1055
1056/* Workspace method: Doxygen documentation will be auto-generated */
1058 const Verbosity&) {
1059 for (Index i = 0; i < propmat_clearsky.NumberOfFrequencies(); i++)
1060 if (propmat_clearsky.Kjj()[i] < 0.0) propmat_clearsky.SetAtPosition(0.0, i);
1061}
1062
1063/* Workspace method: Doxygen documentation will be auto-generated */
1065 SpeciesIsotopologueRatios& isotopologue_ratios, const Verbosity&) {
1066 isotopologue_ratios = Species::isotopologue_ratiosInitFromBuiltin();
1067}
1068
1069/* Workspace method: Doxygen documentation will be auto-generated */
1071 SpeciesIsotopologueRatios& isotopologue_ratios, const Verbosity&) {
1072 isotopologue_ratios = Hitran::isotopologue_ratios();
1073}
1074
1075#ifdef ENABLE_NETCDF
1076/* Workspace method: Doxygen documentation will be auto-generated */
1077/* Included by Claudia Emde, 20100707 */
1078void WriteMolTau( //WS Input
1079 const Vector& f_grid,
1080 const Tensor3& z_field,
1081 const Tensor7& propmat_clearsky_field,
1082 const Index& atmosphere_dim,
1083 //Keyword
1084 const String& filename,
1085 const Verbosity&) {
1086 int retval, ncid;
1087 int nlev_dimid, nlyr_dimid, nwvl_dimid, stokes_dimid, none_dimid;
1088 int dimids[4];
1089 int wvlmin_varid, wvlmax_varid, z_varid, wvl_varid, tau_varid;
1090
1091 ARTS_USER_ERROR_IF(atmosphere_dim != 1,
1092 "WriteMolTau can only be used for atmosphere_dim=1")
1093#pragma omp critical(netcdf__critical_region)
1094 {
1095 // Open file
1096 if ((retval = nc_create(filename.c_str(), NC_CLOBBER, &ncid)))
1097 nca_error(retval, "nc_create");
1098
1099 // Define dimensions
1100 if ((retval = nc_def_dim(ncid, "nlev", (int)z_field.npages(), &nlev_dimid)))
1101 nca_error(retval, "nc_def_dim");
1102
1103 if ((retval =
1104 nc_def_dim(ncid, "nlyr", (int)z_field.npages() - 1, &nlyr_dimid)))
1105 nca_error(retval, "nc_def_dim");
1106
1107 if ((retval = nc_def_dim(ncid, "nwvl", (int)f_grid.nelem(), &nwvl_dimid)))
1108 nca_error(retval, "nc_def_dim");
1109
1110 if ((retval = nc_def_dim(ncid, "none", 1, &none_dimid)))
1111 nca_error(retval, "nc_def_dim");
1112
1113 if ((retval = nc_def_dim(ncid,
1114 "nstk",
1115 (int)propmat_clearsky_field.nbooks(),
1116 &stokes_dimid)))
1117 nca_error(retval, "nc_def_dim");
1118
1119 // Define variables
1120 if ((retval = nc_def_var(
1121 ncid, "wvlmin", NC_DOUBLE, 1, &none_dimid, &wvlmin_varid)))
1122 nca_error(retval, "nc_def_var wvlmin");
1123
1124 if ((retval = nc_def_var(
1125 ncid, "wvlmax", NC_DOUBLE, 1, &none_dimid, &wvlmax_varid)))
1126 nca_error(retval, "nc_def_var wvlmax");
1127
1128 if ((retval = nc_def_var(ncid, "z", NC_DOUBLE, 1, &nlev_dimid, &z_varid)))
1129 nca_error(retval, "nc_def_var z");
1130
1131 if ((retval =
1132 nc_def_var(ncid, "wvl", NC_DOUBLE, 1, &nwvl_dimid, &wvl_varid)))
1133 nca_error(retval, "nc_def_var wvl");
1134
1135 dimids[0] = nlyr_dimid;
1136 dimids[1] = nwvl_dimid;
1137 dimids[2] = stokes_dimid;
1138 dimids[3] = stokes_dimid;
1139
1140 if ((retval =
1141 nc_def_var(ncid, "tau", NC_DOUBLE, 4, &dimids[0], &tau_varid)))
1142 nca_error(retval, "nc_def_var tau");
1143
1144 // Units
1145 if ((retval = nc_put_att_text(ncid, wvlmin_varid, "units", 2, "nm")))
1146 nca_error(retval, "nc_put_att_text");
1147
1148 if ((retval = nc_put_att_text(ncid, wvlmax_varid, "units", 2, "nm")))
1149 nca_error(retval, "nc_put_att_text");
1150
1151 if ((retval = nc_put_att_text(ncid, z_varid, "units", 2, "km")))
1152 nca_error(retval, "nc_put_att_text");
1153
1154 if ((retval = nc_put_att_text(ncid, wvl_varid, "units", 2, "nm")))
1155 nca_error(retval, "nc_put_att_text");
1156
1157 if ((retval = nc_put_att_text(ncid, tau_varid, "units", 1, "-")))
1158 nca_error(retval, "nc_put_att_text");
1159
1160 // End define mode. This tells netCDF we are done defining
1161 // metadata.
1162 if ((retval = nc_enddef(ncid))) nca_error(retval, "nc_enddef");
1163
1164 // Assign data
1165 double wvlmin[1];
1166 wvlmin[0] = SPEED_OF_LIGHT / f_grid[f_grid.nelem() - 1] * 1e9;
1167 if ((retval = nc_put_var_double(ncid, wvlmin_varid, &wvlmin[0])))
1168 nca_error(retval, "nc_put_var");
1169
1170 double wvlmax[1];
1171 wvlmax[0] = SPEED_OF_LIGHT / f_grid[0] * 1e9;
1172 if ((retval = nc_put_var_double(ncid, wvlmax_varid, &wvlmax[0])))
1173 nca_error(retval, "nc_put_var");
1174
1175 double z[z_field.npages()];
1176 for (int iz = 0; iz < z_field.npages(); iz++)
1177 z[iz] = z_field(z_field.npages() - 1 - iz, 0, 0) * 1e-3;
1178
1179 if ((retval = nc_put_var_double(ncid, z_varid, &z[0])))
1180 nca_error(retval, "nc_put_var");
1181
1182 double wvl[f_grid.nelem()];
1183 for (int iv = 0; iv < f_grid.nelem(); iv++)
1184 wvl[iv] = SPEED_OF_LIGHT / f_grid[f_grid.nelem() - 1 - iv] * 1e9;
1185
1186 if ((retval = nc_put_var_double(ncid, wvl_varid, &wvl[0])))
1187 nca_error(retval, "nc_put_var");
1188
1189 const Index zfnp = z_field.npages() - 1;
1190 const Index fgne = f_grid.nelem();
1191 const Index amfnb = propmat_clearsky_field.nbooks();
1192
1193 Tensor4 tau(zfnp, fgne, amfnb, amfnb, 0.);
1194
1195 // Calculate average tau for layers
1196 for (int is = 0; is < propmat_clearsky_field.nlibraries(); is++)
1197 for (int iz = 0; iz < zfnp; iz++)
1198 for (int iv = 0; iv < fgne; iv++)
1199 for (int is1 = 0; is1 < amfnb; is1++)
1200 for (int is2 = 0; is2 < amfnb; is2++)
1201 // sum up all species
1202 tau(iz, iv, is1, is2) +=
1203 0.5 *
1204 (propmat_clearsky_field(is,
1205 f_grid.nelem() - 1 - iv,
1206 is1,
1207 is2,
1208 z_field.npages() - 1 - iz,
1209 0,
1210 0) +
1211 propmat_clearsky_field(is,
1212 f_grid.nelem() - 1 - iv,
1213 is1,
1214 is2,
1215 z_field.npages() - 2 - iz,
1216 0,
1217 0)) *
1218 (z_field(z_field.npages() - 1 - iz, 0, 0) -
1219 z_field(z_field.npages() - 2 - iz, 0, 0));
1220
1221 if ((retval = nc_put_var_double(ncid, tau_varid, tau.get_c_array())))
1222 nca_error(retval, "nc_put_var");
1223
1224 // Close the file
1225 if ((retval = nc_close(ncid))) nca_error(retval, "nc_close");
1226 }
1227}
1228
1229#else
1230
1231void WriteMolTau( //WS Input
1232 const Vector& f_grid _U_,
1233 const Tensor3& z_field _U_,
1234 const Tensor7& propmat_clearsky_field _U_,
1235 const Index& atmosphere_dim _U_,
1236 //Keyword
1237 const String& filename _U_,
1238 const Verbosity&) {
1239 ARTS_USER_ERROR_IF(true,
1240 "The workspace method WriteMolTau is not available"
1241 "because ARTS was compiled without NetCDF support.");
1242}
1243
1244#endif /* ENABLE_NETCDF */
1245
1246void propmat_clearsky_agendaAuto(// Workspace reference:
1247 Workspace& ws,
1248 // WS Output:
1249 Agenda& propmat_clearsky_agenda,
1250 Index& propmat_clearsky_agenda_checked,
1251 // WS Input:
1252 const ArrayOfArrayOfSpeciesTag& abs_species,
1253 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
1254 // WS Generic Input:
1255 const Numeric& H,
1256 const Numeric& T_extrapolfac,
1257 const Numeric& eta,
1258 const Numeric& extpolfac,
1259 const Numeric& force_p,
1260 const Numeric& force_t,
1261 const Index& ignore_errors,
1262 const Numeric& lines_sparse_df,
1263 const Numeric& lines_sparse_lim,
1264 const String& lines_speedup_option,
1265 const Index& manual_mag_field,
1266 const Index& no_negatives,
1267 const Numeric& theta,
1268 const Index& use_abs_as_ext,
1269 const Index& use_abs_lookup_ind,
1270 // Verbosity object:
1271 const Verbosity& verbosity) {
1272 using namespace AgendaManip;
1273
1274 propmat_clearsky_agenda_checked = 0; // In case of crash
1275
1276 AgendaCreator agenda(ws, "propmat_clearsky_agenda");
1277
1278 // Use bool because logic is easier
1279 const bool use_abs_lookup = static_cast<bool>(use_abs_lookup_ind);
1280
1281 const SpeciesTagTypeStatus any_species(abs_species);
1282 const AbsorptionTagTypesStatus any_lines(abs_lines_per_species);
1283
1284 // propmat_clearskyInit
1285 agenda.add("propmat_clearskyInit");
1286
1287 // propmat_clearskyAddFromLookup
1288 if (use_abs_lookup) {
1289 agenda.add("propmat_clearskyAddFromLookup",
1290 SetWsv{"extpolfac", extpolfac},
1291 SetWsv{"no_negatives", no_negatives});
1292 }
1293
1294 // propmat_clearskyAddLines
1295 if (not use_abs_lookup and any_species.Plain and
1296 (any_lines.population.LTE or any_lines.population.NLTE or
1297 any_lines.population.VibTemps)) {
1298 agenda.add("propmat_clearskyAddLines",
1299 SetWsv{"lines_sparse_df", lines_sparse_df},
1300 SetWsv{"lines_sparse_lim", lines_sparse_lim},
1301 SetWsv{"lines_speedup_option", lines_speedup_option},
1302 SetWsv{"no_negatives", no_negatives});
1303 }
1304
1305 // propmat_clearskyAddZeeman
1306 if (any_species.Zeeman and
1307 (any_lines.population.LTE or any_lines.population.NLTE or
1308 any_lines.population.VibTemps)) {
1309 agenda.add("propmat_clearskyAddZeeman",
1310 SetWsv{"manual_mag_field", manual_mag_field},
1311 SetWsv{"H", H},
1312 SetWsv{"theta", theta},
1313 SetWsv{"eta", eta});
1314 }
1315
1316 //propmat_clearskyAddHitranXsec
1317 if (not use_abs_lookup and any_species.XsecFit) {
1318 agenda.add("propmat_clearskyAddXsecFit",
1319 SetWsv{"force_p", force_p},
1320 SetWsv{"force_t", force_t});
1321 }
1322
1323 //propmat_clearskyAddOnTheFlyLineMixing
1324 if (not use_abs_lookup and any_species.Plain and
1325 (any_lines.population.ByMakarovFullRelmat or
1327 agenda.add("propmat_clearskyAddOnTheFlyLineMixing");
1328 }
1329
1330 //propmat_clearskyAddOnTheFlyLineMixingWithZeeman
1331 if (any_species.Zeeman and
1332 (any_lines.population.ByMakarovFullRelmat or
1334 agenda.add("propmat_clearskyAddOnTheFlyLineMixingWithZeeman");
1335 }
1336
1337 //propmat_clearskyAddCIA
1338 if (not use_abs_lookup and any_species.Cia) {
1339 agenda.add("propmat_clearskyAddCIA",
1340 SetWsv{"T_extrapolfac", T_extrapolfac},
1341 SetWsv{"ignore_errors", ignore_errors});
1342 }
1343
1344 //propmat_clearskyAddPredefined
1345 if (not use_abs_lookup and any_species.Predefined) {
1346 agenda.add("propmat_clearskyAddPredefined");
1347 }
1348
1349 //propmat_clearskyAddParticles
1350 if (any_species.Particles) {
1351 agenda.add("propmat_clearskyAddParticles",
1352 SetWsv{"use_abs_as_ext", use_abs_as_ext});
1353 }
1354
1355 //propmat_clearskyAddFaraday
1356 if (any_species.FreeElectrons) {
1357 agenda.add("propmat_clearskyAddFaraday");
1358 }
1359
1360 // propmat_clearskyAddHitranLineMixingLines
1361 if (not use_abs_lookup and any_species.Plain and
1362 (any_lines.population.ByHITRANFullRelmat or
1364 agenda.add("propmat_clearskyAddHitranLineMixingLines");
1365 }
1366
1367 // Extra check (should really never ever fail when species exist)
1368 propmat_clearsky_agenda = agenda.finalize();
1369 propmat_clearsky_agenda_checked = 1;
1370
1372 if (out3.sufficient_priority()) {
1373 out3 << "propmat_clearsky_agendaAuto sets propmat_clearsky_agenda to:\n\n"
1374 << propmat_clearsky_agenda << '\n';
1375 }
1376}
Declarations required for the calculation of absorption coefficients.
AbsorptionSpeciesBandIndex flat_index(Index i, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species)
Get a flat index pair for species and band.
Contains the absorption namespace.
Declarations for agendas.
This file contains the definition of Array.
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:241
The global header file for ARTS.
Constants of physical expressions as constexpr.
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:46
int arts_omp_get_thread_num()
Wrapper for omp_get_thread_num.
Definition: arts_omp.cc:76
bool arts_omp_in_parallel()
Wrapper for omp_in_parallel.
Definition: arts_omp.cc:62
Header file for helper functions for OpenMP.
Stuff related to time in ARTS.
The Agenda class.
Definition: agenda_class.h:69
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
MatrixView real()
Get a view of the real part of the matrix.
VectorView real()
Get a view of the real part of the vector.
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index nlibraries() const noexcept
Definition: matpackVII.h:158
Index nbooks() const noexcept
Definition: matpackVII.h:161
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
void AddAbsorptionVectorAtPosition(const ConstVectorView &x, const Index iv=0, const Index iz=0, const Index ia=0)
Adds as a Stokes vector at position.
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
void SetZero()
Sets all data to zero.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
void AddFaraday(const Numeric &rot, const Index iv=0, const Index iz=0, const Index ia=0)
Adds the Faraday rotation to the PropagationMatrix at required position.
void SetAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Set the At Position object.
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:352
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackIV.cc:342
The Tensor4 class.
Definition: matpackIV.h:435
The Tensor7 class.
Definition: matpackVII.h:2407
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Workspace class.
Definition: workspace_ng.h:53
#define _U_
Definition: config.h:180
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
void find_xml_file(String &filename, const Verbosity &verbosity)
Find an xml file.
Definition: file.cc:355
This file contains basic functions to handle ASCII files.
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1175
bool do_magnetic_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a magnetic derivative.
Definition: jacobian.cc:1179
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1150
Numeric magnetic_field_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the magnetic field perturbation if it exists.
Definition: jacobian.cc:1199
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
Routines for setting up the jacobian.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:218
void sparse_f_gridFromFrequencyGrid(Vector &sparse_f_grid, const Vector &f_grid, const Numeric &sparse_df, const String &speedup_option, const Verbosity &)
WORKSPACE METHOD: sparse_f_gridFromFrequencyGrid.
Definition: m_abs.cc:803
constexpr Numeric ELECTRON_MASS
Definition: m_abs.cc:77
void isotopologue_ratiosInitFromHitran(SpeciesIsotopologueRatios &isotopologue_ratios, const Verbosity &)
WORKSPACE METHOD: isotopologue_ratiosInitFromHitran.
Definition: m_abs.cc:1070
constexpr Numeric SPEED_OF_LIGHT
Definition: m_abs.cc:79
Vector create_sparse_f_grid_internal(const Vector &f_grid, const Numeric &sparse_df, const String &speedup_option, const Verbosity &verbosity)
Definition: m_abs.cc:831
void propmat_clearskyInit(PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Index &stokes_dim, const Index &propmat_clearsky_agenda_checked, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyInit.
Definition: m_abs.cc:352
void nlte_sourceFromTemperatureAndSrcCoefPerSpecies(StokesVector &nlte_source, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfMatrix &src_coef_per_species, const ArrayOfMatrix &dsrc_coef_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Numeric &rtp_temperature, const Verbosity &)
Definition: m_abs.cc:274
void abs_speciesDefineAllInScenario(ArrayOfArrayOfSpeciesTag &tgs, Index &propmat_clearsky_agenda_checked, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesDefineAllInScenario.
Definition: m_abs.cc:177
constexpr Numeric VACUUM_PERMITTIVITY
Definition: m_abs.cc:80
void abs_speciesDefineAll(ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesDefineAll.
Definition: m_abs.cc:225
void AbsInputFromRteScalars(Vector &abs_p, Vector &abs_t, Matrix &abs_vmrs, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const Verbosity &)
Definition: m_abs.cc:83
void propmat_clearskyForceNegativeToZero(PropagationMatrix &propmat_clearsky, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyForceNegativeToZero.
Definition: m_abs.cc:1057
void propmat_clearskyAddLines(PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfSpeciesTag &select_abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const SpeciesIsotopologueRatios &isotopologue_ratios, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Index &nlte_do, const Index &lbl_checked, const Numeric &sparse_df, const Numeric &sparse_lim, const String &speedup_option, const Index &robust, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddLines.
Definition: m_abs.cc:843
void propmat_clearskyZero(PropagationMatrix &propmat_clearsky, const Vector &f_grid, const Index &stokes_dim, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyZero.
Definition: m_abs.cc:1049
void propmat_clearskyAddParticles(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfSpeciesTag &select_abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &rtp_vmr, const Vector &rtp_los, const Numeric &rtp_temperature, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &scat_data_checked, const Index &use_abs_as_ext, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddParticles.
Definition: m_abs.cc:535
constexpr Numeric PI
Definition: m_abs.cc:78
void propmat_clearskyAddFaraday(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfSpeciesTag &select_abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &rtp_vmr, const Vector &rtp_los, const Vector &rtp_mag, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyAddFaraday.
Definition: m_abs.cc:422
void WriteMolTau(const Vector &f_grid, const Tensor3 &z_field, const Tensor7 &propmat_clearsky_field, const Index &atmosphere_dim, const String &filename, const Verbosity &)
WORKSPACE METHOD: WriteMolTau.
Definition: m_abs.cc:1078
void AbsInputFromAtmFields(Vector &abs_p, Vector &abs_t, Matrix &abs_vmrs, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const Tensor4 &vmr_field, const Verbosity &)
WORKSPACE METHOD: AbsInputFromAtmFields.
Definition: m_abs.cc:248
void abs_lines_per_speciesCreateFromLines(ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const ArrayOfAbsorptionLines &abs_lines, const ArrayOfArrayOfSpeciesTag &abs_species, const Verbosity &)
WORKSPACE METHOD: abs_lines_per_speciesCreateFromLines.
Definition: m_abs.cc:106
void isotopologue_ratiosInitFromBuiltin(SpeciesIsotopologueRatios &isotopologue_ratios, const Verbosity &)
WORKSPACE METHOD: isotopologue_ratiosInitFromBuiltin.
Definition: m_abs.cc:1064
constexpr Numeric ELECTRON_CHARGE
Definition: m_abs.cc:76
void propmat_clearsky_agendaAuto(Workspace &ws, Agenda &propmat_clearsky_agenda, Index &propmat_clearsky_agenda_checked, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const Numeric &H, const Numeric &T_extrapolfac, const Numeric &eta, const Numeric &extpolfac, const Numeric &force_p, const Numeric &force_t, const Index &ignore_errors, const Numeric &lines_sparse_df, const Numeric &lines_sparse_lim, const String &lines_speedup_option, const Index &manual_mag_field, const Index &no_negatives, const Numeric &theta, const Index &use_abs_as_ext, const Index &use_abs_lookup_ind, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearsky_agendaAuto.
Definition: m_abs.cc:1246
void abs_speciesSet(ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, const ArrayOfString &names, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesSet.
Workspace methods and template functions for supergeneric XML IO.
constexpr bool any_negative(const MatpackType &var) noexcept
Checks for negative values.
Definition: math_funcs.h:150
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition: messages.h:204
#define CREATE_OUT3
Definition: messages.h:206
#define CREATE_OUT2
Definition: messages.h:205
Declaration of the class MdRecord.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric vacuum_permittivity
Vacuum permittivity [F/m].
constexpr Numeric electron_mass
Mass of resting electron [kg].
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Numeric elementary_charge
Elementary charge [C] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 20...
SpeciesIsotopologueRatios isotopologue_ratios()
void compute(ComputeData &com, ComputeData &sparse_com, const AbsorptionLines &band, const ArrayOfRetrievalQuantity &jacobian_quantities, const EnergyLevelMap &nlte, const Vector &vmrs, const ArrayOfSpeciesTag &self_tag, const Numeric &self_vmr, const Numeric &isot_ratio, const Numeric &P, const Numeric &T, const Numeric &H, const Numeric &sparse_lim, const Zeeman::Polarization zeeman_polarization, const Options::LblSpeedup speedup_type, const bool robust) ARTS_NOEXCEPT
Compute the absorption of an absorption band.
Definition: lineshape.cc:3549
bool good_linear_sparse_f_grid(const Vector &f_grid_dense, const Vector &f_grid_sparse) noexcept
Definition: lineshape.cc:3662
Vector triple_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) noexcept
Definition: lineshape.cc:3677
Vector linear_sparse_f_grid(const Vector &f_grid, const Numeric &sparse_df) ARTS_NOEXCEPT
Definition: lineshape.cc:3639
constexpr IsotopologueRatios isotopologue_ratiosInitFromBuiltin()
Implements Zeeman modeling.
Definition: zeemandata.cc:333
void nca_error(const int e, const std::string_view s)
Throws a runtime error for the given NetCDF error code.
Definition: nc_io.cc:643
This file contains basic functions to handle NetCDF data files.
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
Scattering database structure and functions.
This file contains header information for the dealing with command line parameters.
Numeric dplanck_df(const Numeric &f, const Numeric &t)
dplanck_df
Numeric planck(const Numeric &f, const Numeric &t)
planck
Numeric dplanck_dt(const Numeric &f, const Numeric &t)
dplanck_dt
This file contains declerations of functions of physical character.
Array< StokesVector > ArrayOfStokesVector
Array< PropagationMatrix > ArrayOfPropagationMatrix
bool bad_propmat(const Array< T > &main, const Vector &f_grid, const Index sd=1) noexcept
Checks if a Propagation Matrix or something similar has good grids.
Numeric dotprod_with_los(const ConstVectorView &los, const Numeric &u, const Numeric &v, const Numeric &w, const Index &atmosphere_dim)
Calculates the dot product between a field and a LOS.
Definition: rte.cc:670
void mirror_los(Vector &los_mirrored, const ConstVectorView &los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:1811
Declaration of functions in rte.cc.
void check_abs_species(const ArrayOfArrayOfSpeciesTag &abs_species)
AbsorptionPopulationTagTypeStatus population
SingleLine PopLine(Index) noexcept
Pops a single line.
Array< SingleLine > lines
A list of individual lines.
Index NumLines() const noexcept
Number of lines.
Species::IsotopeRecord Isotopologue() const noexcept
Isotopologue Index.
void ReverseLines() noexcept
Reverses the order of the internal lines.
void AppendSingleLine(SingleLine &&sl)
Appends a single line to the absorption lines.
Helper class to create an agenda.
Definition: agenda_set.h:144
void add(const std::string_view method, Input &&... input)
Add a method with as many inputs as you want. These inputs must be of type Wsv.
Definition: agenda_set.h:155
Agenda finalize()
Check the agenda and return a copy of it (ignoring/touching all agenda input/output not dealt with)
Definition: agenda_set.cc:177
Main computational data for the line shape and strength calculations.
Definition: lineshape.h:840
void interp_add_even(const ComputeData &sparse) ARTS_NOEXCEPT
Add a sparse grid to this grid via linear interpolation.
Definition: lineshape.cc:3703
ComplexVector N
Definition: lineshape.h:841
void interp_add_triplequad(const ComputeData &sparse) ARTS_NOEXCEPT
Add a sparse grid to this grid via square interpolation.
Definition: lineshape.cc:3748
ComplexVector F
Definition: lineshape.h:841
ComplexMatrix dF
Definition: lineshape.h:842
ComplexMatrix dN
Definition: lineshape.h:842
Struct to test of an ArrayOfArrayOfSpeciesTag contains a tagtype.
Definition: species_tags.h:172
#define a
#define b
This file contains basic functions to handle XML data files.