ARTS 2.5.4 (git: 4c0d3b4d)
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
74 Workspace l_ws(ws);
75 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
76
77#pragma omp parallel for if (not arts_omp_in_parallel()) schedule(guided) \
78 firstprivate(l_ws, l_propmat_clearsky_agenda)
79 for (Index i = 0; i < nalt; i++) {
80 for (Index j = 0; j < nlat; j++) {
81 for (Index k = 0; k < nlon; k++) {
82 thread_local Index itmp;
84 l_ws,
85 propmat_field(i, j, k),
86 additional_source_field(i, j, k),
87 itmp,
88 dK_dx,
89 dS_dx,
90 l_propmat_clearsky_agenda,
91 jacobian_quantities,
92 f_grid,
93 mag_field,
94 los,
95 nlte_field(i, j, k),
96 vmr_field(joker, i, j, k),
97 t_field(i, j, k),
98 p_grid[i],
99 false);
100 absorption_field(i, j, k) = propmat_field(i, j, k);
101 }
102 }
103 }
104}
105
107 const FieldOfPropagationMatrix& propmat_field, const Numeric& r)
108{
109 FieldOfTransmissionMatrix transmat_field(
110 propmat_field.npages(), propmat_field.nrows(), propmat_field.ncols());
111 for (size_t ip = 0; ip < propmat_field.npages(); ip++)
112 for (size_t ir = 0; ir < propmat_field.nrows(); ir++)
113 for (size_t ic = 0; ic < propmat_field.ncols(); ic++)
114 transmat_field(ip, ir, ic) =
115 TransmissionMatrix(propmat_field(ip, ir, ic), r);
116 return transmat_field;
117}
118
120 Workspace& ws,
121 ArrayOfRadiationVector& lvl_rad,
122 ArrayOfRadiationVector& src_rad,
125 const FieldOfPropagationMatrix& propmat_field,
126 const FieldOfStokesVector& absorption_field,
127 const FieldOfStokesVector& additional_source_field,
128 const Vector& f_grid,
129 const Tensor3& t_field,
130 const EnergyLevelMap& nlte_field,
131 const Ppath& ppath,
132 const Agenda& iy_main_agenda,
133 const Agenda& iy_space_agenda,
134 const Agenda& iy_surface_agenda,
135 const Agenda& iy_cloudbox_agenda,
136 const Tensor3& surface_props_data,
137 const Verbosity& verbosity)
138{
139 // Size of problem
140 const Index nf = f_grid.nelem();
141 const Index ns = propmat_field(0, 0, 0).StokesDimensions();
142 const Index np = ppath.np;
143
144 // Current limitations
145 ARTS_USER_ERROR_IF (ns not_eq 1, "Only for stokes_dim 1");
146 ARTS_USER_ERROR_IF (ppath.dim not_eq 1, "Only for atmosphere_dim 1");
147
148 // Size of compute variables
149 lvl_rad = ArrayOfRadiationVector(np, RadiationVector(nf, ns));
150 src_rad = ArrayOfRadiationVector(np, RadiationVector(nf, ns));
152
153 // Size radiative variables always used
154 Vector B(nf);
155 PropagationMatrix K_this(nf, ns), K_past(nf, ns);
156
157 // Temporary empty variables to fit available function handles
158 Vector vtmp(0);
159 ArrayOfTensor3 t3tmp;
161 ArrayOfRadiationVector rvtmp(0);
163 ArrayOfStokesVector svtmp(0);
165
166 // Loop ppath points and determine radiative properties
167 for (Index ip = 0; ip < np; ip++) {
169 B, vtmp, f_grid, interp_atmfield_by_gp(1, t_field, ppath.gp_p[ip]), false);
170 K_this = propmat_field(ppath.gp_p[ip]);
171 const StokesVector S(additional_source_field(ppath.gp_p[ip]));
172 const StokesVector a(absorption_field(ppath.gp_p[ip]));
173
174 if (ip)
175 stepwise_transmission(lyr_tra[ip],
176 tmtmp,
177 tmtmp,
178 K_past,
179 K_this,
180 pmtmp,
181 pmtmp,
182 ppath.lstep[ip - 1],
183 0,
184 0,
185 -1);
186
187 stepwise_source(src_rad[ip],
188 rvtmp,
189 K_this,
190 a,
191 S,
192 pmtmp,
193 svtmp,
194 svtmp,
195 B,
196 vtmp,
197 rqtmp,
198 false);
199
200 swap(K_past, K_this);
201 }
202
203 // In case of backwards RT necessary
205 const Tensor3 iy_trans_new = tot_tra[np - 1];
206
207 // Radiative background
208 Matrix iy;
210 iy,
211 t3tmp,
212 iy_trans_new,
213 0,
214 0,
215 rqtmp,
216 ppath,
217 {0},
218 1,
219 nlte_field,
220 0,
221 1,
222 f_grid,
223 "1",
224 surface_props_data,
225 iy_main_agenda,
226 iy_space_agenda,
227 iy_surface_agenda,
228 iy_cloudbox_agenda,
229 1,
230 verbosity);
231 lvl_rad[np - 1] = iy;
232
233 // Radiative transfer calculations
234 for (Index ip = np - 2; ip >= 0; ip--)
235 update_radiation_vector(lvl_rad[ip] = lvl_rad[ip + 1],
236 rvtmp,
237 rvtmp,
238 src_rad[ip],
239 src_rad[ip + 1],
240 rvtmp,
241 rvtmp,
242 lyr_tra[ip + 1],
243 tot_tra[ip],
244 tmtmp,
245 tmtmp,
250 Numeric(),
251 Vector(),
252 Vector(),
253 0,
254 0,
256}
The Agenda class.
Definition: agenda_class.h:51
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
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:1270
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:344
The Tensor4 class.
Definition: matpackIV.h:427
The Vector class.
Definition: matpackI.h:908
Workspace class.
Definition: workspace_ng.h:40
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define ns
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)
Swaps two objects.
const Joker joker
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
#define a
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:1131
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:728
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:1116
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.
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, 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.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
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
Array< TransmissionMatrix > ArrayOfTransmissionMatrix