ARTS 2.5.4 (git: 31ce4f0e)
m_radiation_field.cc
Go to the documentation of this file.
1/* Copyright (C) 2015
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 "absorption.h"
28#include "arts.h"
29#include "arts_omp.h"
30#include "auto_md.h"
31#include "debug.h"
32#include "lineshape.h"
33#include "logic.h"
34#include "physics_funcs.h"
35#include "ppath.h"
36#include "propmat_field.h"
37#include "radiation_field.h"
38
40 Workspace& ws,
41 Matrix& line_irradiance,
42 Tensor3& line_transmission,
43 const ArrayOfArrayOfSpeciesTag& abs_species,
44 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
45 const EnergyLevelMap& nlte_field,
46 const Tensor4& vmr_field,
47 const Tensor3& t_field,
48 const Tensor3& z_field,
49 const Vector& p_grid,
50 const Vector& refellipsoid,
51 const Tensor3& surface_props_data,
52 const Agenda& ppath_agenda,
53 const Agenda& iy_main_agenda,
54 const Agenda& iy_space_agenda,
55 const Agenda& iy_surface_agenda,
56 const Agenda& iy_cloudbox_agenda,
57 const Agenda& propmat_clearsky_agenda,
58 const Numeric& df,
59 const Index& nz,
60 const Index& nf,
61 const Numeric& r,
62 const Verbosity& verbosity)
63{
64 ARTS_USER_ERROR_IF (abs_lines_per_species.nelem() not_eq 1,
65 "Only for one species...");
66 ARTS_USER_ERROR_IF (nf % 2 not_eq 1,
67 "Must hit line center, nf % 2 must be 1.");
68 const Index nl = nelem(abs_lines_per_species);
69 const Index np = p_grid.nelem();
70
71 // Compute variables
72 ArrayOfTensor3 diy_dx;
73 FieldOfPropagationMatrix propmat_field;
74 FieldOfStokesVector absorption_field, additional_source_field;
75 ArrayOfPpath ppath_field;
76
77 // Check that the lines and nf is correct
78 Vector f_grid(nf * nl);
79 Index il=0;
80 for (auto& lines: abs_lines_per_species) {
81 for (auto& band: lines) {
82 for (Index k=0; k<band.NumLines(); k++) {
83 nlinspace(f_grid[Range(il * nf, nf)], band.lines[k].F0 * (1 - df), band.lines[k].F0 * (1 + df), nf);
84 il++;
85 }
86 }
87 }
88
89 ARTS_USER_ERROR_IF(not is_increasing(f_grid), "Frequency grid is not increasing, abs_lines_per_species must be sorted and no overlap is allowed");
90
92 ppath_field,
93 ppath_agenda,
94 -1,
95 1e99,
96 1,
97 z_field,
98 f_grid,
99 0,
100 1,
101 0,
102 Vector(1, 0),
103 Vector(1, 0),
104 Vector(0),
105 refellipsoid,
106 1,
107 nz,
108 verbosity);
109
110 ArrayOfArrayOfIndex sorted_index;
111 ArrayOfVector cos_zenith_angles;
112 sorted_index_of_ppath_field(sorted_index, cos_zenith_angles, ppath_field);
113
114 for (Index ip = 0; ip < np; ip++)
116 "Your lineshape integration does normalize. Increase nf and decrease df until it does.",
117 test_integrate_zenith(cos_zenith_angles[ip], sorted_index[ip]));
118
120 propmat_field,
121 absorption_field,
122 additional_source_field,
123 1,
124 f_grid,
125 p_grid,
126 z_field,
127 t_field,
128 nlte_field,
129 vmr_field,
131 propmat_clearsky_agenda);
132
134 nl, Array<Eigen::VectorXcd>(np, Eigen::VectorXcd(nf * nl)));
135 for (Index ip=0; ip<np; ip++) {
136 il=0;
137 for (auto& lines: abs_lines_per_species) {
138 for (auto& band: lines) {
139 const Numeric DC = band.DopplerConstant(t_field(ip, 0, 0));
140 const Vector vmrs = band.BroadeningSpeciesVMR(vmr_field(joker, ip, 0, 0), abs_species);
141 for (Index k=0; k<band.NumLines(); k++) {
142 const auto X = band.ShapeParameters(k, t_field(ip, 0, 0), p_grid[ip], vmrs);
143 LineShape::Calculator ls(band.lineshapetype, band.lines[k].F0, X, DC, 0, false);
144 for (Index iv=0; iv<f_grid.nelem(); iv++) {
145 lineshapes[il][ip][iv] = ls(f_grid[iv]);
146 }
147 il++;
148 }
149 }
150 }
151 }
152 for (auto& aols : lineshapes)
153 for (auto& ls : aols)
155 "Your lineshape integration does not normalize. Increase nf and decrease df until it does.",
156 test_integrate_convolved(ls, f_grid));
157
158 // Counting the path index so we can make the big loop parallel
159 ArrayOfArrayOfIndex counted_path_index(ppath_field.nelem());
160 ArrayOfIndex counter(np, 0);
161 for (Index i = 0; i < ppath_field.nelem(); i++) {
162 const Ppath& path = ppath_field[i];
163 counted_path_index[i].resize(path.np);
164 for (Index ip_path = 0; ip_path < path.np; ip_path++) {
165 const Index ip_grid = grid_index_from_gp(path.gp_p[ip_path]);
166 counted_path_index[i][ip_path] = counter[ip_grid];
167 counter[ip_grid]++;
168 }
169 }
170
171 line_irradiance = Matrix(nl, np, 0.0);
172 line_transmission = Tensor3(1, nl, np, 0.0);
173
174 ArrayOfMatrix line_radiance(np);
175 for (Index i = 0; i < np; i++)
176 line_radiance[i].resize(sorted_index[i].nelem(), nl);
177
179
180#pragma omp parallel for if (not arts_omp_in_parallel()) \
181 schedule(guided) default(shared) firstprivate(wss, il)
182 for (Index i = 0; i < ppath_field.nelem(); i++) {
183 const Ppath& path = ppath_field[i];
184
185 thread_local ArrayOfRadiationVector lvl_rad;
186 thread_local ArrayOfRadiationVector src_rad;
187 thread_local ArrayOfTransmissionMatrix lyr_tra;
188 thread_local ArrayOfTransmissionMatrix tot_tra;
189
191 lvl_rad,
192 src_rad,
193 lyr_tra,
194 tot_tra,
195 propmat_field,
196 absorption_field,
197 additional_source_field,
198 f_grid,
199 t_field,
200 nlte_field,
201 path,
202 iy_main_agenda,
203 iy_space_agenda,
204 iy_surface_agenda,
205 iy_cloudbox_agenda,
206 surface_props_data,
207 verbosity);
208
209 for (Index ip_path = 0; ip_path < path.np; ip_path++) {
210 const Index ip_grid = grid_index_from_gp(path.gp_p[ip_path]);
211 for (il = 0; il < nl; il++)
212 line_radiance[ip_grid](counted_path_index[i][ip_path], il) =
214 lvl_rad[ip_path], lineshapes[il][ip_grid], f_grid);
215 }
216 }
217
218 for (Index ip = 0; ip < np; ip++) {
219 for (il = 0; il < nl; il++) {
220 line_irradiance(il, ip) = integrate_zenith(line_radiance[ip](joker, il),
221 cos_zenith_angles[ip],
222 sorted_index[ip]);
223 }
224 }
225
226 if (r > 0) {
227 const FieldOfTransmissionMatrix transmat_field =
229 for (Index ip = 0; ip < np; ip++)
230 for (il = 0; il < nl; il++)
231 line_transmission(0, il, ip) = integrate_convolved(
232 transmat_field(ip, 0, 0), lineshapes[il][ip], f_grid);
233 }
234}
235
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
Header file for helper functions for OpenMP.
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
Creates a 3D field of a base unit.
Definition: field.h:33
Line shape calculator.
Definition: lineshape.h:647
The Matrix class.
Definition: matpackI.h:1261
The range class.
Definition: matpackI.h:159
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4 class.
Definition: matpackIV.h:429
The Vector class.
Definition: matpackI.h:899
Workspace class.
Definition: workspace_ng.h:53
Helper macros for debugging.
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:539
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:218
Header file for logic.cc.
void ppath_fieldFromDownUpLimbGeoms(Workspace &ws, ArrayOfPpath &ppath_field, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Tensor3 &z_field, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Vector &refellipsoid, const Index &atmosphere_dim, const Index &zenith_angles_per_position, const Verbosity &verbosity)
WORKSPACE METHOD: ppath_fieldFromDownUpLimbGeoms.
Definition: m_ppath.cc:1554
void line_irradianceCalcForSingleSpeciesNonOverlappingLinesPseudo2D(Workspace &ws, Matrix &line_irradiance, Tensor3 &line_transmission, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const Tensor3 &t_field, const Tensor3 &z_field, const Vector &p_grid, const Vector &refellipsoid, const Tensor3 &surface_props_data, const Agenda &ppath_agenda, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Agenda &propmat_clearsky_agenda, const Numeric &df, const Index &nz, const Index &nf, const Numeric &r, const Verbosity &verbosity)
WORKSPACE METHOD: line_irradianceCalcForSingleSpeciesNonOverlappingLinesPseudo2D.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
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
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric k
Boltzmann constant convenience name [J/K].
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
This file contains declerations of functions of physical character.
Propagation path structure and functions.
void field_of_propagation(Workspace &ws, FieldOfPropagationMatrix &propmat_field, FieldOfStokesVector &absorption_field, FieldOfStokesVector &additional_source_field, const Index &stokes_dim, const Vector &f_grid, const Vector &p_grid, const Tensor3 &z_field, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda)
Creates a field of propagation matrices, absorption vectors, and source vectors.
FieldOfTransmissionMatrix transmat_field_calc_from_propmat_field(const FieldOfPropagationMatrix &propmat_field, const Numeric &r)
Get a field of transmission matrices from the propagation matrix field.
void emission_from_propmat_field(Workspace &ws, ArrayOfRadiationVector &lvl_rad, ArrayOfRadiationVector &src_rad, ArrayOfTransmissionMatrix &lyr_tra, ArrayOfTransmissionMatrix &tot_tra, const FieldOfPropagationMatrix &propmat_field, const FieldOfStokesVector &absorption_field, const FieldOfStokesVector &additional_source_field, const Vector &f_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Ppath &ppath, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Tensor3 &surface_props_data, const Verbosity &verbosity)
Computes the radiation and transmission from fields of atmospheric propagation.
Implements a propagation matrix 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.
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.
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath_struct.h:51