ARTS 2.5.4 (git: 31ce4f0e)
surface.cc
Go to the documentation of this file.
1/* Copyright (C) 2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
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
19/*===========================================================================
20 === File description
21 ===========================================================================*/
22
31/*===========================================================================
32 === External declarations
33 ===========================================================================*/
34
35#include "surface.h"
36#include <cmath>
37#include "auto_md.h"
38#include "check_input.h"
39#include "matpack_complex.h"
40#include "geodetic.h"
41#include "math_funcs.h"
42#include "matpackI.h"
43#include "physics_funcs.h"
44#include "workspace_ng.h"
45
46/*===========================================================================
47 === The functions (in alphabetical order)
48 ===========================================================================*/
49
50// Expression double-checked 210330 (PE)
52 return (180 - abs(rte_los[0]) + abs(specular_los[0])) / 2;
53}
54
56 ConstVectorView z_profile) {
57 Index ip = 0;
58 while (z_surface >= z_profile[ip+1]) {
59 ip++;
60 }
61 return ip;
62}
63
66 ConstMatrixView surface_los,
67 ConstTensor4View surface_rmatrix,
68 ConstMatrixView surface_emission) {
69 // Some sizes
70 const Index nf = I.nrows();
71 const Index stokes_dim = I.ncols();
72 const Index nlos = surface_los.nrows();
73
74 iy = surface_emission;
75
76 // Loop *surface_los*-es. If no such LOS, we are ready.
77 if (nlos > 0) {
78 for (Index ilos = 0; ilos < nlos; ilos++) {
79 Vector rtmp(stokes_dim); // Reflected Stokes vector for 1 frequency
80
81 for (Index iv = 0; iv < nf; iv++) {
82 mult(rtmp, surface_rmatrix(ilos, iv, joker, joker), I(ilos, iv, joker));
83 iy(iv, joker) += rtmp;
84 }
85 }
86 }
87}
88
90 VectorView surface_emission,
91 const Complex& Rv,
92 const Complex& Rh,
93 const Numeric& f,
94 const Index& stokes_dim,
95 const Numeric& surface_skin_t) {
96 ARTS_ASSERT(surface_rmatrix.nrows() == stokes_dim);
97 ARTS_ASSERT(surface_rmatrix.ncols() == stokes_dim);
98 ARTS_ASSERT(surface_emission.nelem() == stokes_dim);
99
100 // Expressions are derived in the surface chapter in the user guide
101
102 surface_rmatrix = 0.0;
103 surface_emission = 0.0;
104
105 Numeric B = planck(f, surface_skin_t);
106
107 const Numeric rv = pow(abs(Rv), 2.0);
108 const Numeric rh = pow(abs(Rh), 2.0);
109 const Numeric rmean = (rv + rh) / 2;
110
111 surface_rmatrix(0, 0) = rmean;
112 surface_emission[0] = B * (1 - rmean);
113
114 if (stokes_dim > 1) {
115 const Numeric rdiff = (rv - rh) / 2;
116
117 surface_rmatrix(1, 0) = rdiff;
118 surface_rmatrix(0, 1) = rdiff;
119 surface_rmatrix(1, 1) = rmean;
120 surface_emission[1] = -B * rdiff;
121
122 if (stokes_dim > 2) {
123 const Complex a = Rh * conj(Rv);
124 const Complex b = Rv * conj(Rh);
125 const Numeric c = real(a + b) / 2.0;
126
127 surface_rmatrix(2, 2) = c;
128
129 if (stokes_dim > 3) {
130 const Numeric d = imag(a - b) / 2.0;
131
132 surface_rmatrix(2, 3) = d;
133 surface_rmatrix(3, 2) = -d;
134 surface_rmatrix(3, 3) = c;
135 }
136 }
137 }
138}
139
140void surface_props_check(const Index& atmosphere_dim,
141 const Vector& lat_grid,
142 const Vector& lon_grid,
143 const Tensor3& surface_props_data,
144 const ArrayOfString& surface_props_names) {
145 // Check sizes
146 ARTS_USER_ERROR_IF (surface_props_data.npages() != surface_props_names.nelem(),
147 "The number of pages in *surface_props_data* and "
148 "length of *surface_props_names* differ.");
149 // If no surface properties, then we are ready
150 if (surface_props_names.nelem() == 0) {
151 return;
152 }
153 ARTS_USER_ERROR_IF (surface_props_data.nrows() !=
154 (atmosphere_dim == 1 ? 1 : lat_grid.nelem()),
155 "Row-size of *surface_props_data* not as expected.");
156 ARTS_USER_ERROR_IF (surface_props_data.ncols() !=
157 (atmosphere_dim <= 2 ? 1 : lon_grid.nelem()),
158 "Column-size of *surface_props_data* not as expected.");
159
160 for (Index i = 0; i < surface_props_names.nelem(); i++) {
161 ARTS_USER_ERROR_IF (surface_props_names[i].nelem() == 0,
162 "Element ", i, " (0-based) of *surface_props_names* is empty.")
163 for (Index j = i + 1; j < surface_props_names.nelem(); j++) {
164 ARTS_USER_ERROR_IF (surface_props_names[j] == surface_props_names[i],
165 "Two surface properties with same name found!\n"
166 "This found for these two properties\n"
167 " index: ", i, '\n',
168 " index: ", j, '\n',
169 " name: ", surface_props_names[i])
170 }
171 }
172}
173
175 const String& vname,
176 const Index& atmosphere_dim,
177 const ArrayOfGridPos& gp_lat,
178 const ArrayOfGridPos& gp_lon,
179 const Matrix& itw,
180 const Tensor3& surface_props_data,
181 const ArrayOfString& surface_props_names) {
182 ARTS_ASSERT(v.nelem() == 1);
183 ARTS_ASSERT(surface_props_data.npages() == surface_props_names.nelem());
184
185 for (Index i = 0; i < surface_props_names.nelem(); i++) {
186 if (surface_props_names[i] == vname) {
188 atmosphere_dim,
189 surface_props_data(i, joker, joker),
190 gp_lat,
191 gp_lon,
192 itw);
193 return;
194 }
195 }
196
198 "The following property was requested\n"
199 " ", vname, '\n',
200 "but it could not be found in *surface_props_names*.")
201}
202
203void dsurface_check(const ArrayOfString& surface_props_names,
204 const ArrayOfString& dsurface_names,
205 const ArrayOfTensor4 dsurface_rmatrix_dx,
206 const ArrayOfMatrix& dsurface_emission_dx) {
207 const Index nq = dsurface_names.nelem();
208
209 ARTS_USER_ERROR_IF (dsurface_rmatrix_dx.nelem() != nq,
210 "The lengths of *dsurface_names* and *dsurface_rmatrix_dx* differ.");
211 ARTS_USER_ERROR_IF (dsurface_emission_dx.nelem() != nq,
212 "The lengths of *dsurface_names* and *dsurface_emission_dx* differ.");
213
214 for (Index i = 0; i < nq; i++) {
215 bool found = false;
216 for (Index j = 0; j < surface_props_names.nelem() && !found; j++) {
217 if (dsurface_names[i] == surface_props_names[j]) {
218 found = true;
219 }
220 }
221 ARTS_USER_ERROR_IF (!found,
222 "String ", i, " (0-based) of *dsurface_names* is \"",
223 dsurface_names[i], "\"\n"
224 "but this string could not be found in *surface_props_names*.\n"
225 "This is likely due to incorrect choice of quantity when\n"
226 " calling *jacobianAddSurfaceQuantity*.")
227 }
228}
229
230
232 Workspace& ws,
233 Matrix& iy_incoming,
234 Index& stars_visible,
235 Vector& specular_los,
236 const Vector& rtp_pos,
237 const Vector& rtp_los,
238 const Index& stokes_dim,
239 const Vector& f_grid,
240 const Index& atmosphere_dim,
241 const Vector& p_grid,
242 const Vector& lat_grid,
243 const Vector& lon_grid,
244 const Tensor3& z_field,
245 const Tensor3& t_field,
246 const EnergyLevelMap& nlte_field,
247 const Tensor4& vmr_field,
248 const ArrayOfArrayOfSpeciesTag& abs_species,
249 const Tensor3& wind_u_field,
250 const Tensor3& wind_v_field,
251 const Tensor3& wind_w_field,
252 const Tensor3& mag_u_field,
253 const Tensor3& mag_v_field,
254 const Tensor3& mag_w_field,
255 const Matrix& z_surface,
256 const Vector& refellipsoid,
257 const Tensor4& pnd_field,
258 const ArrayOfTensor4& dpnd_field_dx,
259 const ArrayOfString& scat_species,
260 const ArrayOfArrayOfSingleScatteringData& scat_data,
261 const Numeric& ppath_lmax,
262 const Numeric& ppath_lraytrace,
263 const Index& ppath_inside_cloudbox_do,
264 const Index& cloudbox_on,
265 const ArrayOfIndex& cloudbox_limits,
266 const Index& gas_scattering_do,
267 const Index& jacobian_do,
268 const ArrayOfRetrievalQuantity& jacobian_quantities,
269 const ArrayOfStar& stars,
270 const Numeric& rte_alonglos_v,
271 const Agenda& propmat_clearsky_agenda,
272 const Agenda& water_p_eq_agenda,
273 const Agenda& gas_scattering_agenda,
274 const Agenda& ppath_step_agenda,
275 const Verbosity& verbosity){
276
277 //Allocate
278 Vector surface_normal;
279 Matrix iy_star_toa;
280
281 //get specular line of sight
282 specular_losCalc(specular_los,
283 surface_normal,
284 rtp_pos,
285 rtp_los,
286 atmosphere_dim,
287 lat_grid,
288 lon_grid,
289 refellipsoid,
290 z_surface,
291 0,
292 verbosity);
293
294 //calculate propagation path from the surface to the space in line of sight
295 Ppath ppath;
296 ppath_calc(ws,
297 ppath,
298 ppath_step_agenda,
299 atmosphere_dim,
300 p_grid,
301 lat_grid,
302 lon_grid,
303 z_field,
304 f_grid,
305 refellipsoid,
306 z_surface,
307 cloudbox_on,
308 cloudbox_limits,
309 rtp_pos,
310 specular_los,
311 ppath_lmax,
312 ppath_lraytrace,
313 ppath_inside_cloudbox_do,
314 verbosity);
315
316
317 //get the incoming spectral radiance of the star at toa. If there is no in
318 //line of sight, then iy_star_toa is simply zero and we are finished. No further
319 //calculations needed.
320 stars_visible=0;
321 get_star_background(iy_star_toa,
322 stars_visible,
323 stars,
324 ppath,
325 f_grid,
326 stokes_dim,
327 atmosphere_dim,
328 refellipsoid);
329
330 if (stars_visible){
331
332 //dummy variables needed for the output and input of iyTransmission
333 ArrayOfMatrix iy_aux_dummy;
334 ArrayOfString iy_aux_vars_dummy;
335 Vector ppvar_p_dummy;
336 Vector ppvar_t_dummy;
337 EnergyLevelMap ppvar_nlte_dummy;
338 Matrix ppvar_vmr_dummy;
339 Matrix ppvar_wind_dummy;
340 Matrix ppvar_mag_dummy;
341 Matrix ppvar_pnd_dummy;
342 Matrix ppvar_f_dummy;
343 Tensor3 ppvar_iy_dummy;
344 Tensor4 ppvar_trans_cumulat_dummy;
345 Tensor4 ppvar_trans_partial_dummy;
346 ArrayOfTensor3 diy_incoming_dummy;
347
348 //Calculate the transmitted radiation from toa to the surface
350 iy_incoming,
351 iy_aux_dummy,
352 diy_incoming_dummy,
353 ppvar_p_dummy,
354 ppvar_t_dummy,
355 ppvar_nlte_dummy,
356 ppvar_vmr_dummy,
357 ppvar_wind_dummy,
358 ppvar_mag_dummy,
359 ppvar_pnd_dummy,
360 ppvar_f_dummy,
361 ppvar_iy_dummy,
362 ppvar_trans_cumulat_dummy,
363 ppvar_trans_partial_dummy,
364 stokes_dim,
365 f_grid,
366 atmosphere_dim,
367 p_grid,
368 t_field,
369 nlte_field,
370 vmr_field,
371 abs_species,
372 wind_u_field,
373 wind_v_field,
374 wind_w_field,
375 mag_u_field,
376 mag_v_field,
377 mag_w_field,
378 cloudbox_on,
379 cloudbox_limits,
380 gas_scattering_do,
381 pnd_field,
382 dpnd_field_dx,
383 scat_species,
384 scat_data,
385 iy_aux_vars_dummy,
386 jacobian_do,
387 jacobian_quantities,
388 ppath,
389 iy_star_toa,
390 propmat_clearsky_agenda,
391 water_p_eq_agenda,
392 gas_scattering_agenda,
393 1,
394 Tensor3(),
395 rte_alonglos_v,
396 verbosity);
397
398 } else {
399 iy_incoming.resize(f_grid.nelem(),stokes_dim);
400 iy_incoming=0;
401 }
402}
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1043
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
A constant view of a Tensor3.
Definition: matpackIII.h:130
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
A constant view of a Tensor4.
Definition: matpackIV.h:131
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The MatrixView class.
Definition: matpackI.h:1164
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4 class.
Definition: matpackIV.h:429
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
Workspace class.
Definition: workspace_ng.h:53
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define abs(x)
void specular_losCalc(Vector &specular_los, Vector &surface_normal, const Vector &rtp_pos, const Vector &rtp_los, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &ignore_surface_slope, const Verbosity &verbosity)
WORKSPACE METHOD: specular_losCalc.
Definition: m_surface.cc:1926
void iyTransmissionStandard(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_pnd, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &gas_scattering_do, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Matrix &iy_transmitter, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Verbosity &)
WORKSPACE METHOD: iyTransmissionStandard.
Implementation of Matrix, Vector, and such stuff.
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
constexpr Numeric imag(Complex c) noexcept
imag
const Joker joker
constexpr Complex conj(Complex c) noexcept
the conjugate of c
std::complex< Numeric > Complex
constexpr Numeric real(Complex c) noexcept
real
Index nelem(const Lines &l)
Number of lines.
Numeric planck(const Numeric &f, const Numeric &t)
planck
This file contains declerations of functions of physical character.
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
This is the core for the WSM ppathStepByStep.
Definition: ppath.cc:5178
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
void interp_atmsurface_by_itw(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates a surface-type variable with pre-calculated weights by interp_atmsurface_gp2itw.
void get_star_background(Matrix &iy, Index &stars_visible, const ArrayOfStar &stars, const Ppath &ppath, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &refellipsoid)
Gets the star background for a given ppath.
Definition: star.cc:116
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
void dsurface_check(const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const ArrayOfTensor4 dsurface_rmatrix_dx, const ArrayOfMatrix &dsurface_emission_dx)
Peforms basic checks of the dsurface variables.
Definition: surface.cc:203
void surface_calc(Matrix &iy, ConstTensor3View I, ConstMatrixView surface_los, ConstTensor4View surface_rmatrix, ConstMatrixView surface_emission)
Weights together downwelling radiation and surface emission.
Definition: surface.cc:64
void surface_props_check(const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names)
Peforms basic checks of surface_props_data and surface_props_names
Definition: surface.cc:140
Numeric calc_incang(ConstVectorView rte_los, ConstVectorView specular_los)
Calculates the incidence angle for a flat surface, based on rte_los and specular_los.
Definition: surface.cc:51
void surface_props_interp(Vector &v, const String &vname, const Index &atmosphere_dim, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, const Matrix &itw, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names)
Peforms an interpolation of surface_props_data
Definition: surface.cc:174
Index index_of_zsurface(const Numeric &z_surface, ConstVectorView z_profile)
Lccates the surface with respect to pressure levels.
Definition: surface.cc:55
void surface_get_incoming_direct(Workspace &ws, Matrix &iy_incoming, Index &stars_visible, Vector &specular_los, const Vector &rtp_pos, const Vector &rtp_los, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Matrix &z_surface, const Vector &refellipsoid, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &ppath_inside_cloudbox_do, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &gas_scattering_do, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfStar &stars, const Numeric &rte_alonglos_v, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
Calculate the incoming direct radiation at the surface for a given line of sight.
Definition: surface.cc:231
void surface_specular_R_and_b(MatrixView surface_rmatrix, VectorView surface_emission, const Complex &Rv, const Complex &Rh, const Numeric &f, const Index &stokes_dim, const Numeric &surface_skin_t)
Sets up the surface reflection matrix and emission vector for the case of specular reflection.
Definition: surface.cc:89
#define d
#define v
#define a
#define c
#define b
This file contains the Workspace class.