ARTS 2.5.9 (git: 825fa5f2)
nlte.cc
Go to the documentation of this file.
1/* Copyright (C) 2018
2 Richard Larsson
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
27#include "nlte.h"
28#include "arts_constants.h"
30
32 const ConstVectorView& Aij,
33 const ConstVectorView& Bij,
34 const ConstVectorView& Bji,
35 const ConstVectorView& Cij,
36 const ConstVectorView& Cji,
37 const ConstVectorView& Jij,
38 const ArrayOfIndex& upper,
39 const ArrayOfIndex& lower) {
40 const Index nlines = Aij.nelem();
41
42 A = 0.0;
43 for (Index iline = 0; iline < nlines; iline++) {
44 const Index i = upper[iline];
45 const Index j = lower[iline];
46
47 A(j, j) -= Bji[iline] * Jij[iline] + Cji[iline];
48 A(i, i) -= Aij[iline] + Bij[iline] * Jij[iline] + Cij[iline];
49
50 A(j, i) += Aij[iline] + Bij[iline] * Jij[iline] + Cij[iline];
51 A(i, j) += Bji[iline] * Jij[iline] + Cji[iline];
52 }
53}
54
56 MatrixView A,
57 const ConstVectorView& x,
58 const ConstVectorView& Aij,
59 const ConstVectorView& Bij,
60 const ConstVectorView& Bji,
61 const ConstVectorView& Cij,
62 const ConstVectorView& Cji,
63 const ConstVectorView& Jij,
64 const ConstVectorView& Lambda,
65 const ArrayOfIndex& upper,
66 const ArrayOfIndex& lower,
67 const Numeric& total_number_count) {
68 const Index nlines = Aij.nelem();
69
70 A = 0.0;
71 for (Index iline = 0; iline < nlines; iline++) {
72 const Index i = upper[iline];
73 const Index j = lower[iline];
74
75 const Numeric Source =
76 total_number_count *
77 (x[i] * Aij[iline] / (x[j] * Bji[iline] - x[i] * Bij[iline]));
78
79 A(j, j) -= Bji[iline] * (Jij[iline] - Lambda[iline] * Source) + Cji[iline];
80 A(i, i) -= Aij[iline] * (1.0 - Lambda[iline]) +
81 Bij[iline] * (Jij[iline] - Lambda[iline] * Source) + Cij[iline];
82
83 A(j, i) += Aij[iline] * (1.0 - Lambda[iline]) +
84 Bij[iline] * (Jij[iline] - Lambda[iline] * Source) + Cij[iline];
85 A(i, j) += Bji[iline] * (Jij[iline] - Lambda[iline] * Source) + Cji[iline];
86 }
87}
88
90 VectorView x,
91 const Numeric& sem_ratio,
92 const Index row) {
93 A(row, joker) = 1.0;
94 x[row] = sem_ratio;
95}
96
98 // Size of problem
99 const Index n = nelem(abs_lines);
100 Vector Aij(n);
101
102 Index i=0;
103 for (auto& lines: abs_lines) {
104 for (auto& band: lines) {
105 for (Index k=0; k<band.NumLines(); k++) {
106 Aij[i] = band.lines[k].A;
107 i++;
108 }
109 }
110 }
111 return Aij;
112}
113
115 constexpr Numeric c0 = 2.0 * Constant::h / Math::pow2(Constant::c);
116
117 // Size of problem
118 const Index n = nelem(abs_lines);
119 Vector Bij(n);
120
121 // Base equation for single state: B21 = A21 c^2 / 2 h f^3 (nb. SI, don't use this without checking your need)
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 Bij[i] = band.lines[k].A / (c0 * Math::pow3(band.lines[k].F0));
127 i++;
128 }
129 }
130 }
131 return Bij;
132}
133
135 // Size of problem
136 const Index n = Bij.nelem();
137 Vector Bji(n);
138
139 // Base equation for single state: B12 = B21 g2 / g1
140 Index i=0;
141 for (auto& lines: abs_lines) {
142 for (auto& band: lines) {
143 for (Index k=0; k<band.NumLines();k++) {
144 Bji[i] = Bij[i] * band.lines[k].gupp / band.lines[k].glow;
145 i++;
146 }
147 }
148 }
149 return Bji;
150}
151
153 const ArrayOfArrayOfAbsorptionLines& abs_lines,
154 const Numeric& T) {
155 const Index n = nelem(abs_lines);
156 Vector Cji(n);
157 setCji(Cji, Cij, abs_lines, T);
158 return Cji;
159}
160
161void setCji(Vector& Cji,
162 const Vector& Cij,
163 const ArrayOfArrayOfAbsorptionLines& abs_lines,
164 const Numeric& T) {
166 const Numeric constant = c0 / T;
167
168 // Base equation for single state: C12 = C21 exp(-hf / kT) g2 / g1
169 Index i=0;
170 for (auto& lines: abs_lines) {
171 for (auto& band: lines) {
172 for (Index k=0; k<band.NumLines(); k++) {
173 Cji[i] = Cij[i] * exp(constant * band.lines[k].F0) * band.lines[k].gupp / band.lines[k].glow;
174 i++;
175 }
176 }
177 }
178}
179
181 Vector& Cij,
182 Vector& Cji,
183 const ArrayOfArrayOfAbsorptionLines& abs_lines,
184 const ArrayOfArrayOfSpeciesTag& abs_species,
185 const ArrayOfArrayOfGriddedField1& collision_coefficients,
186 const ArrayOfQuantumIdentifier& collision_line_identifiers,
187 const SpeciesIsotopologueRatios& isotopologue_ratios,
188 const ConstVectorView& vmr,
189 const Numeric& T,
190 const Numeric& P) {
191 // size of problem
192 const Index nspec = abs_species.nelem();
193 const Index ntrans = collision_line_identifiers.nelem();
194
195 // reset Cij for summing later
196 Cij = 0;
197
198 // For all species
199 for (Index i = 0; i < nspec; i++) {
200 // Compute the number density noting that free_electrons will behave differently
201 const Numeric numden =
202 vmr[i] * (abs_species[i].FreeElectrons() ? 1.0 : P / (Constant::k * T));
203
204 for (Index j = 0; j < ntrans; j++) {
205 Index iline=0;
206 for (auto& lines: abs_lines) {
207 for (auto& band: lines) {
208 const Numeric isot_ratio = isotopologue_ratios[band.Isotopologue()];
209 for (Index k=0; k<band.NumLines(); k++) {
210
211 const auto& transition = collision_line_identifiers[j];
212 const auto& gf1 = collision_coefficients[i][j];
213 const Quantum::Number::StateMatch lt(transition, band.lines[k].localquanta, band.quantumidentity);
214
215 if (lt == Quantum::Number::StateMatchType::Full) {
216 // Standard linear ARTS interpolation
217 const FixedLagrangeInterpolation<1> lag(0, T, gf1.get_numeric_grid(0), false);
218 const auto itw = interpweights(lag);
219
220 Cij[iline] += interp(gf1.data, itw, lag) * numden * isot_ratio;
221 iline++;
222 break;
223 }
224 iline++;
225 }
226 }
227 }
228 }
229 }
230
231 // Compute the reverse
232 setCji(Cji, Cij, abs_lines, T);
233}
234
236 ArrayOfIndex& upper,
237 ArrayOfIndex& lower,
238 const ArrayOfArrayOfAbsorptionLines& abs_lines,
239 const EnergyLevelMap& nlte_field) {
240 const Index nl = nelem(abs_lines), nq = nlte_field.levels.nelem();
241
242 upper = ArrayOfIndex(nl, -1);
243 lower = ArrayOfIndex(nl, -1);
244
245 Index i=0;
246 for (auto& lines: abs_lines) {
247 for (const AbsorptionLines& band: lines) {
248 for (Index k=0; k<band.NumLines(); k++) {
249 for (Index iq = 0; iq < nq; iq++) {
250 const Quantum::Number::StateMatch lt(nlte_field.levels[iq], band.lines[k].localquanta, band.quantumidentity);
251 if (lt == Quantum::Number::StateMatchType::Level and lt.low)
252 lower[i] = iq;
253 if (lt == Quantum::Number::StateMatchType::Level and lt.upp)
254 upper[i] = iq;
255 }
256 i++;
257 }
258 }
259 }
260
261 i = 0;
262 for (Index il = 0; il < nl; il++)
263 if (upper[il] < 0 or lower[il] < 0) i++;
264 ARTS_USER_ERROR_IF (i > 1,
265 "Must set upper and lower levels completely for all but one level");
266}
267
269 const ArrayOfIndex& lower) ARTS_NOEXCEPT {
270 for (const Index& l : lower) {
271 if (std::find(upper.cbegin(), upper.cend(), l) == upper.cend())
272 return l;
273 }
274 return upper.nelem() - 1;
275}
276
277void check_collision_line_identifiers(const ArrayOfQuantumIdentifier& collision_line_identifiers) {
278 auto p =
279 std::find_if(collision_line_identifiers.cbegin(),
280 collision_line_identifiers.cend(),
281 [isot = collision_line_identifiers.front().Isotopologue()](
282 auto& x) { return isot not_eq x.Isotopologue(); });
283 ARTS_USER_ERROR_IF (p not_eq collision_line_identifiers.cend(),
284 *p, "\n"
285 "does not match the requirements for a line identifier\n"
286 "Your list of species is:\n",
287 collision_line_identifiers, "\n"
288 "This contains more than one isotopologue or it contains some non-transition type identifiers.\n"
289 "It will therefore fail in current code. You can only input transitions, and a single isotopologue.\n")
290}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:136
Constants of physical expressions as constexpr.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
#define ARTS_NOEXCEPT
Definition: debug.h:79
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
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
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:114
Vector createBji(const Vector &Bij, const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bji object.
Definition: nlte.cc:134
Vector createAij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Aij object.
Definition: nlte.cc:97
Vector createCji(const Vector &Cij, const ArrayOfArrayOfAbsorptionLines &abs_lines, const Numeric &T)
Create a Cji object.
Definition: nlte.cc:152
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:277
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:89
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:180
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:268
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:235
void setCji(Vector &Cji, const Vector &Cij, const ArrayOfArrayOfAbsorptionLines &abs_lines, const Numeric &T)
Set the Cji object.
Definition: nlte.cc:161
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:55
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:31
Deep calculations for NLTE.
ArrayOfQuantumIdentifier levels
StateMatchType operates so that a check less than a level should be 'better', bar None.