ARTS 2.5.4 (git: 31ce4f0e)
linemixing_jmh.cc
Go to the documentation of this file.
1#include "arts_conversions.h"
2#include "linemixing.h"
3#include "matpackI.h"
4#include "wigner_functions.h"
5#include <cstdlib>
6#include <iostream>
7
8constexpr Numeric FECS(Numeric OMGA, Numeric EC) noexcept {return 1./(1.+EC*OMGA*OMGA)/(1.+EC*OMGA*OMGA);}
9
10Numeric relaxation_matrix_element(const Index li, const Index lf, const Index jji, const Index jjf, const Index jjip, const Index jjfp, const Vector& ECS, const Vector& QL) {
11 constexpr Index idl = 2;
12 Index Ldeb = std::max(std::abs(jji-jjip), std::abs(jjf-jjfp));
13 if ((Ldeb % idl) not_eq 0) Ldeb += 1;
14 const Index Lfin = std::min(jji+jjip, jjf+jjfp);
15
16 Numeric som=0.0;
17 for (Index L=Ldeb; L <= Lfin; L+=idl) {
18 Numeric RJ = wigner3j(jji, jjip, L, li, -li, 0);
19 const Numeric SS=RJ*QL[L]/ECS[L];
20 RJ = wigner3j(jjf, jjfp, L, lf, -lf, 0);
21 som += RJ*SS*wigner6j(jji, jjf, 1, jjfp, jjip, L);
22 }
23 const Numeric ROJ=Numeric(2*jjip+1)*std::sqrt(Numeric((2*jjf+1)*(2*jjfp+1)))*ECS[jji];
24
25 som *= ROJ;
26 return som;
27}
28
29int main() try {
30 constexpr Index Jmax = 128;
31 constexpr Index nRmx = 4 * Jmax;
32 constexpr Index lmax = 4 * Jmax;
33 constexpr std::array typerf{'P', 'Q', 'R'};
34 constexpr Numeric T0{296.0};
35
36 make_wigner_ready(300, 50000000, 6);
37
38 Index nraies = -1;
39
40 for (Index ll=0; ll<=8; ll++) {
41 for (Index Deltal : {0, 1}) {
42 const Index li = ll;
43 const Index lf = ll + Deltal;
44
45 Tensor3 W(nRmx, nRmx, 8, 0.0);
46 ArrayOfIndex ji(nRmx, 0), jf(nRmx, 0);
47 for (Index itemp=0; itemp<8; itemp++) {
48 const auto temp = Numeric(180 + 20 * itemp);
49
50 Tensor3 Wipert(nRmx, nRmx, 2, 0.0);
51
52 for (Index ipp=0; ipp<2; ipp++) {
53 Array<char> typeR(nRmx, '\0');
54
55 const Index iPert = 2 + ipp;
56 const std::array amasse{4., 40., 28., 32.};
57 const std::array aa{35.46e-3, 19.39e-3, 0.0180*std::pow(T0/temp, 0.85), 0.0168*std::pow(T0/temp, 0.5)};
58 const std::array alc{0.53, 3.0, 2.2, 2.4};
59 const std::array alp{1.062, 0.848, 0.81*std::pow(T0/temp, 0.0152), 0.82*std::pow(T0/temp, -0.0910)};
60 const std::array bbet{0., 0.02, 0.008, 0.007};
61 const std::array dx{0.0, aa[iPert], alc[iPert], alp[iPert], bbet[iPert]};
62
63 // QL
64 Vector QL(lmax + 1, 0);
65 Vector ECS(lmax + 1, 0);
66 const Numeric AM = 1./(1./amasse[iPert]+1./44.);
67 const Numeric ECT = 0.0006983*AM*dx[2]*dx[2]/temp;
68 const Numeric AT = dx[1];
69 QL[0] = 0.;
70 ECS[0] = 1;
71 constexpr Index idl = 2;
72 const Numeric brot = 0.39;
73 const Numeric betaa = 1.4388 * brot / temp;
74 for (Index L=1; L <= lmax; L++) {
75 const Numeric QMX = brot * Numeric((L + L + 1 - idl) * idl);
76 ECS[L] = FECS(QMX, ECT);
77 const auto AL2 = Numeric(L*L+L);
78 const Numeric SL0 = AT/std::pow(AL2, dx[3])*std::exp(-betaa*dx[4]*AL2);
79 QL[L] = Numeric(L+L+1) * SL0;
80 }
81
82 // ji,jf
83 nraies = -1;
84 for (Index ibid=3; ibid>=1; ibid--) {
85 const Index itypeR = (ibid % 3) - 1;
86 if (li == 0 and lf == 0 and itypeR == 0) continue;
87 for (Index jji=0; jji<=Jmax; jji++) {
88 Index jjf = jji + itypeR;
89 if (jji >= li and jjf >= lf) {
90 nraies++;
91 ji[nraies] = jji;
92 jf[nraies] = jjf;
93 typeR[nraies] = typerf[itypeR];
94 }
95 }
96 }
97
98 // calculation of W
99 for (Index iR=0; iR<=nraies; iR++) {
100 const Index jji = ji[iR];
101 const Index jjf = jf[iR];
102 for (Index iRp=0; iRp<=nraies; iRp++) {
103 const Index jjip = ji[iRp];
104 const Index jjfp = jf[iRp];
105 if (jjip > jji) continue;
106
107 Wipert(iRp, iR, ipp) = relaxation_matrix_element(li, lf, jji, jjf, jjip, jjfp, ECS, QL);
108 }
109 }
110 }
111
112 for (Index iR=0; iR<=nraies; iR++) {
113 for (Index iRp=0; iRp<=nraies; iRp++) {
114 W(iRp, iR, itemp) = 0.79*Wipert(iRp,iR,0)+0.21*Wipert(iRp,iR,1);
115 }
116 }
117 }
118
119 std::ostringstream os;
120 os << "Wtype" << li << lf << "p.dat";
121 std::ofstream fil(os.str());
122 for (Index iR=0; iR<=nraies; iR++) {
123 for (Index iRp=0; iRp<=nraies; iRp++) {
124 if(ji[iRp] > ji[iR]) continue;
125 if (max(W(iRp, iR, joker)) == min(W(iRp, iR, joker))) continue;
126 char str[200];
127 std::sprintf(
128 str,
129 " %0.12E %0.12E %0.12E %0.12E %0.12E %0.12E %0.12E %0.12E %" PRId64
130 " %" PRId64 " %" PRId64 " %" PRId64 "\n",
131 W(iRp, iR, 0),
132 W(iRp, iR, 1),
133 W(iRp, iR, 2),
134 W(iRp, iR, 3),
135 W(iRp, iR, 4),
136 W(iRp, iR, 5),
137 W(iRp, iR, 6),
138 W(iRp, iR, 7),
139 ji[iR],
140 jf[iR],
141 ji[iRp],
142 jf[iRp]);
143 fil << str;
144 }
145 }
146
147 std::cout << "Done with file: " << os.str() << '\n';
148 }
149 }
150
151 return EXIT_SUCCESS;
152} catch(...) {
153 return EXIT_FAILURE;
154}
Common ARTS conversions.
The Tensor3 class.
Definition: matpackIII.h:346
The Vector class.
Definition: matpackI.h:899
#define abs(x)
#define temp
#define min(a, b)
#define dx
#define ll
#define max(a, b)
Numeric relaxation_matrix_element(const Index li, const Index lf, const Index jji, const Index jjf, const Index jjip, const Index jjfp, const Vector &ECS, const Vector &QL)
constexpr Numeric FECS(Numeric OMGA, Numeric EC) noexcept
int main()
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
Numeric wigner6j(const Rational j1, const Rational j2, const Rational j3, const Rational l1, const Rational l2, const Rational l3)
Wigner 6J symbol.
Index make_wigner_ready(int largest, int fastest, int size)
Ready Wigner.
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
Wigner 3J symbol.
Wigner symbol interactions.