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