ARTS 2.5.0 (git: 9ee3ac6c)
test_propagationmatrix.cc
Go to the documentation of this file.
1/* Copyright (C) 2015 Richard Larsson <ric.larsson@gmail.com>
2 *
3 * This program is free software; you can redistribute it and/or modify it
4 * under the terms of the GNU General Public License as published by the
5 * Free Software Foundation; either version 2, or (at your option) any
6 * later version.
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with this program; if not, write to the Free Software
15 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 * USA. */
17
26#include <random>
27#include "absorption.h"
28#include "arts.h"
29#include "global_data.h"
30#include "hitran_species.h"
31#include "lineshapemodel.h"
32#include "linescaling.h"
33#include "transmissionmatrix.h"
34#include "zeeman.h"
35#include "zeemandata.h"
36#include <Faddeeva/Faddeeva.hh>
37#include "legacy_continua.h"
39#include "wigner_functions.h"
40
41#include "linemixing_hitran.h"
42#include <auto_md.h>
43
45 const Numeric k11 = 1;
46 const Numeric k12 = -0.51;
47 const Numeric k13 = -0.21;
48 const Numeric k14 = 0.31;
49 const Numeric k23 = -0.1;
50 const Numeric k24 = -0.99;
51 const Numeric k34 = 2;
52
53 const Numeric r = 0.5;
54
55 const Numeric a = -k11 * r;
56 const Numeric b = -k12 * r;
57 const Numeric c = -k13 * r;
58 const Numeric d = -k14 * r;
59 const Numeric u = -k23 * r;
60 const Numeric v = -k24 * r;
61 const Numeric w = -k34 * r;
62
63 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
64 w2 = w * w;
65
66 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
67
68 Numeric Const1;
69 Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
70 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
71 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
72 Const1 += u2 * (u2 * 0.5 + v2 + w2);
73 Const1 += v2 * (v2 * 0.5 + w2);
74 Const1 *= 2;
75 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
76 Const1 += w2 * w2;
77
78 if (Const1 > 0.0)
79 Const1 = sqrt(Const1);
80 else
81 Const1 = 0.0;
82
83 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
84 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
85 const Numeric x = sqrt_BpA.real() * sqrt(0.5);
86 const Numeric y = sqrt_BmA.imag() * sqrt(0.5);
87 const Numeric x2 = x * x;
88 const Numeric y2 = y * y;
89 const Numeric cos_y = cos(y);
90 const Numeric sin_y = sin(y);
91 const Numeric cosh_x = cosh(x);
92 const Numeric sinh_x = sinh(x);
93 const Numeric x2y2 = x2 + y2;
94 const Numeric inv_x2y2 = 1.0 / x2y2;
95
96 std::cout << x << " " << y << " " << Const1 << " " << Const2 << "\n";
97
98 Numeric C0, C1, C2, C3;
99 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
100
101 // X and Y cannot both be zero
102 if (x == 0.0) {
103 inv_y = 1.0 / y;
104 C0 = 1.0;
105 C1 = 1.0;
106 C2 = (1.0 - cos_y) * inv_x2y2;
107 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
108 } else if (y == 0.0) {
109 inv_x = 1.0 / x;
110 C0 = 1.0;
111 C1 = 1.0;
112 C2 = (cosh_x - 1.0) * inv_x2y2;
113 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
114 } else {
115 inv_x = 1.0 / x;
116 inv_y = 1.0 / y;
117
118 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
119 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
120 C2 = (cosh_x - cos_y) * inv_x2y2;
121 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
122 }
123
124 std::cout << C0 << " " << C1 << " " << C2 << " " << C3 << "\n";
125
126 Matrix F(4, 4, 0), A(4, 4, 0);
127
128 MatrixViewMap eigF = MapToEigen(F);
129 Eigen::Matrix4d eigA;
130 eigA << 0, b, c, d, b, 0, u, v, c, -u, 0, w, d, -v, -w, 0;
131
132 eigF = C1 * eigA + C2 * eigA * eigA + C3 * eigA * eigA * eigA;
133 eigF(0, 0) += C0;
134 eigF(1, 1) += C0;
135 eigF(2, 2) += C0;
136 eigF(3, 3) += C0;
137 eigF *= exp(a);
138
139 std::cout << F << "\n";
140}
141
143 // Initializes as unity matrices
144 TransmissionMatrix a(2, 4);
145 std::cout << "Initialized TransmissionMatrix(2, 4):\n" << a << "\n";
146
147 // Set a single input
148 Eigen::Matrix4d A;
149 A << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;
150 std::cout << "New Matrix:\n" << A << "\n\n";
151 a.Mat4(0) = A;
152 std::cout << "Updated TransmissionMatrix Position 1 wit New Matrix:\n"
153 << a << "\n";
154
155 // The stream can also set the values
156 String S =
157 "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 125 26 27 28 29 30 31 32";
158 std::cout << "Stream:\n" << S << "\n\n";
159 std::istringstream astream(S);
160 astream >> a;
161 std::cout << "Streamed into TransmissionMatrix:\n" << a << "\n";
162
163 // Initialize empty
164 RadiationVector b(3, 3);
165 std::cout << "Initialized RadiationVector(3, 3)\n" << b << "\n";
166
167 // Set is not defined but add is
168 Eigen::Vector3d B;
169 B << 1, 2, 3;
170 std::cout << "New Vector:\n" << B << "\n\n"; // nb. not transposed
171 b.Vec3(1).noalias() += B;
172 std::cout << "Updated RadiationVector Position 1 with New Vector:\n"
173 << b << "\n";
174
175 // The stream can also set the values
176 String T = "1 2 3 4 5 6 7 8 90";
177 std::cout << "Stream:\n" << T << "\n\n";
178 std::istringstream bstream(T);
179 bstream >> b;
180 std::cout << "Streamed into RadiationVector:\n" << b << "\n";
181}
182
184 auto f = [](Numeric x) { return 0.1 * x; };
185 auto df = [](Numeric x) { return 0.1 + 0 * x; };
186
187 Index nstokes = 4;
188 const Numeric x1 = 30;
189 const Numeric x2 = -0.1;
190 const Numeric r_normal = 1000;
191 const Numeric r_extra1 = r_normal + f(x1);
192 const Numeric r_extra2 = r_normal + f(x2);
193
194 PropagationMatrix a(1, nstokes);
195 a.Kjj() = 10 + Numeric(rand() % 1000) / 100;
196 if (nstokes > 1) {
197 a.K12() = 2 + Numeric(rand() % 1000) / 100;
198 if (nstokes > 2) {
199 a.K13() = 3 + Numeric(rand() % 1000) / 100;
200 a.K23() = -1 + Numeric(rand() % 1000) / 100;
201 if (nstokes > 3) {
202 a.K14() = 5 + Numeric(rand() % 1000) / 100;
203 a.K24() = -3 + Numeric(rand() % 1000) / 100;
204 a.K34() = -2 + Numeric(rand() % 1000) / 100;
205 }
206 }
207 }
208 a.Data() *= 1e-5;
209
210 PropagationMatrix b(1, nstokes);
211 b.Kjj() = 5 + Numeric(rand() % 1000) / 100;
212 if (nstokes > 1) {
213 b.K12() = -1 + Numeric(rand() % 1000) / 100;
214 if (nstokes > 2) {
215 b.K13() = -3 + Numeric(rand() % 1000) / 100;
216 b.K23() = 4 + Numeric(rand() % 1000) / 100;
217 if (nstokes > 3) {
218 b.K14() = 2 + Numeric(rand() % 1000) / 100;
219 b.K24() = -1 + Numeric(rand() % 1000) / 100;
220 b.K34() = 3 + Numeric(rand() % 1000) / 100;
221 }
222 }
223 }
224 b.Data() *= 5e-6;
225
227 da[0].Data() = 0;
228 Tensor3 T_normal(1, nstokes, nstokes), T_extra(1, nstokes, nstokes);
229 Tensor4 dT1(1, 1, nstokes, nstokes), dT2(1, 1, nstokes, nstokes);
230
232 T_normal, dT1, dT2, r_normal, a, b, da, da, df(x1), df(x2), 0);
233
234 std::cout << "Transmission at r=" << r_normal << ":\n"
235 << MapToEigen(T_normal(0, joker, joker)) << "\n"
236 << "\n";
237 std::cout << "First derivative:\n"
238 << MapToEigen(dT1(0, 0, joker, joker)) << "\n"
239 << "\n";
240 std::cout << "Second derivative:\n"
241 << MapToEigen(dT2(0, 0, joker, joker)) << "\n"
242 << "\n";
243
244 compute_transmission_matrix(T_extra, r_extra1, a, b);
245
246 std::cout << "Transmission at perturbed r1=" << r_extra1 << ":\n"
247 << MapToEigen(T_extra(0, joker, joker)) << "\n"
248 << "\n";
249 T_extra -= T_normal;
250 T_extra /= x1;
251 std::cout << "First derivative perturbed:\n"
252 << MapToEigen(T_extra(0, joker, joker)) << "\n"
253 << "\n";
254 T_extra /= dT1(0, joker, joker, joker);
255 std::cout << "First derivative perturbed relative:\n"
256 << MapToEigen(T_extra(0, joker, joker)) << "\n"
257 << "\n";
258
259 compute_transmission_matrix(T_extra, r_extra2, a, b);
260
261 std::cout << "Transmission at perturbed r2=" << r_extra2 << ":\n"
262 << MapToEigen(T_extra(0, joker, joker)) << "\n"
263 << "\n";
264 T_extra -= T_normal;
265 T_extra /= x2;
266 std::cout << "Second derivative perturbed:\n"
267 << MapToEigen(T_extra(0, joker, joker)) << "\n"
268 << "\n";
269 T_extra /= dT2(0, joker, joker, joker);
270 std::cout << "Second derivative perturbed relative:\n"
271 << MapToEigen(T_extra(0, joker, joker)) << "\n"
272 << "\n";
273}
274
276 const Numeric a = 2;
277 const Numeric b = 3;
278 const Numeric c = 4;
279 const Numeric d = 1;
280 const Numeric u = 5;
281 const Numeric v = 1;
282 const Numeric w = 5;
283
288 test1.Kjj() = a;
289 test2.Kjj() = a;
290 test3.Kjj() = a;
291 test4.Kjj() = a;
292 test2.K12() = b;
293 test3.K12() = b;
294 test4.K12() = b;
295 test3.K13() = c;
296 test4.K13() = c;
297 test3.K23() = u;
298 test4.K23() = u;
299 test4.K14() = d;
300 test4.K24() = v;
301 test4.K34() = w;
302
303 std::cout << test1 << "\n\n"
304 << test2 << "\n\n"
305 << test3 << "\n\n"
306 << test4 << "\n\n";
307
308 TransmissionMatrix ans1(1, 1);
309 TransmissionMatrix ans2(1, 2);
310 TransmissionMatrix ans3(1, 3);
311 TransmissionMatrix ans4(1, 4);
312
314
316 empty,
317 empty,
318 test1,
319 test1,
322 1,
323 0,
324 0,
325 -1);
326
328 empty,
329 empty,
330 test2,
331 test2,
334 1,
335 0,
336 0,
337 -1);
338
340 empty,
341 empty,
342 test3,
343 test3,
346 1,
347 0,
348 0,
349 -1);
350
352 empty,
353 empty,
354 test4,
355 test4,
358 1,
359 0,
360 0,
361 -1);
362
363 std::cout << ans1 << "\n\n"
364 << ans2 << "\n\n"
365 << ans3 << "\n\n"
366 << ans4 << "\n\n";
367}
368
370 int i = 1;
372 for (auto& pm : propmats) {
373 pm.K12() = i;
374 i++;
375 pm.K13() = i;
376 i++;
377 pm.K23() = i;
378 i++;
379 pm.K14() = i;
380 i++;
381 pm.K24() = i;
382 i++;
383 pm.K34() = i;
384 i++;
385 pm.Kjj() = 2 * i;
386 i++;
387 pm.Data() *= 1e-2;
388 }
389
390 std::cout << "Propmats:\n" << propmats << "\n\n";
391
394 for (i = 0; i < 4; i++) {
395 stepwise_transmission(layers[i + 1],
396 empty,
397 empty,
398 propmats[i],
399 propmats[i + 1],
402 1,
403 0,
404 0,
405 -1);
406 }
407
408 std::cout << "Layers:\n" << layers << "\n\n";
409
410 ArrayOfTransmissionMatrix cumulative_forward =
412 ArrayOfTransmissionMatrix cumulative_reflect =
414
415 std::cout << "Forward accumulation:\n" << cumulative_forward << "\n\n";
416 std::cout << "Reflect accumulation:\n" << cumulative_reflect << "\n\n";
417}
418
420 Numeric start = 1.0;
421 Numeric end = 1e-7;
422 int n = 10000;
423
424 Vector x(n), sx(n), shx(n), cx(n), chx(n);
425 nlogspace(x, start, end, n);
426
427 for (int i = 0; i < n; i++) {
428 sx[i] = std::sin(x[i]);
429 cx[i] = std::cos(x[i]);
430 shx[i] = std::sinh(x[i]);
431 chx[i] = std::cosh(x[i]);
432 }
433
434 std::cout
435 << std::scientific << std::setprecision(15)
436 << "x\tabs(sx/x-1)\tabs(shx/x-1)\tabs((chx-cx)/2x^2-1/2)\tabs((shx/x-sx/x)/2x^2-1/6)\n";
437
438 for (int i = 0; i < n; i++)
439 std::cout << x[i] << '\t' << std::abs(sx[i] / x[i] - 1.0) << '\t'
440 << std::abs(shx[i] / x[i] - 1.0) << '\t'
441 << std::abs((chx[i] - cx[i]) / (2 * x[i] * x[i]) - 0.5) << '\t'
442 << std::abs((shx[i] / x[i] - sx[i] / x[i]) / (2 * x[i] * x[i]) -
443 1.0 / 6.0)
444 << '\n';
445}
446
448 auto o266 = SpeciesTag("O2-66");
449 auto o268 = SpeciesTag("O2-68");
450
451 Numeric g;
454 qn.Set(QuantumNumberType::Lambda, 0);
455 qn.Set(QuantumNumberType::v1, 0);
456 qn.Set(QuantumNumberType::S, 1);
457
458 std::cout << "Table from Larsson, Lankhaar, Eriksson (2019)\n";
459 for (Index i = 1; i < 51; i++) {
460 qn.Set(QuantumNumberType::J, Rational(i));
461
462 qn.Set(QuantumNumberType::N, qn[QuantumNumberType::J] - 1);
463 std::cout << i << "_" << i - 1;
464
466 QuantumIdentifier(o266.Isotopologue(), qn, qn))
467 .gl();
468 std::cout << '\t' << g;
469
471 QuantumIdentifier(o268.Isotopologue(), qn, qn))
472 .gl();
473 std::cout << '\t' << g;
474
476 QuantumIdentifier(o266.Isotopologue(), qn, qn))
477 .gl();
478 std::cout << '\t' << g;
479
480 qn.Set(QuantumNumberType::N, qn[QuantumNumberType::J]);
481 std::cout << '\t' << i << "_" << i;
482
484 QuantumIdentifier(o266.Isotopologue(), qn, qn))
485 .gl();
486 std::cout << '\t' << g;
487
489 QuantumIdentifier(o268.Isotopologue(), qn, qn))
490 .gl();
491 std::cout << '\t' << g;
492
494 QuantumIdentifier(o266.Isotopologue(), qn, qn))
495 .gl();
496 std::cout << '\t' << g;
497
498 qn.Set(QuantumNumberType::N, qn[QuantumNumberType::J] + 1);
499 std::cout << '\t' << i << "_" << i + 1;
500
502 QuantumIdentifier(o266.Isotopologue(), qn, qn))
503 .gl();
504 std::cout << '\t' << g;
505
507 QuantumIdentifier(o268.Isotopologue(), qn, qn))
508 .gl();
509 std::cout << '\t' << g;
510
512 QuantumIdentifier(o266.Isotopologue(), qn, qn))
513 .gl();
514 std::cout << '\t' << g;
515
516 std::cout << '\n';
517 }
518}
519
520constexpr bool test_quantum_numbers(const QuantumNumbers qns, const Index i)
521{
522 return (i > 0) ? (qns[i].isUndefined() ? test_quantum_numbers(qns, i-1) : false) : true;
523}
524
526{
527 static_assert(test_quantum_numbers(QuantumNumbers(), Index(QuantumNumberType::FINAL) - 1),
528 "Bad last entry in QuantumNumbers. Did you recently expand the list?");
529}
530
531
533{
534 // Test constants
535 constexpr Index nf = 501;
536 constexpr Numeric fstart = 25e9;
537 constexpr Numeric fend = 165e9;
538 constexpr Numeric t = 200;
539 constexpr Numeric p = 1e4;
540 Vector f(nf);
541 nlinspace(f, fstart, fend, nf);
542 PropagationMatrix xsec(nf, 1);
544
545
547 jacs[0].Target(Jacobian::Target(Jacobian::Atm::Temperature));
548 jacs[1].Target(Jacobian::Target(Jacobian::Atm::WindU));
549 Absorption::PredefinedModel::makarov2020_o2_lines_mpm(xsec, dxsec, f, p, t, 1, 0.5, jacs);
550
551 constexpr auto df = 1000;
552 constexpr auto dt = 0.1;
553 Matrix pxsec(nf, 1, 0);
554 Matrix pxsec_dt(nf, 1, 0);
555 Matrix pxsec_df(nf, 1, 0);
556 Vector f_pert = f;
557 f_pert += df;
558 PWR93O2AbsModel(pxsec, 0, 1, 1, 1, "user", "PWR98", f, {p}, {t}, {0.5}, {1}, Verbosity());
559 PWR93O2AbsModel(pxsec_dt, 0, 1, 1, 1, "user", "PWR98", f, {p}, {t+dt}, {0.5}, {1}, Verbosity());
560 PWR93O2AbsModel(pxsec_df, 0, 1, 1, 1, "user", "PWR98", f_pert, {p}, {t}, {0.5}, {1}, Verbosity());
561
562 std::cout<< "xr = np.array([";
563 for (Index i=0; i<nf; i++)
564 std::cout<<'['<<xsec(i, 0)<<','<<' '<<pxsec(i, 0) << ']' << ',' << ' ';
565 std::cout << ']' << ')' << '\n';
566
567 std::cout<< "dxr_dt = np.array([";
568 for (Index i=0; i<nf; i++)
569 std::cout<<'['<<dxsec[0](i, 0)<<','<<' '<<(pxsec_dt(i, 0)-pxsec(i, 0))/dt << ']' << ',' << ' ';
570 std::cout << ']' << ')' << '\n';
571
572 std::cout<< "dxr_df = np.array([";
573 for (Index i=0; i<nf; i++)
574 std::cout<<'['<<dxsec[1](i, 0)<<','<<' '<<(pxsec_df(i, 0)-pxsec(i, 0))/df << ']' << ',' << ' ';
575 std::cout << ']' << ')' << '\n';
576}
577
578void test_hitran2017(bool newtest = true)
579{
580 const Numeric p = 1.0;
581 const Numeric t = 296;
582 const Numeric xco2 = 1.5e-2;
583 const Numeric xh2o = 0;
584 const Numeric sigmin = 600;
585 const Numeric sigmax = 900;
586 const Numeric dsig = 0.005;
587 const Numeric stotmax = 0.1e-21;
588
589 const Index nsig = Index(((sigmax - sigmin) / dsig) + 0.5) + 1;
590 Vector invcm_grid(nsig);
591 Vector f_grid(nsig);
592 Numeric sigc = sigmin-dsig;
593 for (Index isig=0; isig<nsig; isig++) {
594 sigc += dsig;
595 invcm_grid[isig] = sigc;
596 f_grid[isig] = Conversion::kaycm2freq(sigc);
597 }
598
599 constexpr Index n=6;
600 auto types =
601 std::array<std::pair<lm_hitran_2017::ModeOfLineMixing, lm_hitran_2017::calctype>, 6>{
602 std::pair<lm_hitran_2017::ModeOfLineMixing, lm_hitran_2017::calctype>{lm_hitran_2017::ModeOfLineMixing::FullW, lm_hitran_2017::calctype::FullW},
603 {lm_hitran_2017::ModeOfLineMixing::VP_W, lm_hitran_2017::calctype::FullW},
604 {lm_hitran_2017::ModeOfLineMixing::VP, lm_hitran_2017::calctype::NoneVP},
605 {lm_hitran_2017::ModeOfLineMixing::VP_Y, lm_hitran_2017::calctype::NoneRosenkranz},
606 {lm_hitran_2017::ModeOfLineMixing::SDVP, lm_hitran_2017::calctype::SDVP},
607 {lm_hitran_2017::ModeOfLineMixing::SDVP_Y, lm_hitran_2017::calctype::SDRosenkranz}};
608
609 ArrayOfVector absorption(n);
610 make_wigner_ready(int(250), int(20000000), 6);
611
614 for (Index i=0;i<n; i++) {
615 auto type=types[i];
616
617 lm_hitran_2017::read(hitran, bands, Hitran::isotopologue_ratios(Hitran::Type::Newest), "data_new", -1, Conversion::kaycm2freq(sigmin), Conversion::kaycm2freq(sigmax), Conversion::kaycm_per_cmsquared2hz_per_msquared(stotmax), type.first);
618 Vector vmrs = {1-xco2/100-xh2o/100, xh2o/100, xco2/100};
619
620 if (not newtest)
621 absorption[i] = lm_hitran_2017::compute(p, t, xco2, xh2o, invcm_grid, stotmax, type.second);
622 else
623 absorption[i] = lm_hitran_2017::compute(hitran, bands, Hitran::isotopologue_ratios(Hitran::Type::Newest), Conversion::atm2pa(p), t, vmrs, f_grid);
624 }
625
626 for (Index i=0; i<nsig; i++) {
627 for (Index j=0; j<n; j++) {
628 std::cout<<absorption[j][i]<<' ';
629 }
630 std::cout<<'\n';
631 }
632}
633
634
635int main(int n, char **argc) {
636 /*test_speed_of_pressurebroadening();
637 test_transmissionmatrix();
638 test_r_deriv_propagationmatrix();
639 test_transmat_from_propmat();
640 test_transmat_to_cumulativetransmat();
641 test_sinc_likes_0limit();*/
642
643 if (n == 2 and String(argc[1]) == "new") {
644 std::cout<<"new test\n";
645 test_hitran2017(true);
646 }
647 else {
648 std::cout<<"old test\n";
649 test_hitran2017(false);
650 }
651 return 0;
652}
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
This can be used to make arrays out of anything.
Definition: array.h:107
The Matrix class.
Definition: matpackI.h:1225
Container class for Quantum Numbers.
Definition: quantum.h:112
constexpr void Set(Index qn, Rational r)
Set quantum number at position.
Definition: quantum.h:199
Radiation Vector for Stokes dimension 1-4.
Implements rational numbers to work with other ARTS types.
Definition: rational.h:52
The Tensor3 class.
Definition: matpackIII.h:339
The Tensor4 class.
Definition: matpackIV.h:421
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
The Vector class.
Definition: matpackI.h:876
constexpr Numeric & gl() noexcept
Returns the lower state g.
Definition: zeemandata.h:389
#define abs(x)
void PWR93O2AbsModel(MatrixView pxsec, const Numeric CCin, const Numeric CLin, const Numeric CWin, const Numeric COin, const String &model, const String &version, ConstVectorView f_grid, ConstVectorView abs_p, ConstVectorView abs_t, ConstVectorView vmrh2o, ConstVectorView vmr, const Verbosity &verbosity)
Oxygen complex at 60 GHz plus mm O2 lines plus O2 continuum.
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:745
Namespace and functions to deal with HITRAN linemixing.
Constains various line scaling functions.
Contains the line shape namespace.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:225
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Definition: math_funcs.cc:261
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:115
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
std::complex< Numeric > Complex
#define b2
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:287
void makarov2020_o2_lines_mpm(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const Vector &f, const Numeric &p, const Numeric &t, const Numeric &oxygen_vmr, const Numeric &water_vmr, const ArrayOfRetrievalQuantity &jacs)
Adds Makarov MPM2020 O2 absorption lines to the propagation matrix.
constexpr auto kaycm2freq(T x) -> decltype(x *(100 *c))
Conversion from Kayser wavenumber to Hz.
Definition: constants.h:382
constexpr auto atm2pa(T x) -> decltype(x *101 '325.0)
Conversion from Atm to Pa.
Definition: constants.h:430
constexpr auto kaycm_per_cmsquared2hz_per_msquared(T x) -> decltype(x *kaycm2freq(1e-4))
Conversion from cm-1 per molecule per cm^2 to Hz per molecule per m^2.
Definition: constants.h:472
SpeciesIsotopologueRatios isotopologue_ratios(Type type)
char Type type
Temperature
Definition: jacobian.h:57
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
Model GetAdvancedModel(const QuantumIdentifier &qid) ARTS_NOEXCEPT
Returns an advanced Zeeman model.
Definition: zeemandata.cc:103
Model GetSimpleModel(const QuantumIdentifier &qid) ARTS_NOEXCEPT
Returns a simple Zeeman model.
Definition: zeemandata.cc:32
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:109
void read(HitranRelaxationMatrixData &hitran, ArrayOfAbsorptionLines &bands, const SpeciesIsotopologueRatios &isotopologue_ratio, const String &basedir, const Numeric linemixinglimit, const Numeric fmin, const Numeric fmax, const Numeric stot, const ModeOfLineMixing mode)
Read from HITRAN online line mixing file.
Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const ConstVectorView &invcm_grid, const Numeric stotmax, const calctype type)
void compute_transmission_matrix_and_derivative(Tensor3View T, Tensor4View dT_dx_upper_level, Tensor4View dT_dx_lower_level, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const ArrayOfPropagationMatrix &dupper_level_dx, const ArrayOfPropagationMatrix &dlower_level_dx, const Numeric &dr_dTu, const Numeric &dr_dTl, const Index it, const Index iz, const Index ia)
void compute_transmission_matrix(Tensor3View T, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
#define u
#define d
#define v
#define w
#define a
Array< PropagationMatrix > ArrayOfPropagationMatrix
#define c
#define b
Hund
Enum for Hund cases.
Definition: quantum.h:109
Quantum::Identifier QuantumIdentifier
Definition: quantum.h:471
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
#define N
Definition: rng.cc:164
Species::Tag SpeciesTag
Definition: species_tags.h:99
void test2()
void test4()
int test1()
Definition: test_matpack.cc:47
void test_transmat_to_cumulativetransmat()
void test_r_deriv_propagationmatrix()
void test_sinc_likes_0limit()
void test_zeeman()
int main(int n, char **argc)
void test_hitran2017(bool newtest=true)
void test_transmissionmatrix()
constexpr bool test_quantum_numbers(const QuantumNumbers qns, const Index i)
void test_matrix_buildup()
void test_mpm20()
void test_transmat_from_propmat()
void test_quantum()
void test3()
Definition: test_sparse.cc:43
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
Stuff related to the transmission matrix.
Index make_wigner_ready(int largest, int fastest, int size)
Ready Wigner.
Wigner symbol interactions.
Headers and class definition of Zeeman modeling.