ARTS 2.5.10 (git: 2f1c442c)
propmat_field.cc
Go to the documentation of this file.
1/* Copyright (C) 2019 Richard Larsson <ric.larsson@gmail.com>
2 *
3 T his pr*ogram is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
30#include "propmat_field.h"
31#include "rte.h"
32#include "special_interp.h"
33#include "transmissionmatrix.h"
34
36 FieldOfPropagationMatrix& propmat_field,
37 FieldOfStokesVector& absorption_field,
38 FieldOfStokesVector& additional_source_field,
39 const Index& stokes_dim,
40 const Vector& f_grid,
41 const Vector& p_grid,
42 const Tensor3& z_field,
43 const Tensor3& t_field,
44 const EnergyLevelMap& nlte_field,
45 const Tensor4& vmr_field,
46 const ArrayOfRetrievalQuantity& jacobian_quantities,
47 const Agenda& propmat_clearsky_agenda)
48{
49 const Index nalt = z_field.npages();
50 const Index nlat = z_field.nrows();
51 const Index nlon = z_field.ncols();
52 const Index nq = jacobian_quantities.nelem();
53 const Index nf = f_grid.nelem();
54
56 "Does not support Jacobian calculations at this time");
57 ARTS_USER_ERROR_IF (stokes_dim not_eq 1,
58 "Only for stokes_dim 1 at this time.");
59
60 // Compute variables
61 const Vector mag_field = Vector(3, 0);
62 const Vector los = Vector(2, 0);
63 const Vector tmp(0);
64 ArrayOfStokesVector dS_dx(nq);
66
67 propmat_field = FieldOfPropagationMatrix(
68 nalt, nlat, nlon, PropagationMatrix(nf, stokes_dim));
69 absorption_field =
70 FieldOfStokesVector(nalt, nlat, nlon, StokesVector(nf, stokes_dim));
71 additional_source_field =
72 FieldOfStokesVector(nalt, nlat, nlon, StokesVector(nf, stokes_dim));
73
75
76#pragma omp parallel for if (not arts_omp_in_parallel()) collapse(3) \
77 firstprivate(wss)
78 for (Index i = 0; i < nalt; i++) {
79 for (Index j = 0; j < nlat; j++) {
80 for (Index k = 0; k < nlon; k++) {
81 thread_local Index itmp;
83 wss,
84 propmat_field(i, j, k),
85 additional_source_field(i, j, k),
86 itmp,
87 dK_dx,
88 dS_dx,
89 propmat_clearsky_agenda,
90 jacobian_quantities,
91 f_grid,
92 mag_field,
93 los,
94 nlte_field(i, j, k),
95 vmr_field(joker, i, j, k),
96 t_field(i, j, k),
97 p_grid[i],
98 false);
99 absorption_field(i, j, k) = propmat_field(i, j, k);
100 }
101 }
102 }
103}
104
106 const FieldOfPropagationMatrix& propmat_field, const Numeric& r)
107{
108 FieldOfTransmissionMatrix transmat_field(
109 propmat_field.npages(), propmat_field.nrows(), propmat_field.ncols());
110 for (size_t ip = 0; ip < propmat_field.npages(); ip++)
111 for (size_t ir = 0; ir < propmat_field.nrows(); ir++)
112 for (size_t ic = 0; ic < propmat_field.ncols(); ic++)
113 transmat_field(ip, ir, ic) =
114 TransmissionMatrix(propmat_field(ip, ir, ic), r);
115 return transmat_field;
116}
117
119 Workspace& ws,
120 ArrayOfRadiationVector& lvl_rad,
121 ArrayOfRadiationVector& src_rad,
124 const FieldOfPropagationMatrix& propmat_field,
125 const FieldOfStokesVector& absorption_field,
126 const FieldOfStokesVector& additional_source_field,
127 const Vector& f_grid,
128 const Tensor3& t_field,
129 const EnergyLevelMap& nlte_field,
130 const Ppath& ppath,
131 const Agenda& iy_main_agenda,
132 const Agenda& iy_space_agenda,
133 const Agenda& iy_surface_agenda,
134 const Agenda& iy_cloudbox_agenda,
135 const Tensor3& surface_props_data,
136 const Verbosity& verbosity)
137{
138 // Size of problem
139 const Index nf = f_grid.nelem();
140 const Index ns = propmat_field(0, 0, 0).StokesDimensions();
141 const Index np = ppath.np;
142
143 // Current limitations
144 ARTS_USER_ERROR_IF (ns not_eq 1, "Only for stokes_dim 1");
145 ARTS_USER_ERROR_IF (ppath.dim not_eq 1, "Only for atmosphere_dim 1");
146
147 // Size of compute variables
148 lvl_rad = ArrayOfRadiationVector(np, RadiationVector(nf, ns));
149 src_rad = ArrayOfRadiationVector(np, RadiationVector(nf, ns));
150 lyr_tra = ArrayOfTransmissionMatrix(np, TransmissionMatrix(nf, ns));
151
152 RadiationVector J_add_dummy;
153 ArrayOfRadiationVector dJ_add_dummy;
154
155 // Size radiative variables always used
156 Vector B(nf);
157 PropagationMatrix K_this(nf, ns), K_past(nf, ns);
158
159 // Temporary empty variables to fit available function handles
160 Vector vtmp(0);
161 ArrayOfTensor3 t3tmp;
163 ArrayOfRadiationVector rvtmp(0);
165 ArrayOfStokesVector svtmp(0);
167
168 // Loop ppath points and determine radiative properties
169 for (Index ip = 0; ip < np; ip++) {
171 B, vtmp, f_grid, interp_atmfield_by_gp(1, t_field, ppath.gp_p[ip]), false);
172 K_this = propmat_field(ppath.gp_p[ip]);
173 const StokesVector S(additional_source_field(ppath.gp_p[ip]));
174 const StokesVector a(absorption_field(ppath.gp_p[ip]));
175
176 if (ip)
177 stepwise_transmission(lyr_tra[ip],
178 tmtmp,
179 tmtmp,
180 K_past,
181 K_this,
182 pmtmp,
183 pmtmp,
184 ppath.lstep[ip - 1],
185 0,
186 0,
187 -1);
188
189 stepwise_source(src_rad[ip],
190 rvtmp,
191 J_add_dummy,
192 K_this,
193 a,
194 S,
195 pmtmp,
196 svtmp,
197 svtmp,
198 B,
199 vtmp,
200 rqtmp,
201 false);
202
203 swap(K_past, K_this);
204 }
205
206 // In case of backwards RT necessary
208 const Tensor3 iy_trans_new = tot_tra[np - 1];
209
210 // Radiative background
211 Matrix iy;
213 iy,
214 t3tmp,
215 iy_trans_new,
216 0,
217 0,
218 rqtmp,
219 ppath,
220 {0},
221 1,
222 nlte_field,
223 0,
224 1,
225 f_grid,
226 "1",
227 surface_props_data,
228 iy_main_agenda,
229 iy_space_agenda,
230 iy_surface_agenda,
231 iy_cloudbox_agenda,
232 1,
233 verbosity);
234 lvl_rad[np - 1] = iy;
235
236 // Radiative transfer calculations
237 for (Index ip = np - 2; ip >= 0; ip--)
238 update_radiation_vector(lvl_rad[ip] = lvl_rad[ip + 1],
239 rvtmp,
240 rvtmp,
241 src_rad[ip],
242 src_rad[ip + 1],
243 rvtmp,
244 rvtmp,
245 lyr_tra[ip + 1],
246 tot_tra[ip],
247 tmtmp,
248 tmtmp,
253 Numeric(),
254 Vector(),
255 Vector(),
256 0,
257 0,
259}
The Agenda class.
Definition: agenda_class.h:69
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
Creates a 3D field of a base unit.
Definition: field.h:33
size_t ncols() const
Number of columns.
Definition: field.h:176
size_t npages() const
Number of pages.
Definition: field.h:170
size_t nrows() const
Number of rows.
Definition: field.h:173
The Matrix class.
Definition: matpackI.h:1285
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:352
The Tensor4 class.
Definition: matpackIV.h:435
The Vector class.
Definition: matpackI.h:910
Workspace class.
Definition: workspace_ng.h:53
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
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
void swap(ComplexVector &v1, ComplexVector &v2) noexcept
Swaps two objects.
const Joker joker
Array< PropagationMatrix > ArrayOfPropagationMatrix
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.
Field3D< StokesVector > FieldOfStokesVector
Definition: propmat_field.h:43
Field3D< PropagationMatrix > FieldOfPropagationMatrix
Definition: propmat_field.h:42
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_magnetic_field, const ConstVectorView &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const ConstVectorView &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition: rte.cc:1135
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const ConstTensor3View &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const ConstVectorView &rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const ConstVectorView &f_grid, const String &iy_unit, const ConstTensor3View &surface_props_data, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Verbosity &verbosity)
Determines iy of the "background" of a propgation path.
Definition: rte.cc:732
void get_stepwise_blackbody_radiation(VectorView B, VectorView dB_dT, const ConstVectorView &ppath_f_grid, const Numeric &ppath_temperature, const bool &do_temperature_derivative)
Get the blackbody radiation at propagation path point.
Definition: rte.cc:1120
Declaration of functions in rte.cc.
void interp_atmfield_by_gp(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates an atmospheric field given the grid positions.
Header file for special_interp.cc.
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
Vector lstep
The length between ppath points.
Definition: ppath_struct.h:39
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath_struct.h:51
Index dim
Atmospheric dimensionality.
Definition: ppath_struct.h:19
Radiation Vector for Stokes dimension 1-4.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, RadiationVector &J_add, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK, const ArrayOfStokesVector &da, const ArrayOfStokesVector &dS, const ConstVectorView &B, const ConstVectorView &dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Set the stepwise source.
void update_radiation_vector(RadiationVector &I, ArrayOfRadiationVector &dI1, ArrayOfRadiationVector &dI2, const RadiationVector &J1, const RadiationVector &J2, const ArrayOfRadiationVector &dJ1, const ArrayOfRadiationVector &dJ2, const TransmissionMatrix &T, const TransmissionMatrix &PiT, const ArrayOfTransmissionMatrix &dT1, const ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the Radiation Vector.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
Stuff related to the transmission matrix.
Array< RadiationVector > ArrayOfRadiationVector
#define a
Array< TransmissionMatrix > ArrayOfTransmissionMatrix