ARTS  2.4.0(git:4fb77825)
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 "interpolation_poly.h"
29 
31  ConstVectorView Aij,
32  ConstVectorView Bij,
33  ConstVectorView Bji,
34  ConstVectorView Cij,
35  ConstVectorView Cji,
36  ConstVectorView Jij,
37  const ArrayOfIndex& upper,
38  const ArrayOfIndex& lower) {
39  const Index nlines = Aij.nelem();
40 
41  A = 0.0;
42  for (Index iline = 0; iline < nlines; iline++) {
43  const Index i = upper[iline];
44  const Index j = lower[iline];
45 
46  A(j, j) -= Bji[iline] * Jij[iline] + Cji[iline];
47  A(i, i) -= Aij[iline] + Bij[iline] * Jij[iline] + Cij[iline];
48 
49  A(j, i) += Aij[iline] + Bij[iline] * Jij[iline] + Cij[iline];
50  A(i, j) += Bji[iline] * Jij[iline] + Cji[iline];
51  }
52 }
53 
55  MatrixView A,
57  ConstVectorView Aij,
58  ConstVectorView Bij,
59  ConstVectorView Bji,
60  ConstVectorView Cij,
61  ConstVectorView Cji,
62  ConstVectorView Jij,
63  ConstVectorView Lambda,
64  const ArrayOfIndex& upper,
65  const ArrayOfIndex& lower,
66  const Numeric& total_number_count) {
67  const Index nlines = Aij.nelem();
68 
69  A = 0.0;
70  for (Index iline = 0; iline < nlines; iline++) {
71  const Index i = upper[iline];
72  const Index j = lower[iline];
73 
74  const Numeric Source =
75  total_number_count *
76  (x[i] * Aij[iline] / (x[j] * Bji[iline] - x[i] * Bij[iline]));
77 
78  A(j, j) -= Bji[iline] * (Jij[iline] - Lambda[iline] * Source) + Cji[iline];
79  A(i, i) -= Aij[iline] * (1.0 - Lambda[iline]) +
80  Bij[iline] * (Jij[iline] - Lambda[iline] * Source) + Cij[iline];
81 
82  A(j, i) += Aij[iline] * (1.0 - Lambda[iline]) +
83  Bij[iline] * (Jij[iline] - Lambda[iline] * Source) + Cij[iline];
84  A(i, j) += Bji[iline] * (Jij[iline] - Lambda[iline] * Source) + Cji[iline];
85  }
86 }
87 
89  VectorView x,
90  const Numeric& sem_ratio,
91  const Index row) {
92  A(row, joker) = 1.0;
93  x[row] = sem_ratio;
94 }
95 
97  // Size of problem
98  const Index n = nelem(abs_lines);
99  Vector Aij(n);
100 
101  Index i=0;
102  for (auto& lines: abs_lines) {
103  for (auto& band: lines) {
104  for (Index k=0; k<band.NumLines(); k++) {
105  Aij[i] = band.A(k);
106  i++;
107  }
108  }
109  }
110  return Aij;
111 }
112 
114  extern const Numeric PLANCK_CONST, SPEED_OF_LIGHT;
115  const static Numeric c0 =
117 
118  // Size of problem
119  const Index n = nelem(abs_lines);
120  Vector Bij(n);
121 
122  // Base equation for single state: B21 = A21 c^2 / 2 h f^3 (nb. SI, don't use this without checking your need)
123  Index i=0;
124  for (auto& lines: abs_lines) {
125  for (auto& band: lines) {
126  for (Index k=0; k<band.NumLines(); k++) {
127  Bij[i] = band.A(k) / (c0 * band.F0(k) * band.F0(k) * band.F0(k));
128  i++;
129  }
130  }
131  }
132  return Bij;
133 }
134 
136  // Size of problem
137  const Index n = Bij.nelem();
138  Vector Bji(n);
139 
140  // Base equation for single state: B12 = B21 g2 / g1
141  Index i=0;
142  for (auto& lines: abs_lines) {
143  for (auto& band: lines) {
144  for (Index k=0; k<band.NumLines();k++) {
145  Bji[i] = Bij[i] * band.g_upp(k) / band.g_low(k);
146  i++;
147  }
148  }
149  }
150  return Bji;
151 }
152 
155  const Numeric& T) {
156  const Index n = nelem(abs_lines);
157  Vector Cji(n);
158  setCji(Cji, Cij, abs_lines, T);
159  return Cji;
160 }
161 
162 void setCji(Vector& Cji,
163  const Vector& Cij,
165  const Numeric& T) {
166  extern const Numeric PLANCK_CONST, BOLTZMAN_CONST;
167  const static Numeric c0 = -PLANCK_CONST / BOLTZMAN_CONST;
168  const Numeric constant = c0 / T;
169 
170  // Base equation for single state: C12 = C21 exp(-hf / kT) g2 / g1
171  Index i=0;
172  for (auto& lines: abs_lines) {
173  for (auto& band: lines) {
174  for (Index k=0; k<band.NumLines(); k++) {
175  Cji[i] = Cij[i] * exp(constant * band.F0(k)) * band.g_upp(k) / band.g_low(k);
176  i++;
177  }
178  }
179  }
180 }
181 
183  Vector& Cij,
184  Vector& Cji,
190  const ConstVectorView vmr,
191  const Numeric& T,
192  const Numeric& P) {
193  extern const Numeric BOLTZMAN_CONST;
194 
195  // size of problem
196  const Index nspec = abs_species.nelem();
197  const Index ntrans = collision_line_identifiers.nelem();
198 
199  // reset Cij for summing later
200  Cij = 0;
201 
202  // For all species
203  for (Index i = 0; i < nspec; i++) {
204  // Compute the number density noting that free_electrons will behave differently
205  const Numeric numden =
206  vmr[i] * (abs_species[i][0].SpeciesNameMain() == "free_electrons"
207  ? 1.0
208  : P / (BOLTZMAN_CONST * T));
209 
210  for (Index j = 0; j < ntrans; j++) {
211  Index iline=0;
212  for (auto& lines: abs_lines) {
213  for (auto& band: lines) {
214  const Numeric isot_ratio =
215  isotopologue_ratios.getParam(band.Species(), band.Isotopologue())[0]
216  .data[0];
217  for (Index k=0; k<band.NumLines(); k++) {
218 
219  const auto& transition = collision_line_identifiers[j];
220  const auto& gf1 = collision_coefficients[i][j];
221 
222  if (Absorption::id_in_line(band, transition, k)) {
223  // Standard linear ARTS interpolation
224  GridPosPoly gp;
225  gridpos_poly(gp, gf1.get_numeric_grid(0), T, 1, 0.5);
226  Vector itw(gp.idx.nelem());
227  interpweights(itw, gp);
228 
229  Cij[iline] += interp(itw, gf1.data, gp) * numden * isot_ratio;
230  iline++;
231  break;
232  }
233  iline++;
234  }
235  }
236  }
237  }
238  }
239 
240  // Compute the reverse
241  setCji(Cji, Cij, abs_lines, T);
242 }
243 
245  ArrayOfIndex& upper,
246  ArrayOfIndex& lower,
248  const EnergyLevelMap& nlte_field) {
249  const Index nl = nelem(abs_lines), nq = nlte_field.Levels().nelem();
250 
251  upper = ArrayOfIndex(nl, -1);
252  lower = ArrayOfIndex(nl, -1);
253 
254  Index i=0;
255  for (auto& lines: abs_lines) {
256  for (const AbsorptionLines& band: lines) {
257  for (Index k=0; k<band.NumLines(); k++) {
258  for (Index iq = 0; iq < nq; iq++) {
259  if (Absorption::id_in_line_lower(band, nlte_field.Levels()[iq], i))
260  lower[i] = iq;
261  if (Absorption::id_in_line_upper(band, nlte_field.Levels()[iq], i))
262  upper[i] = iq;
263  }
264  i++;
265  }
266  }
267  }
268 
269  i = 0;
270  for (Index il = 0; il < nl; il++)
271  if (upper[il] < 0 or lower[il] < 0) i++;
272  if (i > 1)
273  throw std::runtime_error(
274  "Must set upper and lower levels completely for all but one level");
275 }
276 
278  const ArrayOfIndex& lower) noexcept {
279  for (const Index& l : lower) {
280  if (std::find(upper.cbegin(), upper.cend(), l) == upper.cend())
281  return l;
282  }
283  return upper.nelem() - 1;
284 }
285 
287  auto p = std::find_if(collision_line_identifiers.cbegin(), collision_line_identifiers.cend(),
288  [spec=collision_line_identifiers.front().Species(), isot=collision_line_identifiers.front().Isotopologue()]
289  (auto& x) {
290  return
291  spec not_eq x.Species() or
292  isot not_eq x.Isotopologue() or
293  x.Type() not_eq QuantumIdentifier::TRANSITION;});
294  if (p not_eq collision_line_identifiers.cend()) {
295  std::ostringstream os;
296  os << *p << "\n"
297  << "does not match the requirements for a line identifier\n"
298  << "Your list of species is:\n"
299  << collision_line_identifiers << "\n"
300  << "This contains more than one isotopologue or it contains some non-transition type identifiers.\n"
301  << "It will therefore fail in current code. You can only input transitions, and a single isotopologue.\n";
302 
303  throw std::runtime_error(os.str());
304  }
305 }
ARTS::Var::collision_coefficients
ArrayOfArrayOfGriddedField1 collision_coefficients(Workspace &ws) noexcept
Definition: autoarts.h:2789
check_collision_line_identifiers
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:286
setCji
void setCji(Vector &Cji, const Vector &Cij, const ArrayOfArrayOfAbsorptionLines &abs_lines, const Numeric &T)
Set the Cji object.
Definition: nlte.cc:162
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
nlte_collision_factorsCalcFromCoeffs
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 SpeciesAuxData &isotopologue_ratios, const ConstVectorView vmr, const Numeric &T, const Numeric &P)
Gets collisional factors from coefficients.
Definition: nlte.cc:182
Absorption::id_in_line
bool id_in_line(const Lines &band, const QuantumIdentifier &id, size_t line_index)
Checks if the external quantum identifier match a line's ID.
Definition: absorptionlines.cc:2840
find_first_unique_in_lower
Index find_first_unique_in_lower(const ArrayOfIndex &upper, const ArrayOfIndex &lower) noexcept
Finds a unique lower state if one exists or returns index to last element.
Definition: nlte.cc:277
joker
const Joker joker
PLANCK_CONST
const Numeric PLANCK_CONST
Global constant, the Planck constant [Js].
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
Absorption::id_in_line_lower
bool id_in_line_lower(const Lines &band, const QuantumIdentifier &id, size_t line_index)
Checks if the external quantum identifier match a line's ID.
Definition: absorptionlines.cc:2904
ARTS::Var::collision_line_identifiers
ArrayOfQuantumIdentifier collision_line_identifiers(Workspace &ws) noexcept
Definition: autoarts.h:2797
Absorption::id_in_line_upper
bool id_in_line_upper(const Lines &band, const QuantumIdentifier &id, size_t line_index)
Checks if the external quantum identifier match a line's ID.
Definition: absorptionlines.cc:2877
SpeciesAuxData
Auxiliary data for isotopologues.
Definition: absorption.h:217
createBji
Vector createBji(const Vector &Bij, const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bji object.
Definition: nlte.cc:135
createBij
Vector createBij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bij object.
Definition: nlte.cc:113
ARTS::Var::nlte_field
EnergyLevelMap nlte_field(Workspace &ws) noexcept
Definition: autoarts.h:4604
Array
This can be used to make arrays out of anything.
Definition: array.h:108
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
createAij
Vector createAij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Aij object.
Definition: nlte.cc:96
GridPosPoly
Structure to store a grid position for higher order interpolation.
Definition: interpolation_poly.h:64
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
ArrayOfIndex
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:40
VectorView
The VectorView class.
Definition: matpackI.h:610
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
GridPosPoly::idx
ArrayOfIndex idx
Definition: interpolation_poly.h:67
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Absorption::Lines
Definition: absorptionlines.h:547
EnergyLevelMap
Definition: energylevelmap.h:60
interpolation_poly.h
Header file for interpolation_poly.cc.
createCji
Vector createCji(const Vector &Cij, const ArrayOfArrayOfAbsorptionLines &abs_lines, const Numeric &T)
Create a Cji object.
Definition: nlte.cc:153
nlines
#define nlines
Definition: legacy_continua.cc:22187
ARTS::Var::abs_lines
ArrayOfAbsorptionLines abs_lines(Workspace &ws) noexcept
Definition: autoarts.h:2013
nlte_positions_in_statistical_equilibrium_matrix
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:244
nlte.h
Deep calculations for NLTE.
set_constant_statistical_equilibrium_matrix
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:88
statistical_equilibrium_equation
void statistical_equilibrium_equation(MatrixView A, ConstVectorView Aij, ConstVectorView Bij, ConstVectorView Bji, ConstVectorView Cij, ConstVectorView Cji, ConstVectorView Jij, const ArrayOfIndex &upper, const ArrayOfIndex &lower)
Sets up the solution matrix for linear statistical equilibrium equation.
Definition: nlte.cc:30
ARTS::Var::isotopologue_ratios
SpeciesAuxData isotopologue_ratios(Workspace &ws) noexcept
Definition: autoarts.h:3669
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Definition: interpolation_poly.cc:120
Vector
The Vector class.
Definition: matpackI.h:860
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
BOLTZMAN_CONST
const Numeric BOLTZMAN_CONST
Global constant, the Boltzmann constant [J/K].
dampened_statistical_equilibrium_equation
void dampened_statistical_equilibrium_equation(MatrixView A, ConstVectorView x, ConstVectorView Aij, ConstVectorView Bij, ConstVectorView Bji, ConstVectorView Cij, ConstVectorView Cji, ConstVectorView Jij, 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:54