ARTS  2.4.0(git:4fb77825)
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 "sorting.h"
29 
30 void error_in_integrate(const String& error_msg,
31  const Numeric& value_that_should_be_unity) {
32  if (std::abs(value_that_should_be_unity - 1.0) > 1e-4) {
33  std::ostringstream os;
34  os << "Failure in normalization:\n" << error_msg << "\n";
35  throw std::runtime_error(os.str());
36  }
37 }
38 
39 Numeric test_integrate_convolved(const Eigen::Ref<Eigen::VectorXcd> F,
40  const Vector& f) {
41  Numeric val = 0.0;
42 
43  const Index n = f.nelem();
44  for (Index i = 0; i < n - 1; i++)
45  val += 0.5 * (f[i + 1] - f[i]) * (F[i].real() + F[i + 1].real());
46 
47  return val; // Should return 1.0 for normalized F
48 }
49 
51  const Array<Index>& sorted_index) {
52  Numeric val = 0.0;
53 
54  const Index n = cosza.nelem();
55  for (Index i = 0; i < n - 1; i++)
56  val += 0.5 * (cosza[sorted_index[i]] - cosza[sorted_index[i + 1]]);
57 
58  return val; // Should return 1.0 for normalized cosza
59 }
60 
62  const Eigen::VectorXcd& F,
63  const Vector& f) {
64  Numeric val = 0.0;
65 
66  const Index n = f.nelem();
67  for (Index i = 0; i < n - 1; i++)
68  val += 0.5 * (f[i + 1] - f[i]) *
69  (I(i, 0) * F[i].real() + I(i + 1, 0) * F[i + 1].real());
70 
71  return val;
72 }
73 
75  const Eigen::VectorXcd& F,
76  const Vector& f) {
77  Numeric val = 0.0;
78 
79  const Index n = f.nelem();
80  for (Index i = 0; i < n - 1; i++)
81  val += 0.5 * (f[i + 1] - f[i]) *
82  (T(i, 0, 0) * F[i].real() + T(i + 1, 0, 0) * F[i + 1].real());
83 
84  return 1.0 - val;
85 }
86 
88  const Vector& cosza,
89  const Array<Index>& sorted_index) {
90  Numeric val = 0.0;
91 
92  const Index n = cosza.nelem();
93  for (Index i = 0; i < n - 1; i++)
94  val += 0.25 * (cosza[sorted_index[i]] - cosza[sorted_index[i + 1]]) *
95  (j[sorted_index[i]] + j[sorted_index[i + 1]]);
96 
97  return val;
98 }
99 
101  if (gp.fd[1] == 1.0)
102  return gp.idx;
103  else
104  return gp.idx + 1;
105 }
106 
108  ArrayOfVector& cosza,
109  const ArrayOfPpath& ppath_field) {
110  extern const Numeric DEG2RAD;
111 
112  Index nalt = 0;
113  for (auto& path : ppath_field)
114  for (auto& gp : path.gp_p) nalt = std::max(nalt, grid_index_from_gp(gp));
115  nalt += 1;
116 
117  // Get the data, which is of unknown length
118  Array<ArrayOfNumeric> zeniths_array(nalt, ArrayOfNumeric(0));
119  for (auto& path : ppath_field) {
120  for (Index ip = 0; ip < path.np; ip++) {
121  const Index ind = grid_index_from_gp(path.gp_p[ip]);
122  zeniths_array[ind].push_back(path.los(ip, 0));
123  }
124  }
125 
126  // Finalize sorting
127  sorted_index.resize(nalt);
128  cosza.resize(nalt);
129  for (Index i = 0; i < nalt; i++) {
130  Vector& data = cosza[i];
131  data.resize(zeniths_array[i].nelem());
132 
133  for (Index j = 0; j < data.nelem(); j++) data[j] = zeniths_array[i][j];
134  get_sorted_indexes(sorted_index[i], data);
135 
136  for (Index j = 0; j < data.nelem(); j++)
137  data[j] = std::cos(DEG2RAD * data[j]);
138  }
139 }
ArrayOfNumeric
Array< Numeric > ArrayOfNumeric
An array of Numeric.
Definition: array.h:48
test_integrate_convolved
Numeric test_integrate_convolved(const Eigen::Ref< Eigen::VectorXcd > F, const Vector &f)
Integrate the line shape.
Definition: radiation_field.cc:39
abs
#define abs(x)
Definition: legacy_continua.cc:20626
sorted_index_of_ppath_field
void sorted_index_of_ppath_field(ArrayOfArrayOfIndex &sorted_index, ArrayOfVector &cosza, const ArrayOfPpath &ppath_field)
Get sorting of zenith angles in field of ppath.
Definition: radiation_field.cc:107
data
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Definition: arts_api_classes.cc:232
DEG2RAD
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
GridPos::fd
Numeric fd[2]
Definition: interpolation.h:75
Array< Index >
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
get_sorted_indexes
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:73
integrate_convolved
Numeric integrate_convolved(const RadiationVector &I, const Eigen::VectorXcd &F, const Vector &f)
Convolve intensity and line shape and integrate.
Definition: radiation_field.cc:61
test_integrate_zenith
Numeric test_integrate_zenith(const Vector &cosza, const Array< Index > &sorted_index)
Integrate cos(za) over the angles.
Definition: radiation_field.cc:50
radiation_field.h
Radiation field calculations.
my_basic_string< char >
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
grid_index_from_gp
Index grid_index_from_gp(const GridPos &gp)
Get a discrete position from grid pos.
Definition: radiation_field.cc:100
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ARTS::Var::ppath_field
ArrayOfPpath ppath_field(Workspace &ws) noexcept
Definition: autoarts.h:5158
max
#define max(a, b)
Definition: legacy_continua.cc:20629
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
RadiationVector
Radiation Vector for Stokes dimension 1-4.
Definition: transmissionmatrix.h:395
GridPos::idx
Index idx
Definition: interpolation.h:74
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Vector
The Vector class.
Definition: matpackI.h:860
TransmissionMatrix
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
Definition: transmissionmatrix.h:38
error_in_integrate
void error_in_integrate(const String &error_msg, const Numeric &value_that_should_be_unity)
Throws an error if integration values are bad.
Definition: radiation_field.cc:30
sorting.h
Contains sorting routines.
integrate_zenith
Numeric integrate_zenith(const VectorView j, const Vector &cosza, const Array< Index > &sorted_index)
Convolve source function with 1D sphere and integrate.
Definition: radiation_field.cc:87