ARTS  2.4.0(git:4fb77825)
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 "lineshapemodel.h"
31 #include "linefunctions.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"
38 #include "predefined_absorption_models.h"
39 #include "wigner_functions.h"
40 
41 #include "linemixing_hitran.h"
42 #include <auto_md.h>
43 
44 void test_matrix_buildup() {
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 
142 void test_transmissionmatrix() {
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 
183 void test_r_deriv_propagationmatrix() {
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 
226  ArrayOfPropagationMatrix da(1, PropagationMatrix(1, nstokes));
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 
231  compute_transmission_matrix_and_derivative(
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 
275 void test_transmat_from_propmat() {
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 
284  PropagationMatrix test1(1, 1);
285  PropagationMatrix test2(1, 2);
286  PropagationMatrix test3(1, 3);
287  PropagationMatrix test4(1, 4);
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 
313  ArrayOfTransmissionMatrix empty(0);
314 
315  stepwise_transmission(ans1,
316  empty,
317  empty,
318  test1,
319  test1,
320  ArrayOfPropagationMatrix(0),
321  ArrayOfPropagationMatrix(0),
322  1,
323  0,
324  0,
325  -1);
326 
327  stepwise_transmission(ans2,
328  empty,
329  empty,
330  test2,
331  test2,
332  ArrayOfPropagationMatrix(0),
333  ArrayOfPropagationMatrix(0),
334  1,
335  0,
336  0,
337  -1);
338 
339  stepwise_transmission(ans3,
340  empty,
341  empty,
342  test3,
343  test3,
344  ArrayOfPropagationMatrix(0),
345  ArrayOfPropagationMatrix(0),
346  1,
347  0,
348  0,
349  -1);
350 
351  stepwise_transmission(ans4,
352  empty,
353  empty,
354  test4,
355  test4,
356  ArrayOfPropagationMatrix(0),
357  ArrayOfPropagationMatrix(0),
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 
369 void test_transmat_to_cumulativetransmat() {
370  int i = 1;
371  ArrayOfPropagationMatrix propmats(5, PropagationMatrix(1, 4));
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 
392  ArrayOfTransmissionMatrix layers(5, TransmissionMatrix(1, 4));
393  ArrayOfTransmissionMatrix empty(0);
394  for (i = 0; i < 4; i++) {
395  stepwise_transmission(layers[i + 1],
396  empty,
397  empty,
398  propmats[i],
399  propmats[i + 1],
400  ArrayOfPropagationMatrix(0),
401  ArrayOfPropagationMatrix(0),
402  1,
403  0,
404  0,
405  -1);
406  }
407 
408  std::cout << "Layers:\n" << layers << "\n\n";
409 
410  ArrayOfTransmissionMatrix cumulative_forward =
411  cumulative_transmission(layers, CumulativeTransmission::Forward);
412  ArrayOfTransmissionMatrix cumulative_reflect =
413  cumulative_transmission(layers, CumulativeTransmission::Reverse);
414 
415  std::cout << "Forward accumulation:\n" << cumulative_forward << "\n\n";
416  std::cout << "Reflect accumulation:\n" << cumulative_reflect << "\n\n";
417 }
418 
419 void test_sinc_likes_0limit() {
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 
447 void test_zeeman() {
448  define_species_data();
449  define_species_map();
450 
451  auto o266 = SpeciesTag("O2-66");
452  auto o268 = SpeciesTag("O2-68");
453 
454  Numeric g;
455  QuantumNumbers qn;
456  qn.Set(QuantumNumberType::Hund, Rational(Index(Hund::CaseB)));
457  qn.Set(QuantumNumberType::Lambda, 0_rat);
458  qn.Set(QuantumNumberType::v1, 0_rat);
459  qn.Set(QuantumNumberType::S, 1_rat);
460 
461  std::cout << "Table from Larsson, Lankhaar, Eriksson (2019)\n";
462  for (Index i = 1; i < 51; i++) {
463  qn.Set(QuantumNumberType::J, Rational(i));
464 
465  qn.Set(QuantumNumberType::N, qn[QuantumNumberType::J] - 1);
466  std::cout << i << "_" << i - 1;
467 
468  g = Zeeman::GetAdvancedModel(
469  QuantumIdentifier(o266.Species(), o266.Isotopologue(), qn, qn))
470  .gl();
471  std::cout << '\t' << g;
472 
473  g = Zeeman::GetAdvancedModel(
474  QuantumIdentifier(o268.Species(), o268.Isotopologue(), qn, qn))
475  .gl();
476  std::cout << '\t' << g;
477 
478  g = Zeeman::GetSimpleModel(
479  QuantumIdentifier(o266.Species(), o266.Isotopologue(), qn, qn))
480  .gl();
481  std::cout << '\t' << g;
482 
483  qn.Set(QuantumNumberType::N, qn[QuantumNumberType::J]);
484  std::cout << '\t' << i << "_" << i;
485 
486  g = Zeeman::GetAdvancedModel(
487  QuantumIdentifier(o266.Species(), o266.Isotopologue(), qn, qn))
488  .gl();
489  std::cout << '\t' << g;
490 
491  g = Zeeman::GetAdvancedModel(
492  QuantumIdentifier(o268.Species(), o268.Isotopologue(), qn, qn))
493  .gl();
494  std::cout << '\t' << g;
495 
496  g = Zeeman::GetSimpleModel(
497  QuantumIdentifier(o266.Species(), o266.Isotopologue(), qn, qn))
498  .gl();
499  std::cout << '\t' << g;
500 
501  qn.Set(QuantumNumberType::N, qn[QuantumNumberType::J] + 1);
502  std::cout << '\t' << i << "_" << i + 1;
503 
504  g = Zeeman::GetAdvancedModel(
505  QuantumIdentifier(o266.Species(), o266.Isotopologue(), qn, qn))
506  .gl();
507  std::cout << '\t' << g;
508 
509  g = Zeeman::GetAdvancedModel(
510  QuantumIdentifier(o268.Species(), o268.Isotopologue(), qn, qn))
511  .gl();
512  std::cout << '\t' << g;
513 
514  g = Zeeman::GetSimpleModel(
515  QuantumIdentifier(o266.Species(), o266.Isotopologue(), qn, qn))
516  .gl();
517  std::cout << '\t' << g;
518 
519  std::cout << '\n';
520  }
521 }
522 
523 constexpr bool test_quantum_numbers(const QuantumNumbers qns, const Index i)
524 {
525  return (i > 0) ? (qns[i].isUndefined() ? test_quantum_numbers(qns, i-1) : false) : true;
526 }
527 
528 void test_quantum()
529 {
530  static_assert(test_quantum_numbers(QuantumNumbers(), Index(QuantumNumberType::FINAL_ENTRY) - 1),
531  "Bad last entry in QuantumNumbers. Did you recently expand the list?");
532 }
533 
534 
535 void test_mpm20()
536 {
537  define_species_data();
538  define_species_map();
539 
540  // Test constants
541  constexpr Index nf = 501;
542  constexpr Numeric fstart = 25e9;
543  constexpr Numeric fend = 165e9;
544  constexpr Numeric t = 200;
545  constexpr Numeric p = 1e4;
546  Vector f(nf);
547  nlinspace(f, fstart, fend, nf);
548  Matrix xsec(nf, 1, 0);
549  ArrayOfMatrix dxsec(2, Matrix(nf, 1, 0));
550 
551  ArrayOfRetrievalQuantity jacs(2);
552  jacs[0].PropType(JacPropMatType::Temperature);
553  jacs[1].PropType(JacPropMatType::Frequency);
554  Absorption::PredefinedModel::makarov2020_o2_lines_mpm(xsec, dxsec, f, {p}, {t}, {0.5}, jacs, {0, 1});
555 
556  constexpr auto df = 1000;
557  constexpr auto dt = 0.1;
558  Matrix pxsec(nf, 1, 0);
559  Matrix pxsec_dt(nf, 1, 0);
560  Matrix pxsec_df(nf, 1, 0);
561  Vector f_pert = f;
562  f_pert += df;
563  PWR93O2AbsModel(pxsec, 0, 1, 1, 1, "user", "PWR98", f, {p}, {t}, {0.5}, {1}, Verbosity());
564  PWR93O2AbsModel(pxsec_dt, 0, 1, 1, 1, "user", "PWR98", f, {p}, {t+dt}, {0.5}, {1}, Verbosity());
565  PWR93O2AbsModel(pxsec_df, 0, 1, 1, 1, "user", "PWR98", f_pert, {p}, {t}, {0.5}, {1}, Verbosity());
566 
567  std::cout<< "xr = np.array([";
568  for (Index i=0; i<nf; i++)
569  std::cout<<'['<<xsec(i, 0)<<','<<' '<<pxsec(i, 0) << ']' << ',' << ' ';
570  std::cout << ']' << ')' << '\n';
571 
572  std::cout<< "dxr_dt = np.array([";
573  for (Index i=0; i<nf; i++)
574  std::cout<<'['<<dxsec[0](i, 0)<<','<<' '<<(pxsec_dt(i, 0)-pxsec(i, 0))/dt << ']' << ',' << ' ';
575  std::cout << ']' << ')' << '\n';
576 
577  std::cout<< "dxr_df = np.array([";
578  for (Index i=0; i<nf; i++)
579  std::cout<<'['<<dxsec[1](i, 0)<<','<<' '<<(pxsec_df(i, 0)-pxsec(i, 0))/df << ']' << ',' << ' ';
580  std::cout << ']' << ')' << '\n';
581 }
582 
583 
584 void test_ecs20()
585 {
586  constexpr Index nf = 501;
587  constexpr Numeric fstart = 45e9;
588  constexpr Numeric fend = 85e9;
589  constexpr Numeric t = 296;
590  constexpr Numeric p = 1.00658388e5;
591  Vector f(nf);
592  ComplexVector I(nf);
593  Matrix xsec(nf, 1, 0);
594  ArrayOfMatrix dxsec(0, Matrix(nf, 1, 0));
595  nlinspace(f, fstart, fend, nf);
596 // std::cout << "import numpy as np\n";
597 
598  define_species_data();
599  define_species_map();
600 
601  make_wigner_ready(200, 200, 6);
602 
603  wig_temp_init(200);
604  Absorption::PredefinedModel::makarov2020_o2_lines_ecs(I, f, p, t, 0);
605  ArrayOfRetrievalQuantity jacs(0);
606  Absorption::PredefinedModel::makarov2020_o2_lines_mpm(xsec, dxsec, f, {p}, {t}, {0.0}, jacs, {});
607  wig_temp_free();
608 
609  std::cout<<"I = np.array([";
610  for (Index i=0; i<f.nelem(); i++)
611  std::cout<<I[i].real()<<",\n";
612  std::cout<<"])\n";
613 
614  std::cout<<"I2 = np.array([";
615  for (Index i=0; i<f.nelem(); i++)
616  std::cout<<xsec(i,0)<<",\n";
617  std::cout<<"])\n";
618 }
619 
620 void test_hitran2017(bool newtest = true)
621 {
622  const Numeric p = 1.0;
623  const Numeric t = 296;
624  const Numeric xco2 = 1.5e-2;
625  const Numeric xh2o = 0;
626  const Numeric sigmin = 600;
627  const Numeric sigmax = 900;
628  const Numeric dsig = 0.005;
629  const Numeric stotmax = 0.1e-21;
630 
631  const Index nsig = Index(((sigmax - sigmin) / dsig) + 0.5) + 1;
632  Vector invcm_grid(nsig);
633  Vector f_grid(nsig);
634  Numeric sigc = sigmin-dsig;
635  for (Index isig=0; isig<nsig; isig++) {
636  sigc += dsig;
637  invcm_grid[isig] = sigc;
638  f_grid[isig] = Conversion::kaycm2freq(sigc);
639  }
640 
641  constexpr Index n=6;
642  auto types =
643  std::array<std::pair<lm_hitran_2017::ModeOfLineMixing, lm_hitran_2017::calctype>, 6>{
644  std::pair<lm_hitran_2017::ModeOfLineMixing, lm_hitran_2017::calctype>{lm_hitran_2017::ModeOfLineMixing::FullW, lm_hitran_2017::calctype::FullW},
645  {lm_hitran_2017::ModeOfLineMixing::VP_W, lm_hitran_2017::calctype::FullW},
646  {lm_hitran_2017::ModeOfLineMixing::VP, lm_hitran_2017::calctype::NoneVP},
647  {lm_hitran_2017::ModeOfLineMixing::VP_Y, lm_hitran_2017::calctype::NoneRosenkranz},
648  {lm_hitran_2017::ModeOfLineMixing::SDVP, lm_hitran_2017::calctype::SDVP},
649  {lm_hitran_2017::ModeOfLineMixing::SDVP_Y, lm_hitran_2017::calctype::SDRosenkranz}};
650  define_species_data();
651  define_species_map();
652  ArrayOfVector absorption(n);
653  make_wigner_ready(int(250), int(20000000), 6);
654 
655  ArrayOfAbsorptionLines bands;
656  HitranRelaxationMatrixData hitran;
657  for (Index i=0;i<n; i++) {
658  auto type=types[i];
659 
660  lm_hitran_2017::read(hitran, bands, "data_new", -1, Conversion::kaycm2freq(sigmin), Conversion::kaycm2freq(sigmax), Conversion::hitran2arts_linestrength(stotmax), type.first);
661  Vector vmrs = {1-xco2/100-xh2o/100, xh2o/100, xco2/100};
662  SpeciesAuxData partition_functions;
663  partition_functionsInitFromBuiltin(partition_functions, Verbosity());
664 
665  if (not newtest)
666  absorption[i] = lm_hitran_2017::compute(p, t, xco2, xh2o, invcm_grid, stotmax, type.second);
667  else
668  absorption[i] = lm_hitran_2017::compute(hitran, bands, Conversion::atm2pa(p), t, vmrs, f_grid, partition_functions);
669  }
670 
671  for (Index i=0; i<nsig; i++) {
672  for (Index j=0; j<n; j++) {
673  std::cout<<absorption[j][i]<<' ';
674  }
675  std::cout<<'\n';
676  }
677 }
678 
679 
680 int main(int n, char **argc) {
681  /*test_speed_of_pressurebroadening();
682  test_transmissionmatrix();
683  test_r_deriv_propagationmatrix();
684  test_transmat_from_propmat();
685  test_transmat_to_cumulativetransmat();
686  test_sinc_likes_0limit();*/
687 
688  if (n == 2 and String(argc[1]) == "new") {
689  std::cout<<"new test\n";
690  test_hitran2017(true);
691  }
692  else {
693  std::cout<<"old test\n";
694  test_hitran2017(false);
695  }
696  return 0;
697 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
legacy_continua.h
Zeeman::Model::gl
Numeric & gl() noexcept
Returns the lower state g.
Definition: zeemandata.h:387
test_mpm20
void test_mpm20()
Definition: test_propagationmatrix.cc:535
define_species_map
void define_species_map()
Definition: species_data.cc:1069
QuantumNumberType::N
@ N
ARTS::Group::PropagationMatrix
PropagationMatrix PropagationMatrix
Definition: autoarts.h:93
QuantumNumbers
Container class for Quantum Numbers.
Definition: quantum.h:222
test_quantum_numbers
constexpr bool test_quantum_numbers(const QuantumNumbers qns, const Index i)
Definition: test_propagationmatrix.cc:523
auto_md.h
absorption.h
Declarations required for the calculation of absorption coefficients.
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
Absorption::PredefinedModel::makarov2020_o2_lines_mpm
void makarov2020_o2_lines_mpm(Matrix &xsec, ArrayOfMatrix &dxsec, const Vector &f, const Vector &p, const Vector &t, const Vector &water_vmr, const ArrayOfRetrievalQuantity &jacs, const ArrayOfIndex &jacs_pos)
Adds Makarov MPM2020 O2 absorption lines to the absorption matrix.
Definition: predefined_absorption_models.cc:282
Absorption::PredefinedModel::makarov2020_o2_lines_ecs
void makarov2020_o2_lines_ecs(ComplexVector &I, const Vector &f, Numeric P, Numeric T, Numeric water_vmr)
Definition: predefined_absorption_models.cc:541
Conversion::kaycm2freq
constexpr Numeric kaycm2freq(T x)
Definition: constants.h:383
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
test_transmat_to_cumulativetransmat
void test_transmat_to_cumulativetransmat()
Definition: test_propagationmatrix.cc:369
test_hitran2017
void test_hitran2017(bool newtest=true)
Definition: test_propagationmatrix.cc:620
wigner_functions.h
Wigner symbol interactions.
joker
const Joker joker
linefunctions.h
Stuff related to lineshape functions.
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
ComplexVector
The ComplexVector class.
Definition: complex.h:573
abs
#define abs(x)
Definition: legacy_continua.cc:20626
predefined_absorption_models.h
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
PropagationMatrix
Definition: propagationmatrix.h:87
zeemandata.h
Headers and class definition of Zeeman modeling.
test_matrix_buildup
void test_matrix_buildup()
Definition: test_propagationmatrix.cc:44
test_sinc_likes_0limit
void test_sinc_likes_0limit()
Definition: test_propagationmatrix.cc:419
lm_hitran_2017::read
void read(HitranRelaxationMatrixData &hitran, ArrayOfAbsorptionLines &bands, 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.
Definition: linemixing_hitran.cc:1940
test_transmat_from_propmat
void test_transmat_from_propmat()
Definition: test_propagationmatrix.cc:275
QuantumNumberType::v1
@ v1
CumulativeTransmission::Forward
@ Forward
b2
#define b2
Definition: complex.h:59
PropagationMatrix::K13
VectorView K13(const Index iz=0, const Index ia=0)
Vector view to K(0, 2) elements.
Definition: propagationmatrix.h:819
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
compute_transmission_matrix
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.
Definition: propagationmatrix.cc:33
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
Zeeman::GetSimpleModel
Model GetSimpleModel(const QuantumIdentifier &qid) noexcept
Returns a simple Zeeman model.
Definition: zeemandata.cc:33
Zeeman::GetAdvancedModel
Model GetAdvancedModel(const QuantumIdentifier &qid) noexcept
Returns an advanced Zeeman model.
Definition: zeemandata.cc:105
SpeciesAuxData
Auxiliary data for isotopologues.
Definition: absorption.h:217
make_wigner_ready
Index make_wigner_ready(int largest, [[maybe_unused]] int fastest, int size)
Definition: wigner_functions.cc:276
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
nlogspace
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Definition: math_funcs.cc:267
Array
This can be used to make arrays out of anything.
Definition: array.h:108
SpeciesTag
A tag group can consist of the sum of several of these.
Definition: abs_species_tags.h:44
PropagationMatrix::K24
VectorView K24(const Index iz=0, const Index ia=0)
Vector view to K(1, 3) elements.
Definition: propagationmatrix.h:849
PropagationMatrix::K23
VectorView K23(const Index iz=0, const Index ia=0)
Vector view to K(1, 2) elements.
Definition: propagationmatrix.h:839
test_transmissionmatrix
void test_transmissionmatrix()
Definition: test_propagationmatrix.cc:142
MapToEigen
ComplexConstMatrixViewMap MapToEigen(const ConstComplexMatrixView &A)
Definition: complex.cc:1655
test3
void test3()
Definition: test_sparse.cc:43
HitranRelaxationMatrixData
Definition: linemixing_hitran.h:38
stepwise_transmission
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.
Definition: transmissionmatrix.cc:1353
test2
void test2()
Definition: test_matpack.cc:129
Hund::CaseA
@ CaseA
main
int main(int n, char **argc)
Definition: test_propagationmatrix.cc:680
Conversion::hitran2arts_linestrength
constexpr T hitran2arts_linestrength(T x)
Definition: constants.h:462
PropagationMatrix::K12
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
Definition: propagationmatrix.h:809
my_basic_string< char >
test1
int test1()
Definition: test_matpack.cc:47
lm_hitran_2017::ModeOfLineMixing::VP
@ VP
PropagationMatrix::K34
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
Definition: propagationmatrix.h:859
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
ArrayOfPropagationMatrix
Array< PropagationMatrix > ArrayOfPropagationMatrix
Definition: propagationmatrix.h:978
Hund
Hund
Enum for Hund cases.
Definition: quantum.h:219
ARTS::Group::Rational
Rational Rational
Definition: autoarts.h:96
PWR93O2AbsModel
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.
Definition: legacy_continua.cc:13762
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
RadiationVector::Vec3
const Eigen::Vector3d & Vec3(size_t i) const
Return Vector at position.
Definition: transmissionmatrix.h:503
PropagationMatrix::Data
Tensor4 & Data()
Get full view to data.
Definition: propagationmatrix.h:937
global_data.h
test_ecs20
void test_ecs20()
Definition: test_propagationmatrix.cc:584
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
partition_functionsInitFromBuiltin
void partition_functionsInitFromBuiltin(SpeciesAuxData &partition_functions, const Verbosity &)
WORKSPACE METHOD: partition_functionsInitFromBuiltin.
Definition: m_abs.cc:1631
test_r_deriv_propagationmatrix
void test_r_deriv_propagationmatrix()
Definition: test_propagationmatrix.cc:183
PropagationMatrix::K14
VectorView K14(const Index iz=0, const Index ia=0)
Vector view to K(0, 3) elements.
Definition: propagationmatrix.h:829
Zeeman::end
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:108
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition: transmissionmatrix.cc:1527
test_quantum
void test_quantum()
Definition: test_propagationmatrix.cc:528
RadiationVector
Radiation Vector for Stokes dimension 1-4.
Definition: transmissionmatrix.h:395
MatrixViewMap
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:110
Zeeman::start
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:77
linescaling.h
Constains various line scaling functions.
x1
#define x1
Definition: lineshapemodel.h:414
String
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:280
PropagationMatrix::Kjj
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Definition: propagationmatrix.h:799
ARTS::Var::partition_functions
SpeciesAuxData partition_functions(Workspace &ws) noexcept
Definition: autoarts.h:4860
compute_transmission_matrix_and_derivative
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)
Definition: propagationmatrix.cc:478
Conversion::atm2pa
constexpr Numeric atm2pa(T x)
Definition: constants.h:427
ARTS::Group::TransmissionMatrix
TransmissionMatrix TransmissionMatrix
Definition: autoarts.h:112
lm_hitran_2017::compute
Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const ConstVectorView invcm_grid, const Numeric stotmax, const calctype type)
Definition: linemixing_hitran.cc:1882
transmissionmatrix.h
Stuff related to the transmission matrix.
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
ARTS::Group::QuantumIdentifier
QuantumIdentifier QuantumIdentifier
Definition: autoarts.h:94
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
define_species_data
void define_species_data()
Definition: species_data.cc:158
test4
void test4()
Definition: test_matpack.cc:141
Vector
The Vector class.
Definition: matpackI.h:860
TransmissionMatrix
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
Definition: transmissionmatrix.h:38
QuantumNumbers::Set
void Set(Index qn, Rational r)
Set quantum number at position.
Definition: quantum.h:310
x2
#define x2
Definition: lineshapemodel.h:415
ARTS::Group::Verbosity
Verbosity Verbosity
Definition: autoarts.h:114
test_zeeman
void test_zeeman()
Definition: test_propagationmatrix.cc:447
lineshapemodel.h
Contains the line shape namespace.
zeeman.h
arts.h
The global header file for ARTS.
LineShape::vmrs
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const QuantumIdentifier &self, const ArrayOfSpeciesTag &lineshape_species, bool self_in_list, bool bath_in_list, Type type)
Returns a VMR vector for this model's main calculations.
Definition: lineshapemodel.cc:474
TransmissionMatrix::Mat4
const Eigen::Matrix4d & Mat4(size_t i) const
Get Matrix at position.
Definition: transmissionmatrix.h:129
linemixing_hitran.h
Namespace and functions to deal with HITRAN linemixing.
QuantumNumberType::J
@ J