ARTS  2.4.0(git:4fb77825)
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,
45  const Tensor4& vmr_field,
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 
55  if (nq)
56  throw std::runtime_error(
57  "Does not support Jacobian calculations at this time");
58  if (stokes_dim not_eq 1)
59  throw std::runtime_error("Only for stokes_dim 1 at this time.");
60 
61  // Compute variables
62  const Vector mag_field = Vector(3, 0);
63  const Vector los = Vector(2, 0);
64  const ArrayOfIndex spec_jac(nq);
65  const Vector tmp(0);
66  ArrayOfStokesVector dS_dx(nq);
67  ArrayOfPropagationMatrix dK_dx(nq);
68 
69  propmat_field = FieldOfPropagationMatrix(
70  nalt, nlat, nlon, PropagationMatrix(nf, stokes_dim));
71  absorption_field =
72  FieldOfStokesVector(nalt, nlat, nlon, StokesVector(nf, stokes_dim));
73  additional_source_field =
74  FieldOfStokesVector(nalt, nlat, nlon, StokesVector(nf, stokes_dim));
75 
76  Workspace l_ws(ws);
77  Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
78 
79 #pragma omp parallel for if (not arts_omp_in_parallel()) schedule(guided) \
80  firstprivate(l_ws, l_propmat_clearsky_agenda)
81  for (Index i = 0; i < nalt; i++) {
82  for (Index j = 0; j < nlat; j++) {
83  for (Index k = 0; k < nlon; k++) {
84  thread_local Index itmp;
86  l_ws,
87  propmat_field(i, j, k),
88  additional_source_field(i, j, k),
89  itmp,
90  dK_dx,
91  dS_dx,
92  l_propmat_clearsky_agenda,
94  f_grid,
95  mag_field,
96  los,
97  nlte_field(i, j, k),
98  vmr_field(joker, i, j, k),
99  t_field(i, j, k),
100  p_grid[i],
101  spec_jac,
102  0);
103  absorption_field(i, j, k) = propmat_field(i, j, k);
104  }
105  }
106  }
107 }
108 
110  const FieldOfPropagationMatrix& propmat_field, const Numeric& r)
111 {
112  FieldOfTransmissionMatrix transmat_field(
113  propmat_field.npages(), propmat_field.nrows(), propmat_field.ncols());
114  for (size_t ip = 0; ip < propmat_field.npages(); ip++)
115  for (size_t ir = 0; ir < propmat_field.nrows(); ir++)
116  for (size_t ic = 0; ic < propmat_field.ncols(); ic++)
117  transmat_field(ip, ir, ic) =
118  TransmissionMatrix(propmat_field(ip, ir, ic), r);
119  return transmat_field;
120 }
121 
123  Workspace& ws,
124  ArrayOfRadiationVector& lvl_rad,
125  ArrayOfRadiationVector& src_rad,
126  ArrayOfTransmissionMatrix& lyr_tra,
127  ArrayOfTransmissionMatrix& tot_tra,
128  const FieldOfPropagationMatrix& propmat_field,
129  const FieldOfStokesVector& absorption_field,
130  const FieldOfStokesVector& additional_source_field,
131  const Vector& f_grid,
132  const Tensor3& t_field,
133  const EnergyLevelMap& nlte_field,
134  const Ppath& ppath,
135  const Agenda& iy_main_agenda,
136  const Agenda& iy_space_agenda,
137  const Agenda& iy_surface_agenda,
138  const Agenda& iy_cloudbox_agenda,
140  const Verbosity& verbosity)
141 {
142  // Size of problem
143  const Index nf = f_grid.nelem();
144  const Index ns = propmat_field(0, 0, 0).StokesDimensions();
145  const Index np = ppath.np;
146 
147  // Current limitations
148  if (ns not_eq 1) throw std::runtime_error("Only for stokes_dim 1");
149  if (ppath.dim not_eq 1) throw std::runtime_error("Only for atmosphere_dim 1");
150 
151  // Size of compute variables
152  lvl_rad = ArrayOfRadiationVector(np, RadiationVector(nf, ns));
153  src_rad = ArrayOfRadiationVector(np, RadiationVector(nf, ns));
154  lyr_tra = ArrayOfTransmissionMatrix(np, TransmissionMatrix(nf, ns));
155 
156  // Size radiative variables always used
157  Vector B(nf);
158  PropagationMatrix K_this(nf, ns), K_past(nf, ns);
159 
160  // Temporary empty variables to fit available function handles
161  Vector vtmp(0);
162  ArrayOfTensor3 t3tmp;
163  ArrayOfTransmissionMatrix tmtmp(0);
164  ArrayOfRadiationVector rvtmp(0);
165  ArrayOfPropagationMatrix pmtmp(0);
166  ArrayOfStokesVector svtmp(0);
167  ArrayOfRetrievalQuantity rqtmp(0);
168 
169  // Loop ppath points and determine radiative properties
170  for (Index ip = 0; ip < np; ip++) {
172  B, vtmp, f_grid, interp_atmfield_by_gp(1, t_field, ppath.gp_p[ip]), 0);
173  K_this = propmat_field(ppath.gp_p[ip]);
174  const StokesVector S(additional_source_field(ppath.gp_p[ip]));
175  const StokesVector a(absorption_field(ppath.gp_p[ip]));
176 
177  if (ip)
178  stepwise_transmission(lyr_tra[ip],
179  tmtmp,
180  tmtmp,
181  K_past,
182  K_this,
183  pmtmp,
184  pmtmp,
185  ppath.lstep[ip - 1],
186  0,
187  0,
188  -1);
189 
190  stepwise_source(src_rad[ip],
191  rvtmp,
192  K_this,
193  a,
194  S,
195  pmtmp,
196  svtmp,
197  svtmp,
198  B,
199  vtmp,
200  rqtmp,
201  0);
202 
203  swap(K_past, K_this);
204  }
205 
206  // In case of backwards RT necessary
207  tot_tra = cumulative_transmission(lyr_tra, CumulativeTransmission::Forward);
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",
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,
249  RadiativeTransferSolver::Emission);
250 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
transmat_field_calc_from_propmat_field
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.
Definition: propmat_field.cc:109
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
ARTS::Group::PropagationMatrix
PropagationMatrix PropagationMatrix
Definition: autoarts.h:93
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1075
ARTS::Var::iy_cloudbox_agenda
Agenda iy_cloudbox_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3742
ARTS::Var::jacobian_quantities
ArrayOfRetrievalQuantity jacobian_quantities(Workspace &ws) noexcept
Definition: autoarts.h:3924
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
joker
const Joker joker
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
swap
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Definition: complex.cc:731
PropagationMatrix
Definition: propagationmatrix.h:87
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ARTS::Var::iy
Matrix iy(Workspace &ws) noexcept
Definition: autoarts.h:3690
field_of_propagation
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.
Definition: propmat_field.cc:35
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
ARTS::Var::iy_space_agenda
Agenda iy_space_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3801
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Field3D::npages
size_t npages() const
Number of pages.
Definition: field.h:170
Agenda
The Agenda class.
Definition: agenda_class.h:44
FieldOfStokesVector
Field3D< StokesVector > FieldOfStokesVector
Definition: propmat_field.h:39
emission_from_propmat_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.
Definition: propmat_field.cc:122
ARTS::Var::nlte_field
EnergyLevelMap nlte_field(Workspace &ws) noexcept
Definition: autoarts.h:4604
Array
This can be used to make arrays out of anything.
Definition: array.h:108
get_stepwise_clearsky_propmat
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, ConstVectorView ppath_f_grid, ConstVectorView ppath_magnetic_field, ConstVectorView ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, ConstVectorView ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const ArrayOfIndex &jacobian_species, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition: rte.cc:1322
stepwise_transmission
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.
Definition: transmissionmatrix.cc:1353
stepwise_source
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.
Definition: transmissionmatrix.cc:1371
ArrayOfRadiationVector
Array< RadiationVector > ArrayOfRadiationVector
Definition: transmissionmatrix.h:899
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ARTS::Var::ppath
Ppath ppath(Workspace &ws) noexcept
Definition: autoarts.h:5139
ARTS::Var::surface_props_data
Tensor3 surface_props_data(Workspace &ws) noexcept
Definition: autoarts.h:6745
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
EnergyLevelMap
Definition: energylevelmap.h:60
update_radiation_vector
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 RadiativeTransferSolver solver)
Update the Radiation Vector.
Definition: transmissionmatrix.cc:1484
Field3D::nrows
size_t nrows() const
Number of rows.
Definition: field.h:173
ARTS::Var::iy_surface_agenda
Agenda iy_surface_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3808
ArrayOfTransmissionMatrix
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
Definition: transmissionmatrix.h:896
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition: transmissionmatrix.cc:1527
FieldOfPropagationMatrix
Field3D< PropagationMatrix > FieldOfPropagationMatrix
Definition: propmat_field.h:38
Field3D::ncols
size_t ncols() const
Number of columns.
Definition: field.h:176
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
Field3D
Creates a 3D field of a base unit.
Definition: field.h:33
ARTS::Var::iy_main_agenda
Agenda iy_main_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3794
ARTS::Group::RadiationVector
RadiationVector RadiationVector
Definition: autoarts.h:95
interp_atmfield_by_gp
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.
Definition: special_interp.cc:126
rte.h
Declaration of functions in rte.cc.
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
special_interp.h
Header file for special_interp.cc.
ARTS::Group::TransmissionMatrix
TransmissionMatrix TransmissionMatrix
Definition: autoarts.h:112
ARTS::Var::propmat_clearsky_agenda
Agenda propmat_clearsky_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5405
propmat_field.h
Implements a propagation matrix field.
transmissionmatrix.h
Stuff related to the transmission matrix.
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
get_iy_of_background
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, ConstTensor3View iy_transmission, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, ConstVectorView rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, ConstVectorView f_grid, const String &iy_unit, 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:916
ARTS::Group::StokesVector
StokesVector StokesVector
Definition: autoarts.h:101
get_stepwise_blackbody_radiation
void get_stepwise_blackbody_radiation(VectorView B, VectorView dB_dT, ConstVectorView ppath_f_grid, const Numeric &ppath_temperature, const bool &do_temperature_derivative)
Get the blackbody radiation at propagation path point.
Definition: rte.cc:1307
ns
#define ns
Definition: legacy_continua.cc:22183