ARTS 2.5.4 (git: 31ce4f0e)
radiation_field.cc
Go to the documentation of this file.
1/* Copyright (C) 2019
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 "radiation_field.h"
28#include "arts_conversions.h"
29#include "sorting.h"
30
31void error_in_integrate(const String& error_msg,
32 const Numeric& value_that_should_be_unity) {
33 ARTS_USER_ERROR_IF (std::abs(value_that_should_be_unity - 1.0) > 1e-4,
34 "Failure in normalization:\n", error_msg, "\n")
35}
36
37Numeric test_integrate_convolved(const Eigen::Ref<Eigen::VectorXcd> &F,
38 const Vector& f) {
39 Numeric val = 0.0;
40
41 const Index n = f.nelem();
42 for (Index i = 0; i < n - 1; i++)
43 val += 0.5 * (f[i + 1] - f[i]) * (F[i].real() + F[i + 1].real());
44
45 return val; // Should return 1.0 for normalized F
46}
47
49 const Array<Index>& sorted_index) {
50 Numeric val = 0.0;
51
52 const Index n = cosza.nelem();
53 for (Index i = 0; i < n - 1; i++)
54 val += 0.5 * (cosza[sorted_index[i]] - cosza[sorted_index[i + 1]]);
55
56 return val; // Should return 1.0 for normalized cosza
57}
58
60 const Eigen::VectorXcd& F,
61 const Vector& f) {
62 Numeric val = 0.0;
63
64 const Index n = f.nelem();
65 for (Index i = 0; i < n - 1; i++)
66 val += 0.5 * (f[i + 1] - f[i]) *
67 (I(i, 0) * F[i].real() + I(i + 1, 0) * F[i + 1].real());
68
69 return val;
70}
71
73 const Eigen::VectorXcd& F,
74 const Vector& f) {
75 Numeric val = 0.0;
76
77 const Index n = f.nelem();
78 for (Index i = 0; i < n - 1; i++)
79 val += 0.5 * (f[i + 1] - f[i]) *
80 (T(i, 0, 0) * F[i].real() + T(i + 1, 0, 0) * F[i + 1].real());
81
82 return 1.0 - val;
83}
84
86 const Vector& cosza,
87 const Array<Index>& sorted_index) {
88 Numeric val = 0.0;
89
90 const Index n = cosza.nelem();
91 for (Index i = 0; i < n - 1; i++)
92 val += 0.25 * (cosza[sorted_index[i]] - cosza[sorted_index[i + 1]]) *
93 (j[sorted_index[i]] + j[sorted_index[i + 1]]);
94
95 return val;
96}
97
99 if (gp.fd[1] == 1.0)
100 return gp.idx;
101 return gp.idx + 1;
102}
103
105 ArrayOfVector& cosza,
106 const ArrayOfPpath& ppath_field) {
107 Index nalt = 0;
108 for (auto& path : ppath_field)
109 for (auto& gp : path.gp_p) nalt = std::max(nalt, grid_index_from_gp(gp));
110 nalt += 1;
111
112 // Get the data, which is of unknown length
113 Array<ArrayOfNumeric> zeniths_array(nalt, ArrayOfNumeric(0));
114 for (auto& path : ppath_field) {
115 for (Index ip = 0; ip < path.np; ip++) {
116 const Index ind = grid_index_from_gp(path.gp_p[ip]);
117 zeniths_array[ind].push_back(path.los(ip, 0));
118 }
119 }
120
121 // Finalize sorting
122 sorted_index.resize(nalt);
123 cosza.resize(nalt);
124 for (Index i = 0; i < nalt; i++) {
125 Vector& data = cosza[i];
126 data.resize(zeniths_array[i].nelem());
127
128 for (Index j = 0; j < data.nelem(); j++) data[j] = zeniths_array[i][j];
129 get_sorted_indexes(sorted_index[i], data);
130
131 for (Index j = 0; j < data.nelem(); j++)
132 data[j] = Conversion::cosd(data[j]);
133 }
134}
Array< Numeric > ArrayOfNumeric
An array of Numeric.
Definition: array.h:141
Common ARTS conversions.
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define abs(x)
#define max(a, b)
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
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric e
Elementary charge convenience name [C].
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
void sorted_index_of_ppath_field(ArrayOfArrayOfIndex &sorted_index, ArrayOfVector &cosza, const ArrayOfPpath &ppath_field)
Get sorting of zenith angles in field of ppath.
Numeric test_integrate_convolved(const Eigen::Ref< Eigen::VectorXcd > &F, const Vector &f)
Integrate the line shape.
void error_in_integrate(const String &error_msg, const Numeric &value_that_should_be_unity)
Throws an error if integration values are bad.
Numeric integrate_zenith(const ConstVectorView &j, const Vector &cosza, const Array< Index > &sorted_index)
Convolve source function with 1D sphere and integrate.
Index grid_index_from_gp(const GridPos &gp)
Get a discrete position from grid pos.
Numeric test_integrate_zenith(const Vector &cosza, const Array< Index > &sorted_index)
Integrate cos(za) over the angles.
Numeric integrate_convolved(const RadiationVector &I, const Eigen::VectorXcd &F, const Vector &f)
Convolve intensity and line shape and integrate.
Radiation field calculations.
Contains sorting routines.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:57
Structure to store a grid position.
Definition: interpolation.h:73
std::array< Numeric, 2 > fd
Definition: interpolation.h:75
Index idx
Definition: interpolation.h:74
Radiation Vector for Stokes dimension 1-4.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.