ARTS 2.5.11 (git: 6827797f)
nlte.cc
Go to the documentation of this file.
1
9#include "nlte.h"
10#include "arts_constants.h"
11#include "interp.h"
12
14 const ConstVectorView& Aij,
15 const ConstVectorView& Bij,
16 const ConstVectorView& Bji,
17 const ConstVectorView& Cij,
18 const ConstVectorView& Cji,
19 const ConstVectorView& Jij,
20 const ArrayOfIndex& upper,
21 const ArrayOfIndex& lower) {
22 const Index nlines = Aij.nelem();
23
24 A = 0.0;
25 for (Index iline = 0; iline < nlines; iline++) {
26 const Index i = upper[iline];
27 const Index j = lower[iline];
28
29 A(j, j) -= Bji[iline] * Jij[iline] + Cji[iline];
30 A(i, i) -= Aij[iline] + Bij[iline] * Jij[iline] + Cij[iline];
31
32 A(j, i) += Aij[iline] + Bij[iline] * Jij[iline] + Cij[iline];
33 A(i, j) += Bji[iline] * Jij[iline] + Cji[iline];
34 }
35}
36
38 MatrixView A,
39 const ConstVectorView& x,
40 const ConstVectorView& Aij,
41 const ConstVectorView& Bij,
42 const ConstVectorView& Bji,
43 const ConstVectorView& Cij,
44 const ConstVectorView& Cji,
45 const ConstVectorView& Jij,
46 const ConstVectorView& Lambda,
47 const ArrayOfIndex& upper,
48 const ArrayOfIndex& lower,
49 const Numeric& total_number_count) {
50 const Index nlines = Aij.nelem();
51
52 A = 0.0;
53 for (Index iline = 0; iline < nlines; iline++) {
54 const Index i = upper[iline];
55 const Index j = lower[iline];
56
57 const Numeric Source =
58 total_number_count *
59 (x[i] * Aij[iline] / (x[j] * Bji[iline] - x[i] * Bij[iline]));
60
61 A(j, j) -= Bji[iline] * (Jij[iline] - Lambda[iline] * Source) + Cji[iline];
62 A(i, i) -= Aij[iline] * (1.0 - Lambda[iline]) +
63 Bij[iline] * (Jij[iline] - Lambda[iline] * Source) + Cij[iline];
64
65 A(j, i) += Aij[iline] * (1.0 - Lambda[iline]) +
66 Bij[iline] * (Jij[iline] - Lambda[iline] * Source) + Cij[iline];
67 A(i, j) += Bji[iline] * (Jij[iline] - Lambda[iline] * Source) + Cji[iline];
68 }
69}
70
72 VectorView x,
73 const Numeric& sem_ratio,
74 const Index row) {
75 A(row, joker) = 1.0;
76 x[row] = sem_ratio;
77}
78
79Vector createAij(const ArrayOfArrayOfAbsorptionLines& abs_lines) {
80 // Size of problem
81 const Index n = nelem(abs_lines);
82 Vector Aij(n);
83
84 Index i=0;
85 for (auto& lines: abs_lines) {
86 for (auto& band: lines) {
87 for (Index k=0; k<band.NumLines(); k++) {
88 Aij[i] = band.lines[k].A;
89 i++;
90 }
91 }
92 }
93 return Aij;
94}
95
96Vector createBij(const ArrayOfArrayOfAbsorptionLines& abs_lines) {
97 constexpr Numeric c0 = 2.0 * Constant::h / Math::pow2(Constant::c);
98
99 // Size of problem
100 const Index n = nelem(abs_lines);
101 Vector Bij(n);
102
103 // Base equation for single state: B21 = A21 c^2 / 2 h f^3 (nb. SI, don't use this without checking your need)
104 Index i=0;
105 for (auto& lines: abs_lines) {
106 for (auto& band: lines) {
107 for (Index k=0; k<band.NumLines(); k++) {
108 Bij[i] = band.lines[k].A / (c0 * Math::pow3(band.lines[k].F0));
109 i++;
110 }
111 }
112 }
113 return Bij;
114}
115
116Vector createBji(const Vector& Bij, const ArrayOfArrayOfAbsorptionLines& abs_lines) {
117 // Size of problem
118 const Index n = Bij.nelem();
119 Vector Bji(n);
120
121 // Base equation for single state: B12 = B21 g2 / g1
122 Index i=0;
123 for (auto& lines: abs_lines) {
124 for (auto& band: lines) {
125 for (Index k=0; k<band.NumLines();k++) {
126 Bji[i] = Bij[i] * band.lines[k].gupp / band.lines[k].glow;
127 i++;
128 }
129 }
130 }
131 return Bji;
132}
133
134Vector createCji(const Vector& Cij,
135 const ArrayOfArrayOfAbsorptionLines& abs_lines,
136 const Numeric& T) {
137 const Index n = nelem(abs_lines);
138 Vector Cji(n);
139 setCji(Cji, Cij, abs_lines, T);
140 return Cji;
141}
142
143void setCji(Vector& Cji,
144 const Vector& Cij,
145 const ArrayOfArrayOfAbsorptionLines& abs_lines,
146 const Numeric& T) {
147 static constexpr Numeric c0 = -Constant::planck_constant / Constant::boltzmann_constant;
148 const Numeric constant = c0 / T;
149
150 // Base equation for single state: C12 = C21 exp(-hf / kT) g2 / g1
151 Index i=0;
152 for (auto& lines: abs_lines) {
153 for (auto& band: lines) {
154 for (Index k=0; k<band.NumLines(); k++) {
155 Cji[i] = Cij[i] * exp(constant * band.lines[k].F0) * band.lines[k].gupp / band.lines[k].glow;
156 i++;
157 }
158 }
159 }
160}
161
163 Vector& Cij,
164 Vector& Cji,
165 const ArrayOfArrayOfAbsorptionLines& abs_lines,
166 const ArrayOfArrayOfSpeciesTag& abs_species,
167 const ArrayOfArrayOfGriddedField1& collision_coefficients,
168 const ArrayOfQuantumIdentifier& collision_line_identifiers,
169 const SpeciesIsotopologueRatios& isotopologue_ratios,
170 const ConstVectorView& vmr,
171 const Numeric& T,
172 const Numeric& P) {
173 // size of problem
174 const Index nspec = abs_species.nelem();
175 const Index ntrans = collision_line_identifiers.nelem();
176
177 // reset Cij for summing later
178 Cij = 0;
179
180 // For all species
181 for (Index i = 0; i < nspec; i++) {
182 // Compute the number density noting that free_electrons will behave differently
183 const Numeric numden =
184 vmr[i] * (abs_species[i].FreeElectrons() ? 1.0 : P / (Constant::k * T));
185
186 for (Index j = 0; j < ntrans; j++) {
187 Index iline=0;
188 for (auto& lines: abs_lines) {
189 for (auto& band: lines) {
190 const Numeric isot_ratio = isotopologue_ratios[band.Isotopologue()];
191 for (Index k=0; k<band.NumLines(); k++) {
192
193 const auto& transition = collision_line_identifiers[j];
194 const auto& gf1 = collision_coefficients[i][j];
195 const Quantum::Number::StateMatch lt(transition, band.lines[k].localquanta, band.quantumidentity);
196
197 if (lt == Quantum::Number::StateMatchType::Full) {
198 // Standard linear ARTS interpolation
199 const FixedLagrangeInterpolation<1> lag(0, T, gf1.get_numeric_grid(0));
200 const auto itw = interpweights(lag);
201
202 Cij[iline] += interp(gf1.data, itw, lag) * numden * isot_ratio;
203 iline++;
204 break;
205 }
206 iline++;
207 }
208 }
209 }
210 }
211 }
212
213 // Compute the reverse
214 setCji(Cji, Cij, abs_lines, T);
215}
216
218 ArrayOfIndex& upper,
219 ArrayOfIndex& lower,
220 const ArrayOfArrayOfAbsorptionLines& abs_lines,
221 const EnergyLevelMap& nlte_field) {
222 const Index nl = nelem(abs_lines), nq = nlte_field.levels.nelem();
223
224 upper = ArrayOfIndex(nl, -1);
225 lower = ArrayOfIndex(nl, -1);
226
227 Index i=0;
228 for (auto& lines: abs_lines) {
229 for (const AbsorptionLines& band: lines) {
230 for (Index k=0; k<band.NumLines(); k++) {
231 for (Index iq = 0; iq < nq; iq++) {
232 const Quantum::Number::StateMatch lt(nlte_field.levels[iq], band.lines[k].localquanta, band.quantumidentity);
233 if (lt == Quantum::Number::StateMatchType::Level and lt.low)
234 lower[i] = iq;
235 if (lt == Quantum::Number::StateMatchType::Level and lt.upp)
236 upper[i] = iq;
237 }
238 i++;
239 }
240 }
241 }
242
243 i = 0;
244 for (Index il = 0; il < nl; il++)
245 if (upper[il] < 0 or lower[il] < 0) i++;
246 ARTS_USER_ERROR_IF (i > 1,
247 "Must set upper and lower levels completely for all but one level");
248}
249
251 const ArrayOfIndex& lower) ARTS_NOEXCEPT {
252 for (const Index& l : lower) {
253 if (std::find(upper.cbegin(), upper.cend(), l) == upper.cend())
254 return l;
255 }
256 return upper.nelem() - 1;
257}
258
259void check_collision_line_identifiers(const ArrayOfQuantumIdentifier& collision_line_identifiers) {
260 auto p =
261 std::find_if(collision_line_identifiers.cbegin(),
262 collision_line_identifiers.cend(),
263 [isot = collision_line_identifiers.front().Isotopologue()](
264 auto& x) { return isot not_eq x.Isotopologue(); });
265 ARTS_USER_ERROR_IF (p not_eq collision_line_identifiers.cend(),
266 *p, "\n"
267 "does not match the requirements for a line identifier\n"
268 "Your list of species is:\n",
269 collision_line_identifiers, "\n"
270 "This contains more than one isotopologue or it contains some non-transition type identifiers.\n"
271 "It will therefore fail in current code. You can only input transitions, and a single isotopologue.\n")
272}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:119
Constants of physical expressions as constexpr.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
#define ARTS_NOEXCEPT
Definition: debug.h:81
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:135
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
constexpr Numeric boltzmann_constant
Boltzmann constant [J/K] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date:...
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric planck_constant
Planck constant [J s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 20...
constexpr Numeric c
Speed of light convenience name [m/s].
constexpr Numeric h
Planck constant convenience name [J s].
constexpr auto pow3(auto x) noexcept
power of three
constexpr auto pow2(auto x) noexcept
power of two
Vector createBij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bij object.
Definition: nlte.cc:96
Vector createBji(const Vector &Bij, const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bji object.
Definition: nlte.cc:116
Vector createAij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Aij object.
Definition: nlte.cc:79
Vector createCji(const Vector &Cij, const ArrayOfArrayOfAbsorptionLines &abs_lines, const Numeric &T)
Create a Cji object.
Definition: nlte.cc:134
void check_collision_line_identifiers(const ArrayOfQuantumIdentifier &collision_line_identifiers)
Checks that a WSV is OK or throws a run-time error.
Definition: nlte.cc:259
void set_constant_statistical_equilibrium_matrix(MatrixView A, VectorView x, const Numeric &sem_ratio, const Index row)
Set a row of the SEE matrix and level distribution vector to constant.
Definition: nlte.cc:71
void nlte_collision_factorsCalcFromCoeffs(Vector &Cij, Vector &Cji, const ArrayOfArrayOfAbsorptionLines &abs_lines, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfGriddedField1 &collision_coefficients, const ArrayOfQuantumIdentifier &collision_line_identifiers, const SpeciesIsotopologueRatios &isotopologue_ratios, const ConstVectorView &vmr, const Numeric &T, const Numeric &P)
Gets collisional factors from coefficients.
Definition: nlte.cc:162
Index find_first_unique_in_lower(const ArrayOfIndex &upper, const ArrayOfIndex &lower) ARTS_NOEXCEPT
Finds a unique lower state if one exists or returns index to last element.
Definition: nlte.cc:250
void nlte_positions_in_statistical_equilibrium_matrix(ArrayOfIndex &upper, ArrayOfIndex &lower, const ArrayOfArrayOfAbsorptionLines &abs_lines, const EnergyLevelMap &nlte_field)
Finds upper and lower states in SEE Matrix.
Definition: nlte.cc:217
void setCji(Vector &Cji, const Vector &Cij, const ArrayOfArrayOfAbsorptionLines &abs_lines, const Numeric &T)
Set the Cji object.
Definition: nlte.cc:143
void dampened_statistical_equilibrium_equation(MatrixView A, const ConstVectorView &x, const ConstVectorView &Aij, const ConstVectorView &Bij, const ConstVectorView &Bji, const ConstVectorView &Cij, const ConstVectorView &Cji, const ConstVectorView &Jij, const ConstVectorView &Lambda, const ArrayOfIndex &upper, const ArrayOfIndex &lower, const Numeric &total_number_count)
Sets up the solution matrix for linear dampened statistical equilibrium equation.
Definition: nlte.cc:37
void statistical_equilibrium_equation(MatrixView A, const ConstVectorView &Aij, const ConstVectorView &Bij, const ConstVectorView &Bji, const ConstVectorView &Cij, const ConstVectorView &Cji, const ConstVectorView &Jij, const ArrayOfIndex &upper, const ArrayOfIndex &lower)
Sets up the solution matrix for linear statistical equilibrium equation.
Definition: nlte.cc:13
Deep calculations for NLTE.
ArrayOfQuantumIdentifier levels
StateMatchType operates so that a check less than a level should be 'better', bar None.