ARTS 2.5.4 (git: 31ce4f0e)
m_surface.cc
Go to the documentation of this file.
1/* Copyright (C) 2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3 Stefan Buehler <sbuehler@ltu.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
20/*===========================================================================
21 === File description
22 ===========================================================================*/
23
35/*===========================================================================
36 === External declarations
37 ===========================================================================*/
38
39#include <cmath>
40#include "array.h"
41#include "arts.h"
42#include "arts_constants.h"
43#include "auto_md.h"
44#include "check_input.h"
45#include "matpack_complex.h"
46#include "fastem.h"
47#include "geodetic.h"
48#include "interpolation.h"
50#include "math_funcs.h"
51#include "messages.h"
52#include "physics_funcs.h"
53#include "ppath.h"
54#include "rte.h"
55#include "special_interp.h"
56#include "surface.h"
57#include "tessem.h"
58#include "arts_conversions.h"
59#include "gas_scattering.h"
60
61using Constant::pi;
63
65inline constexpr Numeric DEG2RAD = Conversion::deg2rad(1);
66
67/*===========================================================================
68 === The functions (in alphabetical order)
69 ===========================================================================*/
70
71/* Workspace method: Doxygen documentation will be auto-generated */
72void FastemStandAlone(Matrix& emissivity,
73 Matrix& reflectivity,
74 const Vector& f_grid,
75 const Numeric& surface_skin_t,
76 const Numeric& za,
77 const Numeric& salinity,
78 const Numeric& wind_speed,
79 const Numeric& rel_aa,
80 const Vector& transmittance,
81 const Index& fastem_version,
82 const Verbosity&) {
83 const Index nf = f_grid.nelem();
84
85 chk_if_in_range("zenith angle", za, 90, 180);
86 chk_if_in_range_exclude("surface skin temperature", surface_skin_t, 260, 373);
87 chk_if_in_range_exclude_high("salinity", salinity, 0, 1);
88 chk_if_in_range_exclude_high("wind speed", wind_speed, 0, 100);
89 chk_if_in_range("azimuth angle", rel_aa, -180, 180);
90 chk_vector_length("transmittance", "f_grid", transmittance, f_grid);
91 if (fastem_version < 3 || fastem_version > 6)
92 throw std::runtime_error(
93 "Invalid fastem version: 3 <= fastem_version <= 6");
94
95 emissivity.resize(nf, 4);
96 reflectivity.resize(nf, 4);
97
98 const Numeric t = max(surface_skin_t, Numeric(270));
99
100 for (Index i = 0; i < nf; i++) {
101 if (f_grid[i] > 250e9)
102 throw std::runtime_error("Only frequency <= 250 GHz are allowed");
103 chk_if_in_range("transmittance", transmittance[i], 0, 1);
104
105 Vector e, r;
106
107 fastem(e,
108 r,
109 f_grid[i],
110 za,
111 t,
112 salinity,
113 wind_speed,
114 transmittance[i],
115 rel_aa,
116 fastem_version);
117
118 emissivity(i, joker) = e;
119 reflectivity(i, joker) = r;
120 }
121
122 // FASTEM does not work close to the horizon (at least v6). Make sure values
123 // are inside [0,1]. Then seems best to make sure that e+r=1.
124 for (Index i = 0; i < nf; i++) {
125 for (Index s = 0; s < 2; s++) {
126 if (emissivity(i, s) > 1) {
127 emissivity(i, s) = 1;
128 reflectivity(i, s) = 0;
129 }
130 if (emissivity(i, s) < 0) {
131 emissivity(i, s) = 0;
132 reflectivity(i, s) = 1;
133 }
134 if (reflectivity(i, s) > 1) {
135 emissivity(i, s) = 0;
136 reflectivity(i, s) = 1;
137 }
138 if (reflectivity(i, s) < 0) {
139 emissivity(i, s) = 1;
140 reflectivity(i, s) = 0;
141 }
142 }
143 }
144}
145
146/* Workspace method: Doxygen documentation will be auto-generated */
148 const Index& atmosphere_dim,
149 const Vector& lat_grid,
150 const Vector& lat_true,
151 const Vector& lon_true,
152 const Vector& rtp_pos,
153 const GriddedField2& gfield2,
154 const Verbosity&) {
155 // Set expected order of grids
156 Index gfield_latID = 0;
157 Index gfield_lonID = 1;
158
159 // Basic checks and sizes
160 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
161 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
162 chk_rte_pos(atmosphere_dim, rtp_pos);
163 gfield2.checksize_strict();
164 //
165 chk_griddedfield_gridname(gfield2, gfield_latID, "Latitude");
166 chk_griddedfield_gridname(gfield2, gfield_lonID, "Longitude");
167 //
168 const Index nlat = gfield2.data.nrows();
169 const Index nlon = gfield2.data.ncols();
170 //
171 if (nlat < 2 || nlon < 2) {
172 ostringstream os;
173 os << "The data in *gfield2* must span a geographical region. That is,\n"
174 << "the latitude and longitude grids must have a length >= 2.";
175 throw runtime_error(os.str());
176 }
177
178 const Vector& GFlat = gfield2.get_numeric_grid(gfield_latID);
179 const Vector& GFlon = gfield2.get_numeric_grid(gfield_lonID);
180
181 // Determine true geographical position
182 Vector lat(1), lon(1);
184 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
185
186 // Ensure correct coverage of lon grid
187 Vector lon_shifted;
188 lon_shiftgrid(lon_shifted, GFlon, lon[0]);
189
190 // Check if lat/lon we need are actually covered
191 chk_if_in_range("rtp_pos.lat", lat[0], GFlat[0], GFlat[nlat - 1]);
192 chk_if_in_range("rtp_pos.lon", lon[0], lon_shifted[0], lon_shifted[nlon - 1]);
193
194 // Interpolate in lat and lon
195 //
196 GridPos gp_lat, gp_lon;
197 gridpos(gp_lat, GFlat, lat[0]);
198 gridpos(gp_lon, lon_shifted, lon[0]);
199 Vector itw(4);
200 interpweights(itw, gp_lat, gp_lon);
201 outvalue = interp(itw, gfield2.data, gp_lat, gp_lon);
202}
203
204
205/* Workspace method: Doxygen documentation will be auto-generated */
207 const Index& atmosphere_dim,
208 const Vector& lat_grid,
209 const Vector& lon_grid,
210 const Vector& rtp_pos,
211 const Matrix& z_surface,
212 const Matrix& field,
213 const Verbosity& verbosity) {
214 // Input checks (dummy p_grid)
215 chk_atm_grids(atmosphere_dim, Vector(2, 2, -1), lat_grid, lon_grid);
217 "input argument *field*", field, atmosphere_dim, lat_grid, lon_grid);
218 chk_rte_pos(atmosphere_dim, rtp_pos);
219 //
220 const Numeric zmax = max(z_surface);
221 const Numeric zmin = min(z_surface);
222 const Numeric dzok = 1;
223 if (rtp_pos[0] < zmin - dzok || rtp_pos[0] > zmax + dzok) {
224 ostringstream os;
225 os << "The given position does not match *z_surface*.\nThe altitude in "
226 << "*rtp_pos* is " << rtp_pos[0] / 1e3 << " km.\n"
227 << "The altitude range covered by *z_surface* is [" << zmin / 1e3 << ","
228 << zmax / 1e3 << "] km.\n"
229 << "One possible mistake is to mix up *rtp_pos* and *rte_los*.";
230 throw runtime_error(os.str());
231 }
232
233 if (atmosphere_dim == 1) {
234 outvalue = field(0, 0);
235 } else {
236 chk_interpolation_grids("Latitude interpolation", lat_grid, rtp_pos[1]);
237 GridPos gp_lat, gp_lon;
238 gridpos(gp_lat, lat_grid, rtp_pos[1]);
239 if (atmosphere_dim == 3) {
240 chk_interpolation_grids("Longitude interpolation", lon_grid, rtp_pos[2]);
241 gridpos(gp_lon, lon_grid, rtp_pos[2]);
242 }
243 //
244 outvalue = interp_atmsurface_by_gp(atmosphere_dim, field, gp_lat, gp_lon);
245 }
246
247 // Interpolate
249 out3 << " Result = " << outvalue << "\n";
250}
251
252/* Workspace method: Doxygen documentation will be auto-generated */
254 Matrix& iy,
255 ArrayOfTensor3& diy_dx,
256 const String& iy_unit,
257 const Tensor3& iy_transmittance,
258 const Index& iy_id,
259 const Index& cloudbox_on,
260 const Index& jacobian_do,
261 const Vector& f_grid,
262 const Agenda& iy_main_agenda,
263 const Vector& rtp_pos,
264 const Vector& rtp_los,
265 const Vector& rte_pos2,
266 const ArrayOfAgenda& iy_surface_agenda_array,
267 const ArrayOfIndex& surface_types,
268 const Vector& surface_types_aux,
269 const Vector& surface_types_weights,
270 const Verbosity&) {
271
272 const Index ntypes = surface_types.nelem();
273 if (ntypes < 1)
274 throw runtime_error("*surface_types* is empty!");
275
276 // Loop surface types and sum up
277 //
278 Numeric wtot = 0;
279 //
280 for (Index t=0; t<ntypes; t++ ) {
281
282 if (surface_types[t] < 0)
283 throw runtime_error(
284 "No element in *surface_types* is allowed to be negative.");
285 if (surface_types[t] >= iy_surface_agenda_array.nelem()) {
286 ostringstream os;
287 os << "*iy_surface_agenda_array* has only "
288 << iy_surface_agenda_array.nelem()
289 << " elements,\n while you have selected element " << surface_types[t];
290 throw runtime_error(os.str());
291 }
292
293 Matrix iy1;
294 ArrayOfTensor3 diy_dx1;
295
297 iy1,
298 diy_dx1,
299 surface_types[t],
300 iy_unit,
301 iy_transmittance,
302 iy_id,
303 cloudbox_on,
304 jacobian_do,
305 iy_main_agenda,
306 f_grid,
307 rtp_pos,
308 rtp_los,
309 rte_pos2,
310 surface_types_aux[t],
311 iy_surface_agenda_array);
312
313 iy1 *= surface_types_weights[t];
314 for (Index i=0; diy_dx1.nelem(); i++ )
315 diy_dx1[i] *= surface_types_weights[t];
316 wtot += surface_types_weights[t];;
317
318 if (t==0) {
319 iy = iy1;
320 diy_dx = diy_dx1;
321 } else {
322 iy += iy1;
323 for (Index i=0; diy_dx1.nelem(); i++ )
324 diy_dx[i] += diy_dx1[i];
325 }
326 }
327 if (abs(wtot-1)>1e-4)
328 throw runtime_error("Sum of *surface_types_weights* deviates from 1.");
329}
330
331/* Workspace method: Doxygen documentation will be auto-generated */
333 Matrix& iy,
334 ArrayOfTensor3& diy_dx,
335 const Tensor3& iy_transmittance,
336 const Index& iy_id,
337 const Index& jacobian_do,
338 const Index& atmosphere_dim,
339 const EnergyLevelMap& nlte_field,
340 const Index& cloudbox_on,
341 const Index& stokes_dim,
342 const Vector& f_grid,
343 const Vector& rtp_pos,
344 const Vector& rtp_los,
345 const Vector& rte_pos2,
346 const String& iy_unit,
347 const Agenda& iy_main_agenda,
348 const Numeric& surface_skin_t,
349 const Numeric& salinity,
350 const Numeric& wind_speed,
351 const Numeric& wind_direction,
352 const Index& fastem_version,
353 const Verbosity& verbosity) {
354 // Most obvious input checks are performed in specular_losCalc and surfaceFastem
355
356 // Obtain radiance and transmission for specular direction
357
358 // Determine specular direction
359 Vector specular_los, surface_normal;
360 specular_losCalcNoTopography(specular_los,
361 surface_normal,
362 rtp_pos,
363 rtp_los,
364 atmosphere_dim,
365 verbosity);
366
367 // Use iy_aux to get optical depth for downwelling radiation.
368 ArrayOfString iy_aux_vars(1);
369 iy_aux_vars[0] = "Optical depth";
370
371 // Calculate iy for downwelling radiation
372 // Note that iy_transmittance used here lacks surface R. Fixed below.
373 //
374 const Index nf = f_grid.nelem();
375 Vector transmittance(nf);
376 Vector geo_pos;
377 ArrayOfMatrix iy_aux;
378 Ppath ppath;
379 //
381 iy,
382 iy_aux,
383 ppath,
384 diy_dx,
385 geo_pos,
386 0,
387 iy_transmittance,
388 iy_aux_vars,
389 iy_id,
390 iy_unit,
391 cloudbox_on,
392 jacobian_do,
393 f_grid,
394 nlte_field,
395 rtp_pos,
396 specular_los,
397 rte_pos2,
398 iy_main_agenda);
399
400 // Convert tau to transmissions
401 for (Index i = 0; i < nf; i++) {
402 transmittance[i] = exp(-iy_aux[0](i, 0));
403 }
404
405 // Call Fastem by surface_RTprop version
406 //
407 Matrix surface_los;
408 Tensor4 surface_rmatrix;
409 Matrix surface_emission;
410 //
411 surfaceFastem(surface_los,
412 surface_rmatrix,
413 surface_emission,
414 atmosphere_dim,
415 stokes_dim,
416 f_grid,
417 rtp_pos,
418 rtp_los,
419 surface_skin_t,
420 salinity,
421 wind_speed,
422 wind_direction,
423 transmittance,
424 fastem_version,
425 verbosity);
426
427 // Add up
428 //
429 Tensor3 I(1, nf, stokes_dim);
430 I(0, joker, joker) = iy;
431 Matrix sensor_los_dummy(1, 1, 0);
432 //
433 surface_calc(iy, I, sensor_los_dummy, surface_rmatrix, surface_emission);
434
435 // Adjust diy_dx, if necessary.
436 // For vector cases this is a slight approximation, as the order of the
437 // different transmission and reflectivities matters.
438 if (iy_transmittance.npages()) {
439 for (Index q = 0; q < diy_dx.nelem(); q++) {
440 for (Index p = 0; p < diy_dx[q].npages(); p++) {
441 for (Index i = 0; i < nf; i++) {
442 Vector x = diy_dx[q](p, i, joker);
443 mult(diy_dx[q](p, i, joker), surface_rmatrix(0, i, joker, joker), x);
444 }
445 }
446 }
447 }
448}
449
450/* Workspace method: Doxygen documentation will be auto-generated */
452 Matrix& iy,
453 ArrayOfTensor3& diy_dx,
454 ArrayOfTensor4& dsurface_rmatrix_dx,
455 ArrayOfMatrix& dsurface_emission_dx,
456 const Tensor3& iy_transmittance,
457 const Index& iy_id,
458 const Index& jacobian_do,
459 const Index& stars_do,
460 const Index& atmosphere_dim,
461 const EnergyLevelMap& nlte_field,
462 const Index& cloudbox_on,
463 const Index& stokes_dim,
464 const Vector& f_grid,
465 const Vector& lat_grid,
466 const Vector& lon_grid,
467 const Matrix& z_surface,
468 const Vector& refellipsoid,
469 const Vector& rtp_pos,
470 const Vector& rtp_los,
471 const Vector& rte_pos2,
472 const String& iy_unit,
473 const Tensor3& surface_reflectivity,
474 const Tensor3& surface_props_data,
475 const ArrayOfString& surface_props_names,
476 const ArrayOfString& dsurface_names,
477 const ArrayOfRetrievalQuantity& jacobian_quantities,
478 const Agenda& iy_main_agenda,
479 const Verbosity& verbosity) {
480
481 // Input checks
482 ARTS_USER_ERROR_IF(atmosphere_dim==2, "This method does not work for 2d atmospheres.")
483 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
484 chk_rte_pos(atmosphere_dim, rtp_pos);
485 chk_rte_los(atmosphere_dim, rtp_los);
486 chk_size("iy",iy,f_grid.nelem(),stokes_dim);
487
488 // Check surface_data
489 surface_props_check(atmosphere_dim,
490 lat_grid,
491 lon_grid,
492 surface_props_data,
493 surface_props_names);
494
495 // Interplation grid positions and weights
496 ArrayOfGridPos gp_lat(1), gp_lon(1);
497 Matrix itw;
499 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
500 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
501
502 // Skin temperature
503 Vector surface_skin_t(1);
504 surface_props_interp(surface_skin_t,
505 "Skin temperature",
506 atmosphere_dim,
507 gp_lat,
508 gp_lon,
509 itw,
510 surface_props_data,
511 surface_props_names);
512
513 chk_not_negative("surface_skin_t", surface_skin_t[0]);
514
515
516
517 //Allocate
518 Tensor3 iy_trans_new;
519 Matrix iy_incoming;
520 Vector specular_los, surface_normal;
521 ArrayOfMatrix iy_aux;
522 Ppath ppath;
523 ArrayOfTensor3 diy_dx_dumb;
524 Matrix surface_los;
525 Tensor4 surface_rmatrix;
526 Matrix surface_emission;
527 Tensor3 I(1,f_grid.nelem(),stokes_dim);
528 if (iy_transmittance.npages()) {
529 iy_trans_new=iy_transmittance;
530 }
531
532 //get specular line of sight
533 specular_losCalc(specular_los,
534 surface_normal,
535 rtp_pos,
536 rtp_los,
537 atmosphere_dim,
538 lat_grid,
539 lon_grid,
540 refellipsoid,
541 z_surface,
542 0,
543 verbosity);
544
545 // Calculate incoming radiation directions, surface reflection matrix and
546 // emission vector
547 surfaceFlatReflectivity(surface_los,
548 surface_rmatrix,
549 surface_emission,
550 f_grid,
551 stokes_dim,
552 atmosphere_dim,
553 rtp_pos,
554 rtp_los,
555 specular_los,
556 surface_skin_t[0],
557 surface_reflectivity,
558 verbosity);
559
560
561
562 // Jacobian part
563 if (jacobian_do) {
564 dsurface_check(surface_props_names,
565 dsurface_names,
566 dsurface_rmatrix_dx,
567 dsurface_emission_dx);
568
569 Index irq;
570
571 // Skin temperature
572 irq = find_first(dsurface_names, String("Skin temperature"));
573 if (irq >= 0) {
574 const Numeric dd = 0.1;
575 Matrix surface_los2;
576 surfaceFlatReflectivity(surface_los2,
577 dsurface_rmatrix_dx[irq],
578 dsurface_emission_dx[irq],
579 f_grid,
580 stokes_dim,
581 atmosphere_dim,
582 rtp_pos,
583 rtp_los,
584 specular_los,
585 surface_skin_t[0]+dd,
586 surface_reflectivity,
587 verbosity);
588 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
589 dsurface_rmatrix_dx[irq] /= dd;
590 dsurface_emission_dx[irq] -= surface_emission;
591 dsurface_emission_dx[irq] /= dd;
592 }
593 }
594
596 iy,
597 diy_dx,
598 surface_los,
599 surface_rmatrix,
600 surface_emission,
601 dsurface_names,
602 dsurface_rmatrix_dx,
603 dsurface_emission_dx,
604 iy_transmittance,
605 iy_id,
606 jacobian_do,
607 stars_do,
608 jacobian_quantities,
609 atmosphere_dim,
610 nlte_field,
611 cloudbox_on,
612 stokes_dim,
613 f_grid,
614 rtp_pos,
615 rtp_los,
616 rte_pos2,
617 iy_unit,
618 iy_main_agenda,
619 verbosity);
620
621
622}
623
624/* Workspace method: Doxygen documentation will be auto-generated */
626 Workspace& ws,
627 Matrix& iy,
628 const Vector& rtp_pos,
629 const Vector& rtp_los,
630 const Index& stokes_dim,
631 const Vector& f_grid,
632 const Index& atmosphere_dim,
633 const Vector& p_grid,
634 const Vector& lat_grid,
635 const Vector& lon_grid,
636 const Tensor3& z_field,
637 const Tensor3& t_field,
638 const EnergyLevelMap& nlte_field,
639 const Tensor4& vmr_field,
640 const ArrayOfArrayOfSpeciesTag& abs_species,
641 const Tensor3& wind_u_field,
642 const Tensor3& wind_v_field,
643 const Tensor3& wind_w_field,
644 const Tensor3& mag_u_field,
645 const Tensor3& mag_v_field,
646 const Tensor3& mag_w_field,
647 const Matrix& z_surface,
648 const Tensor3& surface_reflectivity,
649 const Vector& refellipsoid,
650 const Tensor4& pnd_field,
651 const ArrayOfTensor4& dpnd_field_dx,
652 const ArrayOfString& scat_species,
653 const ArrayOfArrayOfSingleScatteringData& scat_data,
654 const Numeric& ppath_lmax,
655 const Numeric& ppath_lraytrace,
656 const Index& ppath_inside_cloudbox_do,
657 const Index& cloudbox_on,
658 const ArrayOfIndex& cloudbox_limits,
659 const Index& stars_do,
660 const Index& gas_scattering_do,
661 const Index& jacobian_do,
662 const ArrayOfRetrievalQuantity& jacobian_quantities,
663 const ArrayOfStar& stars,
664 const Numeric& rte_alonglos_v,
665 const String& iy_unit,
666 const Agenda& propmat_clearsky_agenda,
667 const Agenda& water_p_eq_agenda,
668 const Agenda& gas_scattering_agenda,
669 const Agenda& ppath_step_agenda,
670 const Verbosity& verbosity) {
671 //Check for correct unit
672 ARTS_USER_ERROR_IF(iy_unit != "1" && stars_do,
673 "If stars are present only iy_unit=\"1\" can be used.");
674
675 chk_size("iy", iy, f_grid.nelem(), stokes_dim);
676
677 if (stars_do) {
678 Matrix iy_incoming;
679 Index stars_visible;
680 Vector specular_los;
681
682 //Get incoming direct radiation, if no star is in line of sight, then
683 //iy_incoming is zero.
685 iy_incoming,
686 stars_visible,
687 specular_los,
688 rtp_pos,
689 rtp_los,
690 stokes_dim,
691 f_grid,
692 atmosphere_dim,
693 p_grid,
694 lat_grid,
695 lon_grid,
696 z_field,
697 t_field,
698 nlte_field,
699 vmr_field,
700 abs_species,
701 wind_u_field,
702 wind_v_field,
703 wind_w_field,
704 mag_u_field,
705 mag_v_field,
706 mag_w_field,
707 z_surface,
708 refellipsoid,
709 pnd_field,
710 dpnd_field_dx,
711 scat_species,
712 scat_data,
713 ppath_lmax,
714 ppath_lraytrace,
715 ppath_inside_cloudbox_do,
716 cloudbox_on,
717 cloudbox_limits,
718 gas_scattering_do,
719 jacobian_do,
720 jacobian_quantities,
721 stars,
722 rte_alonglos_v,
723 propmat_clearsky_agenda,
724 water_p_eq_agenda,
725 gas_scattering_agenda,
726 ppath_step_agenda,
727 verbosity);
728
729 if (stars_visible) {
730 Matrix surface_los;
731 Tensor4 surface_rmatrix;
732 Matrix surface_emission;
733 Numeric surface_skin_t_dummy = 1.;
734
735 //As we only consider the reflection of the direct radiation we do not consider
736 //any type of surface emission, therefore we set the surface skin temperature
737 //to a dummy value and later on set the surface emission vector (matrix) to zero.
738 surfaceFlatReflectivity(surface_los,
739 surface_rmatrix,
740 surface_emission,
741 f_grid,
742 stokes_dim,
743 atmosphere_dim,
744 rtp_pos,
745 rtp_los,
746 specular_los,
747 surface_skin_t_dummy,
748 surface_reflectivity,
749 verbosity);
750
751 surface_emission *= 0.;
752
753 Tensor3 I(1, f_grid.nelem(), stokes_dim);
754 I(0, joker, joker) = iy_incoming;
755
756 surface_calc(iy, I, surface_los, surface_rmatrix, surface_emission);
757
758 //TODO: add reflectivity jacobian
759 }
760 }
761}
762
763/* Workspace method: Doxygen documentation will be auto-generated */
765 Matrix& iy,
766 ArrayOfTensor3& diy_dx,
767 ArrayOfTensor4& dsurface_rmatrix_dx,
768 ArrayOfMatrix& dsurface_emission_dx,
769 const Tensor3& iy_transmittance,
770 const Index& iy_id,
771 const Index& jacobian_do,
772 const Index& stars_do,
773 const Index& atmosphere_dim,
774 const EnergyLevelMap& nlte_field,
775 const Index& cloudbox_on,
776 const Index& stokes_dim,
777 const Vector& f_grid,
778 const Vector& lat_grid,
779 const Vector& lon_grid,
780 const Matrix& z_surface,
781 const Vector& refellipsoid,
782 const Vector& rtp_pos,
783 const Vector& rtp_los,
784 const Vector& rte_pos2,
785 const String& iy_unit,
786 const GriddedField3& surface_complex_refr_index,
787 const Tensor3& surface_props_data,
788 const ArrayOfString& surface_props_names,
789 const ArrayOfString& dsurface_names,
790 const ArrayOfRetrievalQuantity& jacobian_quantities,
791 const Agenda& iy_main_agenda,
792 const Verbosity& verbosity) {
793
794 // Input checks
795 ARTS_USER_ERROR_IF(atmosphere_dim==2, "This method does not work for 2d atmospheres.")
796 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
797 chk_rte_pos(atmosphere_dim, rtp_pos);
798 chk_rte_los(atmosphere_dim, rtp_los);
799 chk_size("iy",iy,f_grid.nelem(),stokes_dim);
800
801 // Check surface_data
802 surface_props_check(atmosphere_dim,
803 lat_grid,
804 lon_grid,
805 surface_props_data,
806 surface_props_names);
807
808 // Interplation grid positions and weights
809 ArrayOfGridPos gp_lat(1), gp_lon(1);
810 Matrix itw;
812 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
813 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
814
815 // Skin temperature
816 Vector surface_skin_t(1);
817 surface_props_interp(surface_skin_t,
818 "Skin temperature",
819 atmosphere_dim,
820 gp_lat,
821 gp_lon,
822 itw,
823 surface_props_data,
824 surface_props_names);
825
826 chk_not_negative("surface_skin_t", surface_skin_t[0]);
827
828
829 //Allocate
830 Tensor3 iy_trans_new;
831 Matrix iy_incoming;
832 Vector specular_los, surface_normal;
833 ArrayOfMatrix iy_aux;
834 Ppath ppath;
835 ArrayOfTensor3 diy_dx_dumb;
836 Matrix surface_los;
837 Tensor4 surface_rmatrix;
838 Matrix surface_emission;
839 Tensor3 I(1,f_grid.nelem(),stokes_dim);
840 if (iy_transmittance.npages()) {
841 iy_trans_new=iy_transmittance;
842 }
843
844 //get specular line of sight
845 specular_losCalc(specular_los,
846 surface_normal,
847 rtp_pos,
848 rtp_los,
849 atmosphere_dim,
850 lat_grid,
851 lon_grid,
852 refellipsoid,
853 z_surface,
854 0,
855 verbosity);
856
857 // Calculate incoming radiation directions, surface reflection matrix and
858 // emission vector
859 surfaceFlatRefractiveIndex(surface_los,
860 surface_rmatrix,
861 surface_emission,
862 f_grid,
863 stokes_dim,
864 atmosphere_dim,
865 rtp_pos,
866 rtp_los,
867 specular_los,
868 surface_skin_t[0],
869 surface_complex_refr_index,
870 verbosity);
871
872
873
874 // Jacobian part
875 if (jacobian_do) {
876 dsurface_check(surface_props_names,
877 dsurface_names,
878 dsurface_rmatrix_dx,
879 dsurface_emission_dx);
880
881 Index irq;
882
883 // Skin temperature
884 irq = find_first(dsurface_names, String("Skin temperature"));
885 if (irq >= 0) {
886 const Numeric dd = 0.1;
887 Matrix surface_los2;
888 surfaceFlatRefractiveIndex(surface_los2,
889 dsurface_rmatrix_dx[irq],
890 dsurface_emission_dx[irq],
891 f_grid,
892 stokes_dim,
893 atmosphere_dim,
894 rtp_pos,
895 rtp_los,
896 specular_los,
897 surface_skin_t[0]+dd,
898 surface_complex_refr_index,
899 verbosity);
900 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
901 dsurface_rmatrix_dx[irq] /= dd;
902 dsurface_emission_dx[irq] -= surface_emission;
903 dsurface_emission_dx[irq] /= dd;
904 }
905 }
906
908 iy,
909 diy_dx,
910 surface_los,
911 surface_rmatrix,
912 surface_emission,
913 dsurface_names,
914 dsurface_rmatrix_dx,
915 dsurface_emission_dx,
916 iy_transmittance,
917 iy_id,
918 jacobian_do,
919 stars_do,
920 jacobian_quantities,
921 atmosphere_dim,
922 nlte_field,
923 cloudbox_on,
924 stokes_dim,
925 f_grid,
926 rtp_pos,
927 rtp_los,
928 rte_pos2,
929 iy_unit,
930 iy_main_agenda,
931 verbosity);
932
933
934}
935
936/* Workspace method: Doxygen documentation will be auto-generated */
938 Workspace& ws,
939 Matrix& iy,
940 const Vector& rtp_pos,
941 const Vector& rtp_los,
942 const Index& stokes_dim,
943 const Vector& f_grid,
944 const Index& atmosphere_dim,
945 const Vector& p_grid,
946 const Vector& lat_grid,
947 const Vector& lon_grid,
948 const Tensor3& z_field,
949 const Tensor3& t_field,
950 const EnergyLevelMap& nlte_field,
951 const Tensor4& vmr_field,
952 const ArrayOfArrayOfSpeciesTag& abs_species,
953 const Tensor3& wind_u_field,
954 const Tensor3& wind_v_field,
955 const Tensor3& wind_w_field,
956 const Tensor3& mag_u_field,
957 const Tensor3& mag_v_field,
958 const Tensor3& mag_w_field,
959 const Matrix& z_surface,
960 const GriddedField3& surface_complex_refr_index,
961 const Vector& refellipsoid,
962 const Tensor4& pnd_field,
963 const ArrayOfTensor4& dpnd_field_dx,
964 const ArrayOfString& scat_species,
965 const ArrayOfArrayOfSingleScatteringData& scat_data,
966 const Numeric& ppath_lmax,
967 const Numeric& ppath_lraytrace,
968 const Index& ppath_inside_cloudbox_do,
969 const Index& cloudbox_on,
970 const ArrayOfIndex& cloudbox_limits,
971 const Index& stars_do,
972 const Index& gas_scattering_do,
973 const Index& jacobian_do,
974 const ArrayOfRetrievalQuantity& jacobian_quantities,
975 const ArrayOfStar& stars,
976 const Numeric& rte_alonglos_v,
977 const String& iy_unit,
978 const Agenda& propmat_clearsky_agenda,
979 const Agenda& water_p_eq_agenda,
980 const Agenda& gas_scattering_agenda,
981 const Agenda& ppath_step_agenda,
982 const Verbosity& verbosity) {
983 //Check for correct unit
984 ARTS_USER_ERROR_IF(iy_unit != "1" && stars_do,
985 "If stars are present only iy_unit=\"1\" can be used.");
986
987 chk_size("iy", iy, f_grid.nelem(), stokes_dim);
988
989 if (stars_do) {
990 Matrix iy_incoming;
991 Index stars_visible;
992 Vector specular_los;
993
994 //Get incoming direct radiation, if no star is in line of sight, then
995 //iy_incoming is zero.
997 iy_incoming,
998 stars_visible,
999 specular_los,
1000 rtp_pos,
1001 rtp_los,
1002 stokes_dim,
1003 f_grid,
1004 atmosphere_dim,
1005 p_grid,
1006 lat_grid,
1007 lon_grid,
1008 z_field,
1009 t_field,
1010 nlte_field,
1011 vmr_field,
1012 abs_species,
1013 wind_u_field,
1014 wind_v_field,
1015 wind_w_field,
1016 mag_u_field,
1017 mag_v_field,
1018 mag_w_field,
1019 z_surface,
1020 refellipsoid,
1021 pnd_field,
1022 dpnd_field_dx,
1023 scat_species,
1024 scat_data,
1025 ppath_lmax,
1026 ppath_lraytrace,
1027 ppath_inside_cloudbox_do,
1028 cloudbox_on,
1029 cloudbox_limits,
1030 gas_scattering_do,
1031 jacobian_do,
1032 jacobian_quantities,
1033 stars,
1034 rte_alonglos_v,
1035 propmat_clearsky_agenda,
1036 water_p_eq_agenda,
1037 gas_scattering_agenda,
1038 ppath_step_agenda,
1039 verbosity);
1040
1041 if (stars_visible) {
1042 Matrix surface_los;
1043 Tensor4 surface_rmatrix;
1044 Matrix surface_emission;
1045 Numeric surface_skin_t_dummy = 1.;
1046
1047 //As we only consider the reflection of the direct radiation we do not consider
1048 //any type of surface emission, therefore we set the surface skin temperature
1049 //to a dummy value and later on set the surface emission vector (matrix) to zero.
1050 surfaceFlatRefractiveIndex(surface_los,
1051 surface_rmatrix,
1052 surface_emission,
1053 f_grid,
1054 stokes_dim,
1055 atmosphere_dim,
1056 rtp_pos,
1057 rtp_los,
1058 specular_los,
1059 surface_skin_t_dummy,
1060 surface_complex_refr_index,
1061 verbosity);
1062
1063 surface_emission *= 0.;
1064
1065 Tensor3 I(1, f_grid.nelem(), stokes_dim);
1066 I(0, joker, joker) = iy_incoming;
1067
1068 surface_calc(iy, I, surface_los, surface_rmatrix, surface_emission);
1069
1070 //TODO: add reflectivity jacobian
1071 }
1072 }
1073}
1074
1075/* Workspace method: Doxygen documentation will be auto-generated */
1077 const Vector& f_grid,
1078 const Index& stokes_dim,
1079 const Verbosity&) {
1080 iy.resize(f_grid.nelem(), stokes_dim);
1081 iy = 0.;
1082}
1083
1084/* Workspace method: Doxygen documentation will be auto-generated */
1086 Matrix& iy,
1087 ArrayOfTensor3& diy_dx,
1088 const Tensor3& iy_transmittance,
1089 const Index& iy_id,
1090 const Index& jacobian_do,
1091 const Index& stars_do,
1092 const Index& atmosphere_dim,
1093 const EnergyLevelMap& nlte_field,
1094 const Index& cloudbox_on,
1095 const Index& stokes_dim,
1096 const Vector& f_grid,
1097 const Vector& lat_grid,
1098 const Vector& lon_grid,
1099 const Matrix& z_surface,
1100 const Vector& refellipsoid,
1101 const Vector& rtp_pos,
1102 const Vector& rtp_los,
1103 const Vector& rte_pos2,
1104 const String& iy_unit,
1105 const Vector& surface_scalar_reflectivity,
1106 const Tensor3& surface_props_data,
1107 const ArrayOfString& surface_props_names,
1108 const ArrayOfString& dsurface_names,
1109 const ArrayOfRetrievalQuantity& jacobian_quantities,
1110 const Agenda& iy_main_agenda,
1111 const Index& N_za,
1112 const Index& N_aa,
1113 const Verbosity& verbosity) {
1114 // Input checks
1115 ARTS_USER_ERROR_IF(surface_scalar_reflectivity.nelem() != f_grid.nelem() &&
1116 surface_scalar_reflectivity.nelem() != 1,
1117 "The number of elements in *surface_scalar_reflectivity* should\n",
1118 "match length of *f_grid* or be 1.",
1119 "\n length of *f_grid* : ", f_grid.nelem(),
1120 "\n length of *surface_scalar_reflectivity* : ",
1121 surface_scalar_reflectivity.nelem());
1122 ARTS_USER_ERROR_IF(min(surface_scalar_reflectivity) < 0 ||
1123 max(surface_scalar_reflectivity) > 1,
1124 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
1125 ARTS_USER_ERROR_IF(atmosphere_dim==2, "This method does not work for 2d atmospheres.");
1126 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1127 chk_rte_pos(atmosphere_dim, rtp_pos);
1128 chk_rte_los(atmosphere_dim, rtp_los);
1129 chk_size("iy",iy,f_grid.nelem(),stokes_dim);
1130
1131 // Check surface_data
1132 surface_props_check(atmosphere_dim,
1133 lat_grid,
1134 lon_grid,
1135 surface_props_data,
1136 surface_props_names);
1137
1138 // Interplation grid positions and weights
1139 ArrayOfGridPos gp_lat(1), gp_lon(1);
1140 Matrix itw;
1142 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
1143 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
1144
1145 // Skin temperature
1146 Vector surface_skin_t(1);
1147 surface_props_interp(surface_skin_t,
1148 "Skin temperature",
1149 atmosphere_dim,
1150 gp_lat,
1151 gp_lon,
1152 itw,
1153 surface_props_data,
1154 surface_props_names);
1155
1156 chk_not_negative("surface_skin_t", surface_skin_t[0]);
1157
1158
1159
1160
1161 const Index nf = f_grid.nelem();
1162
1163 Vector za_grid, aa_grid, za_grid_weights;
1164
1165 //get angular grids. here we use gauss quadrature. As the function is not intended for
1166 //a hemisphere but for full sphere, we have to tweak it a little bit and use only
1167 //the first half of the angles.
1169 aa_grid,
1170 za_grid_weights,
1171 N_za * 2,
1172 N_aa,
1173 "double_gauss",
1174 verbosity);
1175
1176 Tensor3 iy_trans_new;
1177
1178
1179 ArrayOfString iy_aux_var(0);
1180 if (stars_do) {
1181 iy_aux_var.resize(1);
1182 iy_aux_var[0] = "Direct radiation";
1183 }
1184
1185 Matrix iy_temp;
1186 Matrix Flx(nf, 1, 0.);
1187
1188 if (iy_transmittance.npages()) {
1189 iy_trans_new=iy_transmittance;
1190 }
1191
1192
1193 ArrayOfTensor3 diy_dx_dumb;
1194 iy.resize(nf, stokes_dim);
1195 iy = 0;
1196
1197 Vector los;
1198
1199 //No azimuth dependency
1200 if (N_aa == 1) {
1201 //Set azimuth angle index
1202 Index i_aa = 0;
1203
1204 for (Index i_za = 0; i_za < N_za; i_za++) {
1205
1206 if (atmosphere_dim==1){
1207 los.resize(1);
1208 }
1209 else{
1210 los.resize(2);
1211 // For clearsky, the skylight is in general approximately independent of the azimuth direction.
1212 // So, we consider only the line of sight plane.
1213 los[1] = 0.;
1214 }
1215 los[0] = za_grid[i_za];
1216
1217
1218 // Calculate downwelling radiation for a los
1219 ArrayOfMatrix iy_aux;
1220 Ppath ppath;
1221 Vector geo_pos_dumb;
1222 Index iy_id_new = iy_id + i_za + i_aa * N_za + 1;
1224 iy_temp,
1225 iy_aux,
1226 ppath,
1227 diy_dx_dumb,
1228 geo_pos_dumb,
1229 0,
1230 iy_trans_new,
1231 iy_aux_var,
1232 iy_id_new,
1233 iy_unit,
1234 cloudbox_on,
1235 jacobian_do,
1236 f_grid,
1237 nlte_field,
1238 rtp_pos,
1239 los,
1240 rte_pos2,
1241 iy_main_agenda);
1242
1243 //For the case that a star is present and the los is towards a star, we
1244 //subtract the direct radiation, so that only the diffuse radiation is considered here.
1245 //If star is within los iy_aux[0] is the incoming and attenuated direct (star)
1246 //radiation at the surface. Otherwise, it is zero.
1247 if (stars_do) {
1248 iy_temp -= iy_aux[0];
1249 }
1250 iy_temp *= abs(cos(deg2rad(los[0])) * sin(deg2rad(los[0]))) *
1251 za_grid_weights[i_za];
1252 Flx(joker, 0) += iy_temp(joker, 0);
1253 }
1254 //Azimuthal integration just gives a factor of 2*pi
1255 Flx *= 2 * pi;
1256
1257 //Azimuth dependency
1258 } else {
1259 // before we start with the actual RT calculation, we calculate the integration
1260 // weights for a trapezoidal integration over the azimuth.
1261 Vector aa_grid_weights;
1262 Vector aa_grid_rad = aa_grid;
1263
1264 //For the integration we need the aa grid in rad
1265 for (Index i = 0; i < aa_grid_rad.nelem(); i++) {
1266 aa_grid_rad[i] = deg2rad(aa_grid_rad[i]);
1267 }
1268
1269 calculate_int_weights_arbitrary_grid(aa_grid_weights, aa_grid_rad);
1270
1271 //Allocate
1272 los.resize(2);
1273
1274 for (Index i_aa = 0; i_aa < N_aa; i_aa++) {
1275 for (Index i_za = 0; i_za < N_za; i_za++) {
1276
1277
1278 los[0] = za_grid[i_za];
1279 los[1] = aa_grid[i_aa];
1280
1281 // Calculate downwelling radiation for LOS ilos
1282 ArrayOfMatrix iy_aux;
1283 Ppath ppath;
1284 Vector geo_pos_dumb;
1285 Index iy_id_new = iy_id + i_za + i_aa * N_za + 1;
1287 iy_temp,
1288 iy_aux,
1289 ppath,
1290 diy_dx_dumb,
1291 geo_pos_dumb,
1292 0,
1293 iy_trans_new,
1294 iy_aux_var,
1295 iy_id_new,
1296 iy_unit,
1297 cloudbox_on,
1298 jacobian_do,
1299 f_grid,
1300 nlte_field,
1301 rtp_pos,
1302 los,
1303 rte_pos2,
1304 iy_main_agenda);
1305
1306 //For the case that a star is present and the los is towards a star, we
1307 //subtract the direct radiation, so that only the diffuse radiation is considered here.
1308 //If star is within los iy_aux[0] is the incoming and attenuated direct (star)
1309 //radiation at the surface. Otherwise it is zero.
1310 if (stars_do) {
1311 iy_temp -= iy_aux[0];
1312 }
1313 iy_temp *= abs(cos(deg2rad(los[0])) * sin(deg2rad(los[0]))) *
1314 za_grid_weights[i_za] * aa_grid_weights[i_aa];
1315 Flx(joker, 0) += iy_temp(joker, 0);
1316 }
1317 }
1318 }
1319
1320 Vector specular_los, surface_normal;
1321 specular_losCalc(specular_los,
1322 surface_normal,
1323 rtp_pos,
1324 rtp_los,
1325 atmosphere_dim,
1326 lat_grid,
1327 lon_grid,
1328 refellipsoid,
1329 z_surface,
1330 0,
1331 verbosity);
1332
1333 //Surface emission
1334 Vector b(nf);
1335 planck(b, f_grid, surface_skin_t[0]);
1336
1337 //Upwelling radiation
1338 Numeric r;
1339 for (Index i_freq = 0; i_freq < f_grid.nelem(); i_freq++) {
1340 if (i_freq == 0 || surface_scalar_reflectivity.nelem() > 1) {
1341 r = surface_scalar_reflectivity[i_freq];
1342 }
1343 iy(i_freq, 0) = r *
1344 cos(deg2rad(specular_los[0])) / pi * Flx(i_freq, 0) +
1345 (1 - r) * b[i_freq];
1346 }
1347
1348 // Surface Jacobians
1349 if (jacobian_do && dsurface_names.nelem()) {
1350
1351 Index irq;
1352 irq = find_first(dsurface_names, String("Skin temperature"));
1353
1354 if (irq >= 0){
1355 // Find index among jacobian quantities
1356 Index ihit = -1;
1357 for (Index j = 0; j < jacobian_quantities.nelem() && ihit < 0; j++) {
1358 if (dsurface_names[irq] == jacobian_quantities[j].Subtag()) {
1359 ihit = j;
1360 }
1361 }
1362
1363 // Derivative of surface skin temperature, as observed at the surface
1364 Matrix diy_dpos0(f_grid.nelem(), stokes_dim, 0.), diy_dpos;
1365 dplanck_dt(diy_dpos0(joker, 0), f_grid, surface_skin_t[0]);
1366
1367 Vector emissivity=surface_scalar_reflectivity;
1368 emissivity*=-1;
1369 emissivity+=1;
1370 diy_dpos0*=emissivity;
1371
1372 // Weight with transmission to sensor
1373 iy_transmittance_mult(diy_dpos, iy_transmittance, diy_dpos0);
1374
1375 // Put into diy_dx
1376 diy_from_pos_to_rgrids(diy_dx[ihit],
1377 jacobian_quantities[ihit],
1378 diy_dpos,
1379 atmosphere_dim,
1380 rtp_pos);
1381
1382 }
1383 }
1384}
1385
1386/* Workspace method: Doxygen documentation will be auto-generated */
1388 Workspace& ws,
1389 Matrix& iy,
1390 const Vector& rtp_pos,
1391 const Index& stokes_dim,
1392 const Vector& f_grid,
1393 const Index& atmosphere_dim,
1394 const Vector& p_grid,
1395 const Vector& lat_grid,
1396 const Vector& lon_grid,
1397 const Tensor3& z_field,
1398 const Tensor3& t_field,
1399 const EnergyLevelMap& nlte_field,
1400 const Tensor4& vmr_field,
1401 const ArrayOfArrayOfSpeciesTag& abs_species,
1402 const Tensor3& wind_u_field,
1403 const Tensor3& wind_v_field,
1404 const Tensor3& wind_w_field,
1405 const Tensor3& mag_u_field,
1406 const Tensor3& mag_v_field,
1407 const Tensor3& mag_w_field,
1408 const Matrix& z_surface,
1409 const Vector& surface_scalar_reflectivity,
1410 const Vector& refellipsoid,
1411 const Tensor4& pnd_field,
1412 const ArrayOfTensor4& dpnd_field_dx,
1413 const ArrayOfString& scat_species,
1414 const ArrayOfArrayOfSingleScatteringData& scat_data,
1415 const Numeric& ppath_lmax,
1416 const Numeric& ppath_lraytrace,
1417 const Index& cloudbox_on,
1418 const ArrayOfIndex& cloudbox_limits,
1419 const Index& stars_do,
1420 const Index& gas_scattering_do,
1421 const Index& jacobian_do,
1422 const ArrayOfRetrievalQuantity& jacobian_quantities,
1423 const ArrayOfStar& stars,
1424 const Numeric& rte_alonglos_v,
1425 const String& iy_unit,
1426 const Agenda& propmat_clearsky_agenda,
1427 const Agenda& water_p_eq_agenda,
1428 const Agenda& gas_scattering_agenda,
1429 const Agenda& ppath_step_agenda,
1430 const Verbosity& verbosity) {
1431 //Allocate
1432 ArrayOfVector star_rte_los(stars.nelem());
1433 ArrayOfMatrix transmitted_starlight;
1434 ArrayOfArrayOfTensor3 dtransmitted_starlight;
1435 ArrayOfPpath star_ppaths(stars.nelem());
1436 ArrayOfIndex stars_visible(stars.nelem());
1437
1438 //Check for correct unit
1439 ARTS_USER_ERROR_IF(iy_unit != "1" && stars_do,
1440 "If stars are present only iy_unit=\"1\" can be used.");
1441 //Check size of iy
1442 chk_size("iy",iy,f_grid.nelem(),stokes_dim);
1443
1444 ARTS_USER_ERROR_IF(surface_scalar_reflectivity.nelem() != f_grid.nelem() &&
1445 surface_scalar_reflectivity.nelem() != 1,
1446 "The number of elements in *surface_scalar_reflectivity* should\n",
1447 "match length of *f_grid* or be 1.",
1448 "\n length of *f_grid* : ", f_grid.nelem(),
1449 "\n length of *surface_scalar_reflectivity* : ",
1450 surface_scalar_reflectivity.nelem());
1451 ARTS_USER_ERROR_IF(min(surface_scalar_reflectivity) < 0 ||
1452 max(surface_scalar_reflectivity) > 1,
1453 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
1454 // Input checks
1455 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1456 chk_rte_pos(atmosphere_dim, rtp_pos);
1457
1458
1459 //do something only if there is a star
1460 if (stars_do) {
1461 //get star ppaths
1462 get_star_ppaths(ws,
1463 star_ppaths,
1464 stars_visible,
1465 star_rte_los,
1466 rtp_pos,
1467 stars,
1468 f_grid,
1469 atmosphere_dim,
1470 p_grid,
1471 lat_grid,
1472 lon_grid,
1473 z_field,
1474 z_surface,
1475 refellipsoid,
1476 ppath_lmax,
1477 ppath_lraytrace,
1478 ppath_step_agenda,
1479 verbosity);
1480
1482 transmitted_starlight,
1483 dtransmitted_starlight,
1484 stokes_dim,
1485 f_grid,
1486 atmosphere_dim,
1487 p_grid,
1488 lat_grid,
1489 lon_grid,
1490 t_field,
1491 nlte_field,
1492 vmr_field,
1493 abs_species,
1494 wind_u_field,
1495 wind_v_field,
1496 wind_w_field,
1497 mag_u_field,
1498 mag_v_field,
1499 mag_w_field,
1500 cloudbox_on,
1501 cloudbox_limits,
1502 gas_scattering_do,
1503 1,
1504 star_ppaths,
1505 stars,
1506 stars_visible,
1507 refellipsoid,
1508 pnd_field,
1509 dpnd_field_dx,
1510 scat_species,
1511 scat_data,
1512 jacobian_do,
1513 jacobian_quantities,
1514 propmat_clearsky_agenda,
1515 water_p_eq_agenda,
1516 gas_scattering_agenda,
1517 rte_alonglos_v,
1518 verbosity);
1519
1520 //loop over stars
1521 Vector specular_los, surface_normal;
1522
1523 for (Index i_star = 0; i_star < stars.nelem(); i_star++) {
1524 if (stars_visible[i_star]) {
1525 //star_rte_los is the direction toward the star, but we need from the star
1526
1527 Vector incoming_los;
1528 mirror_los(incoming_los,star_rte_los[i_star], atmosphere_dim);
1529
1530 specular_losCalc(specular_los,
1531 surface_normal,
1532 rtp_pos,
1533 incoming_los,
1534 atmosphere_dim,
1535 lat_grid,
1536 lon_grid,
1537 refellipsoid,
1538 z_surface,
1539 0,
1540 verbosity);
1541
1542 // Only the first component of transmitted_starlight is relevant.
1543 // Comment taken from surfaceLambertianSimple.
1544 // This follows VDISORT that refers to: K. L. Coulson, Polarization and Intensity of Light in
1545 // the Atmosphere (1989), page 229
1546 // (Thanks to Michael Kahnert for providing this information!)
1547 // Update: Is the above for a later edition? We have not found a copy of
1548 // that edition. In a 1988 version of the book, the relevant page seems
1549 // to be 232.
1550
1551 Vector iy_surface_direct(f_grid.nelem());
1552 Numeric r;
1553 for (Index i_freq = 0; i_freq < f_grid.nelem(); i_freq++) {
1554 if (i_freq == 0 || surface_scalar_reflectivity.nelem() > 1) {
1555 r = surface_scalar_reflectivity[i_freq];
1556 }
1557 iy_surface_direct[i_freq] = r / pi *
1558 transmitted_starlight[i_star](i_freq, 0) *
1559 cos(deg2rad(specular_los[0]));
1560 }
1561
1562 //Add the surface contribution of the direct part
1563 iy(joker, 0) += iy_surface_direct;
1564
1565
1566 }
1567 }
1568 }
1569}
1570
1571/* Workspace method: Doxygen documentation will be auto-generated */
1573 Matrix& iy,
1574 ArrayOfTensor3& diy_dx,
1575 Numeric& surface_skin_t,
1576 Matrix& surface_los,
1577 Tensor4& surface_rmatrix,
1578 Matrix& surface_emission,
1579 const Tensor3& iy_transmittance,
1580 const Index& iy_id,
1581 const Index& jacobian_do,
1582 const Index& stars_do,
1583 const Index& atmosphere_dim,
1584 const EnergyLevelMap& nlte_field,
1585 const Index& cloudbox_on,
1586 const Index& stokes_dim,
1587 const Vector& f_grid,
1588 const Vector& rtp_pos,
1589 const Vector& rtp_los,
1590 const Vector& rte_pos2,
1591 const String& iy_unit,
1592 const Agenda& iy_main_agenda,
1593 const Agenda& surface_rtprop_agenda,
1594 const Verbosity&) {
1595 // Input checks
1596 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1597 chk_rte_pos(atmosphere_dim, rtp_pos);
1598 chk_rte_los(atmosphere_dim, rtp_los);
1599
1600 // Call *surface_rtprop_agenda*
1602 surface_skin_t,
1603 surface_emission,
1604 surface_los,
1605 surface_rmatrix,
1606 f_grid,
1607 rtp_pos,
1608 rtp_los,
1609 surface_rtprop_agenda);
1610
1611 // Check output of *surface_rtprop_agenda*
1612 const Index nlos = surface_los.nrows();
1613 const Index nf = f_grid.nelem();
1614 //
1615 if (nlos) // if 0, blackbody ground and not all checks are needed
1616 {
1617 if (surface_los.ncols() != rtp_los.nelem())
1618 throw runtime_error("Number of columns in *surface_los* is not correct.");
1619 if (nlos != surface_rmatrix.nbooks())
1620 throw runtime_error(
1621 "Mismatch in size of *surface_los* and *surface_rmatrix*.");
1622 if (surface_rmatrix.npages() != nf)
1623 throw runtime_error(
1624 "Mismatch in size of *surface_rmatrix* and *f_grid*.");
1625 if (surface_rmatrix.nrows() != stokes_dim ||
1626 surface_rmatrix.ncols() != stokes_dim)
1627 throw runtime_error(
1628 "Mismatch between size of *surface_rmatrix* and *stokes_dim*.");
1629 }
1630 if (surface_emission.ncols() != stokes_dim)
1631 throw runtime_error(
1632 "Mismatch between size of *surface_emission* and *stokes_dim*.");
1633 if (surface_emission.nrows() != nf)
1634 throw runtime_error("Mismatch in size of *surface_emission* and f_grid*.");
1635
1636 // Variable to hold down-welling radiation
1637 Tensor3 I(nlos, nf, stokes_dim);
1638
1639 ArrayOfString iy_aux_var(0);
1640 if (stars_do) iy_aux_var.emplace_back("Direct radiation");
1641
1642 // Loop *surface_los*-es. If no such LOS, we are ready.
1643 if (nlos > 0) {
1644 for (Index ilos = 0; ilos < nlos; ilos++) {
1645 Vector los = surface_los(ilos, joker);
1646
1647 // Include surface reflection matrix in *iy_transmittance*
1648 // If iy_transmittance is empty, this is interpreted as the
1649 // variable is not needed.
1650 //
1651 Tensor3 iy_trans_new;
1652 //
1653 if (iy_transmittance.npages()) {
1654 iy_transmittance_mult(iy_trans_new,
1655 iy_transmittance,
1656 surface_rmatrix(ilos, joker, joker, joker));
1657 }
1658
1659 // Calculate downwelling radiation for LOS ilos
1660 //
1661 {
1662 ArrayOfMatrix iy_aux;
1663 Ppath ppath;
1664 Vector geo_pos;
1665 Index iy_id_new = iy_id + ilos + 1;
1667 iy,
1668 iy_aux,
1669 ppath,
1670 diy_dx,
1671 geo_pos,
1672 0,
1673 iy_trans_new,
1674 iy_aux_var,
1675 iy_id_new,
1676 iy_unit,
1677 cloudbox_on,
1678 jacobian_do,
1679 f_grid,
1680 nlte_field,
1681 rtp_pos,
1682 los,
1683 rte_pos2,
1684 iy_main_agenda);
1685
1686 //For the case that a star is present and the los is towards a star, we
1687 //subtract the direct radiation, so that only the diffuse radiation is considered here.
1688 //If star is within los iy_aux[0] is the incoming and attenuated direct (star)
1689 //radiation at the surface. Otherwise it is zero.
1690 if (stars_do){
1691 iy-=iy_aux[0];
1692 }
1693
1694 }
1695
1696 if (iy.ncols() != stokes_dim || iy.nrows() != nf) {
1697 ostringstream os;
1698 os << "The size of *iy* returned from *" << iy_main_agenda.name()
1699 << "* is\n"
1700 << "not correct:\n"
1701 << " expected size = [" << nf << "," << stokes_dim << "]\n"
1702 << " size of iy = [" << iy.nrows() << "," << iy.ncols() << "]\n";
1703 throw runtime_error(os.str());
1704 }
1705
1706 I(ilos, joker, joker) = iy;
1707 }
1708 }
1709
1710 // Add up
1711 surface_calc(iy, I, surface_los, surface_rmatrix, surface_emission);
1712}
1713
1714/* Workspace method: Doxygen documentation will be auto-generated */
1716 Matrix& iy,
1717 ArrayOfTensor3& diy_dx,
1718 const Matrix& surface_los,
1719 const Tensor4& surface_rmatrix,
1720 const Matrix& surface_emission,
1721 const ArrayOfString& dsurface_names,
1722 const ArrayOfTensor4& dsurface_rmatrix_dx,
1723 const ArrayOfMatrix& dsurface_emission_dx,
1724 const Tensor3& iy_transmittance,
1725 const Index& iy_id,
1726 const Index& jacobian_do,
1727 const Index& stars_do,
1728 const ArrayOfRetrievalQuantity& jacobian_quantities,
1729 const Index& atmosphere_dim,
1730 const EnergyLevelMap& nlte_field,
1731 const Index& cloudbox_on,
1732 const Index& stokes_dim,
1733 const Vector& f_grid,
1734 const Vector& rtp_pos,
1735 const Vector& rtp_los,
1736 const Vector& rte_pos2,
1737 const String& iy_unit,
1738 const Agenda& iy_main_agenda,
1739 const Verbosity&) {
1740 // Input checks
1741 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1742 chk_rte_pos(atmosphere_dim, rtp_pos);
1743 chk_rte_los(atmosphere_dim, rtp_los);
1744
1745 // Check provided surface rtprop variables
1746 const Index nlos = surface_los.nrows();
1747 const Index nf = f_grid.nelem();
1748 //
1749 if (nlos) // if 0, blackbody ground and not all checks are needed
1750 {
1751 if (surface_los.ncols() != rtp_los.nelem())
1752 throw runtime_error("Number of columns in *surface_los* is not correct.");
1753 if (nlos != surface_rmatrix.nbooks())
1754 throw runtime_error(
1755 "Mismatch in size of *surface_los* and *surface_rmatrix*.");
1756 if (surface_rmatrix.npages() != nf)
1757 throw runtime_error(
1758 "Mismatch in size of *surface_rmatrix* and *f_grid*.");
1759 if (surface_rmatrix.nrows() != stokes_dim ||
1760 surface_rmatrix.ncols() != stokes_dim)
1761 throw runtime_error(
1762 "Mismatch between size of *surface_rmatrix* and *stokes_dim*.");
1763 }
1764 if (surface_emission.ncols() != stokes_dim)
1765 throw runtime_error(
1766 "Mismatch between size of *surface_emission* and *stokes_dim*.");
1767 if (surface_emission.nrows() != nf)
1768 throw runtime_error("Mismatch in size of *surface_emission* and f_grid*.");
1769
1770 // Variable to hold down-welling radiation
1771 Tensor3 I(nlos, nf, stokes_dim);
1772
1773 ArrayOfString iy_aux_var(0);
1774 if (stars_do) iy_aux_var.emplace_back("Direct radiation");
1775
1776 // Loop *surface_los*-es.
1777 if (nlos > 0) {
1778 for (Index ilos = 0; ilos < nlos; ilos++) {
1779 Vector los = surface_los(ilos, joker);
1780
1781 // Include surface reflection matrix in *iy_transmittance*
1782 // If iy_transmittance is empty, this is interpreted as the
1783 // variable is not needed.
1784 //
1785 Tensor3 iy_trans_new;
1786 //
1787 if (iy_transmittance.npages()) {
1788 iy_transmittance_mult(iy_trans_new,
1789 iy_transmittance,
1790 surface_rmatrix(ilos, joker, joker, joker));
1791 }
1792
1793 // Calculate downwelling radiation for LOS ilos
1794 //
1795 {
1796 ArrayOfMatrix iy_aux;
1797 Ppath ppath;
1798 Vector geo_pos;
1800 iy,
1801 iy_aux,
1802 ppath,
1803 diy_dx,
1804 geo_pos,
1805 0,
1806 iy_trans_new,
1807 iy_aux_var,
1808 iy_id,
1809 iy_unit,
1810 cloudbox_on,
1811 jacobian_do,
1812 f_grid,
1813 nlte_field,
1814 rtp_pos,
1815 los,
1816 rte_pos2,
1817 iy_main_agenda);
1818
1819 //For the case that a star is present and the los is towards a star, we
1820 //subtract the direct radiation, so that only the diffuse radiation is considered here.
1821 //If star is within los iy_aux[0] is the incoming and attenuated direct (star)
1822 //radiation at the surface. Otherwise it is zero.
1823 if (stars_do){
1824 iy-=iy_aux[0];
1825 }
1826
1827 }
1828
1829 if (iy.ncols() != stokes_dim || iy.nrows() != nf) {
1830 ostringstream os;
1831 os << "The size of *iy* returned from *" << iy_main_agenda.name()
1832 << "* is\n"
1833 << "not correct:\n"
1834 << " expected size = [" << nf << "," << stokes_dim << "]\n"
1835 << " size of iy = [" << iy.nrows() << "," << iy.ncols() << "]\n";
1836 throw runtime_error(os.str());
1837 }
1838
1839 I(ilos, joker, joker) = iy;
1840 }
1841 }
1842
1843 // Add up
1844 surface_calc(iy, I, surface_los, surface_rmatrix, surface_emission);
1845
1846 // Surface Jacobians
1847 if (jacobian_do && dsurface_names.nelem()) {
1848 // Loop dsurface_names
1849 for (Index i = 0; i < dsurface_names.nelem(); i++) {
1850 // Error if derivatives not calculated
1851 // Or should we accept this?
1852 if (dsurface_emission_dx[i].empty() || dsurface_rmatrix_dx[i].empty()) {
1853 ostringstream os;
1854 os << "The derivatives for surface quantity: " << dsurface_names[i]
1855 << "\nwere not calculated by *iy_surface_agenda*.\n"
1856 << "That is, *dsurface_emission_dx* and/or *dsurface_rmatrix_dx*\n"
1857 << "are empty.";
1858 throw runtime_error(os.str());
1859 } else {
1860 // Find index among jacobian quantities
1861 Index ihit = -1;
1862 for (Index j = 0; j < jacobian_quantities.nelem() && ihit < 0; j++) {
1863 if (dsurface_names[i] == jacobian_quantities[j].Subtag()) {
1864 ihit = j;
1865 }
1866 }
1867 ARTS_ASSERT(ihit >= 0);
1868 // Derivative, as observed at the surface
1869 Matrix diy_dpos0, diy_dpos;
1870 surface_calc(diy_dpos0,
1871 I,
1872 surface_los,
1873 dsurface_rmatrix_dx[i],
1874 dsurface_emission_dx[i]);
1875 // Weight with transmission to sensor
1876 iy_transmittance_mult(diy_dpos, iy_transmittance, diy_dpos0);
1877 // Put into diy_dx
1878 diy_from_pos_to_rgrids(diy_dx[ihit],
1879 jacobian_quantities[ihit],
1880 diy_dpos,
1881 atmosphere_dim,
1882 rtp_pos);
1883 }
1884 }
1885 }
1886}
1887
1888/* Workspace method: Doxygen documentation will be auto-generated */
1890 Vector& surface_normal,
1891 const Vector& rtp_pos,
1892 const Vector& rtp_los,
1893 const Index& atmosphere_dim,
1894 const Verbosity&) {
1895 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1896 chk_rte_pos(atmosphere_dim, rtp_pos);
1897 chk_rte_los(atmosphere_dim, rtp_los);
1898
1899 surface_normal.resize(max(Index(1), atmosphere_dim - 1));
1900 specular_los.resize(max(Index(1), atmosphere_dim - 1));
1901
1902 if (atmosphere_dim == 1) {
1903 surface_normal[0] = 0;
1904 if (rtp_los[0] < 90) {
1905 throw runtime_error(
1906 "Invalid zenith angle. The zenith angle corresponds "
1907 "to observe the surface from below.");
1908 }
1909 specular_los[0] = 180 - rtp_los[0];
1910 }
1911
1912 else if (atmosphere_dim == 2) {
1913 specular_los[0] = sign(rtp_los[0]) * 180 - rtp_los[0];
1914 surface_normal[0] = 0;
1915 }
1916
1917 else if (atmosphere_dim == 3) {
1918 specular_los[0] = 180 - rtp_los[0];
1919 specular_los[1] = rtp_los[1];
1920 surface_normal[0] = 0;
1921 surface_normal[1] = 0;
1922 }
1923}
1924
1925/* Workspace method: Doxygen documentation will be auto-generated */
1926void specular_losCalc(Vector& specular_los,
1927 Vector& surface_normal,
1928 const Vector& rtp_pos,
1929 const Vector& rtp_los,
1930 const Index& atmosphere_dim,
1931 const Vector& lat_grid,
1932 const Vector& lon_grid,
1933 const Vector& refellipsoid,
1934 const Matrix& z_surface,
1935 const Index& ignore_surface_slope,
1936 const Verbosity& verbosity) {
1937 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1938 chk_rte_pos(atmosphere_dim, rtp_pos);
1939 chk_rte_los(atmosphere_dim, rtp_los);
1940 chk_if_in_range("ignore_surface_slope", ignore_surface_slope, 0, 1);
1941
1942 // Use special function if there is no slope, or it is ignored
1943 if (atmosphere_dim == 1 || ignore_surface_slope) {
1944 specular_losCalcNoTopography(specular_los,
1945 surface_normal,
1946 rtp_pos,
1947 rtp_los,
1948 atmosphere_dim,
1949 verbosity);
1950 return;
1951 }
1952
1953 surface_normal.resize(max(Index(1), atmosphere_dim - 1));
1954 specular_los.resize(max(Index(1), atmosphere_dim - 1));
1955
1956 if (atmosphere_dim == 2) {
1957 chk_interpolation_grids("Latitude interpolation", lat_grid, rtp_pos[1]);
1958 GridPos gp_lat;
1959 gridpos(gp_lat, lat_grid, rtp_pos[1]);
1960 Numeric c1; // Radial slope of the surface
1962 c1, lat_grid, refellipsoid, z_surface(joker, 0), gp_lat, rtp_los[0]);
1963 Vector itw(2);
1964 interpweights(itw, gp_lat);
1965 const Numeric rv_surface = refell2d(refellipsoid, lat_grid, gp_lat) +
1966 interp(itw, z_surface(joker, 0), gp_lat);
1967 surface_normal[0] = -plevel_angletilt(rv_surface, c1);
1968 if (abs(rtp_los[0] - surface_normal[0]) < 90) {
1969 throw runtime_error(
1970 "Invalid zenith angle. The zenith angle corresponds "
1971 "to observe the surface from below.");
1972 }
1973 specular_los[0] =
1974 sign(rtp_los[0]) * 180 - rtp_los[0] + 2 * surface_normal[0];
1975 }
1976
1977 else {
1978 // Calculate surface normal in South-North direction
1979 chk_interpolation_grids("Latitude interpolation", lat_grid, rtp_pos[1]);
1980 chk_interpolation_grids("Longitude interpolation", lon_grid, rtp_pos[2]);
1981 GridPos gp_lat, gp_lon;
1982 gridpos(gp_lat, lat_grid, rtp_pos[1]);
1983 gridpos(gp_lon, lon_grid, rtp_pos[2]);
1984 Numeric c1, c2;
1986 c1, c2, lat_grid, lon_grid, refellipsoid, z_surface, gp_lat, gp_lon, 0);
1987 Vector itw(4);
1988 interpweights(itw, gp_lat, gp_lon);
1989 const Numeric rv_surface = refell2d(refellipsoid, lat_grid, gp_lat) +
1990 interp(itw, z_surface, gp_lat, gp_lon);
1991 const Numeric zaSN = 90 - plevel_angletilt(rv_surface, c1);
1992 // The same for East-West
1993 plevel_slope_3d(c1,
1994 c2,
1995 lat_grid,
1996 lon_grid,
1997 refellipsoid,
1998 z_surface,
1999 gp_lat,
2000 gp_lon,
2001 90);
2002 const Numeric zaEW = 90 - plevel_angletilt(rv_surface, c1);
2003 // Convert to Cartesian, and determine normal by cross-product
2004 Vector tangentSN(3), tangentEW(3), normal(3);
2005 zaaa2cart(tangentSN[0], tangentSN[1], tangentSN[2], zaSN, 0);
2006 zaaa2cart(tangentEW[0], tangentEW[1], tangentEW[2], zaEW, 90);
2007 cross3(normal, tangentSN, tangentEW);
2008 // Convert rtp_los to cartesian and flip direction
2009 Vector di(3);
2010 zaaa2cart(di[0], di[1], di[2], rtp_los[0], rtp_los[1]);
2011 di *= -1;
2012 // Set LOS normal vector
2013 cart2zaaa(
2014 surface_normal[0], surface_normal[1], normal[0], normal[1], normal[2]);
2015 if (abs(rtp_los[0] - surface_normal[0]) < 90) {
2016 throw runtime_error(
2017 "Invalid zenith angle. The zenith angle corresponds "
2018 "to observe the surface from below.");
2019 }
2020 // Specular direction is 2(dn*di)dn-di, where dn is the normal vector
2021 Vector speccart(3);
2022 const Numeric fac = 2 * (normal * di);
2023 for (Index i = 0; i < 3; i++) {
2024 speccart[i] = fac * normal[i] - di[i];
2025 }
2026 cart2zaaa(specular_los[0],
2027 specular_los[1],
2028 speccart[0],
2029 speccart[1],
2030 speccart[2]);
2031 }
2032}
2033
2034/* Workspace method: Doxygen documentation will be auto-generated */
2035void surfaceBlackbody(Matrix& surface_los,
2036 Tensor4& surface_rmatrix,
2037 Matrix& surface_emission,
2038 const Index& atmosphere_dim,
2039 const Vector& f_grid,
2040 const Index& stokes_dim,
2041 const Vector& rtp_pos,
2042 const Vector& rtp_los,
2043 const Numeric& surface_skin_t,
2044 const Verbosity& verbosity) {
2045 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2046 chk_rte_pos(atmosphere_dim, rtp_pos);
2047 chk_rte_los(atmosphere_dim, rtp_los);
2048 chk_not_negative("surface_skin_t", surface_skin_t);
2049
2051 out2 << " Sets variables to model a blackbody surface with a temperature "
2052 << " of " << surface_skin_t << " K.\n";
2053
2054 surface_los.resize(0, 0);
2055 surface_rmatrix.resize(0, 0, 0, 0);
2056
2057 const Index nf = f_grid.nelem();
2058
2059 Vector b(nf);
2060 planck(b, f_grid, surface_skin_t);
2061
2062 surface_emission.resize(nf, stokes_dim);
2063 surface_emission = 0.0;
2064
2065 for (Index iv = 0; iv < nf; iv++) {
2066 surface_emission(iv, 0) = b[iv];
2067 for (Index is = 1; is < stokes_dim; is++) {
2068 surface_emission(iv, is) = 0;
2069 }
2070 }
2071}
2072
2073/* Workspace method: Doxygen documentation will be auto-generated */
2074void surfaceFastem(Matrix& surface_los,
2075 Tensor4& surface_rmatrix,
2076 Matrix& surface_emission,
2077 const Index& atmosphere_dim,
2078 const Index& stokes_dim,
2079 const Vector& f_grid,
2080 const Vector& rtp_pos,
2081 const Vector& rtp_los,
2082 const Numeric& surface_skin_t,
2083 const Numeric& salinity,
2084 const Numeric& wind_speed,
2085 const Numeric& wind_direction,
2086 const Vector& transmittance,
2087 const Index& fastem_version,
2088 const Verbosity& verbosity) {
2089 // Input checks
2090 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2091 chk_rte_pos(atmosphere_dim, rtp_pos);
2092 chk_rte_los(atmosphere_dim, rtp_los);
2093 chk_if_in_range("wind_direction", wind_direction, -180, 180);
2094
2095 const Index nf = f_grid.nelem();
2096
2097 // Determine specular direction
2098 Vector specular_los, surface_normal;
2099 specular_losCalcNoTopography(specular_los,
2100 surface_normal,
2101 rtp_pos,
2102 rtp_los,
2103 atmosphere_dim,
2104 verbosity);
2105
2106 // Determine relative azimuth
2107 //
2108 // According to email from James Hocking, UkMet:
2109 // All azimuth angles are defined as being measured clockwise from north
2110 // (i.e. if the projection onto the Earth's surface of the view path lies due
2111 // north the azimuth angle is 0 and if it lies due east the azimuth angle is
2112 // 90 degrees). The relative azimuth is the wind direction (azimuth) angle
2113 // minus the satellite azimuth angle.
2114 Numeric rel_azimuth = wind_direction; // Always true for 1D
2115 if (atmosphere_dim == 2 && rtp_los[0] < 0) {
2116 rel_azimuth -= 180;
2117 resolve_lon(rel_azimuth, -180, 180);
2118 } else if (atmosphere_dim == 3) {
2119 rel_azimuth -= rtp_los[1];
2120 resolve_lon(rel_azimuth, -180, 180);
2121 }
2122
2123 // Call FASTEM
2124 Matrix emissivity, reflectivity;
2125 FastemStandAlone(emissivity,
2126 reflectivity,
2127 f_grid,
2128 surface_skin_t,
2129 abs(rtp_los[0]),
2130 salinity,
2131 wind_speed,
2132 rel_azimuth,
2133 transmittance,
2134 fastem_version,
2135 verbosity);
2136
2137 // Set surface_los
2138 surface_los.resize(1, specular_los.nelem());
2139 surface_los(0, joker) = specular_los;
2140
2141 // Surface emission
2142 //
2143 Vector b(nf);
2144 planck(b, f_grid, surface_skin_t);
2145 //
2146 surface_emission.resize(nf, stokes_dim);
2147 for (Index i = 0; i < nf; i++) {
2148 // I
2149 surface_emission(i, 0) = b[i] * 0.5 * (emissivity(i, 0) + emissivity(i, 1));
2150 // Q
2151 if (stokes_dim >= 2) {
2152 surface_emission(i, 1) =
2153 b[i] * 0.5 * (emissivity(i, 0) - emissivity(i, 1));
2154 }
2155 // U and V
2156 for (Index j = 2; j < stokes_dim; j++) {
2157 surface_emission(i, j) = b[i] * emissivity(i, j);
2158 }
2159 }
2160
2161 // Surface reflectivity matrix
2162 //
2163 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
2164 surface_rmatrix = 0.0;
2165 for (Index iv = 0; iv < nf; iv++) {
2166 surface_rmatrix(0, iv, 0, 0) =
2167 0.5 * (reflectivity(iv, 0) + reflectivity(iv, 1));
2168 if (stokes_dim >= 2) {
2169 surface_rmatrix(0, iv, 0, 1) =
2170 0.5 * (reflectivity(iv, 0) - reflectivity(iv, 1));
2171 ;
2172 surface_rmatrix(0, iv, 1, 0) = surface_rmatrix(0, iv, 0, 1);
2173 surface_rmatrix(0, iv, 1, 1) = surface_rmatrix(0, iv, 0, 0);
2174
2175 for (Index i = 2; i < stokes_dim; i++) {
2176 surface_rmatrix(0, iv, i, i) = surface_rmatrix(0, iv, 0, 0);
2177 }
2178 }
2179 }
2180}
2181
2182/* Workspace method: Doxygen documentation will be auto-generated */
2184 Tensor4& surface_rmatrix,
2185 const Index& stokes_dim,
2186 const Numeric& pol_angle,
2187 const Verbosity& ) {
2188 ARTS_USER_ERROR_IF (stokes_dim != 1,
2189 "You should only use this method where the main calculations are "
2190 "done with *stokes_dim* set to 1.");
2191
2192 Index local_stokes = surface_emission.ncols();
2193 ARTS_USER_ERROR_IF (local_stokes < 2,
2194 "This method requires that the input surface proporties match a Stokes "
2195 "dimension of 2-4. Incoming *surface_emission* matches stokes_dim=1.");
2196
2197 const Index nf = surface_emission.nrows();
2198 const Index nlos = surface_rmatrix.nbooks();
2199 Matrix se = surface_emission;
2200 Tensor4 sr = surface_rmatrix;
2201 surface_emission.resize(nf, 1);
2202 surface_rmatrix.resize(nlos, nf, 1, 1);
2203
2204 if (local_stokes == 2) {
2205 const Numeric alpha = DEG2RAD * pol_angle;
2206 const Numeric c2 = pow( cos(alpha), 2 );
2207 const Numeric s2 = pow( sin(alpha), 2 );
2208 for (Index f=0; f<nf; ++f) {
2209 const Numeric bv = se(f,0) + se(f,1); // Note that we have to multiply with 2
2210 const Numeric bh = se(f,0) - se(f,1); // as we place a single pol as I
2211 surface_emission(f,0) = c2*bv + s2*bh;
2212 for (Index l=0; l<nlos; ++l) {
2213 const Numeric rv = sr(l,f,0,0) + sr(l,f,1,0);
2214 const Numeric rh = sr(l,f,0,0) - sr(l,f,1,0);
2215 surface_rmatrix(l,f,0,0) = c2*rv + s2*rh;
2216 }
2217 }
2218 } else {
2219 Matrix Cm, Cs, Lp, Lm;
2220 mueller_stokes2modif(Cm, local_stokes);
2221 mueller_modif2stokes(Cs, local_stokes);
2222 mueller_rotation(Lp, local_stokes, pol_angle);
2223 mueller_rotation(Lm, local_stokes, -pol_angle);
2224 Matrix Mleft(local_stokes,local_stokes), Mright(local_stokes,local_stokes);
2225 mult(Mleft, Cm, Lp);
2226 mult(Mright, Lm, Cs);
2227 //
2228 Vector Vr(local_stokes);
2229 Matrix Tmp(local_stokes,local_stokes), Mr(local_stokes,local_stokes);
2230 //
2231 for (Index f=0; f<nf; ++f) {
2232 mult(Vr, Mleft, se(f,joker) ); // Note that we have to multiply with 2
2233 surface_emission(f,0) = 2.0 * Vr[0]; // as we place a single pol as I
2234 for (Index l=0; l<nlos; ++l) {
2235 mult(Tmp, sr(l,f,joker,joker), Mright);
2236 mult(Mr, Mleft, Tmp);
2237 surface_rmatrix(l,f,0,0) = Tmp(0,0);
2238 }
2239 }
2240 }
2241}
2242
2243/* Workspace method: Doxygen documentation will be auto-generated */
2244void surfaceTelsem(Matrix& surface_los,
2245 Tensor4& surface_rmatrix,
2246 Matrix& surface_emission,
2247 const Index& atmosphere_dim,
2248 const Index& stokes_dim,
2249 const Vector& f_grid,
2250 const Vector& lat_grid,
2251 const Vector& lat_true,
2252 const Vector& lon_true,
2253 const Vector& rtp_pos,
2254 const Vector& rtp_los,
2255 const Vector& specular_los,
2256 const Numeric& surface_skin_t,
2257 const TelsemAtlas& atlas,
2258 const Numeric& r_min,
2259 const Numeric& r_max,
2260 const Numeric& d_max,
2261 const Verbosity& verbosity) {
2262 // Input checks
2263 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2264 chk_rte_pos(atmosphere_dim, rtp_pos);
2265 chk_rte_los(atmosphere_dim, rtp_los);
2266 chk_rte_los(atmosphere_dim, specular_los);
2268 "surface skin temperature", surface_skin_t, 190.0, 373.0);
2269
2270 const Index nf = f_grid.nelem();
2271 Matrix surface_rv_rh(nf, 2);
2272
2273 // Lookup cell in atlas.
2274
2275 Numeric lat, lon;
2277 lat, lon, atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
2278 chk_if_in_range("Latitude input to TELSEM2", lat, -90.0, 90.0);
2279
2280 lon = fmod(lon, 360.0);
2281 if (lon < 0) {
2282 lon += 360.0;
2283 }
2284 chk_if_in_range("Longitude input to TELSEM2", lon, 0.0, 360.0);
2285
2286 Index cellnumber = atlas.calc_cellnum(lat, lon);
2287 // Check if cell is in atlas.
2288 if (!atlas.contains(cellnumber)) {
2289 if (d_max <= 0.0) {
2290 throw std::runtime_error(
2291 "Given coordinates are not contained in "
2292 " TELSEM atlas. To enable nearest neighbour"
2293 "interpolation set *d_max* to a positive "
2294 "value.");
2295 } else {
2296 cellnumber = atlas.calc_cellnum_nearest_neighbor(lat, lon);
2297 Numeric lat_nn, lon_nn;
2298 std::tie(lat_nn, lon_nn) = atlas.get_coordinates(cellnumber);
2299 Numeric d = sphdist(lat, lon, lat_nn, lon_nn);
2300 if (d > d_max) {
2301 std::ostringstream out{};
2302 out << "Distance of nearest neighbour exceeds provided limit (";
2303 out << d << " > " << d_max << ").";
2304 throw std::runtime_error(out.str());
2305 }
2306 }
2307 }
2308
2309 Index class1 = atlas.get_class1(cellnumber);
2310 Index class2 = atlas.get_class2(cellnumber);
2311
2312 Vector emis_h = atlas.get_emis_h(cellnumber);
2313 Vector emis_v = atlas.get_emis_v(cellnumber);
2314 Numeric e_v, e_h;
2315
2316 // Incidence angle
2317 const Numeric theta = calc_incang(rtp_los, specular_los);
2318 ARTS_ASSERT(theta <= 90);
2319
2320 for (Index i = 0; i < nf; ++i) {
2321 if (f_grid[i] < 5e9)
2322 throw std::runtime_error("Only frequency >= 5 GHz are allowed");
2323 if (f_grid[i] > 900e9)
2324 throw std::runtime_error("Only frequency <= 900 GHz are allowed");
2325
2326 // For frequencies 700 <= f <= 900 GHz use 700 GHz to avoid extrapolation
2327 // outside of TELSEM specifications.
2328 Numeric f = std::min(f_grid[i], 700e9) * 1e-9;
2329 std::tie(e_v, e_h) =
2330 atlas.emis_interp(theta, f, class1, class2, emis_v, emis_h);
2331
2332 surface_rv_rh(i, 0) = min(max(1.0 - e_v, r_min), r_max);
2333 surface_rv_rh(i, 1) = min(max(1.0 - e_h, r_min), r_max);
2334 }
2335
2336 surfaceFlatRvRh(surface_los,
2337 surface_rmatrix,
2338 surface_emission,
2339 f_grid,
2340 stokes_dim,
2341 atmosphere_dim,
2342 rtp_pos,
2343 rtp_los,
2344 specular_los,
2345 surface_skin_t,
2346 surface_rv_rh,
2347 verbosity);
2348}
2349
2350/* Workspace method: Doxygen documentation will be auto-generated */
2351void surfaceTessem(Matrix& surface_los,
2352 Tensor4& surface_rmatrix,
2353 Matrix& surface_emission,
2354 const Index& atmosphere_dim,
2355 const Index& stokes_dim,
2356 const Vector& f_grid,
2357 const Vector& rtp_pos,
2358 const Vector& rtp_los,
2359 const Numeric& surface_skin_t,
2360 const TessemNN& net_h,
2361 const TessemNN& net_v,
2362 const Numeric& salinity,
2363 const Numeric& wind_speed,
2364 const Verbosity& verbosity) {
2365 // Input checks
2366 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2367 chk_rte_pos(atmosphere_dim, rtp_pos);
2368 chk_rte_los(atmosphere_dim, rtp_los);
2370 "surface skin temperature", surface_skin_t, 260.0, 373.0);
2371 chk_if_in_range_exclude_high("salinity", salinity, 0, 1);
2372 chk_if_in_range_exclude_high("wind speed", wind_speed, 0, 100);
2373
2374 // Determine specular direction
2375 Vector specular_los, surface_normal;
2376 specular_losCalcNoTopography(specular_los,
2377 surface_normal,
2378 rtp_pos,
2379 rtp_los,
2380 atmosphere_dim,
2381 verbosity);
2382
2383 // TESSEM in and out
2384 //
2385 Vector out(2);
2386 VectorView e_h = out[Range(0, 1)];
2387 VectorView e_v = out[Range(1, 1)];
2388 //
2389 Vector in(5);
2390 in[1] = 180.0 - abs(rtp_los[0]);
2391 in[2] = wind_speed;
2392 in[3] = surface_skin_t;
2393 in[4] = salinity;
2394
2395 // Get Rv and Rh
2396 //
2397 const Index nf = f_grid.nelem();
2398 Matrix surface_rv_rh(nf, 2);
2399 //
2400 for (Index i = 0; i < nf; ++i) {
2401 if (f_grid[i] < 5e9)
2402 throw std::runtime_error("Only frequency >= 5 GHz are allowed");
2403 if (f_grid[i] > 900e9)
2404 throw std::runtime_error("Only frequency <= 900 GHz are allowed");
2405
2406 in[0] = f_grid[i];
2407
2408 tessem_prop_nn(e_h, net_h, in);
2409 tessem_prop_nn(e_v, net_v, in);
2410
2411 surface_rv_rh(i, 0) = min(max(1 - e_v[0], (Numeric)0), (Numeric)1);
2412 surface_rv_rh(i, 1) = min(max(1 - e_h[0], (Numeric)0), (Numeric)1);
2413 }
2414
2415 surfaceFlatRvRh(surface_los,
2416 surface_rmatrix,
2417 surface_emission,
2418 f_grid,
2419 stokes_dim,
2420 atmosphere_dim,
2421 rtp_pos,
2422 rtp_los,
2423 specular_los,
2424 surface_skin_t,
2425 surface_rv_rh,
2426 verbosity);
2427}
2428
2429/* Workspace method: Doxygen documentation will be auto-generated */
2431 Tensor4& surface_rmatrix,
2432 Matrix& surface_emission,
2433 const Vector& f_grid,
2434 const Index& stokes_dim,
2435 const Index& atmosphere_dim,
2436 const Vector& rtp_pos,
2437 const Vector& rtp_los,
2438 const Vector& specular_los,
2439 const Numeric& surface_skin_t,
2440 const GriddedField3& surface_complex_refr_index,
2441 const Verbosity& verbosity) {
2444
2445 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2446 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2447 chk_rte_pos(atmosphere_dim, rtp_pos);
2448 chk_rte_los(atmosphere_dim, rtp_los);
2449 chk_rte_los(atmosphere_dim, specular_los);
2450 chk_not_negative("surface_skin_t", surface_skin_t);
2451
2452 // Interpolate *surface_complex_refr_index*
2453 //
2454 const Index nf = f_grid.nelem();
2455 //
2456 Matrix n_real(nf, 1), n_imag(nf, 1);
2457 //
2458 complex_n_interp(n_real,
2459 n_imag,
2460 surface_complex_refr_index,
2461 "surface_complex_refr_index",
2462 f_grid,
2463 Vector(1, surface_skin_t));
2464
2465 out2 << " Sets variables to model a flat surface\n";
2466 out3 << " surface temperature: " << surface_skin_t << " K.\n";
2467
2468 surface_los.resize(1, specular_los.nelem());
2469 surface_los(0, joker) = specular_los;
2470
2471 surface_emission.resize(nf, stokes_dim);
2472 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
2473
2474 // Incidence angle
2475 const Numeric incang = calc_incang(rtp_los, specular_los);
2476 ARTS_ASSERT(incang <= 90);
2477
2478 // Complex (amplitude) reflection coefficients
2479 Complex Rv, Rh;
2480
2481 for (Index iv = 0; iv < nf; iv++) {
2482 // Set n2 (refractive index of surface medium)
2483 Complex n2(n_real(iv, 0), n_imag(iv, 0));
2484
2485 // Amplitude reflection coefficients
2486 fresnel(Rv, Rh, Numeric(1.0), n2, incang);
2487
2488 // Fill reflection matrix and emission vector
2489 surface_specular_R_and_b(surface_rmatrix(0, iv, joker, joker),
2490 surface_emission(iv, joker),
2491 Rv,
2492 Rh,
2493 f_grid[iv],
2494 stokes_dim,
2495 surface_skin_t);
2496 }
2497}
2498
2499/* Workspace method: Doxygen documentation will be auto-generated */
2501 Tensor4& surface_rmatrix,
2502 Matrix& surface_emission,
2503 const Vector& f_grid,
2504 const Index& stokes_dim,
2505 const Index& atmosphere_dim,
2506 const Vector& rtp_pos,
2507 const Vector& rtp_los,
2508 const Vector& specular_los,
2509 const Numeric& surface_skin_t,
2510 const Tensor3& surface_reflectivity,
2511 const Verbosity& verbosity) {
2513
2514 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2515 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2516 chk_rte_pos(atmosphere_dim, rtp_pos);
2517 chk_rte_los(atmosphere_dim, rtp_los);
2518 chk_rte_los(atmosphere_dim, specular_los);
2519 chk_not_negative("surface_skin_t", surface_skin_t);
2520
2521 const Index nf = f_grid.nelem();
2522
2523 if (surface_reflectivity.nrows() != stokes_dim &&
2524 surface_reflectivity.ncols() != stokes_dim) {
2525 ostringstream os;
2526 os << "The number of rows and columnss in *surface_reflectivity* must\n"
2527 << "match *stokes_dim*."
2528 << "\n stokes_dim : " << stokes_dim
2529 << "\n number of rows in *surface_reflectivity* : "
2530 << surface_reflectivity.nrows()
2531 << "\n number of columns in *surface_reflectivity* : "
2532 << surface_reflectivity.ncols() << "\n";
2533 throw runtime_error(os.str());
2534 }
2535
2536 if (surface_reflectivity.npages() != nf &&
2537 surface_reflectivity.npages() != 1) {
2538 ostringstream os;
2539 os << "The number of pages in *surface_reflectivity* should\n"
2540 << "match length of *f_grid* or be 1."
2541 << "\n length of *f_grid* : " << nf
2542 << "\n dimension of *surface_reflectivity* : "
2543 << surface_reflectivity.npages() << "\n";
2544 throw runtime_error(os.str());
2545 }
2546
2547 out2 << " Sets variables to model a flat surface\n";
2548
2549 surface_los.resize(1, specular_los.nelem());
2550 surface_los(0, joker) = specular_los;
2551
2552 surface_emission.resize(nf, stokes_dim);
2553 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
2554
2555 Matrix R, IR(stokes_dim, stokes_dim);
2556
2557 Vector b(nf);
2558 planck(b, f_grid, surface_skin_t);
2559
2560 Vector B(stokes_dim, 0);
2561
2562 for (Index iv = 0; iv < nf; iv++) {
2563 if (iv == 0 || surface_reflectivity.npages() > 1) {
2564 R = surface_reflectivity(iv, joker, joker);
2565 for (Index i = 0; i < stokes_dim; i++) {
2566 for (Index j = 0; j < stokes_dim; j++) {
2567 if (i == j) {
2568 IR(i, j) = 1 - R(i, j);
2569 } else {
2570 IR(i, j) = -R(i, j);
2571 }
2572 }
2573 }
2574 }
2575
2576 surface_rmatrix(0, iv, joker, joker) = R;
2577
2578 B[0] = b[iv];
2579 mult(surface_emission(iv, joker), IR, B);
2580 }
2581}
2582
2583/* Workspace method: Doxygen documentation will be auto-generated */
2584void surfaceFlatRvRh(Matrix& surface_los,
2585 Tensor4& surface_rmatrix,
2586 Matrix& surface_emission,
2587 const Vector& f_grid,
2588 const Index& stokes_dim,
2589 const Index& atmosphere_dim,
2590 const Vector& rtp_pos,
2591 const Vector& rtp_los,
2592 const Vector& specular_los,
2593 const Numeric& surface_skin_t,
2594 const Matrix& surface_rv_rh,
2595 const Verbosity&) {
2596 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2597 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2598 chk_rte_pos(atmosphere_dim, rtp_pos);
2599 chk_rte_los(atmosphere_dim, rtp_los);
2600 chk_rte_los(atmosphere_dim, specular_los);
2601 chk_not_negative("surface_skin_t", surface_skin_t);
2602
2603 const Index nf = f_grid.nelem();
2604
2605 if (surface_rv_rh.ncols() != 2) {
2606 ostringstream os;
2607 os << "The number of columns in *surface_rv_rh* must be two,\n"
2608 << "but the actual number of columns is " << surface_rv_rh.ncols()
2609 << "\n";
2610 throw runtime_error(os.str());
2611 }
2612
2613 if (surface_rv_rh.nrows() != nf && surface_rv_rh.nrows() != 1) {
2614 ostringstream os;
2615 os << "The number of rows in *surface_rv_rh* should\n"
2616 << "match length of *f_grid* or be 1."
2617 << "\n length of *f_grid* : " << nf
2618 << "\n rows in *surface_rv_rh* : " << surface_rv_rh.nrows() << "\n";
2619 throw runtime_error(os.str());
2620 }
2621
2622 if (min(surface_rv_rh) < 0 || max(surface_rv_rh) > 1) {
2623 throw runtime_error("All values in *surface_rv_rh* must be inside [0,1].");
2624 }
2625
2626 surface_los.resize(1, specular_los.nelem());
2627 surface_los(0, joker) = specular_los;
2628
2629 surface_emission.resize(nf, stokes_dim);
2630 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
2631
2632 surface_emission = 0;
2633 surface_rmatrix = 0;
2634
2635 Vector b(nf);
2636 planck(b, f_grid, surface_skin_t);
2637
2638 Numeric rmean = 0.0, rdiff = 0.0;
2639
2640 for (Index iv = 0; iv < nf; iv++) {
2641 if (iv == 0 || surface_rv_rh.nrows() > 1) {
2642 rmean = 0.5 * (surface_rv_rh(iv, 0) + surface_rv_rh(iv, 1));
2643 rdiff = 0.5 * (surface_rv_rh(iv, 0) - surface_rv_rh(iv, 1));
2644 }
2645
2646 surface_emission(iv, 0) = (1.0 - rmean) * b[iv];
2647 surface_rmatrix(0, iv, 0, 0) = rmean;
2648
2649 if (stokes_dim > 1) {
2650 surface_emission(iv, 1) = -rdiff * b[iv];
2651
2652 surface_rmatrix(0, iv, 0, 1) = rdiff;
2653 surface_rmatrix(0, iv, 1, 0) = rdiff;
2654 surface_rmatrix(0, iv, 1, 1) = rmean;
2655
2656 for (Index i = 2; i < stokes_dim; i++) {
2657 surface_rmatrix(0, iv, i, i) = rmean;
2658 }
2659 }
2660 }
2661}
2662
2663/* Workspace method: Doxygen documentation will be auto-generated */
2665 Tensor4& surface_rmatrix,
2666 Matrix& surface_emission,
2667 const Vector& f_grid,
2668 const Index& stokes_dim,
2669 const Index& atmosphere_dim,
2670 const Vector& rtp_pos,
2671 const Vector& rtp_los,
2672 const Vector& specular_los,
2673 const Numeric& surface_skin_t,
2674 const Vector& surface_scalar_reflectivity,
2675 const Verbosity&) {
2676 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2677 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2678 chk_rte_pos(atmosphere_dim, rtp_pos);
2679 chk_rte_los(atmosphere_dim, rtp_los);
2680 chk_rte_los(atmosphere_dim, specular_los);
2681 chk_not_negative("surface_skin_t", surface_skin_t);
2682
2683 const Index nf = f_grid.nelem();
2684
2685 if (surface_scalar_reflectivity.nelem() != nf &&
2686 surface_scalar_reflectivity.nelem() != 1) {
2687 ostringstream os;
2688 os << "The number of elements in *surface_scalar_reflectivity* should\n"
2689 << "match length of *f_grid* or be 1."
2690 << "\n length of *f_grid* : " << nf
2691 << "\n length of *surface_scalar_reflectivity* : "
2692 << surface_scalar_reflectivity.nelem() << "\n";
2693 throw runtime_error(os.str());
2694 }
2695
2696 if (min(surface_scalar_reflectivity) < 0 ||
2697 max(surface_scalar_reflectivity) > 1) {
2698 throw runtime_error(
2699 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
2700 }
2701
2702 surface_los.resize(1, specular_los.nelem());
2703 surface_los(0, joker) = specular_los;
2704
2705 surface_emission.resize(nf, stokes_dim);
2706 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
2707
2708 surface_emission = 0;
2709 surface_rmatrix = 0;
2710
2711 Vector b(nf);
2712 planck(b, f_grid, surface_skin_t);
2713
2714 Numeric r = 0.0;
2715
2716 for (Index iv = 0; iv < nf; iv++) {
2717 if (iv == 0 || surface_scalar_reflectivity.nelem() > 1) {
2718 r = surface_scalar_reflectivity[iv];
2719 }
2720
2721 surface_emission(iv, 0) = (1.0 - r) * b[iv];
2722 surface_rmatrix(0, iv, 0, 0) = r;
2723 for (Index i = 1; i < stokes_dim; i++) {
2724 surface_rmatrix(0, iv, i, i) = r;
2725 }
2726 }
2727}
2728
2729/* Workspace method: Doxygen documentation will be auto-generated */
2731 Tensor4& surface_rmatrix,
2732 Matrix& surface_emission,
2733 const Vector& f_grid,
2734 const Index& stokes_dim,
2735 const Index& atmosphere_dim,
2736 const Vector& rtp_pos,
2737 const Vector& rtp_los,
2738 const Vector& surface_normal,
2739 const Numeric& surface_skin_t,
2740 const Vector& surface_scalar_reflectivity,
2741 const Index& lambertian_nza,
2742 const Numeric& za_pos,
2743 const Verbosity&) {
2744 const Index nf = f_grid.nelem();
2745
2746 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2747 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2748 chk_rte_pos(atmosphere_dim, rtp_pos);
2749 chk_rte_los(atmosphere_dim, rtp_los);
2750 chk_not_negative("surface_skin_t", surface_skin_t);
2751 chk_if_in_range("za_pos", za_pos, 0, 1);
2752
2753 if (surface_scalar_reflectivity.nelem() != nf &&
2754 surface_scalar_reflectivity.nelem() != 1) {
2755 ostringstream os;
2756 os << "The number of elements in *surface_scalar_reflectivity* should\n"
2757 << "match length of *f_grid* or be 1."
2758 << "\n length of *f_grid* : " << nf
2759 << "\n length of *surface_scalar_reflectivity* : "
2760 << surface_scalar_reflectivity.nelem() << "\n";
2761 throw runtime_error(os.str());
2762 }
2763
2764 if (min(surface_scalar_reflectivity) < 0 ||
2765 max(surface_scalar_reflectivity) > 1) {
2766 throw runtime_error(
2767 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
2768 }
2769
2770 // Allocate and init everything to zero
2771 //
2772 surface_los.resize(lambertian_nza, rtp_los.nelem());
2773 surface_rmatrix.resize(lambertian_nza, nf, stokes_dim, stokes_dim);
2774 surface_emission.resize(nf, stokes_dim);
2775 //
2776 surface_los = 0.0;
2777 surface_rmatrix = 0.0;
2778 surface_emission = 0.0;
2779
2780 // Help variables
2781 //
2782 const Numeric dza = (90.0 - abs(surface_normal[0])) / (Numeric)lambertian_nza;
2783 const Vector za_lims(0.0, lambertian_nza + 1, dza);
2784
2785 // surface_los
2786 for (Index ip = 0; ip < lambertian_nza; ip++) {
2787 surface_los(ip, 0) = za_lims[ip] + za_pos * dza;
2788 if (atmosphere_dim == 2) {
2789 if (rtp_los[0] < 0) {
2790 surface_los(ip, 0) *= -1.0;
2791 }
2792
2793 } else if (atmosphere_dim == 3) {
2794 surface_los(ip, 1) = rtp_los[1];
2795 }
2796 }
2797
2798 Vector b(nf);
2799 planck(b, f_grid, surface_skin_t);
2800
2801 // Loop frequencies and set remaining values
2802 //
2803 Numeric r = 0.0;
2804 //
2805 for (Index iv = 0; iv < nf; iv++) {
2806 // Get reflectivity
2807 if (iv == 0 || surface_scalar_reflectivity.nelem() > 1) {
2808 r = surface_scalar_reflectivity[iv];
2809 }
2810
2811 // surface_rmatrix:
2812 // Only element (0,0) is set to be non-zero. This follows VDISORT
2813 // that refers to: K. L. Coulson, Polarization and Intensity of Light in
2814 // the Atmosphere (1989), page 229
2815 // (Thanks to Michael Kahnert for providing this information!)
2816 // Update: Is the above for a later edition? We have not found a copy of
2817 // that edition. In a 1988 version of the book, the relevant page seems
2818 // to be 232.
2819 for (Index ip = 0; ip < lambertian_nza; ip++) {
2820 const Numeric w =
2821 r * 0.5 *
2822 (cos(2 * DEG2RAD * za_lims[ip]) - cos(2 * DEG2RAD * za_lims[ip + 1]));
2823 surface_rmatrix(ip, iv, 0, 0) = w;
2824 }
2825
2826 // surface_emission
2827 surface_emission(iv, 0) = (1 - r) * b[iv];
2828 }
2829}
2830
2831/* Workspace method: Doxygen documentation will be auto-generated */
2833 Numeric& surface_skin_t,
2834 Matrix& surface_los,
2835 Tensor4& surface_rmatrix,
2836 Matrix& surface_emission,
2837 const Index& atmosphere_dim,
2838 const Vector& f_grid,
2839 const Vector& rtp_pos,
2840 const Vector& rtp_los,
2841 const Agenda& surface_rtprop_sub_agenda,
2842 const Numeric& specular_factor,
2843 const Numeric& dza,
2844 const Verbosity&) {
2845 chk_rte_pos(atmosphere_dim, rtp_pos);
2846 chk_rte_los(atmosphere_dim, rtp_los);
2847
2848 // Checks of GIN variables
2849 if (specular_factor > 1 || specular_factor < 1.0 / 3.0)
2850 throw runtime_error("The valid range for *specular_factor* is [1/3,1].");
2851 if (dza > 45 || dza <= 0)
2852 throw runtime_error("The valid range for *dza* is ]0,45].");
2853
2854 // Obtain data for specular direction
2855 //
2856 Matrix los1, emission1;
2857 Tensor4 rmatrix1;
2858 //
2860 surface_skin_t,
2861 emission1,
2862 los1,
2863 rmatrix1,
2864 f_grid,
2865 rtp_pos,
2866 rtp_los,
2867 surface_rtprop_sub_agenda);
2868 if (los1.nrows() != 1)
2869 throw runtime_error(
2870 "*surface_rtprop_sub_agenda* must return data "
2871 "describing a specular surface.");
2872
2873 // Standard number of beams. Set to 2 if try/catch below fails
2874 Index nbeams = 3;
2875
2876 // Test if some lower zenith angle works.
2877 // It will fail if a higher za results in looking at the surface from below.
2878 //
2879 Matrix los2, emission2;
2880 Tensor4 rmatrix2;
2881 //
2882 Numeric skin_t_dummy;
2883 Numeric dza_try = dza;
2884 bool failed = true;
2885 while (failed && dza_try > 0) {
2886 try {
2887 Vector los_new = rtp_los;
2888 los_new[0] -=
2889 sign(rtp_los[0]) * dza_try; // Sign to also handle 2D negative za
2890 adjust_los(los_new, atmosphere_dim);
2892 skin_t_dummy,
2893 emission2,
2894 los2,
2895 rmatrix2,
2896 f_grid,
2897 rtp_pos,
2898 los_new,
2899 surface_rtprop_sub_agenda);
2900 failed = false;
2901 } catch (const runtime_error& e) {
2902 dza_try -= 1.0;
2903 }
2904 }
2905 if (failed) {
2906 nbeams = 2;
2907 }
2908
2909 // Allocate output WSVs
2910 //
2911 surface_emission.resize(emission1.nrows(), emission1.ncols());
2912 surface_emission = 0;
2913 surface_los.resize(nbeams, los1.ncols());
2914 surface_rmatrix.resize(
2915 nbeams, rmatrix1.npages(), rmatrix1.nrows(), rmatrix1.ncols());
2916
2917 // Put in specular direction at index 1
2918 //
2919 Numeric w;
2920 if (nbeams == 3) {
2921 w = specular_factor;
2922 } else {
2923 w = specular_factor + (1.0 - specular_factor) / 2.0;
2924 }
2925 //
2926 surface_los(1, joker) = los1(0, joker);
2927 for (Index p = 0; p < rmatrix1.npages(); p++) {
2928 for (Index r = 0; r < rmatrix1.nrows(); r++) {
2929 surface_emission(p, r) += w * emission1(p, r);
2930 for (Index c = 0; c < rmatrix1.ncols(); c++) {
2931 surface_rmatrix(1, p, r, c) = w * rmatrix1(0, p, r, c);
2932 }
2933 }
2934 }
2935
2936 // Put in lower za as index 2, if worked
2937 //
2938 w = (1.0 - specular_factor) / 2.0;
2939 //
2940 if (nbeams == 3) {
2941 surface_los(2, joker) = los2(0, joker);
2942 for (Index p = 0; p < rmatrix2.npages(); p++) {
2943 for (Index r = 0; r < rmatrix2.nrows(); r++) {
2944 surface_emission(p, r) += w * emission2(p, r);
2945 for (Index c = 0; c < rmatrix1.ncols(); c++) {
2946 surface_rmatrix(2, p, r, c) = w * rmatrix2(0, p, r, c);
2947 }
2948 }
2949 }
2950 }
2951
2952 // Do higher za and put in as index 0 (reusing variables for beam 2)
2953 //
2954 Vector los_new = rtp_los;
2955 los_new[0] += sign(rtp_los[0]) * dza; // Sign to also handle 2D negative za
2956 adjust_los(los_new, atmosphere_dim);
2958 skin_t_dummy,
2959 emission2,
2960 los2,
2961 rmatrix2,
2962 f_grid,
2963 rtp_pos,
2964 los_new,
2965 surface_rtprop_sub_agenda);
2966 //
2967 surface_los(0, joker) = los2(0, joker);
2968 for (Index p = 0; p < rmatrix2.npages(); p++) {
2969 for (Index r = 0; r < rmatrix2.nrows(); r++) {
2970 surface_emission(p, r) += w * emission2(p, r);
2971 for (Index c = 0; c < rmatrix1.ncols(); c++) {
2972 surface_rmatrix(0, p, r, c) = w * rmatrix2(0, p, r, c);
2973 }
2974 }
2975 }
2976}
2977
2978/* Workspace method: Doxygen documentation will be auto-generated */
2980 Tensor4& surface_rmatrix,
2981 const Index& atmosphere_dim,
2982 const Vector& rtp_pos,
2983 const Vector& rtp_los,
2984 const Numeric& specular_factor,
2985 const Numeric& dza,
2986 const Verbosity&) {
2987 chk_rte_pos(atmosphere_dim, rtp_pos);
2988 chk_rte_los(atmosphere_dim, rtp_los);
2989
2990 // Check that input surface data are of specular type
2991 if (surface_los.nrows() != 1)
2992 throw runtime_error(
2993 "Input surface data must be of specular type. That is, "
2994 "*surface_los* must contain a single direction.");
2995 if (surface_rmatrix.nbooks() != 1)
2996 throw runtime_error(
2997 "*surface_rmatrix* describes a different number of "
2998 "directions than *surface_los*.");
2999
3000 // Checks of GIN variables
3001 if (specular_factor > 1 || specular_factor < 1.0 / 3.0)
3002 throw runtime_error("The valid range for *specular_factor* is [1/3,1].");
3003 if (dza > 45 || dza <= 0)
3004 throw runtime_error("The valid range for *dza* is ]0,45].");
3005
3006 // Make copies of input data
3007 const Matrix los1 = surface_los;
3008 const Tensor4 rmatrix1 = surface_rmatrix;
3009
3010 // Use abs(za) in all expressions below, to also handle 2D
3011
3012 // Calculate highest possible za for downwelling radiation, with 1 degree
3013 // margin to the surface. This can be derived from za in rtp_los and surface_los.
3014 // (The directions in surface_los are not allowed to point into the surface)
3015 const Numeric za_max = 89 + (180 - abs(los1(0, 0)) - abs(rtp_los[0])) / 2.0;
3016
3017 // Number of downwelling beams
3018 Index nbeams = 3;
3019 if (abs(los1(0, 0)) > za_max) {
3020 nbeams = 2;
3021 }
3022
3023 // New los-s
3024 //
3025 surface_los.resize(nbeams, los1.ncols());
3026 //
3027 for (Index r = 0; r < nbeams; r++) {
3028 surface_los(r, 0) = ((Numeric)r - 1.0) * dza + abs(los1(0, 0));
3029 if (r == 2 && surface_los(r, 0) > za_max) {
3030 surface_los(r, 0) = za_max;
3031 }
3032 for (Index c = 1; c < los1.ncols(); c++) {
3033 surface_los(r, c) = los1(0, c);
3034 }
3035 }
3036
3037 // New rmatrix
3038 //
3039 surface_rmatrix.resize(
3040 nbeams, rmatrix1.npages(), rmatrix1.nrows(), rmatrix1.ncols());
3041 //
3042 for (Index b = 0; b < nbeams; b++) {
3043 Numeric w;
3044 if (b == 1 && nbeams == 3) // Specular direction with nbeams==3
3045 {
3046 w = specular_factor;
3047 } else if (b == 1) // Specular direction with nbeams==2
3048 {
3049 w = specular_factor + (1.0 - specular_factor) / 2.0;
3050 } else // Side directions
3051 {
3052 w = (1.0 - specular_factor) / 2.0;
3053 }
3054
3055 for (Index p = 0; p < rmatrix1.npages(); p++) {
3056 for (Index r = 0; r < rmatrix1.nrows(); r++) {
3057 for (Index c = 0; c < rmatrix1.ncols(); c++) {
3058 surface_rmatrix(b, p, r, c) = w * rmatrix1(0, p, r, c);
3059 }
3060 }
3061 }
3062 }
3063
3064 // Handle sign of za
3065 if (atmosphere_dim == 1) {
3066 // We only need to make sure that first direction has positive za
3067 surface_los(0, 0) = abs(surface_los(0, 0));
3068 } else if (atmosphere_dim == 2) {
3069 // Change sign if specular direction has za < 0
3070 if (los1(0, 0) < 0) {
3071 for (Index r = 0; r < rmatrix1.nrows(); r++) {
3072 surface_los(r, 0) = -surface_los(r, 0);
3073 }
3074 }
3075 } else if (atmosphere_dim == 1) {
3076 // We only need to make sure that first direction has positive za
3077 if (surface_los(0, 0) < 0) {
3078 surface_los(0, 0) = -surface_los(0, 0);
3079 surface_los(0, 1) += 180;
3080 if (surface_los(0, 1) > 180) {
3081 surface_los(0, 1) -= 360;
3082 }
3083 }
3084 }
3085}
3086
3087/* Workspace method: Doxygen documentation will be auto-generated */
3089 GriddedField3& surface_complex_refr_index,
3090 const Index& atmosphere_dim,
3091 const Vector& lat_grid,
3092 const Vector& lat_true,
3093 const Vector& lon_true,
3094 const Vector& rtp_pos,
3095 const GriddedField5& complex_n_field,
3096 const Verbosity&) {
3097 // Set expected order of grids
3098 Index gfield_fID = 0;
3099 Index gfield_tID = 1;
3100 Index gfield_compID = 2;
3101 Index gfield_latID = 3;
3102 Index gfield_lonID = 4;
3103
3104 // Basic checks and sizes
3105 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3106 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
3107 chk_rte_pos(atmosphere_dim, rtp_pos);
3108 complex_n_field.checksize_strict();
3109 //
3110 chk_griddedfield_gridname(complex_n_field, gfield_fID, "Frequency");
3111 chk_griddedfield_gridname(complex_n_field, gfield_tID, "Temperature");
3112 chk_griddedfield_gridname(complex_n_field, gfield_compID, "Complex");
3113 chk_griddedfield_gridname(complex_n_field, gfield_latID, "Latitude");
3114 chk_griddedfield_gridname(complex_n_field, gfield_lonID, "Longitude");
3115 //
3116 const Index nf = complex_n_field.data.nshelves();
3117 const Index nt = complex_n_field.data.nbooks();
3118 const Index nn = complex_n_field.data.npages();
3119 const Index nlat = complex_n_field.data.nrows();
3120 const Index nlon = complex_n_field.data.ncols();
3121 //
3122 if (nlat < 2 || nlon < 2) {
3123 ostringstream os;
3124 os << "The data in *complex_refr_index_field* must span a geographical "
3125 << "region. That is,\nthe latitude and longitude grids must have a "
3126 << "length >= 2.";
3127 throw runtime_error(os.str());
3128 }
3129 //
3130 if (nn != 2) {
3131 ostringstream os;
3132 os << "The data in *complex_refr_index_field* must have exactly two "
3133 << "pages. One page each\nfor the real and imaginary part of the "
3134 << "complex refractive index.";
3135 throw runtime_error(os.str());
3136 }
3137
3138 const Vector& GFlat = complex_n_field.get_numeric_grid(gfield_latID);
3139 const Vector& GFlon = complex_n_field.get_numeric_grid(gfield_lonID);
3140
3141 // Determine true geographical position
3142 Vector lat(1), lon(1);
3144 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
3145
3146 // Ensure correct coverage of lon grid
3147 Vector lon_shifted;
3148 lon_shiftgrid(lon_shifted, GFlon, lon[0]);
3149
3150 // Check if lat/lon we need are actually covered
3151 chk_if_in_range("rtp_pos.lat", lat[0], GFlat[0], GFlat[nlat - 1]);
3152 chk_if_in_range("rtp_pos.lon", lon[0], lon_shifted[0], lon_shifted[nlon - 1]);
3153
3154 // Size and fills grids of *surface_complex_refr_index*
3155 surface_complex_refr_index.resize(nf, nt, 2);
3156 surface_complex_refr_index.set_grid_name(0, "Frequency");
3157 surface_complex_refr_index.set_grid(
3158 0, complex_n_field.get_numeric_grid(gfield_fID));
3159 surface_complex_refr_index.set_grid_name(1, "Temperature");
3160 surface_complex_refr_index.set_grid(
3161 1, complex_n_field.get_numeric_grid(gfield_tID));
3162 surface_complex_refr_index.set_grid_name(2, "Complex");
3163 surface_complex_refr_index.set_grid(2, {"real", "imaginary"});
3164
3165 // Interpolate in lat and lon
3166 //
3167 GridPos gp_lat, gp_lon;
3168 gridpos(gp_lat, GFlat, lat[0]);
3169 gridpos(gp_lon, lon_shifted, lon[0]);
3170 Vector itw(4);
3171 interpweights(itw, gp_lat, gp_lon);
3172 //
3173 for (Index iv = 0; iv < nf; iv++) {
3174 for (Index it = 0; it < nt; it++) {
3175 surface_complex_refr_index.data(iv, it, 0) = interp(
3176 itw, complex_n_field.data(iv, it, 0, joker, joker), gp_lat, gp_lon);
3177 surface_complex_refr_index.data(iv, it, 1) = interp(
3178 itw, complex_n_field.data(iv, it, 1, joker, joker), gp_lat, gp_lon);
3179 }
3180 }
3181}
3182
3183/* Workspace method: Doxygen documentation will be auto-generated */
3185 const Index& stokes_dim,
3186 const Vector& f_grid,
3187 const Index& atmosphere_dim,
3188 const Vector& lat_grid,
3189 const Vector& lat_true,
3190 const Vector& lon_true,
3191 const Vector& rtp_pos,
3192 const Vector& rtp_los,
3193 const GriddedField6& r_field,
3194 const Verbosity&) {
3195 // Basic checks and sizes
3196 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3197 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
3198 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
3199 chk_rte_pos(atmosphere_dim, rtp_pos);
3200 chk_rte_los(atmosphere_dim, rtp_los);
3201 r_field.checksize_strict();
3202 chk_griddedfield_gridname(r_field, 0, "Frequency");
3203 chk_griddedfield_gridname(r_field, 1, "Stokes element");
3204 chk_griddedfield_gridname(r_field, 2, "Stokes element");
3205 chk_griddedfield_gridname(r_field, 3, "Incidence angle");
3206 chk_griddedfield_gridname(r_field, 4, "Latitude");
3207 chk_griddedfield_gridname(r_field, 5, "Longitude");
3208 //
3209 const Index nf_in = r_field.data.nvitrines();
3210 const Index ns2 = r_field.data.nshelves();
3211 const Index ns1 = r_field.data.nbooks();
3212 const Index nza = r_field.data.npages();
3213 const Index nlat = r_field.data.nrows();
3214 const Index nlon = r_field.data.ncols();
3215 //
3216 if (nlat < 2 || nlon < 2) {
3217 ostringstream os;
3218 os << "The data in *r_field* must span a geographical region. That is,\n"
3219 << "the latitude and longitude grids must have a length >= 2.";
3220 throw runtime_error(os.str());
3221 }
3222 //
3223 if (nza < 2) {
3224 ostringstream os;
3225 os << "The data in *r_field* must span a range of zenith angles. That\n"
3226 << "is the zenith angle grid must have a length >= 2.";
3227 throw runtime_error(os.str());
3228 }
3229 if (ns1 < stokes_dim || ns2 < stokes_dim || ns1 > 4 || ns2 > 4) {
3230 ostringstream os;
3231 os << "The \"Stokes dimensions\" must have a size that is >= "
3232 << "*stokes_dim* (but not exceeding 4).";
3233 throw runtime_error(os.str());
3234 }
3235
3236 // Determine true geographical position
3237 Vector lat(1), lon(1);
3239 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
3240
3241 // Ensure correct coverage of lon grid
3242 Vector lon_shifted;
3243 lon_shiftgrid(lon_shifted, r_field.get_numeric_grid(5), lon[0]);
3244
3245 // Interpolate in lat and lon
3246 Tensor4 r_f_za(nf_in, stokes_dim, stokes_dim, nza);
3247 {
3249 "Latitude interpolation", r_field.get_numeric_grid(4), lat[0]);
3250 chk_interpolation_grids("Longitude interpolation", lon_shifted, lon[0]);
3251 GridPos gp_lat, gp_lon;
3252 gridpos(gp_lat, r_field.get_numeric_grid(4), lat[0]);
3253 gridpos(gp_lon, lon_shifted, lon[0]);
3254 Vector itw(4);
3255 interpweights(itw, gp_lat, gp_lon);
3256 for (Index iv = 0; iv < nf_in; iv++) {
3257 for (Index iz = 0; iz < nza; iz++) {
3258 for (Index is1 = 0; is1 < stokes_dim; is1++) {
3259 for (Index is2 = 0; is2 < stokes_dim; is2++) {
3260 r_f_za(iv, is1, is2, iz) =
3261 interp(itw,
3262 r_field.data(iv, is1, is2, iz, joker, joker),
3263 gp_lat,
3264 gp_lon);
3265 }
3266 }
3267 }
3268 }
3269 }
3270
3271 // Interpolate in incidence angle, cubic if possible
3272 Tensor3 r_f(nf_in, stokes_dim, stokes_dim);
3273 Index order = 3;
3274 if (nza < 4) {
3275 order = 1;
3276 }
3277 {
3279 "Incidence angle interpolation", r_field.get_numeric_grid(3), 180 - rtp_los[0]);
3280 const LagrangeInterpolation lag(0, 180 - rtp_los[0], r_field.get_numeric_grid(3), order);
3281 const auto itw = interpweights(lag);
3282 //
3283 for (Index i = 0; i < nf_in; i++) {
3284 for (Index is1 = 0; is1 < stokes_dim; is1++) {
3285 for (Index is2 = 0; is2 < stokes_dim; is2++) {
3286 r_f(i, is1, is2) = interp(r_f_za(i, is1, is2, joker), itw, lag);
3287 }
3288 }
3289 }
3290 }
3291
3292 // Extract or interpolate in frequency
3293 //
3294 if (nf_in == 1) {
3295 surface_reflectivity = r_f;
3296 } else {
3298 "Frequency interpolation", r_field.get_numeric_grid(0), f_grid);
3299 const Index nf_out = f_grid.nelem();
3300 surface_reflectivity.resize(nf_out, stokes_dim, stokes_dim);
3301 //
3302 ArrayOfGridPos gp(nf_out);
3303 Matrix itw(nf_out, 2);
3304 gridpos(gp, r_field.get_numeric_grid(0), f_grid);
3305 interpweights(itw, gp);
3306 for (Index is1 = 0; is1 < stokes_dim; is1++) {
3307 for (Index is2 = 0; is2 < stokes_dim; is2++) {
3308 interp(surface_reflectivity(joker, is1, is2),
3309 itw,
3310 r_f(joker, is1, is2),
3311 gp);
3312 }
3313 }
3314 }
3315}
3316
3317/* Workspace method: Doxygen documentation will be auto-generated */
3319 Vector& surface_scalar_reflectivity,
3320 const Index& stokes_dim,
3321 const Vector& f_grid,
3322 const Index& atmosphere_dim,
3323 const Vector& lat_grid,
3324 const Vector& lat_true,
3325 const Vector& lon_true,
3326 const Vector& rtp_pos,
3327 const Vector& rtp_los,
3328 const GriddedField4& r_field,
3329 const Verbosity&) {
3330 // Basic checks and sizes
3331 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3332 chk_if_in_range("stokes_dim", stokes_dim, 1, 1);
3333 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
3334 chk_rte_pos(atmosphere_dim, rtp_pos);
3335 chk_rte_los(atmosphere_dim, rtp_los);
3336 r_field.checksize_strict();
3337 chk_griddedfield_gridname(r_field, 0, "Frequency");
3338 chk_griddedfield_gridname(r_field, 1, "Incidence angle");
3339 chk_griddedfield_gridname(r_field, 2, "Latitude");
3340 chk_griddedfield_gridname(r_field, 3, "Longitude");
3341 //
3342 const Index nf_in = r_field.data.nbooks();
3343 const Index nza = r_field.data.npages();
3344 const Index nlat = r_field.data.nrows();
3345 const Index nlon = r_field.data.ncols();
3346 //
3347 if (nlat < 2 || nlon < 2) {
3348 ostringstream os;
3349 os << "The data in *r_field* must span a geographical region. That is,\n"
3350 << "the latitude and longitude grids must have a length >= 2.";
3351 throw runtime_error(os.str());
3352 }
3353 //
3354 if (nza < 2) {
3355 ostringstream os;
3356 os << "The data in *r_field* must span a range of zenith angles. That\n"
3357 << "is the zenith angle grid must have a length >= 2.";
3358 throw runtime_error(os.str());
3359 }
3360
3361 // Determine true geographical position
3362 Vector lat(1), lon(1);
3364 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
3365
3366 // Ensure correct coverage of lon grid
3367 Vector lon_shifted;
3368 lon_shiftgrid(lon_shifted, r_field.get_numeric_grid(3), lon[0]);
3369
3370 // Interpolate in lat and lon
3371 Matrix r_f_za(nf_in, nza);
3372 {
3374 "Latitude interpolation", r_field.get_numeric_grid(2), lat[0]);
3375 chk_interpolation_grids("Longitude interpolation", lon_shifted, lon[0]);
3376 GridPos gp_lat, gp_lon;
3377 gridpos(gp_lat, r_field.get_numeric_grid(2), lat[0]);
3378 gridpos(gp_lon, lon_shifted, lon[0]);
3379 Vector itw(4);
3380 interpweights(itw, gp_lat, gp_lon);
3381 for (Index iv = 0; iv < nf_in; iv++) {
3382 for (Index iz = 0; iz < nza; iz++) {
3383 r_f_za(iv, iz) =
3384 interp(itw, r_field.data(iv, iz, joker, joker), gp_lat, gp_lon);
3385 }
3386 }
3387 }
3388
3389 // Interpolate in incidence angle, cubic if possible
3390 Vector r_f(nf_in);
3391 Index order = 3;
3392 if (nza < 4) {
3393 order = 1;
3394 }
3395 {
3397 "Incidence angle interpolation", r_field.get_numeric_grid(1), 180 - rtp_los[0]);
3398 const LagrangeInterpolation lag(0, 180 - rtp_los[0], r_field.get_numeric_grid(1), order);
3399 const auto itw=interpweights(lag);
3400 for (Index i = 0; i < nf_in; i++) {
3401 r_f[i] = interp(r_f_za(i, joker), itw, lag);
3402 }
3403 }
3404
3405 // Extract or interpolate in frequency
3406 //
3407 if (nf_in == 1) {
3408 surface_scalar_reflectivity.resize(1);
3409 surface_scalar_reflectivity[0] = r_f[0];
3410 } else {
3412 "Frequency interpolation", r_field.get_numeric_grid(0), f_grid);
3413 const Index nf_out = f_grid.nelem();
3414 surface_scalar_reflectivity.resize(nf_out);
3415 //
3416 ArrayOfGridPos gp(nf_out);
3417 Matrix itw(nf_out, 2);
3418 gridpos(gp, r_field.get_numeric_grid(0), f_grid);
3419 interpweights(itw, gp);
3420 interp(surface_scalar_reflectivity, itw, r_f, gp);
3421 }
3422}
3423
3424/* Workspace method: Doxygen documentation will be auto-generated */
3426 Vector& surface_scalar_reflectivity,
3427 const Tensor4& surface_rmatrix,
3428 const Verbosity&) {
3429 const Index nf = surface_rmatrix.npages();
3430 const Index nlos = surface_rmatrix.nbooks();
3431
3432 surface_scalar_reflectivity.resize(nf);
3433 surface_scalar_reflectivity = 0;
3434
3435 for (Index i = 0; i < nf; i++) {
3436 for (Index l = 0; l < nlos; l++) {
3437 surface_scalar_reflectivity[i] += surface_rmatrix(l, i, 0, 0);
3438 }
3439 }
3440}
3441
3442/* Workspace method: Doxygen documentation will be auto-generated */
3444 Vector& surface_types_aux,
3445 Vector& surface_types_weights,
3446 const Index& atmosphere_dim,
3447 const Vector& lat_grid,
3448 const Vector& lat_true,
3449 const Vector& lon_true,
3450 const Vector& rtp_pos,
3451 const GriddedField2& surface_type_mask,
3452 const String& method,
3453 const Verbosity&) {
3454 // Set expected order of grids
3455 Index gfield_latID = 0;
3456 Index gfield_lonID = 1;
3457
3458 // Basic checks and sizes
3459 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3460 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
3461 chk_rte_pos(atmosphere_dim, rtp_pos);
3462 surface_type_mask.checksize_strict();
3463 //
3464 chk_griddedfield_gridname(surface_type_mask, gfield_latID, "Latitude");
3465 chk_griddedfield_gridname(surface_type_mask, gfield_lonID, "Longitude");
3466 //
3467 const Index nlat = surface_type_mask.data.nrows();
3468 const Index nlon = surface_type_mask.data.ncols();
3469 //
3470 if (nlat < 2 || nlon < 2) {
3471 ostringstream os;
3472 os << "The data in *surface_type_mask* must span a geographical "
3473 << "region. Accordingly,\nthe latitude and longitude grids must "
3474 << "both have a length >= 2.";
3475 throw runtime_error(os.str());
3476 }
3477
3478 const Vector& GFlat = surface_type_mask.get_numeric_grid(gfield_latID);
3479 const Vector& GFlon = surface_type_mask.get_numeric_grid(gfield_lonID);
3480
3481 // Determine true geographical position
3482 Vector lat(1), lon(1);
3484 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
3485
3486 // Ensure correct coverage of lon grid
3487 Vector lon_shifted;
3488 lon_shiftgrid(lon_shifted, GFlon, lon[0]);
3489
3490 // Check if lat/lon we need are actually covered
3491 chk_if_in_range("rtp_pos.lat", lat[0], GFlat[0], GFlat[nlat - 1]);
3492 chk_if_in_range("rtp_pos.lon", lon[0], lon_shifted[0], lon_shifted[nlon - 1]);
3493
3494 // Grid positions
3495 GridPos gp_lat, gp_lon;
3496 gridpos(gp_lat, GFlat, lat[0]);
3497 gridpos(gp_lon, lon_shifted, lon[0]);
3498
3499 if (method == "nearest" ) {
3500 // Extract closest point
3501 Index ilat, ilon;
3502 if (gp_lat.fd[0] < 0.5) {
3503 ilat = gp_lat.idx;
3504 } else {
3505 ilat = gp_lat.idx + 1;
3506 }
3507 if (gp_lon.fd[0] < 0.5) {
3508 ilon = gp_lon.idx;
3509 } else {
3510 ilon = gp_lon.idx + 1;
3511 }
3512 //
3513 surface_types.resize(1);
3514 surface_types_aux.resize(1);
3515 surface_types_weights.resize(1);
3516 surface_types[0] = (Index)floor(surface_type_mask.data(ilat, ilon));
3517 surface_types_aux[0] = surface_type_mask.data(ilat, ilon) -
3518 Numeric(surface_types[0]);
3519 surface_types_weights[0] = 1.0;
3520 }
3521
3522 else if (method == "linear" ) {
3523 // Determine types involved
3524 ArrayOfIndex types0(4), types(4);
3525 Index ntypes = 1;
3526 //
3527 types0[0] = (Index)floor(surface_type_mask.data(gp_lat.idx,gp_lon.idx));
3528 types0[1] = (Index)floor(surface_type_mask.data(gp_lat.idx,gp_lon.idx+1));
3529 types0[2] = (Index)floor(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx));
3530 types0[3] = (Index)floor(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx+1));
3531 //
3532 types[0] = types0[0];
3533 for (Index t=1; t<4; t++) {
3534 bool unique = true;
3535 for (Index n=0; n<ntypes && unique; n++)
3536 if (types0[t] == types[n] ) { unique = false; }
3537 if (unique) {
3538 types[ntypes] = types0[t];
3539 ntypes += 1;
3540 }
3541 }
3542 //
3543 surface_types.resize(ntypes);
3544 surface_types_aux.resize(ntypes);
3545 surface_types_weights.resize(ntypes);
3546
3547 // Interpolation weights
3548 Vector itw(4);
3549 interpweights(itw, gp_lat, gp_lon);
3550
3551 // Determine weight for each type, and make an interpolation of aux value
3552 for (Index t=0; t<ntypes; t++ ) {
3553 Numeric wtot = 0, auxtot = 0;
3554 Numeric ntype = (Numeric)types[t];
3555 if (types[t] == types0[0]) { wtot += itw[0];
3556 auxtot += itw[0]*(surface_type_mask.data(gp_lat.idx,gp_lon.idx)-ntype);
3557 };
3558 if (types[t] == types0[1]) { wtot += itw[1];
3559 auxtot += itw[1]*(surface_type_mask.data(gp_lat.idx,gp_lon.idx+1)-ntype);
3560 };
3561 if (types[t] == types0[2]) { wtot += itw[2];
3562 auxtot += itw[2]*(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx)-ntype);
3563 };
3564 if (types[t] == types0[3]) { wtot += itw[3];
3565 auxtot += itw[3]*(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx+1)-ntype);
3566 };
3567 //
3568 surface_types[t] = types[t];
3569 surface_types_aux[t] = auxtot / wtot;
3570 surface_types_weights[t] = wtot;
3571 }
3572 }
3573
3574 else {
3575 ostringstream os;
3576 os << "The allowed options for *method are: \"nearest\" and \"linear\",\n"
3577 << "but you have selected: " << method;
3578 throw runtime_error(os.str());
3579 }
3580}
3581
3582/* Workspace method: Doxygen documentation will be auto-generated */
3584 Numeric& surface_skin_t,
3585 Matrix& surface_los,
3586 Tensor4& surface_rmatrix,
3587 Matrix& surface_emission,
3588 const Vector& f_grid,
3589 const Vector& rtp_pos,
3590 const Vector& rtp_los,
3591 const ArrayOfAgenda& surface_rtprop_agenda_array,
3592 const ArrayOfIndex& surface_types,
3593 const Vector& surface_types_aux,
3594 const Vector& surface_types_weights,
3595 const Verbosity&) {
3596 if (surface_types.nelem() != 1) {
3597 ostringstream os;
3598 os << "This method requires that *surface_types* have length 1.\n"
3599 << "If you are trying to apply a mixture model, you have to use\n"
3600 << "a setup based on *iySurfaceCallAgendaX* instead.";
3601 throw runtime_error(os.str());
3602 }
3603 if (surface_types[0] < 0)
3604 throw runtime_error("*surface_type* is not allowed to be negative.");
3605 if (surface_types[0] >= surface_rtprop_agenda_array.nelem()) {
3606 ostringstream os;
3607 os << "*surface_rtprop_agenda_array* has only "
3608 << surface_rtprop_agenda_array.nelem()
3609 << " elements,\n while you have selected element " << surface_types[0];
3610 throw runtime_error(os.str());
3611 }
3612 if (abs(surface_types_weights[0]-1)>1e-4)
3613 throw runtime_error("Sum of *surface_types_weights* deviates from 1.");
3614
3616 surface_skin_t,
3617 surface_emission,
3618 surface_los,
3619 surface_rmatrix,
3620 surface_types[0],
3621 f_grid,
3622 rtp_pos,
3623 rtp_los,
3624 surface_types_aux[0],
3625 surface_rtprop_agenda_array);
3626}
3627
3628/* Workspace method: Doxygen documentation will be auto-generated */
3629void SurfaceBlackbody(Matrix& surface_los,
3630 Tensor4& surface_rmatrix,
3631 ArrayOfTensor4& dsurface_rmatrix_dx,
3632 Matrix& surface_emission,
3633 ArrayOfMatrix& dsurface_emission_dx,
3634 const Index& stokes_dim,
3635 const Index& atmosphere_dim,
3636 const Vector& lat_grid,
3637 const Vector& lon_grid,
3638 const Vector& f_grid,
3639 const Vector& rtp_pos,
3640 const Vector& rtp_los,
3641 const Tensor3& surface_props_data,
3642 const ArrayOfString& surface_props_names,
3643 const ArrayOfString& dsurface_names,
3644 const Index& jacobian_do,
3645 const Verbosity& verbosity) {
3646 // Check surface_data
3647 surface_props_check(atmosphere_dim,
3648 lat_grid,
3649 lon_grid,
3650 surface_props_data,
3651 surface_props_names);
3652
3653 // Interplation grid positions and weights
3654 ArrayOfGridPos gp_lat(1), gp_lon(1);
3655 Matrix itw;
3657 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
3658 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
3659
3660 // Skin temperature
3661 Vector skin_t(1);
3662 surface_props_interp(skin_t,
3663 "Skin temperature",
3664 atmosphere_dim,
3665 gp_lat,
3666 gp_lon,
3667 itw,
3668 surface_props_data,
3669 surface_props_names);
3670
3671 surfaceBlackbody(surface_los,
3672 surface_rmatrix,
3673 surface_emission,
3674 atmosphere_dim,
3675 f_grid,
3676 stokes_dim,
3677 rtp_pos,
3678 rtp_los,
3679 skin_t[0],
3680 verbosity);
3681
3682 surface_rmatrix.resize(1, f_grid.nelem(), stokes_dim, stokes_dim);
3683 surface_rmatrix = 0.0;
3684
3685 // Jacobian part
3686 if (jacobian_do) {
3687 dsurface_check(surface_props_names,
3688 dsurface_names,
3689 dsurface_rmatrix_dx,
3690 dsurface_emission_dx);
3691
3692 Index irq;
3693
3694 // Skin temperature
3695 irq = find_first(dsurface_names, String("Skin temperature"));
3696 if (irq >= 0) {
3697 Matrix dbdt(f_grid.nelem(), 1);
3698 dplanck_dt(dbdt(joker, 0), f_grid, skin_t[0]);
3699 dsurface_emission_dx[irq] = dbdt;
3700 dsurface_rmatrix_dx[irq].resize(surface_rmatrix.nbooks(),
3701 surface_rmatrix.npages(),
3702 surface_rmatrix.nrows(),
3703 surface_rmatrix.ncols());
3704 dsurface_rmatrix_dx[irq] = 0;
3705 }
3706 }
3707}
3708
3709/* Workspace method: Doxygen documentation will be auto-generated */
3710void SurfaceDummy(ArrayOfTensor4& dsurface_rmatrix_dx,
3711 ArrayOfMatrix& dsurface_emission_dx,
3712 const Index& atmosphere_dim,
3713 const Vector& lat_grid,
3714 const Vector& lon_grid,
3715 const Tensor3& surface_props_data,
3716 const ArrayOfString& surface_props_names,
3717 const ArrayOfString& dsurface_names,
3718 const Index& jacobian_do,
3719 const Verbosity&) {
3720 if (surface_props_names.nelem())
3721 throw runtime_error(
3722 "When calling this method, *surface_props_names* should be empty.");
3723
3724 surface_props_check(atmosphere_dim,
3725 lat_grid,
3726 lon_grid,
3727 surface_props_data,
3728 surface_props_names);
3729
3730 if (jacobian_do) {
3731 dsurface_check(surface_props_names,
3732 dsurface_names,
3733 dsurface_rmatrix_dx,
3734 dsurface_emission_dx);
3735 }
3736}
3737
3738/* Workspace method: Doxygen documentation will be auto-generated */
3739void SurfaceFastem(Matrix& surface_los,
3740 Tensor4& surface_rmatrix,
3741 ArrayOfTensor4& dsurface_rmatrix_dx,
3742 Matrix& surface_emission,
3743 ArrayOfMatrix& dsurface_emission_dx,
3744 const Index& stokes_dim,
3745 const Index& atmosphere_dim,
3746 const Vector& lat_grid,
3747 const Vector& lon_grid,
3748 const Vector& f_grid,
3749 const Vector& rtp_pos,
3750 const Vector& rtp_los,
3751 const Tensor3& surface_props_data,
3752 const ArrayOfString& surface_props_names,
3753 const ArrayOfString& dsurface_names,
3754 const Index& jacobian_do,
3755 const Vector& transmittance,
3756 const Index& fastem_version,
3757 const Verbosity& verbosity) {
3758 // Check surface_data
3759 surface_props_check(atmosphere_dim,
3760 lat_grid,
3761 lon_grid,
3762 surface_props_data,
3763 surface_props_names);
3764
3765 // Interplation grid positions and weights
3766 ArrayOfGridPos gp_lat(1), gp_lon(1);
3767 Matrix itw;
3769 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
3770 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
3771
3772 // Skin temperature
3773 Vector skin_t(1);
3774 surface_props_interp(skin_t,
3775 "Water skin temperature",
3776 atmosphere_dim,
3777 gp_lat,
3778 gp_lon,
3779 itw,
3780 surface_props_data,
3781 surface_props_names);
3782
3783 // Wind speed
3784 Vector wind_speed(1);
3785 surface_props_interp(wind_speed,
3786 "Wind speed",
3787 atmosphere_dim,
3788 gp_lat,
3789 gp_lon,
3790 itw,
3791 surface_props_data,
3792 surface_props_names);
3793
3794 // Wind direction
3795 Vector wind_direction(1);
3796 surface_props_interp(wind_direction,
3797 "Wind direction",
3798 atmosphere_dim,
3799 gp_lat,
3800 gp_lon,
3801 itw,
3802 surface_props_data,
3803 surface_props_names);
3804
3805 // Salinity
3806 Vector salinity(1);
3807 surface_props_interp(salinity,
3808 "Salinity",
3809 atmosphere_dim,
3810 gp_lat,
3811 gp_lon,
3812 itw,
3813 surface_props_data,
3814 surface_props_names);
3815
3816 // Call FASTEM
3817 surfaceFastem(surface_los,
3818 surface_rmatrix,
3819 surface_emission,
3820 atmosphere_dim,
3821 stokes_dim,
3822 f_grid,
3823 rtp_pos,
3824 rtp_los,
3825 skin_t[0],
3826 salinity[0],
3827 wind_speed[0],
3828 wind_direction[0],
3829 transmittance,
3830 fastem_version,
3831 verbosity);
3832
3833 // Jacobian part
3834 if (jacobian_do) {
3835 dsurface_check(surface_props_names,
3836 dsurface_names,
3837 dsurface_rmatrix_dx,
3838 dsurface_emission_dx);
3839
3840 Index irq;
3841
3842 // Skin temperature
3843 irq = find_first(dsurface_names, String("Water skin temperature"));
3844 if (irq >= 0) {
3845 const Numeric dd = 0.1;
3846 Matrix surface_los2;
3847 surfaceFastem(surface_los2,
3848 dsurface_rmatrix_dx[irq],
3849 dsurface_emission_dx[irq],
3850 atmosphere_dim,
3851 stokes_dim,
3852 f_grid,
3853 rtp_pos,
3854 rtp_los,
3855 skin_t[0] + dd,
3856 salinity[0],
3857 wind_speed[0],
3858 wind_direction[0],
3859 transmittance,
3860 fastem_version,
3861 verbosity);
3862 //
3863 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3864 dsurface_rmatrix_dx[irq] /= dd;
3865 dsurface_emission_dx[irq] -= surface_emission;
3866 dsurface_emission_dx[irq] /= dd;
3867 }
3868
3869 // Wind speed
3870 irq = find_first(dsurface_names, String("Wind speed"));
3871 if (irq >= 0) {
3872 const Numeric dd = 0.1;
3873 Matrix surface_los2;
3874 surfaceFastem(surface_los2,
3875 dsurface_rmatrix_dx[irq],
3876 dsurface_emission_dx[irq],
3877 atmosphere_dim,
3878 stokes_dim,
3879 f_grid,
3880 rtp_pos,
3881 rtp_los,
3882 skin_t[0],
3883 salinity[0],
3884 wind_speed[0] + dd,
3885 wind_direction[0],
3886 transmittance,
3887 fastem_version,
3888 verbosity);
3889 //
3890 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3891 dsurface_rmatrix_dx[irq] /= dd;
3892 dsurface_emission_dx[irq] -= surface_emission;
3893 dsurface_emission_dx[irq] /= dd;
3894 }
3895
3896 // Wind direction
3897 irq = find_first(dsurface_names, String("Wind direction"));
3898 if (irq >= 0) {
3899 const Numeric dd = 1;
3900 Matrix surface_los2;
3901 surfaceFastem(surface_los2,
3902 dsurface_rmatrix_dx[irq],
3903 dsurface_emission_dx[irq],
3904 atmosphere_dim,
3905 stokes_dim,
3906 f_grid,
3907 rtp_pos,
3908 rtp_los,
3909 skin_t[0],
3910 salinity[0],
3911 wind_speed[0],
3912 wind_direction[0] + dd,
3913 transmittance,
3914 fastem_version,
3915 verbosity);
3916 //
3917 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3918 dsurface_rmatrix_dx[irq] /= dd;
3919 dsurface_emission_dx[irq] -= surface_emission;
3920 dsurface_emission_dx[irq] /= dd;
3921 }
3922
3923 // Salinity
3924 irq = find_first(dsurface_names, String("Salinity"));
3925 if (irq >= 0) {
3926 const Numeric dd = 0.0005;
3927 Matrix surface_los2;
3928 surfaceFastem(surface_los2,
3929 dsurface_rmatrix_dx[irq],
3930 dsurface_emission_dx[irq],
3931 atmosphere_dim,
3932 stokes_dim,
3933 f_grid,
3934 rtp_pos,
3935 rtp_los,
3936 skin_t[0],
3937 salinity[0] + dd,
3938 wind_speed[0],
3939 wind_direction[0],
3940 transmittance,
3941 fastem_version,
3942 verbosity);
3943 //
3944 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3945 dsurface_rmatrix_dx[irq] /= dd;
3946 dsurface_emission_dx[irq] -= surface_emission;
3947 dsurface_emission_dx[irq] /= dd;
3948 }
3949 }
3950}
3951
3952
3953/* Workspace method: Doxygen documentation will be auto-generated */
3955 Tensor4& surface_rmatrix,
3956 ArrayOfTensor4& dsurface_rmatrix_dx,
3957 Matrix& surface_emission,
3958 ArrayOfMatrix& dsurface_emission_dx,
3959 const Index& stokes_dim,
3960 const Index& atmosphere_dim,
3961 const Vector& lat_grid,
3962 const Vector& lon_grid,
3963 const Vector& f_grid,
3964 const Vector& rtp_pos,
3965 const Vector& rtp_los,
3966 const Vector& specular_los,
3967 const Tensor3& surface_props_data,
3968 const ArrayOfString& surface_props_names,
3969 const ArrayOfString& dsurface_names,
3970 const Index& jacobian_do,
3971 const Vector& f_reflectivities,
3972 const Verbosity& verbosity) {
3973 // Check surface_data
3974 surface_props_check(atmosphere_dim,
3975 lat_grid,
3976 lon_grid,
3977 surface_props_data,
3978 surface_props_names);
3979
3980 // Check of GINs
3981 chk_if_increasing("GIN *f_reflectivities*", f_reflectivities);
3982
3983 // Sizes
3984 const Index nr = f_reflectivities.nelem();
3985 const Index nf = f_grid.nelem();
3986
3987 // Interplation grid positions and weights
3988 ArrayOfGridPos gp_lat(1), gp_lon(1);
3989 Matrix itw;
3991 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
3992 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
3993
3994 // Skin temperature
3995 Vector surface_skin_t(1);
3996 surface_props_interp(surface_skin_t,
3997 "Skin temperature",
3998 atmosphere_dim,
3999 gp_lat,
4000 gp_lon,
4001 itw,
4002 surface_props_data,
4003 surface_props_names);
4004
4005 // Reflectivities
4006 Vector reflectivities(nr);
4007 //
4008 for (Index i=0; i<nr; i++) {
4009 Vector rv(1);
4010 ostringstream sstr;
4011 sstr << "Scalar reflectivity " << i;
4013 sstr.str(),
4014 atmosphere_dim,
4015 gp_lat,
4016 gp_lon,
4017 itw,
4018 surface_props_data,
4019 surface_props_names);
4020 reflectivities[i] = rv[0];
4021 }
4022
4023 // Create surface_scalar_reflectivity
4024 Vector surface_scalar_reflectivity(nf);
4025 ArrayOfGridPos gp(nf);
4026 if (nr==1) {
4027 gp4length1grid(gp);
4028 } else {
4029 gridpos(gp, f_reflectivities, f_grid, 1e99);
4031 }
4032 Matrix itw2(nf, 2);
4033 interpweights(itw2, gp);
4034 interp( surface_scalar_reflectivity, itw2, reflectivities, gp);
4035
4036 // Call non-Jacobian method
4038 surface_rmatrix,
4039 surface_emission,
4040 f_grid,
4041 stokes_dim,
4042 atmosphere_dim,
4043 rtp_pos,
4044 rtp_los,
4045 specular_los,
4046 surface_skin_t[0],
4047 surface_scalar_reflectivity,
4048 verbosity);
4049
4050 // Jacobian part
4051 if (jacobian_do) {
4052 dsurface_check(surface_props_names,
4053 dsurface_names,
4054 dsurface_rmatrix_dx,
4055 dsurface_emission_dx);
4056
4057 Index irq;
4058
4059 // Skin temperature
4060 irq = find_first(dsurface_names, String("Skin temperature"));
4061 if (irq >= 0) {
4062 const Numeric dd = 0.1;
4063 Matrix surface_los2;
4064 surfaceFlatScalarReflectivity(surface_los2,
4065 dsurface_rmatrix_dx[irq],
4066 dsurface_emission_dx[irq],
4067 f_grid,
4068 stokes_dim,
4069 atmosphere_dim,
4070 rtp_pos,
4071 rtp_los,
4072 specular_los,
4073 surface_skin_t[0]+dd,
4074 surface_scalar_reflectivity,
4075 verbosity);
4076 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
4077 dsurface_rmatrix_dx[irq] /= dd;
4078 dsurface_emission_dx[irq] -= surface_emission;
4079 dsurface_emission_dx[irq] /= dd;
4080 }
4081
4082 // Reflectivities
4083 for (Index i=0; i<nr; i++) {
4084 ostringstream sstr;
4085 sstr << "Scalar reflectivity " << i;
4086 irq = find_first(dsurface_names, String(sstr.str()));
4087 if (irq >= 0) {
4088 const Numeric dd = 1e-4;
4089 Vector rpert = reflectivities;
4090 rpert[i] += dd;
4091 interp( surface_scalar_reflectivity, itw2, rpert, gp);
4092 //
4093 Matrix surface_los2;
4094 surfaceFlatScalarReflectivity(surface_los2,
4095 dsurface_rmatrix_dx[irq],
4096 dsurface_emission_dx[irq],
4097 f_grid,
4098 stokes_dim,
4099 atmosphere_dim,
4100 rtp_pos,
4101 rtp_los,
4102 specular_los,
4103 surface_skin_t[0],
4104 surface_scalar_reflectivity,
4105 verbosity);
4106 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
4107 dsurface_rmatrix_dx[irq] /= dd;
4108 dsurface_emission_dx[irq] -= surface_emission;
4109 dsurface_emission_dx[irq] /= dd;
4110 }
4111 }
4112 }
4113}
4114
4115/* Workspace method: Doxygen documentation will be auto-generated */
4116void SurfaceTessem(Matrix& surface_los,
4117 Tensor4& surface_rmatrix,
4118 ArrayOfTensor4& dsurface_rmatrix_dx,
4119 Matrix& surface_emission,
4120 ArrayOfMatrix& dsurface_emission_dx,
4121 const Index& stokes_dim,
4122 const Index& atmosphere_dim,
4123 const Vector& lat_grid,
4124 const Vector& lon_grid,
4125 const Vector& f_grid,
4126 const Vector& rtp_pos,
4127 const Vector& rtp_los,
4128 const TessemNN& net_h,
4129 const TessemNN& net_v,
4130 const Tensor3& surface_props_data,
4131 const ArrayOfString& surface_props_names,
4132 const ArrayOfString& dsurface_names,
4133 const Index& jacobian_do,
4134 const Verbosity& verbosity) {
4135 // Check surface_data
4136 surface_props_check(atmosphere_dim,
4137 lat_grid,
4138 lon_grid,
4139 surface_props_data,
4140 surface_props_names);
4141
4142 // Interplation grid positions and weights
4143 ArrayOfGridPos gp_lat(1), gp_lon(1);
4144 Matrix itw;
4146 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
4147 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
4148
4149 // Skin temperature
4150 Vector skin_t(1);
4151 surface_props_interp(skin_t,
4152 "Water skin temperature",
4153 atmosphere_dim,
4154 gp_lat,
4155 gp_lon,
4156 itw,
4157 surface_props_data,
4158 surface_props_names);
4159
4160 // Wind speed
4161 Vector wind_speed(1);
4162 surface_props_interp(wind_speed,
4163 "Wind speed",
4164 atmosphere_dim,
4165 gp_lat,
4166 gp_lon,
4167 itw,
4168 surface_props_data,
4169 surface_props_names);
4170
4171 // Salinity
4172 Vector salinity(1);
4173 surface_props_interp(salinity,
4174 "Salinity",
4175 atmosphere_dim,
4176 gp_lat,
4177 gp_lon,
4178 itw,
4179 surface_props_data,
4180 surface_props_names);
4181
4182 // Call TESSEM
4183 surfaceTessem(surface_los,
4184 surface_rmatrix,
4185 surface_emission,
4186 atmosphere_dim,
4187 stokes_dim,
4188 f_grid,
4189 rtp_pos,
4190 rtp_los,
4191 skin_t[0],
4192 net_h,
4193 net_v,
4194 salinity[0],
4195 wind_speed[0],
4196 verbosity);
4197
4198 // Jacobian part
4199 if (jacobian_do) {
4200 dsurface_check(surface_props_names,
4201 dsurface_names,
4202 dsurface_rmatrix_dx,
4203 dsurface_emission_dx);
4204
4205 Index irq;
4206
4207 // Skin temperature
4208 irq = find_first(dsurface_names, String("Water skin temperature"));
4209 if (irq >= 0) {
4210 const Numeric dd = 0.1;
4211 Matrix surface_los2;
4212 surfaceTessem(surface_los2,
4213 dsurface_rmatrix_dx[irq],
4214 dsurface_emission_dx[irq],
4215 atmosphere_dim,
4216 stokes_dim,
4217 f_grid,
4218 rtp_pos,
4219 rtp_los,
4220 skin_t[0] + dd,
4221 net_h,
4222 net_v,
4223 salinity[0],
4224 wind_speed[0],
4225 verbosity);
4226 //
4227 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
4228 dsurface_rmatrix_dx[irq] /= dd;
4229 dsurface_emission_dx[irq] -= surface_emission;
4230 dsurface_emission_dx[irq] /= dd;
4231 }
4232
4233 // Wind speed
4234 irq = find_first(dsurface_names, String("Wind speed"));
4235 if (irq >= 0) {
4236 const Numeric dd = 0.1;
4237 Matrix surface_los2;
4238 surfaceTessem(surface_los2,
4239 dsurface_rmatrix_dx[irq],
4240 dsurface_emission_dx[irq],
4241 atmosphere_dim,
4242 stokes_dim,
4243 f_grid,
4244 rtp_pos,
4245 rtp_los,
4246 skin_t[0],
4247 net_h,
4248 net_v,
4249 salinity[0],
4250 wind_speed[0] + dd,
4251 verbosity);
4252 //
4253 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
4254 dsurface_rmatrix_dx[irq] /= dd;
4255 dsurface_emission_dx[irq] -= surface_emission;
4256 dsurface_emission_dx[irq] /= dd;
4257 }
4258
4259 // Salinity
4260 irq = find_first(dsurface_names, String("Salinity"));
4261 if (irq >= 0) {
4262 const Numeric dd = 0.0005;
4263 Matrix surface_los2;
4264 surfaceTessem(surface_los2,
4265 dsurface_rmatrix_dx[irq],
4266 dsurface_emission_dx[irq],
4267 atmosphere_dim,
4268 stokes_dim,
4269 f_grid,
4270 rtp_pos,
4271 rtp_los,
4272 skin_t[0],
4273 net_h,
4274 net_v,
4275 salinity[0] + dd,
4276 wind_speed[0],
4277 verbosity);
4278 //
4279 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
4280 dsurface_rmatrix_dx[irq] /= dd;
4281 dsurface_emission_dx[irq] -= surface_emission;
4282 dsurface_emission_dx[irq] /= dd;
4283 }
4284 }
4285}
4286
4287/* Workspace method: Doxygen documentation will be auto-generated */
4289 const ArrayOfString& iy_aux_vars,
4290 const ArrayOfMatrix& iy_aux,
4291 const Verbosity&) {
4292 Index ihit = -1;
4293
4294 for (Index i = 0; i < iy_aux_vars.nelem(); i++) {
4295 if (iy_aux_vars[i] == "Optical depth") {
4296 ihit = i;
4297 break;
4298 }
4299 }
4300
4301 if (ihit < 0)
4302 throw runtime_error("No element in *iy_aux* holds optical depths.");
4303
4304 const Index n = iy_aux[ihit].nrows();
4305
4306 transmittance.resize(n);
4307
4308 for (Index i = 0; i < n; i++) {
4309 transmittance[i] = exp(-iy_aux[ihit](i, 0));
4310 }
4311}
This file contains the definition of Array.
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:190
The global header file for ARTS.
Constants of physical expressions as constexpr.
Common ARTS conversions.
void surface_rtprop_agenda_arrayExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Index agenda_array_index, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Numeric surface_type_aux, const ArrayOfAgenda &input_agenda_array)
Definition: auto_md.cc:25881
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, Vector &geo_pos, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const String &iy_unit, const Index cloudbox_on, const Index jacobian_do, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition: auto_md.cc:25104
void surface_rtprop_sub_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:25937
void iy_surface_agenda_arrayExecute(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Index agenda_array_index, const String &iy_unit, const Tensor3 &iy_transmittance, const Index iy_id, const Index cloudbox_on, const Index jacobian_do, const Agenda &iy_main_agenda, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const Numeric surface_type_aux, const ArrayOfAgenda &input_agenda_array)
Definition: auto_md.cc:25314
void surface_rtprop_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:25839
void chk_atm_surface(const String &x_name, const Matrix &x, const Index &dim, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_surface
void chk_atm_grids(const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_grids
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:906
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
void chk_rte_los(const Index &atmosphere_dim, ConstVectorView rte_los)
chk_rte_los
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:421
void chk_griddedfield_gridname(const GriddedField &gf, const Index gridindex, const String &gridname)
Check name of grid in GriddedField.
void chk_latlon_true(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true)
chk_latlon_true
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:86
void chk_if_in_range_exclude(const String &x_name, const Numeric &x, const Numeric &x_low, const Numeric &x_high)
chk_if_in_range_exclude
Definition: check_input.cc:230
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:111
void chk_if_in_range_exclude_high(const String &x_name, const Numeric &x, const Numeric &x_low, const Numeric &x_high)
chk_if_in_range_exclude_high
Definition: check_input.cc:205
void chk_vector_length(const String &x_name, ConstVectorView x, const Index &l)
chk_vector_length
Definition: check_input.cc:257
void chk_not_negative(const String &x_name, const Numeric &x)
chk_not_negative
Definition: check_input.cc:134
The Agenda class.
Definition: agenda_class.h:69
String name() const
Agenda name.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
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 ncols() const noexcept
Definition: matpackIV.h:142
Index nrows() const noexcept
Definition: matpackIV.h:141
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index npages() const noexcept
Definition: matpackIV.h:140
Index nrows() const noexcept
Definition: matpackV.h:152
Index ncols() const noexcept
Definition: matpackV.h:153
Index npages() const noexcept
Definition: matpackV.h:151
Index nbooks() const noexcept
Definition: matpackV.h:150
Index nshelves() const noexcept
Definition: matpackV.h:149
Index nbooks() const noexcept
Definition: matpackVI.h:157
Index nvitrines() const noexcept
Definition: matpackVI.h:155
Index ncols() const noexcept
Definition: matpackVI.h:160
Index npages() const noexcept
Definition: matpackVI.h:158
Index nshelves() const noexcept
Definition: matpackVI.h:156
Index nrows() const noexcept
Definition: matpackVI.h:159
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
void checksize_strict() const final
Strict consistency check.
void resize(const GriddedField3 &gf)
Make this GriddedField3 the same size as the given one.
void checksize_strict() const final
Strict consistency check.
void checksize_strict() const final
Strict consistency check.
void checksize_strict() const final
Strict consistency check.
void set_grid_name(Index i, const String &s)
Set grid name.
void set_grid(Index i, const Vector &g)
Set a numeric grid.
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
The range class.
Definition: matpackI.h:159
A telsem atlas.
Definition: telsem.h:59
Index calc_cellnum(Numeric lat, Numeric lon) const
Definition: telsem.cc:144
bool contains(Index cellnumber) const
Definition: telsem.h:85
Index get_class2(Index cellnumber) const
Definition: telsem.h:120
Index calc_cellnum_nearest_neighbor(Numeric lat, Numeric lon) const
Definition: telsem.cc:172
Vector get_emis_h(Index cellnum) const
Definition: telsem.h:159
std::pair< Numeric, Numeric > get_coordinates(Index cellnum) const
Definition: telsem.cc:227
Index get_class1(Index cellnumber) const
Definition: telsem.h:102
std::pair< Numeric, Numeric > emis_interp(Numeric theta, Numeric freq, Index class1, Index class2, const ConstVectorView &ev, const ConstVectorView &eh) const
Definition: telsem.cc:289
Vector get_emis_v(Index i) const
Definition: telsem.h:137
The Tensor3 class.
Definition: matpackIII.h:346
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:429
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
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_IF(condition,...)
Definition: debug.h:134
void fastem(Vector &emissivity, Vector &reflectivity, const Numeric frequency, const Numeric za, const Numeric temperature, const Numeric salinity, const Numeric wind_speed, const Numeric transmittance, const Numeric rel_azimuth, const Index fastem_version)
Calculate the surface emissivity using FASTEM.
Definition: fastem.cc:109
This file contains functions that are adapted from FASTEM code which is used to calculate surface emi...
Header file for functions related to gas scattering.
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1336
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
Definition: geodetic.cc:1409
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1305
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
void diy_from_pos_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstMatrixView diy_dpos, const Index &atmosphere_dim, ConstVectorView rtp_pos)
diy_from_pos_to_rgrids
Definition: jacobian.cc:475
void jacobian_type_extrapol(ArrayOfGridPos &gp)
Adopts grid positions to extrapolation used for jacobians.
Definition: jacobian.cc:836
#define abs(x)
#define q
#define min(a, b)
#define max(a, b)
void AngularGridsSetFluxCalc(Vector &za_grid, Vector &aa_grid, Vector &za_grid_weights, const Index &N_za_grid, const Index &N_aa_grid, const String &za_grid_type, const Verbosity &)
WORKSPACE METHOD: AngularGridsSetFluxCalc.
Definition: m_fluxes.cc:60
void iySurfaceRtpropAgenda(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, Numeric &surface_skin_t, Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const Index &stars_do, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &surface_rtprop_agenda, const Verbosity &)
WORKSPACE METHOD: iySurfaceRtpropAgenda.
Definition: m_surface.cc:1572
void surfaceTelsem(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &f_grid, const Vector &lat_grid, const Vector &lat_true, const Vector &lon_true, const Vector &rtp_pos, const Vector &rtp_los, const Vector &specular_los, const Numeric &surface_skin_t, const TelsemAtlas &atlas, const Numeric &r_min, const Numeric &r_max, const Numeric &d_max, const Verbosity &verbosity)
WORKSPACE METHOD: surfaceTelsem.
Definition: m_surface.cc:2244
void surfaceSplitSpecularTo3beams(Matrix &surface_los, Tensor4 &surface_rmatrix, const Index &atmosphere_dim, const Vector &rtp_pos, const Vector &rtp_los, const Numeric &specular_factor, const Numeric &dza, const Verbosity &)
WORKSPACE METHOD: surfaceSplitSpecularTo3beams.
Definition: m_surface.cc:2979
void iySurfaceLambertianDirect(Workspace &ws, Matrix &iy, const Vector &rtp_pos, 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 &surface_scalar_reflectivity, 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 &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stars_do, const Index &gas_scattering_do, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfStar &stars, const Numeric &rte_alonglos_v, const String &iy_unit, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceLambertianDirect.
Definition: m_surface.cc:1387
void surface_scalar_reflectivityFromGriddedField4(Vector &surface_scalar_reflectivity, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lat_true, const Vector &lon_true, const Vector &rtp_pos, const Vector &rtp_los, const GriddedField4 &r_field, const Verbosity &)
WORKSPACE METHOD: surface_scalar_reflectivityFromGriddedField4.
Definition: m_surface.cc:3318
void surface_complex_refr_indexFromGriddedField5(GriddedField3 &surface_complex_refr_index, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lat_true, const Vector &lon_true, const Vector &rtp_pos, const GriddedField5 &complex_n_field, const Verbosity &)
WORKSPACE METHOD: surface_complex_refr_indexFromGriddedField5.
Definition: m_surface.cc:3088
void InterpSurfaceFieldToPosition(Numeric &outvalue, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &rtp_pos, const Matrix &z_surface, const Matrix &field, const Verbosity &verbosity)
WORKSPACE METHOD: InterpSurfaceFieldToPosition.
Definition: m_surface.cc:206
void iySurfaceFlatRefractiveIndex(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, ArrayOfTensor4 &dsurface_rmatrix_dx, ArrayOfMatrix &dsurface_emission_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const Index &stars_do, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Vector &lat_grid, const Vector &lon_grid, const Matrix &z_surface, const Vector &refellipsoid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const String &iy_unit, const GriddedField3 &surface_complex_refr_index, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &iy_main_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceFlatRefractiveIndex.
Definition: m_surface.cc:764
void surfaceTessem(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Numeric &surface_skin_t, const TessemNN &net_h, const TessemNN &net_v, const Numeric &salinity, const Numeric &wind_speed, const Verbosity &verbosity)
WORKSPACE METHOD: surfaceTessem.
Definition: m_surface.cc:2351
void surface_rtpropCallAgendaX(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const ArrayOfAgenda &surface_rtprop_agenda_array, const ArrayOfIndex &surface_types, const Vector &surface_types_aux, const Vector &surface_types_weights, const Verbosity &)
WORKSPACE METHOD: surface_rtpropCallAgendaX.
Definition: m_surface.cc:3583
constexpr Numeric EARTH_RADIUS
Definition: m_surface.cc:64
void surfaceSemiSpecularBy3beams(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Index &atmosphere_dim, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &surface_rtprop_sub_agenda, const Numeric &specular_factor, const Numeric &dza, const Verbosity &)
WORKSPACE METHOD: surfaceSemiSpecularBy3beams.
Definition: m_surface.cc:2832
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
constexpr Numeric DEG2RAD
Definition: m_surface.cc:65
void surfaceFastem(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Numeric &surface_skin_t, const Numeric &salinity, const Numeric &wind_speed, const Numeric &wind_direction, const Vector &transmittance, const Index &fastem_version, const Verbosity &verbosity)
WORKSPACE METHOD: surfaceFastem.
Definition: m_surface.cc:2074
void surfaceMapToLinearPolarisation(Matrix &surface_emission, Tensor4 &surface_rmatrix, const Index &stokes_dim, const Numeric &pol_angle, const Verbosity &)
WORKSPACE METHOD: surfaceMapToLinearPolarisation.
Definition: m_surface.cc:2183
void surfaceFlatRefractiveIndex(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &rtp_pos, const Vector &rtp_los, const Vector &specular_los, const Numeric &surface_skin_t, const GriddedField3 &surface_complex_refr_index, const Verbosity &verbosity)
WORKSPACE METHOD: surfaceFlatRefractiveIndex.
Definition: m_surface.cc:2430
void SurfaceBlackbody(Matrix &surface_los, Tensor4 &surface_rmatrix, ArrayOfTensor4 &dsurface_rmatrix_dx, Matrix &surface_emission, ArrayOfMatrix &dsurface_emission_dx, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const Index &jacobian_do, const Verbosity &verbosity)
WORKSPACE METHOD: SurfaceBlackbody.
Definition: m_surface.cc:3629
void specular_losCalcNoTopography(Vector &specular_los, Vector &surface_normal, const Vector &rtp_pos, const Vector &rtp_los, const Index &atmosphere_dim, const Verbosity &)
WORKSPACE METHOD: specular_losCalcNoTopography.
Definition: m_surface.cc:1889
void surfaceBlackbody(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Index &atmosphere_dim, const Vector &f_grid, const Index &stokes_dim, const Vector &rtp_pos, const Vector &rtp_los, const Numeric &surface_skin_t, const Verbosity &verbosity)
WORKSPACE METHOD: surfaceBlackbody.
Definition: m_surface.cc:2035
void iySurfaceFlatRefractiveIndexDirect(Workspace &ws, Matrix &iy, 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 GriddedField3 &surface_complex_refr_index, 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 &stars_do, const Index &gas_scattering_do, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfStar &stars, const Numeric &rte_alonglos_v, const String &iy_unit, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceFlatRefractiveIndexDirect.
Definition: m_surface.cc:937
void iySurfaceLambertian(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const Index &stars_do, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Vector &lat_grid, const Vector &lon_grid, const Matrix &z_surface, const Vector &refellipsoid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const String &iy_unit, const Vector &surface_scalar_reflectivity, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &iy_main_agenda, const Index &N_za, const Index &N_aa, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceLambertian.
Definition: m_surface.cc:1085
void SurfaceFastem(Matrix &surface_los, Tensor4 &surface_rmatrix, ArrayOfTensor4 &dsurface_rmatrix_dx, Matrix &surface_emission, ArrayOfMatrix &dsurface_emission_dx, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const Index &jacobian_do, const Vector &transmittance, const Index &fastem_version, const Verbosity &verbosity)
WORKSPACE METHOD: SurfaceFastem.
Definition: m_surface.cc:3739
void surfaceFlatScalarReflectivity(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &rtp_pos, const Vector &rtp_los, const Vector &specular_los, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, const Verbosity &)
WORKSPACE METHOD: surfaceFlatScalarReflectivity.
Definition: m_surface.cc:2664
void surface_typeInterpTypeMask(ArrayOfIndex &surface_types, Vector &surface_types_aux, Vector &surface_types_weights, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lat_true, const Vector &lon_true, const Vector &rtp_pos, const GriddedField2 &surface_type_mask, const String &method, const Verbosity &)
WORKSPACE METHOD: surface_typeInterpTypeMask.
Definition: m_surface.cc:3443
void surface_reflectivityFromGriddedField6(Tensor3 &surface_reflectivity, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lat_true, const Vector &lon_true, const Vector &rtp_pos, const Vector &rtp_los, const GriddedField6 &r_field, const Verbosity &)
WORKSPACE METHOD: surface_reflectivityFromGriddedField6.
Definition: m_surface.cc:3184
void surfaceFlatReflectivity(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &rtp_pos, const Vector &rtp_los, const Vector &specular_los, const Numeric &surface_skin_t, const Tensor3 &surface_reflectivity, const Verbosity &verbosity)
WORKSPACE METHOD: surfaceFlatReflectivity.
Definition: m_surface.cc:2500
void surfaceFlatRvRh(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &rtp_pos, const Vector &rtp_los, const Vector &specular_los, const Numeric &surface_skin_t, const Matrix &surface_rv_rh, const Verbosity &)
WORKSPACE METHOD: surfaceFlatRvRh.
Definition: m_surface.cc:2584
void iySurfaceFlatReflectivityDirect(Workspace &ws, Matrix &iy, 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 Tensor3 &surface_reflectivity, 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 &stars_do, const Index &gas_scattering_do, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfStar &stars, const Numeric &rte_alonglos_v, const String &iy_unit, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceFlatReflectivityDirect.
Definition: m_surface.cc:625
void iySurfaceCallAgendaX(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const String &iy_unit, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &cloudbox_on, const Index &jacobian_do, const Vector &f_grid, const Agenda &iy_main_agenda, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const ArrayOfAgenda &iy_surface_agenda_array, const ArrayOfIndex &surface_types, const Vector &surface_types_aux, const Vector &surface_types_weights, const Verbosity &)
WORKSPACE METHOD: iySurfaceCallAgendaX.
Definition: m_surface.cc:253
void transmittanceFromIy_aux(Vector &transmittance, const ArrayOfString &iy_aux_vars, const ArrayOfMatrix &iy_aux, const Verbosity &)
WORKSPACE METHOD: transmittanceFromIy_aux.
Definition: m_surface.cc:4288
void iySurfaceFlatReflectivity(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, ArrayOfTensor4 &dsurface_rmatrix_dx, ArrayOfMatrix &dsurface_emission_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const Index &stars_do, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Vector &lat_grid, const Vector &lon_grid, const Matrix &z_surface, const Vector &refellipsoid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const String &iy_unit, const Tensor3 &surface_reflectivity, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &iy_main_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceFlatReflectivity.
Definition: m_surface.cc:451
void surfaceLambertianSimple(Matrix &surface_los, Tensor4 &surface_rmatrix, Matrix &surface_emission, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &rtp_pos, const Vector &rtp_los, const Vector &surface_normal, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, const Index &lambertian_nza, const Numeric &za_pos, const Verbosity &)
WORKSPACE METHOD: surfaceLambertianSimple.
Definition: m_surface.cc:2730
void iySurfaceFastem(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda, const Numeric &surface_skin_t, const Numeric &salinity, const Numeric &wind_speed, const Numeric &wind_direction, const Index &fastem_version, const Verbosity &verbosity)
WORKSPACE METHOD: iySurfaceFastem.
Definition: m_surface.cc:332
void surface_scalar_reflectivityFromSurface_rmatrix(Vector &surface_scalar_reflectivity, const Tensor4 &surface_rmatrix, const Verbosity &)
WORKSPACE METHOD: surface_scalar_reflectivityFromSurface_rmatrix.
Definition: m_surface.cc:3425
void SurfaceTessem(Matrix &surface_los, Tensor4 &surface_rmatrix, ArrayOfTensor4 &dsurface_rmatrix_dx, Matrix &surface_emission, ArrayOfMatrix &dsurface_emission_dx, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const TessemNN &net_h, const TessemNN &net_v, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const Index &jacobian_do, const Verbosity &verbosity)
WORKSPACE METHOD: SurfaceTessem.
Definition: m_surface.cc:4116
void SurfaceDummy(ArrayOfTensor4 &dsurface_rmatrix_dx, ArrayOfMatrix &dsurface_emission_dx, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const Index &jacobian_do, const Verbosity &)
WORKSPACE METHOD: SurfaceDummy.
Definition: m_surface.cc:3710
void iySurfaceInit(Matrix &iy, const Vector &f_grid, const Index &stokes_dim, const Verbosity &)
WORKSPACE METHOD: iySurfaceInit.
Definition: m_surface.cc:1076
void InterpGriddedField2ToPosition(Numeric &outvalue, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lat_true, const Vector &lon_true, const Vector &rtp_pos, const GriddedField2 &gfield2, const Verbosity &)
WORKSPACE METHOD: InterpGriddedField2ToPosition.
Definition: m_surface.cc:147
void iySurfaceRtpropCalc(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Matrix &surface_los, const Tensor4 &surface_rmatrix, const Matrix &surface_emission, const ArrayOfString &dsurface_names, const ArrayOfTensor4 &dsurface_rmatrix_dx, const ArrayOfMatrix &dsurface_emission_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const Index &stars_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda, const Verbosity &)
WORKSPACE METHOD: iySurfaceRtpropCalc.
Definition: m_surface.cc:1715
void FastemStandAlone(Matrix &emissivity, Matrix &reflectivity, const Vector &f_grid, const Numeric &surface_skin_t, const Numeric &za, const Numeric &salinity, const Numeric &wind_speed, const Numeric &rel_aa, const Vector &transmittance, const Index &fastem_version, const Verbosity &)
WORKSPACE METHOD: FastemStandAlone.
Definition: m_surface.cc:72
void SurfaceFlatScalarReflectivity(Matrix &surface_los, Tensor4 &surface_rmatrix, ArrayOfTensor4 &dsurface_rmatrix_dx, Matrix &surface_emission, ArrayOfMatrix &dsurface_emission_dx, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Vector &specular_los, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const Index &jacobian_do, const Vector &f_reflectivities, const Verbosity &verbosity)
WORKSPACE METHOD: SurfaceFlatScalarReflectivity.
Definition: m_surface.cc:3954
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:65
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:442
void calculate_int_weights_arbitrary_grid(Vector &w, const Vector &x)
Calculates trapezoidal integration weights for arbitray grid.
Definition: math_funcs.cc:830
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1348
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
const Joker joker
std::complex< Numeric > Complex
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:206
#define CREATE_OUT2
Definition: messages.h:205
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:216
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric R
Ideal gas constant convenience name [J/mol K].
constexpr Numeric earth_radius
Global constant, the radius of the Earth [m].
constexpr Numeric e
Elementary charge convenience name [C].
constexpr Numeric alpha
Fine structure constant convenience name [-].
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
Numeric planck(const Numeric &f, const Numeric &t)
planck
Numeric dplanck_dt(const Numeric &f, const Numeric &t)
dplanck_dt
void fresnel(Complex &Rv, Complex &Rh, const Complex &n1, const Complex &n2, const Numeric &theta)
fresnel
This file contains declerations of functions of physical character.
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to a cartesian unit vector.
Definition: ppath.cc:313
void cart2zaaa(Numeric &za, Numeric &aa, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Converts a cartesian directional vector to zenith and azimuth.
Definition: ppath.cc:300
void plevel_slope_3d(Numeric &c1, Numeric &c2, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon, const Numeric &aa)
Definition: ppath.cc:1068
void resolve_lon(Numeric &lon, const Numeric &lon5, const Numeric &lon6)
Resolves which longitude angle that shall be used.
Definition: ppath.cc:501
Numeric plevel_angletilt(const Numeric &r, const Numeric &c1)
Calculates the angular tilt of the surface or a pressure level.
Definition: ppath.cc:617
void plevel_slope_2d(Numeric &c1, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstVectorView z_surf, const GridPos &gp, const Numeric &za)
Calculates the radial slope of the surface or a pressure level for 2D.
Definition: ppath.cc:580
Propagation path structure and functions.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
void mueller_stokes2modif(Matrix &Cm, const Index &stokes_dim)
mueller_stokes2modif
Definition: rte.cc:1913
void pos2true_latlon(Numeric &lat, Numeric &lon, const Index &atmosphere_dim, const ConstVectorView &lat_grid, const ConstVectorView &lat_true, const ConstVectorView &lon_true, const ConstVectorView &pos)
Determines the true alt and lon for an "ARTS position".
Definition: rte.cc:1934
void iy_transmittance_mult(Tensor3 &iy_trans_total, const ConstTensor3View &iy_trans_old, const ConstTensor3View &iy_trans_new)
Multiplicates iy_transmittance with transmissions.
Definition: rte.cc:1775
void adjust_los(VectorView los, const Index &atmosphere_dim)
Ensures that the zenith and azimuth angles of a line-of-sight vector are inside defined ranges.
Definition: rte.cc:86
void mueller_modif2stokes(Matrix &Cs, const Index &stokes_dim)
mueller_modif2stokes
Definition: rte.cc:1865
void mueller_rotation(Matrix &L, const Index &stokes_dim, const Numeric &rotangle)
mueller_rotation
Definition: rte.cc:1886
void mirror_los(Vector &los_mirrored, const ConstVectorView &los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:1812
Declaration of functions in rte.cc.
void complex_n_interp(MatrixView n_real, MatrixView n_imag, const GriddedField3 &complex_n, const String &varname, ConstVectorView f_grid, ConstVectorView t_grid)
General function for interpolating data of complex n type.
void interp_atmsurface_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of a surface-type variable.
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
Converts a geographical position (rte_pos) to grid positions for p, lat and lon.
void interp_atmsurface_by_gp(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates a surface-type variable given the grid positions.
Header file for special_interp.cc.
void get_star_ppaths(Workspace &ws, ArrayOfPpath &star_ppaths, ArrayOfIndex &stars_visible, ArrayOfVector &star_rte_los, const Vector &rte_pos, const ArrayOfStar &stars, 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 Matrix &z_surface, const Vector &refellipsoid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
Calculates the ppath towards the stars from a given position and indicates if star is visible or not.
Definition: star.cc:391
void get_direct_radiation(Workspace &ws, ArrayOfMatrix &direct_radiation, ArrayOfArrayOfTensor3 &ddirect_radiation_dx, 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 &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 Index &irradiance_flag, const ArrayOfPpath &star_ppaths, const ArrayOfStar &stars, const ArrayOfIndex &stars_visible, const Vector &refellipsoid, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Numeric &rte_alonglos_v, const Verbosity &verbosity)
Calculates the transmitted star radiation at the end position of the ppath.
Definition: star.cc:215
Structure to store a grid position.
Definition: interpolation.h:73
std::array< Numeric, 2 > fd
Definition: interpolation.h:75
Index idx
Definition: interpolation.h:74
A Lagrange interpolation computer.
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
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
void tessem_prop_nn(VectorView &ny, const TessemNN &net, ConstVectorView nx)
Definition: tessem.cc:87
This file contains functions that are adapted from TESSEM code which is used to calculate surface emi...
#define d
#define w
#define c
#define b