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