ARTS 2.5.4 (git: 4c0d3b4d)
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 "auto_md.h"
43#include "check_input.h"
44#include "matpack_complex.h"
45#include "fastem.h"
46#include "geodetic.h"
47#include "interpolation.h"
49#include "math_funcs.h"
50#include "messages.h"
51#include "physics_funcs.h"
52#include "ppath.h"
53#include "rte.h"
54#include "special_interp.h"
55#include "surface.h"
56#include "tessem.h"
57
59extern const Numeric DEG2RAD;
60
61/*===========================================================================
62 === The functions (in alphabetical order)
63 ===========================================================================*/
64
65/* Workspace method: Doxygen documentation will be auto-generated */
66void FastemStandAlone(Matrix& emissivity,
67 Matrix& reflectivity,
68 const Vector& f_grid,
69 const Numeric& surface_skin_t,
70 const Numeric& za,
71 const Numeric& salinity,
72 const Numeric& wind_speed,
73 const Numeric& rel_aa,
74 const Vector& transmittance,
75 const Index& fastem_version,
76 const Verbosity&) {
77 const Index nf = f_grid.nelem();
78
79 chk_if_in_range("zenith angle", za, 90, 180);
80 chk_if_in_range_exclude("surface skin temperature", surface_skin_t, 260, 373);
81 chk_if_in_range_exclude_high("salinity", salinity, 0, 1);
82 chk_if_in_range_exclude_high("wind speed", wind_speed, 0, 100);
83 chk_if_in_range("azimuth angle", rel_aa, -180, 180);
84 chk_vector_length("transmittance", "f_grid", transmittance, f_grid);
85 if (fastem_version < 3 || fastem_version > 6)
86 throw std::runtime_error(
87 "Invalid fastem version: 3 <= fastem_version <= 6");
88
89 emissivity.resize(nf, 4);
90 reflectivity.resize(nf, 4);
91
92 const Numeric t = max(surface_skin_t, Numeric(270));
93
94 for (Index i = 0; i < nf; i++) {
95 if (f_grid[i] > 250e9)
96 throw std::runtime_error("Only frequency <= 250 GHz are allowed");
97 chk_if_in_range("transmittance", transmittance[i], 0, 1);
98
99 Vector e, r;
100
101 fastem(e,
102 r,
103 f_grid[i],
104 za,
105 t,
106 salinity,
107 wind_speed,
108 transmittance[i],
109 rel_aa,
110 fastem_version);
111
112 emissivity(i, joker) = e;
113 reflectivity(i, joker) = r;
114 }
115
116 // FASTEM does not work close to the horizon (at least v6). Make sure values
117 // are inside [0,1]. Then seems best to make sure that e+r=1.
118 for (Index i = 0; i < nf; i++) {
119 for (Index s = 0; s < 2; s++) {
120 if (emissivity(i, s) > 1) {
121 emissivity(i, s) = 1;
122 reflectivity(i, s) = 0;
123 }
124 if (emissivity(i, s) < 0) {
125 emissivity(i, s) = 0;
126 reflectivity(i, s) = 1;
127 }
128 if (reflectivity(i, s) > 1) {
129 emissivity(i, s) = 0;
130 reflectivity(i, s) = 1;
131 }
132 if (reflectivity(i, s) < 0) {
133 emissivity(i, s) = 1;
134 reflectivity(i, s) = 0;
135 }
136 }
137 }
138}
139
140/* Workspace method: Doxygen documentation will be auto-generated */
142 const Index& atmosphere_dim,
143 const Vector& lat_grid,
144 const Vector& lat_true,
145 const Vector& lon_true,
146 const Vector& rtp_pos,
147 const GriddedField2& gfield2,
148 const Verbosity&) {
149 // Set expected order of grids
150 Index gfield_latID = 0;
151 Index gfield_lonID = 1;
152
153 // Basic checks and sizes
154 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
155 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
156 chk_rte_pos(atmosphere_dim, rtp_pos);
157 gfield2.checksize_strict();
158 //
159 chk_griddedfield_gridname(gfield2, gfield_latID, "Latitude");
160 chk_griddedfield_gridname(gfield2, gfield_lonID, "Longitude");
161 //
162 const Index nlat = gfield2.data.nrows();
163 const Index nlon = gfield2.data.ncols();
164 //
165 if (nlat < 2 || nlon < 2) {
166 ostringstream os;
167 os << "The data in *gfield2* must span a geographical region. That is,\n"
168 << "the latitude and longitude grids must have a length >= 2.";
169 throw runtime_error(os.str());
170 }
171
172 const Vector& GFlat = gfield2.get_numeric_grid(gfield_latID);
173 const Vector& GFlon = gfield2.get_numeric_grid(gfield_lonID);
174
175 // Determine true geographical position
176 Vector lat(1), lon(1);
178 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
179
180 // Ensure correct coverage of lon grid
181 Vector lon_shifted;
182 lon_shiftgrid(lon_shifted, GFlon, lon[0]);
183
184 // Check if lat/lon we need are actually covered
185 chk_if_in_range("rtp_pos.lat", lat[0], GFlat[0], GFlat[nlat - 1]);
186 chk_if_in_range("rtp_pos.lon", lon[0], lon_shifted[0], lon_shifted[nlon - 1]);
187
188 // Interpolate in lat and lon
189 //
190 GridPos gp_lat, gp_lon;
191 gridpos(gp_lat, GFlat, lat[0]);
192 gridpos(gp_lon, lon_shifted, lon[0]);
193 Vector itw(4);
194 interpweights(itw, gp_lat, gp_lon);
195 outvalue = interp(itw, gfield2.data, gp_lat, gp_lon);
196}
197
198
199/* Workspace method: Doxygen documentation will be auto-generated */
201 const Index& atmosphere_dim,
202 const Vector& lat_grid,
203 const Vector& lon_grid,
204 const Vector& rtp_pos,
205 const Matrix& z_surface,
206 const Matrix& field,
207 const Verbosity& verbosity) {
208 // Input checks (dummy p_grid)
209 chk_atm_grids(atmosphere_dim, Vector(2, 2, -1), lat_grid, lon_grid);
211 "input argument *field*", field, atmosphere_dim, lat_grid, lon_grid);
212 chk_rte_pos(atmosphere_dim, rtp_pos);
213 //
214 const Numeric zmax = max(z_surface);
215 const Numeric zmin = min(z_surface);
216 const Numeric dzok = 1;
217 if (rtp_pos[0] < zmin - dzok || rtp_pos[0] > zmax + dzok) {
218 ostringstream os;
219 os << "The given position does not match *z_surface*.\nThe altitude in "
220 << "*rtp_pos* is " << rtp_pos[0] / 1e3 << " km.\n"
221 << "The altitude range covered by *z_surface* is [" << zmin / 1e3 << ","
222 << zmax / 1e3 << "] km.\n"
223 << "One possible mistake is to mix up *rtp_pos* and *rte_los*.";
224 throw runtime_error(os.str());
225 }
226
227 if (atmosphere_dim == 1) {
228 outvalue = field(0, 0);
229 } else {
230 chk_interpolation_grids("Latitude interpolation", lat_grid, rtp_pos[1]);
231 GridPos gp_lat, gp_lon;
232 gridpos(gp_lat, lat_grid, rtp_pos[1]);
233 if (atmosphere_dim == 3) {
234 chk_interpolation_grids("Longitude interpolation", lon_grid, rtp_pos[2]);
235 gridpos(gp_lon, lon_grid, rtp_pos[2]);
236 }
237 //
238 outvalue = interp_atmsurface_by_gp(atmosphere_dim, field, gp_lat, gp_lon);
239 }
240
241 // Interpolate
243 out3 << " Result = " << outvalue << "\n";
244}
245
246/* Workspace method: Doxygen documentation will be auto-generated */
248 Matrix& iy,
249 ArrayOfTensor3& diy_dx,
250 const String& iy_unit,
251 const Tensor3& iy_transmittance,
252 const Index& iy_id,
253 const Index& cloudbox_on,
254 const Index& jacobian_do,
255 const Vector& f_grid,
256 const Agenda& iy_main_agenda,
257 const Vector& rtp_pos,
258 const Vector& rtp_los,
259 const Vector& rte_pos2,
260 const ArrayOfAgenda& iy_surface_agenda_array,
261 const ArrayOfIndex& surface_types,
262 const Vector& surface_types_aux,
263 const Vector& surface_types_weights,
264 const Verbosity&) {
265
266 const Index ntypes = surface_types.nelem();
267 if (ntypes < 1)
268 throw runtime_error("*surface_types* is empty!");
269
270 // Loop surface types and sum up
271 //
272 Numeric wtot = 0;
273 //
274 for (Index t=0; t<ntypes; t++ ) {
275
276 if (surface_types[t] < 0)
277 throw runtime_error(
278 "No element in *surface_types* is allowed to be negative.");
279 if (surface_types[t] >= iy_surface_agenda_array.nelem()) {
280 ostringstream os;
281 os << "*iy_surface_agenda_array* has only "
282 << iy_surface_agenda_array.nelem()
283 << " elements,\n while you have selected element " << surface_types[t];
284 throw runtime_error(os.str());
285 }
286
287 Matrix iy1;
288 ArrayOfTensor3 diy_dx1;
289
291 iy1,
292 diy_dx1,
293 surface_types[t],
294 iy_unit,
295 iy_transmittance,
296 iy_id,
297 cloudbox_on,
298 jacobian_do,
299 iy_main_agenda,
300 f_grid,
301 rtp_pos,
302 rtp_los,
303 rte_pos2,
304 surface_types_aux[t],
305 iy_surface_agenda_array);
306
307 iy1 *= surface_types_weights[t];
308 for (Index i=0; diy_dx1.nelem(); i++ )
309 diy_dx1[i] *= surface_types_weights[t];
310 wtot += surface_types_weights[t];;
311
312 if (t==0) {
313 iy = iy1;
314 diy_dx = diy_dx1;
315 } else {
316 iy += iy1;
317 for (Index i=0; diy_dx1.nelem(); i++ )
318 diy_dx[i] += diy_dx1[i];
319 }
320 }
321 if (abs(wtot-1)>1e-4)
322 throw runtime_error("Sum of *surface_types_weights* deviates from 1.");
323}
324
325/* Workspace method: Doxygen documentation will be auto-generated */
327 Matrix& iy,
328 ArrayOfTensor3& diy_dx,
329 const Tensor3& iy_transmittance,
330 const Index& iy_id,
331 const Index& jacobian_do,
332 const Index& atmosphere_dim,
333 const EnergyLevelMap& nlte_field,
334 const Index& cloudbox_on,
335 const Index& stokes_dim,
336 const Vector& f_grid,
337 const Vector& rtp_pos,
338 const Vector& rtp_los,
339 const Vector& rte_pos2,
340 const String& iy_unit,
341 const Agenda& iy_main_agenda,
342 const Numeric& surface_skin_t,
343 const Numeric& salinity,
344 const Numeric& wind_speed,
345 const Numeric& wind_direction,
346 const Index& fastem_version,
347 const Verbosity& verbosity) {
348 // Most obvious input checks are performed in specular_losCalc and surfaceFastem
349
350 // Obtain radiance and transmission for specular direction
351
352 // Determine specular direction
353 Vector specular_los, surface_normal;
354 specular_losCalcNoTopography(specular_los,
355 surface_normal,
356 rtp_pos,
357 rtp_los,
358 atmosphere_dim,
359 verbosity);
360
361 // Use iy_aux to get optical depth for downwelling radiation.
362 ArrayOfString iy_aux_vars(1);
363 iy_aux_vars[0] = "Optical depth";
364
365 // Calculate iy for downwelling radiation
366 // Note that iy_transmittance used here lacks surface R. Fixed below.
367 //
368 const Index nf = f_grid.nelem();
369 Vector transmittance(nf);
370 ArrayOfMatrix iy_aux;
371 Ppath ppath;
372 //
374 iy,
375 iy_aux,
376 ppath,
377 diy_dx,
378 0,
379 iy_transmittance,
380 iy_aux_vars,
381 iy_id,
382 iy_unit,
383 cloudbox_on,
384 jacobian_do,
385 f_grid,
386 nlte_field,
387 rtp_pos,
388 specular_los,
389 rte_pos2,
390 iy_main_agenda);
391
392 // Convert tau to transmissions
393 for (Index i = 0; i < nf; i++) {
394 transmittance[i] = exp(-iy_aux[0](i, 0));
395 }
396
397 // Call Fastem by surface_RTprop version
398 //
399 Matrix surface_los;
400 Tensor4 surface_rmatrix;
401 Matrix surface_emission;
402 //
403 surfaceFastem(surface_los,
404 surface_rmatrix,
405 surface_emission,
406 atmosphere_dim,
407 stokes_dim,
408 f_grid,
409 rtp_pos,
410 rtp_los,
411 surface_skin_t,
412 salinity,
413 wind_speed,
414 wind_direction,
415 transmittance,
416 fastem_version,
417 verbosity);
418
419 // Add up
420 //
421 Tensor3 I(1, nf, stokes_dim);
422 I(0, joker, joker) = iy;
423 Matrix sensor_los_dummy(1, 1, 0);
424 //
425 surface_calc(iy, I, sensor_los_dummy, surface_rmatrix, surface_emission);
426
427 // Adjust diy_dx, if necessary.
428 // For vector cases this is a slight approximation, as the order of the
429 // different transmission and reflectivities matters.
430 if (iy_transmittance.npages()) {
431 for (Index q = 0; q < diy_dx.nelem(); q++) {
432 for (Index p = 0; p < diy_dx[q].npages(); p++) {
433 for (Index i = 0; i < nf; i++) {
434 Vector x = diy_dx[q](p, i, joker);
435 mult(diy_dx[q](p, i, joker), surface_rmatrix(0, i, joker, joker), x);
436 }
437 }
438 }
439 }
440}
441
442/* Workspace method: Doxygen documentation will be auto-generated */
444 Matrix& iy,
445 ArrayOfTensor3& diy_dx,
446 Numeric& surface_skin_t,
447 Matrix& surface_los,
448 Tensor4& surface_rmatrix,
449 Matrix& surface_emission,
450 const Tensor3& iy_transmittance,
451 const Index& iy_id,
452 const Index& jacobian_do,
453 const Index& atmosphere_dim,
454 const EnergyLevelMap& nlte_field,
455 const Index& cloudbox_on,
456 const Index& stokes_dim,
457 const Vector& f_grid,
458 const Vector& rtp_pos,
459 const Vector& rtp_los,
460 const Vector& rte_pos2,
461 const String& iy_unit,
462 const Agenda& iy_main_agenda,
463 const Agenda& surface_rtprop_agenda,
464 const Verbosity&) {
465 // Input checks
466 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
467 chk_rte_pos(atmosphere_dim, rtp_pos);
468 chk_rte_los(atmosphere_dim, rtp_los);
469
470 // Call *surface_rtprop_agenda*
472 surface_skin_t,
473 surface_emission,
474 surface_los,
475 surface_rmatrix,
476 f_grid,
477 rtp_pos,
478 rtp_los,
479 surface_rtprop_agenda);
480
481 // Check output of *surface_rtprop_agenda*
482 const Index nlos = surface_los.nrows();
483 const Index nf = f_grid.nelem();
484 //
485 if (nlos) // if 0, blackbody ground and not all checks are needed
486 {
487 if (surface_los.ncols() != rtp_los.nelem())
488 throw runtime_error("Number of columns in *surface_los* is not correct.");
489 if (nlos != surface_rmatrix.nbooks())
490 throw runtime_error(
491 "Mismatch in size of *surface_los* and *surface_rmatrix*.");
492 if (surface_rmatrix.npages() != nf)
493 throw runtime_error(
494 "Mismatch in size of *surface_rmatrix* and *f_grid*.");
495 if (surface_rmatrix.nrows() != stokes_dim ||
496 surface_rmatrix.ncols() != stokes_dim)
497 throw runtime_error(
498 "Mismatch between size of *surface_rmatrix* and *stokes_dim*.");
499 }
500 if (surface_emission.ncols() != stokes_dim)
501 throw runtime_error(
502 "Mismatch between size of *surface_emission* and *stokes_dim*.");
503 if (surface_emission.nrows() != nf)
504 throw runtime_error("Mismatch in size of *surface_emission* and f_grid*.");
505
506 // Variable to hold down-welling radiation
507 Tensor3 I(nlos, nf, stokes_dim);
508
509 // Loop *surface_los*-es. If no such LOS, we are ready.
510 if (nlos > 0) {
511 for (Index ilos = 0; ilos < nlos; ilos++) {
512 Vector los = surface_los(ilos, joker);
513
514 // Include surface reflection matrix in *iy_transmittance*
515 // If iy_transmittance is empty, this is interpreted as the
516 // variable is not needed.
517 //
518 Tensor3 iy_trans_new;
519 //
520 if (iy_transmittance.npages()) {
521 iy_transmittance_mult(iy_trans_new,
522 iy_transmittance,
523 surface_rmatrix(ilos, joker, joker, joker));
524 }
525
526 // Calculate downwelling radiation for LOS ilos
527 //
528 {
529 ArrayOfMatrix iy_aux;
530 Ppath ppath;
531 Index iy_id_new = iy_id + ilos + 1;
533 iy,
534 iy_aux,
535 ppath,
536 diy_dx,
537 0,
538 iy_trans_new,
539 ArrayOfString(0),
540 iy_id_new,
541 iy_unit,
542 cloudbox_on,
543 jacobian_do,
544 f_grid,
545 nlte_field,
546 rtp_pos,
547 los,
548 rte_pos2,
549 iy_main_agenda);
550 }
551
552 if (iy.ncols() != stokes_dim || iy.nrows() != nf) {
553 ostringstream os;
554 os << "The size of *iy* returned from *" << iy_main_agenda.name()
555 << "* is\n"
556 << "not correct:\n"
557 << " expected size = [" << nf << "," << stokes_dim << "]\n"
558 << " size of iy = [" << iy.nrows() << "," << iy.ncols() << "]\n";
559 throw runtime_error(os.str());
560 }
561
562 I(ilos, joker, joker) = iy;
563 }
564 }
565
566 // Add up
567 surface_calc(iy, I, surface_los, surface_rmatrix, surface_emission);
568}
569
570/* Workspace method: Doxygen documentation will be auto-generated */
572 Matrix& iy,
573 ArrayOfTensor3& diy_dx,
574 const Matrix& surface_los,
575 const Tensor4& surface_rmatrix,
576 const Matrix& surface_emission,
577 const ArrayOfString& dsurface_names,
578 const ArrayOfTensor4& dsurface_rmatrix_dx,
579 const ArrayOfMatrix& dsurface_emission_dx,
580 const Tensor3& iy_transmittance,
581 const Index& iy_id,
582 const Index& jacobian_do,
583 const ArrayOfRetrievalQuantity& jacobian_quantities,
584 const Index& atmosphere_dim,
585 const EnergyLevelMap& nlte_field,
586 const Index& cloudbox_on,
587 const Index& stokes_dim,
588 const Vector& f_grid,
589 const Vector& rtp_pos,
590 const Vector& rtp_los,
591 const Vector& rte_pos2,
592 const String& iy_unit,
593 const Agenda& iy_main_agenda,
594 const Verbosity&) {
595 // Input checks
596 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
597 chk_rte_pos(atmosphere_dim, rtp_pos);
598 chk_rte_los(atmosphere_dim, rtp_los);
599
600 // Check provided surface rtprop variables
601 const Index nlos = surface_los.nrows();
602 const Index nf = f_grid.nelem();
603 //
604 if (nlos) // if 0, blackbody ground and not all checks are needed
605 {
606 if (surface_los.ncols() != rtp_los.nelem())
607 throw runtime_error("Number of columns in *surface_los* is not correct.");
608 if (nlos != surface_rmatrix.nbooks())
609 throw runtime_error(
610 "Mismatch in size of *surface_los* and *surface_rmatrix*.");
611 if (surface_rmatrix.npages() != nf)
612 throw runtime_error(
613 "Mismatch in size of *surface_rmatrix* and *f_grid*.");
614 if (surface_rmatrix.nrows() != stokes_dim ||
615 surface_rmatrix.ncols() != stokes_dim)
616 throw runtime_error(
617 "Mismatch between size of *surface_rmatrix* and *stokes_dim*.");
618 }
619 if (surface_emission.ncols() != stokes_dim)
620 throw runtime_error(
621 "Mismatch between size of *surface_emission* and *stokes_dim*.");
622 if (surface_emission.nrows() != nf)
623 throw runtime_error("Mismatch in size of *surface_emission* and f_grid*.");
624
625 // Variable to hold down-welling radiation
626 Tensor3 I(nlos, nf, stokes_dim);
627
628 // Loop *surface_los*-es.
629 if (nlos > 0) {
630 for (Index ilos = 0; ilos < nlos; ilos++) {
631 Vector los = surface_los(ilos, joker);
632
633 // Include surface reflection matrix in *iy_transmittance*
634 // If iy_transmittance is empty, this is interpreted as the
635 // variable is not needed.
636 //
637 Tensor3 iy_trans_new;
638 //
639 if (iy_transmittance.npages()) {
640 iy_transmittance_mult(iy_trans_new,
641 iy_transmittance,
642 surface_rmatrix(ilos, joker, joker, joker));
643 }
644
645 // Calculate downwelling radiation for LOS ilos
646 //
647 {
648 ArrayOfMatrix iy_aux;
649 Ppath ppath;
651 iy,
652 iy_aux,
653 ppath,
654 diy_dx,
655 0,
656 iy_trans_new,
657 ArrayOfString(0),
658 iy_id,
659 iy_unit,
660 cloudbox_on,
661 jacobian_do,
662 f_grid,
663 nlte_field,
664 rtp_pos,
665 los,
666 rte_pos2,
667 iy_main_agenda);
668 }
669
670 if (iy.ncols() != stokes_dim || iy.nrows() != nf) {
671 ostringstream os;
672 os << "The size of *iy* returned from *" << iy_main_agenda.name()
673 << "* is\n"
674 << "not correct:\n"
675 << " expected size = [" << nf << "," << stokes_dim << "]\n"
676 << " size of iy = [" << iy.nrows() << "," << iy.ncols() << "]\n";
677 throw runtime_error(os.str());
678 }
679
680 I(ilos, joker, joker) = iy;
681 }
682 }
683
684 // Add up
685 surface_calc(iy, I, surface_los, surface_rmatrix, surface_emission);
686
687 // Surface Jacobians
688 if (jacobian_do && dsurface_names.nelem()) {
689 // Loop dsurface_names
690 for (Index i = 0; i < dsurface_names.nelem(); i++) {
691 // Error if derivatives not calculated
692 // Or should we accept this?
693 if (dsurface_emission_dx[i].empty() || dsurface_rmatrix_dx[i].empty()) {
694 ostringstream os;
695 os << "The derivatives for surface quantity: " << dsurface_names[i]
696 << "\nwere not calculated by *iy_surface_agenda*.\n"
697 << "That is, *dsurface_emission_dx* and/or *dsurface_rmatrix_dx*\n"
698 << "are empty.";
699 throw runtime_error(os.str());
700 } else {
701 // Find index among jacobian quantities
702 Index ihit = -1;
703 for (Index j = 0; j < jacobian_quantities.nelem() && ihit < 0; j++) {
704 if (dsurface_names[i] == jacobian_quantities[j].Subtag()) {
705 ihit = j;
706 }
707 }
708 ARTS_ASSERT(ihit >= 0);
709 // Derivative, as observed at the surface
710 Matrix diy_dpos0, diy_dpos;
711 surface_calc(diy_dpos0,
712 I,
713 surface_los,
714 dsurface_rmatrix_dx[i],
715 dsurface_emission_dx[i]);
716 // Weight with transmission to sensor
717 iy_transmittance_mult(diy_dpos, iy_transmittance, diy_dpos0);
718 // Put into diy_dx
719 diy_from_pos_to_rgrids(diy_dx[ihit],
720 jacobian_quantities[ihit],
721 diy_dpos,
722 atmosphere_dim,
723 rtp_pos);
724 }
725 }
726 }
727}
728
729/* Workspace method: Doxygen documentation will be auto-generated */
731 Vector& surface_normal,
732 const Vector& rtp_pos,
733 const Vector& rtp_los,
734 const Index& atmosphere_dim,
735 const Verbosity&) {
736 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
737 chk_rte_pos(atmosphere_dim, rtp_pos);
738 chk_rte_los(atmosphere_dim, rtp_los);
739
740 surface_normal.resize(max(Index(1), atmosphere_dim - 1));
741 specular_los.resize(max(Index(1), atmosphere_dim - 1));
742
743 if (atmosphere_dim == 1) {
744 surface_normal[0] = 0;
745 if (rtp_los[0] < 90) {
746 throw runtime_error(
747 "Invalid zenith angle. The zenith angle corresponds "
748 "to observe the surface from below.");
749 }
750 specular_los[0] = 180 - rtp_los[0];
751 }
752
753 else if (atmosphere_dim == 2) {
754 specular_los[0] = sign(rtp_los[0]) * 180 - rtp_los[0];
755 surface_normal[0] = 0;
756 }
757
758 else if (atmosphere_dim == 3) {
759 specular_los[0] = 180 - rtp_los[0];
760 specular_los[1] = rtp_los[1];
761 surface_normal[0] = 0;
762 surface_normal[1] = 0;
763 }
764}
765
766/* Workspace method: Doxygen documentation will be auto-generated */
767void specular_losCalc(Vector& specular_los,
768 Vector& surface_normal,
769 const Vector& rtp_pos,
770 const Vector& rtp_los,
771 const Index& atmosphere_dim,
772 const Vector& lat_grid,
773 const Vector& lon_grid,
774 const Vector& refellipsoid,
775 const Matrix& z_surface,
776 const Index& ignore_surface_slope,
777 const Verbosity& verbosity) {
778 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
779 chk_rte_pos(atmosphere_dim, rtp_pos);
780 chk_rte_los(atmosphere_dim, rtp_los);
781 chk_if_in_range("ignore_surface_slope", ignore_surface_slope, 0, 1);
782
783 // Use special function if there is no slope, or it is ignored
784 if (atmosphere_dim == 1 || ignore_surface_slope) {
785 specular_losCalcNoTopography(specular_los,
786 surface_normal,
787 rtp_pos,
788 rtp_los,
789 atmosphere_dim,
790 verbosity);
791 return;
792 }
793
794 surface_normal.resize(max(Index(1), atmosphere_dim - 1));
795 specular_los.resize(max(Index(1), atmosphere_dim - 1));
796
797 if (atmosphere_dim == 2) {
798 chk_interpolation_grids("Latitude interpolation", lat_grid, rtp_pos[1]);
799 GridPos gp_lat;
800 gridpos(gp_lat, lat_grid, rtp_pos[1]);
801 Numeric c1; // Radial slope of the surface
803 c1, lat_grid, refellipsoid, z_surface(joker, 0), gp_lat, rtp_los[0]);
804 Vector itw(2);
805 interpweights(itw, gp_lat);
806 const Numeric rv_surface = refell2d(refellipsoid, lat_grid, gp_lat) +
807 interp(itw, z_surface(joker, 0), gp_lat);
808 surface_normal[0] = -plevel_angletilt(rv_surface, c1);
809 if (abs(rtp_los[0] - surface_normal[0]) < 90) {
810 throw runtime_error(
811 "Invalid zenith angle. The zenith angle corresponds "
812 "to observe the surface from below.");
813 }
814 specular_los[0] =
815 sign(rtp_los[0]) * 180 - rtp_los[0] + 2 * surface_normal[0];
816 }
817
818 else {
819 // Calculate surface normal in South-North direction
820 chk_interpolation_grids("Latitude interpolation", lat_grid, rtp_pos[1]);
821 chk_interpolation_grids("Longitude interpolation", lon_grid, rtp_pos[2]);
822 GridPos gp_lat, gp_lon;
823 gridpos(gp_lat, lat_grid, rtp_pos[1]);
824 gridpos(gp_lon, lon_grid, rtp_pos[2]);
825 Numeric c1, c2;
827 c1, c2, lat_grid, lon_grid, refellipsoid, z_surface, gp_lat, gp_lon, 0);
828 Vector itw(4);
829 interpweights(itw, gp_lat, gp_lon);
830 const Numeric rv_surface = refell2d(refellipsoid, lat_grid, gp_lat) +
831 interp(itw, z_surface, gp_lat, gp_lon);
832 const Numeric zaSN = 90 - plevel_angletilt(rv_surface, c1);
833 // The same for East-West
835 c2,
836 lat_grid,
837 lon_grid,
838 refellipsoid,
839 z_surface,
840 gp_lat,
841 gp_lon,
842 90);
843 const Numeric zaEW = 90 - plevel_angletilt(rv_surface, c1);
844 // Convert to Cartesian, and determine normal by cross-product
845 Vector tangentSN(3), tangentEW(3), normal(3);
846 zaaa2cart(tangentSN[0], tangentSN[1], tangentSN[2], zaSN, 0);
847 zaaa2cart(tangentEW[0], tangentEW[1], tangentEW[2], zaEW, 90);
848 cross3(normal, tangentSN, tangentEW);
849 // Convert rtp_los to cartesian and flip direction
850 Vector di(3);
851 zaaa2cart(di[0], di[1], di[2], rtp_los[0], rtp_los[1]);
852 di *= -1;
853 // Set LOS normal vector
854 cart2zaaa(
855 surface_normal[0], surface_normal[1], normal[0], normal[1], normal[2]);
856 if (abs(rtp_los[0] - surface_normal[0]) < 90) {
857 throw runtime_error(
858 "Invalid zenith angle. The zenith angle corresponds "
859 "to observe the surface from below.");
860 }
861 // Specular direction is 2(dn*di)dn-di, where dn is the normal vector
862 Vector speccart(3);
863 const Numeric fac = 2 * (normal * di);
864 for (Index i = 0; i < 3; i++) {
865 speccart[i] = fac * normal[i] - di[i];
866 }
867 cart2zaaa(specular_los[0],
868 specular_los[1],
869 speccart[0],
870 speccart[1],
871 speccart[2]);
872 }
873}
874
875/* Workspace method: Doxygen documentation will be auto-generated */
876void surfaceBlackbody(Matrix& surface_los,
877 Tensor4& surface_rmatrix,
878 Matrix& surface_emission,
879 const Index& atmosphere_dim,
880 const Vector& f_grid,
881 const Index& stokes_dim,
882 const Vector& rtp_pos,
883 const Vector& rtp_los,
884 const Numeric& surface_skin_t,
885 const Verbosity& verbosity) {
886 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
887 chk_rte_pos(atmosphere_dim, rtp_pos);
888 chk_rte_los(atmosphere_dim, rtp_los);
889 chk_not_negative("surface_skin_t", surface_skin_t);
890
892 out2 << " Sets variables to model a blackbody surface with a temperature "
893 << " of " << surface_skin_t << " K.\n";
894
895 surface_los.resize(0, 0);
896 surface_rmatrix.resize(0, 0, 0, 0);
897
898 const Index nf = f_grid.nelem();
899
900 Vector b(nf);
901 planck(b, f_grid, surface_skin_t);
902
903 surface_emission.resize(nf, stokes_dim);
904 surface_emission = 0.0;
905
906 for (Index iv = 0; iv < nf; iv++) {
907 surface_emission(iv, 0) = b[iv];
908 for (Index is = 1; is < stokes_dim; is++) {
909 surface_emission(iv, is) = 0;
910 }
911 }
912}
913
914/* Workspace method: Doxygen documentation will be auto-generated */
915void surfaceFastem(Matrix& surface_los,
916 Tensor4& surface_rmatrix,
917 Matrix& surface_emission,
918 const Index& atmosphere_dim,
919 const Index& stokes_dim,
920 const Vector& f_grid,
921 const Vector& rtp_pos,
922 const Vector& rtp_los,
923 const Numeric& surface_skin_t,
924 const Numeric& salinity,
925 const Numeric& wind_speed,
926 const Numeric& wind_direction,
927 const Vector& transmittance,
928 const Index& fastem_version,
929 const Verbosity& verbosity) {
930 // Input checks
931 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
932 chk_rte_pos(atmosphere_dim, rtp_pos);
933 chk_rte_los(atmosphere_dim, rtp_los);
934 chk_if_in_range("wind_direction", wind_direction, -180, 180);
935
936 const Index nf = f_grid.nelem();
937
938 // Determine specular direction
939 Vector specular_los, surface_normal;
940 specular_losCalcNoTopography(specular_los,
941 surface_normal,
942 rtp_pos,
943 rtp_los,
944 atmosphere_dim,
945 verbosity);
946
947 // Determine relative azimuth
948 //
949 // According to email from James Hocking, UkMet:
950 // All azimuth angles are defined as being measured clockwise from north
951 // (i.e. if the projection onto the Earth's surface of the view path lies due
952 // north the azimuth angle is 0 and if it lies due east the azimuth angle is
953 // 90 degrees). The relative azimuth is the wind direction (azimuth) angle
954 // minus the satellite azimuth angle.
955 Numeric rel_azimuth = wind_direction; // Always true for 1D
956 if (atmosphere_dim == 2 && rtp_los[0] < 0) {
957 rel_azimuth -= 180;
958 resolve_lon(rel_azimuth, -180, 180);
959 } else if (atmosphere_dim == 3) {
960 rel_azimuth -= rtp_los[1];
961 resolve_lon(rel_azimuth, -180, 180);
962 }
963
964 // Call FASTEM
965 Matrix emissivity, reflectivity;
966 FastemStandAlone(emissivity,
967 reflectivity,
968 f_grid,
969 surface_skin_t,
970 abs(rtp_los[0]),
971 salinity,
972 wind_speed,
973 rel_azimuth,
974 transmittance,
975 fastem_version,
976 verbosity);
977
978 // Set surface_los
979 surface_los.resize(1, specular_los.nelem());
980 surface_los(0, joker) = specular_los;
981
982 // Surface emission
983 //
984 Vector b(nf);
985 planck(b, f_grid, surface_skin_t);
986 //
987 surface_emission.resize(nf, stokes_dim);
988 for (Index i = 0; i < nf; i++) {
989 // I
990 surface_emission(i, 0) = b[i] * 0.5 * (emissivity(i, 0) + emissivity(i, 1));
991 // Q
992 if (stokes_dim >= 2) {
993 surface_emission(i, 1) =
994 b[i] * 0.5 * (emissivity(i, 0) - emissivity(i, 1));
995 }
996 // U and V
997 for (Index j = 2; j < stokes_dim; j++) {
998 surface_emission(i, j) = b[i] * emissivity(i, j);
999 }
1000 }
1001
1002 // Surface reflectivity matrix
1003 //
1004 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
1005 surface_rmatrix = 0.0;
1006 for (Index iv = 0; iv < nf; iv++) {
1007 surface_rmatrix(0, iv, 0, 0) =
1008 0.5 * (reflectivity(iv, 0) + reflectivity(iv, 1));
1009 if (stokes_dim >= 2) {
1010 surface_rmatrix(0, iv, 0, 1) =
1011 0.5 * (reflectivity(iv, 0) - reflectivity(iv, 1));
1012 ;
1013 surface_rmatrix(0, iv, 1, 0) = surface_rmatrix(0, iv, 0, 1);
1014 surface_rmatrix(0, iv, 1, 1) = surface_rmatrix(0, iv, 0, 0);
1015
1016 for (Index i = 2; i < stokes_dim; i++) {
1017 surface_rmatrix(0, iv, i, i) = surface_rmatrix(0, iv, 0, 0);
1018 }
1019 }
1020 }
1021}
1022
1023/* Workspace method: Doxygen documentation will be auto-generated */
1025 Tensor4& surface_rmatrix,
1026 const Index& stokes_dim,
1027 const Numeric& pol_angle,
1028 const Verbosity& ) {
1029 ARTS_USER_ERROR_IF (stokes_dim != 1,
1030 "You should only use this method where the main calculations are "
1031 "done with *stokes_dim* set to 1.");
1032
1033 Index local_stokes = surface_emission.ncols();
1034 ARTS_USER_ERROR_IF (local_stokes < 2,
1035 "This method requires that the input surface proporties match a Stokes "
1036 "dimension of 2-4. Incoming *surface_emission* matches stokes_dim=1.");
1037
1038 const Index nf = surface_emission.nrows();
1039 const Index nlos = surface_rmatrix.nbooks();
1040 Matrix se = surface_emission;
1041 Tensor4 sr = surface_rmatrix;
1042 surface_emission.resize(nf, 1);
1043 surface_rmatrix.resize(nlos, nf, 1, 1);
1044
1045 if (local_stokes == 2) {
1046 const Numeric alpha = DEG2RAD * pol_angle;
1047 const Numeric c2 = pow( cos(alpha), 2 );
1048 const Numeric s2 = pow( sin(alpha), 2 );
1049 for (Index f=0; f<nf; ++f) {
1050 const Numeric bv = se(f,0) + se(f,1); // Note that we have to multiply with 2
1051 const Numeric bh = se(f,0) - se(f,1); // as we place a single pol as I
1052 surface_emission(f,0) = c2*bv + s2*bh;
1053 for (Index l=0; l<nlos; ++l) {
1054 const Numeric rv = sr(l,f,0,0) + sr(l,f,1,0);
1055 const Numeric rh = sr(l,f,0,0) - sr(l,f,1,0);
1056 surface_rmatrix(l,f,0,0) = c2*rv + s2*rh;
1057 }
1058 }
1059 } else {
1060 Matrix Cm, Cs, Lp, Lm;
1061 mueller_stokes2modif(Cm, local_stokes);
1062 mueller_modif2stokes(Cs, local_stokes);
1063 mueller_rotation(Lp, local_stokes, pol_angle);
1064 mueller_rotation(Lm, local_stokes, -pol_angle);
1065 Matrix Mleft(local_stokes,local_stokes), Mright(local_stokes,local_stokes);
1066 mult(Mleft, Cm, Lp);
1067 mult(Mright, Lm, Cs);
1068 //
1069 Vector Vr(local_stokes);
1070 Matrix Tmp(local_stokes,local_stokes), Mr(local_stokes,local_stokes);
1071 //
1072 for (Index f=0; f<nf; ++f) {
1073 mult(Vr, Mleft, se(f,joker) ); // Note that we have to multiply with 2
1074 surface_emission(f,0) = 2.0 * Vr[0]; // as we place a single pol as I
1075 for (Index l=0; l<nlos; ++l) {
1076 mult(Tmp, sr(l,f,joker,joker), Mright);
1077 mult(Mr, Mleft, Tmp);
1078 surface_rmatrix(l,f,0,0) = Tmp(0,0);
1079 }
1080 }
1081 }
1082}
1083
1084/* Workspace method: Doxygen documentation will be auto-generated */
1085void surfaceTelsem(Matrix& surface_los,
1086 Tensor4& surface_rmatrix,
1087 Matrix& surface_emission,
1088 const Index& atmosphere_dim,
1089 const Index& stokes_dim,
1090 const Vector& f_grid,
1091 const Vector& lat_grid,
1092 const Vector& lat_true,
1093 const Vector& lon_true,
1094 const Vector& rtp_pos,
1095 const Vector& rtp_los,
1096 const Vector& specular_los,
1097 const Numeric& surface_skin_t,
1098 const TelsemAtlas& atlas,
1099 const Numeric& r_min,
1100 const Numeric& r_max,
1101 const Numeric& d_max,
1102 const Verbosity& verbosity) {
1103 // Input checks
1104 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1105 chk_rte_pos(atmosphere_dim, rtp_pos);
1106 chk_rte_los(atmosphere_dim, rtp_los);
1107 chk_rte_los(atmosphere_dim, specular_los);
1109 "surface skin temperature", surface_skin_t, 190.0, 373.0);
1110
1111 const Index nf = f_grid.nelem();
1112 Matrix surface_rv_rh(nf, 2);
1113
1114 // Lookup cell in atlas.
1115
1116 Numeric lat, lon;
1118 lat, lon, atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
1119 chk_if_in_range("Latitude input to TELSEM2", lat, -90.0, 90.0);
1120
1121 lon = fmod(lon, 360.0);
1122 if (lon < 0) {
1123 lon += 360.0;
1124 }
1125 chk_if_in_range("Longitude input to TELSEM2", lon, 0.0, 360.0);
1126
1127 Index cellnumber = atlas.calc_cellnum(lat, lon);
1128 // Check if cell is in atlas.
1129 if (!atlas.contains(cellnumber)) {
1130 if (d_max <= 0.0) {
1131 throw std::runtime_error(
1132 "Given coordinates are not contained in "
1133 " TELSEM atlas. To enable nearest neighbour"
1134 "interpolation set *d_max* to a positive "
1135 "value.");
1136 } else {
1137 cellnumber = atlas.calc_cellnum_nearest_neighbor(lat, lon);
1138 Numeric lat_nn, lon_nn;
1139 std::tie(lat_nn, lon_nn) = atlas.get_coordinates(cellnumber);
1140 Numeric d = sphdist(lat, lon, lat_nn, lon_nn);
1141 if (d > d_max) {
1142 std::ostringstream out{};
1143 out << "Distance of nearest neighbour exceeds provided limit (";
1144 out << d << " > " << d_max << ").";
1145 throw std::runtime_error(out.str());
1146 }
1147 }
1148 }
1149
1150 Index class1 = atlas.get_class1(cellnumber);
1151 Index class2 = atlas.get_class2(cellnumber);
1152
1153 Vector emis_h = atlas.get_emis_h(cellnumber);
1154 Vector emis_v = atlas.get_emis_v(cellnumber);
1155 Numeric e_v, e_h;
1156
1157 // Incidence angle
1158 const Numeric theta = calc_incang(rtp_los, specular_los);
1159 ARTS_ASSERT(theta <= 90);
1160
1161 for (Index i = 0; i < nf; ++i) {
1162 if (f_grid[i] < 5e9)
1163 throw std::runtime_error("Only frequency >= 5 GHz are allowed");
1164 if (f_grid[i] > 900e9)
1165 throw std::runtime_error("Only frequency <= 900 GHz are allowed");
1166
1167 // For frequencies 700 <= f <= 900 GHz use 700 GHz to avoid extrapolation
1168 // outside of TELSEM specifications.
1169 Numeric f = std::min(f_grid[i], 700e9) * 1e-9;
1170 std::tie(e_v, e_h) =
1171 atlas.emis_interp(theta, f, class1, class2, emis_v, emis_h);
1172
1173 surface_rv_rh(i, 0) = min(max(1.0 - e_v, r_min), r_max);
1174 surface_rv_rh(i, 1) = min(max(1.0 - e_h, r_min), r_max);
1175 }
1176
1177 surfaceFlatRvRh(surface_los,
1178 surface_rmatrix,
1179 surface_emission,
1180 f_grid,
1181 stokes_dim,
1182 atmosphere_dim,
1183 rtp_pos,
1184 rtp_los,
1185 specular_los,
1186 surface_skin_t,
1187 surface_rv_rh,
1188 verbosity);
1189}
1190
1191/* Workspace method: Doxygen documentation will be auto-generated */
1192void surfaceTessem(Matrix& surface_los,
1193 Tensor4& surface_rmatrix,
1194 Matrix& surface_emission,
1195 const Index& atmosphere_dim,
1196 const Index& stokes_dim,
1197 const Vector& f_grid,
1198 const Vector& rtp_pos,
1199 const Vector& rtp_los,
1200 const Numeric& surface_skin_t,
1201 const TessemNN& net_h,
1202 const TessemNN& net_v,
1203 const Numeric& salinity,
1204 const Numeric& wind_speed,
1205 const Verbosity& verbosity) {
1206 // Input checks
1207 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1208 chk_rte_pos(atmosphere_dim, rtp_pos);
1209 chk_rte_los(atmosphere_dim, rtp_los);
1211 "surface skin temperature", surface_skin_t, 260.0, 373.0);
1212 chk_if_in_range_exclude_high("salinity", salinity, 0, 1);
1213 chk_if_in_range_exclude_high("wind speed", wind_speed, 0, 100);
1214
1215 // Determine specular direction
1216 Vector specular_los, surface_normal;
1217 specular_losCalcNoTopography(specular_los,
1218 surface_normal,
1219 rtp_pos,
1220 rtp_los,
1221 atmosphere_dim,
1222 verbosity);
1223
1224 // TESSEM in and out
1225 //
1226 Vector out(2);
1227 VectorView e_h = out[Range(0, 1)];
1228 VectorView e_v = out[Range(1, 1)];
1229 //
1230 Vector in(5);
1231 in[1] = 180.0 - abs(rtp_los[0]);
1232 in[2] = wind_speed;
1233 in[3] = surface_skin_t;
1234 in[4] = salinity;
1235
1236 // Get Rv and Rh
1237 //
1238 const Index nf = f_grid.nelem();
1239 Matrix surface_rv_rh(nf, 2);
1240 //
1241 for (Index i = 0; i < nf; ++i) {
1242 if (f_grid[i] < 5e9)
1243 throw std::runtime_error("Only frequency >= 5 GHz are allowed");
1244 if (f_grid[i] > 900e9)
1245 throw std::runtime_error("Only frequency <= 900 GHz are allowed");
1246
1247 in[0] = f_grid[i];
1248
1249 tessem_prop_nn(e_h, net_h, in);
1250 tessem_prop_nn(e_v, net_v, in);
1251
1252 surface_rv_rh(i, 0) = min(max(1 - e_v[0], (Numeric)0), (Numeric)1);
1253 surface_rv_rh(i, 1) = min(max(1 - e_h[0], (Numeric)0), (Numeric)1);
1254 }
1255
1256 surfaceFlatRvRh(surface_los,
1257 surface_rmatrix,
1258 surface_emission,
1259 f_grid,
1260 stokes_dim,
1261 atmosphere_dim,
1262 rtp_pos,
1263 rtp_los,
1264 specular_los,
1265 surface_skin_t,
1266 surface_rv_rh,
1267 verbosity);
1268}
1269
1270/* Workspace method: Doxygen documentation will be auto-generated */
1272 Tensor4& surface_rmatrix,
1273 Matrix& surface_emission,
1274 const Vector& f_grid,
1275 const Index& stokes_dim,
1276 const Index& atmosphere_dim,
1277 const Vector& rtp_pos,
1278 const Vector& rtp_los,
1279 const Vector& specular_los,
1280 const Numeric& surface_skin_t,
1281 const GriddedField3& surface_complex_refr_index,
1282 const Verbosity& verbosity) {
1285
1286 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1287 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1288 chk_rte_pos(atmosphere_dim, rtp_pos);
1289 chk_rte_los(atmosphere_dim, rtp_los);
1290 chk_rte_los(atmosphere_dim, specular_los);
1291 chk_not_negative("surface_skin_t", surface_skin_t);
1292
1293 // Interpolate *surface_complex_refr_index*
1294 //
1295 const Index nf = f_grid.nelem();
1296 //
1297 Matrix n_real(nf, 1), n_imag(nf, 1);
1298 //
1299 complex_n_interp(n_real,
1300 n_imag,
1301 surface_complex_refr_index,
1302 "surface_complex_refr_index",
1303 f_grid,
1304 Vector(1, surface_skin_t));
1305
1306 out2 << " Sets variables to model a flat surface\n";
1307 out3 << " surface temperature: " << surface_skin_t << " K.\n";
1308
1309 surface_los.resize(1, specular_los.nelem());
1310 surface_los(0, joker) = specular_los;
1311
1312 surface_emission.resize(nf, stokes_dim);
1313 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
1314
1315 // Incidence angle
1316 const Numeric incang = calc_incang(rtp_los, specular_los);
1317 ARTS_ASSERT(incang <= 90);
1318
1319 // Complex (amplitude) reflection coefficients
1320 Complex Rv, Rh;
1321
1322 for (Index iv = 0; iv < nf; iv++) {
1323 // Set n2 (refractive index of surface medium)
1324 Complex n2(n_real(iv, 0), n_imag(iv, 0));
1325
1326 // Amplitude reflection coefficients
1327 fresnel(Rv, Rh, Numeric(1.0), n2, incang);
1328
1329 // Fill reflection matrix and emission vector
1330 surface_specular_R_and_b(surface_rmatrix(0, iv, joker, joker),
1331 surface_emission(iv, joker),
1332 Rv,
1333 Rh,
1334 f_grid[iv],
1335 stokes_dim,
1336 surface_skin_t);
1337 }
1338}
1339
1340/* Workspace method: Doxygen documentation will be auto-generated */
1342 Tensor4& surface_rmatrix,
1343 Matrix& surface_emission,
1344 const Vector& f_grid,
1345 const Index& stokes_dim,
1346 const Index& atmosphere_dim,
1347 const Vector& rtp_pos,
1348 const Vector& rtp_los,
1349 const Vector& specular_los,
1350 const Numeric& surface_skin_t,
1351 const Tensor3& surface_reflectivity,
1352 const Verbosity& verbosity) {
1354
1355 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1356 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1357 chk_rte_pos(atmosphere_dim, rtp_pos);
1358 chk_rte_los(atmosphere_dim, rtp_los);
1359 chk_rte_los(atmosphere_dim, specular_los);
1360 chk_not_negative("surface_skin_t", surface_skin_t);
1361
1362 const Index nf = f_grid.nelem();
1363
1364 if (surface_reflectivity.nrows() != stokes_dim &&
1365 surface_reflectivity.ncols() != stokes_dim) {
1366 ostringstream os;
1367 os << "The number of rows and columnss in *surface_reflectivity* must\n"
1368 << "match *stokes_dim*."
1369 << "\n stokes_dim : " << stokes_dim
1370 << "\n number of rows in *surface_reflectivity* : "
1371 << surface_reflectivity.nrows()
1372 << "\n number of columns in *surface_reflectivity* : "
1373 << surface_reflectivity.ncols() << "\n";
1374 throw runtime_error(os.str());
1375 }
1376
1377 if (surface_reflectivity.npages() != nf &&
1378 surface_reflectivity.npages() != 1) {
1379 ostringstream os;
1380 os << "The number of pages in *surface_reflectivity* should\n"
1381 << "match length of *f_grid* or be 1."
1382 << "\n length of *f_grid* : " << nf
1383 << "\n dimension of *surface_reflectivity* : "
1384 << surface_reflectivity.npages() << "\n";
1385 throw runtime_error(os.str());
1386 }
1387
1388 out2 << " Sets variables to model a flat surface\n";
1389
1390 surface_los.resize(1, specular_los.nelem());
1391 surface_los(0, joker) = specular_los;
1392
1393 surface_emission.resize(nf, stokes_dim);
1394 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
1395
1396 Matrix R, IR(stokes_dim, stokes_dim);
1397
1398 Vector b(nf);
1399 planck(b, f_grid, surface_skin_t);
1400
1401 Vector B(stokes_dim, 0);
1402
1403 for (Index iv = 0; iv < nf; iv++) {
1404 if (iv == 0 || surface_reflectivity.npages() > 1) {
1405 R = surface_reflectivity(iv, joker, joker);
1406 for (Index i = 0; i < stokes_dim; i++) {
1407 for (Index j = 0; j < stokes_dim; j++) {
1408 if (i == j) {
1409 IR(i, j) = 1 - R(i, j);
1410 } else {
1411 IR(i, j) = -R(i, j);
1412 }
1413 }
1414 }
1415 }
1416
1417 surface_rmatrix(0, iv, joker, joker) = R;
1418
1419 B[0] = b[iv];
1420 mult(surface_emission(iv, joker), IR, B);
1421 }
1422}
1423
1424/* Workspace method: Doxygen documentation will be auto-generated */
1425void surfaceFlatRvRh(Matrix& surface_los,
1426 Tensor4& surface_rmatrix,
1427 Matrix& surface_emission,
1428 const Vector& f_grid,
1429 const Index& stokes_dim,
1430 const Index& atmosphere_dim,
1431 const Vector& rtp_pos,
1432 const Vector& rtp_los,
1433 const Vector& specular_los,
1434 const Numeric& surface_skin_t,
1435 const Matrix& surface_rv_rh,
1436 const Verbosity&) {
1437 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1438 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1439 chk_rte_pos(atmosphere_dim, rtp_pos);
1440 chk_rte_los(atmosphere_dim, rtp_los);
1441 chk_rte_los(atmosphere_dim, specular_los);
1442 chk_not_negative("surface_skin_t", surface_skin_t);
1443
1444 const Index nf = f_grid.nelem();
1445
1446 if (surface_rv_rh.ncols() != 2) {
1447 ostringstream os;
1448 os << "The number of columns in *surface_rv_rh* must be two,\n"
1449 << "but the actual number of columns is " << surface_rv_rh.ncols()
1450 << "\n";
1451 throw runtime_error(os.str());
1452 }
1453
1454 if (surface_rv_rh.nrows() != nf && surface_rv_rh.nrows() != 1) {
1455 ostringstream os;
1456 os << "The number of rows in *surface_rv_rh* should\n"
1457 << "match length of *f_grid* or be 1."
1458 << "\n length of *f_grid* : " << nf
1459 << "\n rows in *surface_rv_rh* : " << surface_rv_rh.nrows() << "\n";
1460 throw runtime_error(os.str());
1461 }
1462
1463 if (min(surface_rv_rh) < 0 || max(surface_rv_rh) > 1) {
1464 throw runtime_error("All values in *surface_rv_rh* must be inside [0,1].");
1465 }
1466
1467 surface_los.resize(1, specular_los.nelem());
1468 surface_los(0, joker) = specular_los;
1469
1470 surface_emission.resize(nf, stokes_dim);
1471 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
1472
1473 surface_emission = 0;
1474 surface_rmatrix = 0;
1475
1476 Vector b(nf);
1477 planck(b, f_grid, surface_skin_t);
1478
1479 Numeric rmean = 0.0, rdiff = 0.0;
1480
1481 for (Index iv = 0; iv < nf; iv++) {
1482 if (iv == 0 || surface_rv_rh.nrows() > 1) {
1483 rmean = 0.5 * (surface_rv_rh(iv, 0) + surface_rv_rh(iv, 1));
1484 rdiff = 0.5 * (surface_rv_rh(iv, 0) - surface_rv_rh(iv, 1));
1485 }
1486
1487 surface_emission(iv, 0) = (1.0 - rmean) * b[iv];
1488 surface_rmatrix(0, iv, 0, 0) = rmean;
1489
1490 if (stokes_dim > 1) {
1491 surface_emission(iv, 1) = -rdiff * b[iv];
1492
1493 surface_rmatrix(0, iv, 0, 1) = rdiff;
1494 surface_rmatrix(0, iv, 1, 0) = rdiff;
1495 surface_rmatrix(0, iv, 1, 1) = rmean;
1496
1497 for (Index i = 2; i < stokes_dim; i++) {
1498 surface_rmatrix(0, iv, i, i) = rmean;
1499 }
1500 }
1501 }
1502}
1503
1504/* Workspace method: Doxygen documentation will be auto-generated */
1506 Tensor4& surface_rmatrix,
1507 Matrix& surface_emission,
1508 const Vector& f_grid,
1509 const Index& stokes_dim,
1510 const Index& atmosphere_dim,
1511 const Vector& rtp_pos,
1512 const Vector& rtp_los,
1513 const Vector& specular_los,
1514 const Numeric& surface_skin_t,
1515 const Vector& surface_scalar_reflectivity,
1516 const Verbosity&) {
1517 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1518 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1519 chk_rte_pos(atmosphere_dim, rtp_pos);
1520 chk_rte_los(atmosphere_dim, rtp_los);
1521 chk_rte_los(atmosphere_dim, specular_los);
1522 chk_not_negative("surface_skin_t", surface_skin_t);
1523
1524 const Index nf = f_grid.nelem();
1525
1526 if (surface_scalar_reflectivity.nelem() != nf &&
1527 surface_scalar_reflectivity.nelem() != 1) {
1528 ostringstream os;
1529 os << "The number of elements in *surface_scalar_reflectivity* should\n"
1530 << "match length of *f_grid* or be 1."
1531 << "\n length of *f_grid* : " << nf
1532 << "\n length of *surface_scalar_reflectivity* : "
1533 << surface_scalar_reflectivity.nelem() << "\n";
1534 throw runtime_error(os.str());
1535 }
1536
1537 if (min(surface_scalar_reflectivity) < 0 ||
1538 max(surface_scalar_reflectivity) > 1) {
1539 throw runtime_error(
1540 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
1541 }
1542
1543 surface_los.resize(1, specular_los.nelem());
1544 surface_los(0, joker) = specular_los;
1545
1546 surface_emission.resize(nf, stokes_dim);
1547 surface_rmatrix.resize(1, nf, stokes_dim, stokes_dim);
1548
1549 surface_emission = 0;
1550 surface_rmatrix = 0;
1551
1552 Vector b(nf);
1553 planck(b, f_grid, surface_skin_t);
1554
1555 Numeric r = 0.0;
1556
1557 for (Index iv = 0; iv < nf; iv++) {
1558 if (iv == 0 || surface_scalar_reflectivity.nelem() > 1) {
1559 r = surface_scalar_reflectivity[iv];
1560 }
1561
1562 surface_emission(iv, 0) = (1.0 - r) * b[iv];
1563 surface_rmatrix(0, iv, 0, 0) = r;
1564 for (Index i = 1; i < stokes_dim; i++) {
1565 surface_rmatrix(0, iv, i, i) = r;
1566 }
1567 }
1568}
1569
1570/* Workspace method: Doxygen documentation will be auto-generated */
1572 Tensor4& surface_rmatrix,
1573 Matrix& surface_emission,
1574 const Vector& f_grid,
1575 const Index& stokes_dim,
1576 const Index& atmosphere_dim,
1577 const Vector& rtp_pos,
1578 const Vector& rtp_los,
1579 const Vector& surface_normal,
1580 const Numeric& surface_skin_t,
1581 const Vector& surface_scalar_reflectivity,
1582 const Index& lambertian_nza,
1583 const Numeric& za_pos,
1584 const Verbosity&) {
1585 const Index nf = f_grid.nelem();
1586
1587 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1588 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1589 chk_rte_pos(atmosphere_dim, rtp_pos);
1590 chk_rte_los(atmosphere_dim, rtp_los);
1591 chk_not_negative("surface_skin_t", surface_skin_t);
1592 chk_if_in_range("za_pos", za_pos, 0, 1);
1593
1594 if (surface_scalar_reflectivity.nelem() != nf &&
1595 surface_scalar_reflectivity.nelem() != 1) {
1596 ostringstream os;
1597 os << "The number of elements in *surface_scalar_reflectivity* should\n"
1598 << "match length of *f_grid* or be 1."
1599 << "\n length of *f_grid* : " << nf
1600 << "\n length of *surface_scalar_reflectivity* : "
1601 << surface_scalar_reflectivity.nelem() << "\n";
1602 throw runtime_error(os.str());
1603 }
1604
1605 if (min(surface_scalar_reflectivity) < 0 ||
1606 max(surface_scalar_reflectivity) > 1) {
1607 throw runtime_error(
1608 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
1609 }
1610
1611 // Allocate and init everything to zero
1612 //
1613 surface_los.resize(lambertian_nza, rtp_los.nelem());
1614 surface_rmatrix.resize(lambertian_nza, nf, stokes_dim, stokes_dim);
1615 surface_emission.resize(nf, stokes_dim);
1616 //
1617 surface_los = 0.0;
1618 surface_rmatrix = 0.0;
1619 surface_emission = 0.0;
1620
1621 // Help variables
1622 //
1623 const Numeric dza = (90.0 - abs(surface_normal[0])) / (Numeric)lambertian_nza;
1624 const Vector za_lims(0.0, lambertian_nza + 1, dza);
1625
1626 // surface_los
1627 for (Index ip = 0; ip < lambertian_nza; ip++) {
1628 surface_los(ip, 0) = za_lims[ip] + za_pos * dza;
1629 if (atmosphere_dim == 2) {
1630 if (rtp_los[0] < 0) {
1631 surface_los(ip, 0) *= -1.0;
1632 }
1633
1634 } else if (atmosphere_dim == 3) {
1635 surface_los(ip, 1) = rtp_los[1];
1636 }
1637 }
1638
1639 Vector b(nf);
1640 planck(b, f_grid, surface_skin_t);
1641
1642 // Loop frequencies and set remaining values
1643 //
1644 Numeric r = 0.0;
1645 //
1646 for (Index iv = 0; iv < nf; iv++) {
1647 // Get reflectivity
1648 if (iv == 0 || surface_scalar_reflectivity.nelem() > 1) {
1649 r = surface_scalar_reflectivity[iv];
1650 }
1651
1652 // surface_rmatrix:
1653 // Only element (0,0) is set to be non-zero. This follows VDISORT
1654 // that refers to: K. L. Coulson, Polarization and Intensity of Light in
1655 // the Atmosphere (1989), page 229
1656 // (Thanks to Michael Kahnert for providing this information!)
1657 // Update: Is the above for a later edition? We have not found a copy of
1658 // that edition. In a 1988 version of the book, the relevant page seems
1659 // to be 232.
1660 for (Index ip = 0; ip < lambertian_nza; ip++) {
1661 const Numeric w =
1662 r * 0.5 *
1663 (cos(2 * DEG2RAD * za_lims[ip]) - cos(2 * DEG2RAD * za_lims[ip + 1]));
1664 surface_rmatrix(ip, iv, 0, 0) = w;
1665 }
1666
1667 // surface_emission
1668 surface_emission(iv, 0) = (1 - r) * b[iv];
1669 }
1670}
1671
1672/* Workspace method: Doxygen documentation will be auto-generated */
1674 Numeric& surface_skin_t,
1675 Matrix& surface_los,
1676 Tensor4& surface_rmatrix,
1677 Matrix& surface_emission,
1678 const Index& atmosphere_dim,
1679 const Vector& f_grid,
1680 const Vector& rtp_pos,
1681 const Vector& rtp_los,
1682 const Agenda& surface_rtprop_sub_agenda,
1683 const Numeric& specular_factor,
1684 const Numeric& dza,
1685 const Verbosity&) {
1686 chk_rte_pos(atmosphere_dim, rtp_pos);
1687 chk_rte_los(atmosphere_dim, rtp_los);
1688
1689 // Checks of GIN variables
1690 if (specular_factor > 1 || specular_factor < 1.0 / 3.0)
1691 throw runtime_error("The valid range for *specular_factor* is [1/3,1].");
1692 if (dza > 45 || dza <= 0)
1693 throw runtime_error("The valid range for *dza* is ]0,45].");
1694
1695 // Obtain data for specular direction
1696 //
1697 Matrix los1, emission1;
1698 Tensor4 rmatrix1;
1699 //
1701 surface_skin_t,
1702 emission1,
1703 los1,
1704 rmatrix1,
1705 f_grid,
1706 rtp_pos,
1707 rtp_los,
1708 surface_rtprop_sub_agenda);
1709 if (los1.nrows() != 1)
1710 throw runtime_error(
1711 "*surface_rtprop_sub_agenda* must return data "
1712 "describing a specular surface.");
1713
1714 // Standard number of beams. Set to 2 if try/catch below fails
1715 Index nbeams = 3;
1716
1717 // Test if some lower zenith angle works.
1718 // It will fail if a higher za results in looking at the surface from below.
1719 //
1720 Matrix los2, emission2;
1721 Tensor4 rmatrix2;
1722 //
1723 Numeric skin_t_dummy;
1724 Numeric dza_try = dza;
1725 bool failed = true;
1726 while (failed && dza_try > 0) {
1727 try {
1728 Vector los_new = rtp_los;
1729 los_new[0] -=
1730 sign(rtp_los[0]) * dza_try; // Sign to also handle 2D negative za
1731 adjust_los(los_new, atmosphere_dim);
1733 skin_t_dummy,
1734 emission2,
1735 los2,
1736 rmatrix2,
1737 f_grid,
1738 rtp_pos,
1739 los_new,
1740 surface_rtprop_sub_agenda);
1741 failed = false;
1742 } catch (const runtime_error& e) {
1743 dza_try -= 1.0;
1744 }
1745 }
1746 if (failed) {
1747 nbeams = 2;
1748 }
1749
1750 // Allocate output WSVs
1751 //
1752 surface_emission.resize(emission1.nrows(), emission1.ncols());
1753 surface_emission = 0;
1754 surface_los.resize(nbeams, los1.ncols());
1755 surface_rmatrix.resize(
1756 nbeams, rmatrix1.npages(), rmatrix1.nrows(), rmatrix1.ncols());
1757
1758 // Put in specular direction at index 1
1759 //
1760 Numeric w;
1761 if (nbeams == 3) {
1762 w = specular_factor;
1763 } else {
1764 w = specular_factor + (1.0 - specular_factor) / 2.0;
1765 }
1766 //
1767 surface_los(1, joker) = los1(0, joker);
1768 for (Index p = 0; p < rmatrix1.npages(); p++) {
1769 for (Index r = 0; r < rmatrix1.nrows(); r++) {
1770 surface_emission(p, r) += w * emission1(p, r);
1771 for (Index c = 0; c < rmatrix1.ncols(); c++) {
1772 surface_rmatrix(1, p, r, c) = w * rmatrix1(0, p, r, c);
1773 }
1774 }
1775 }
1776
1777 // Put in lower za as index 2, if worked
1778 //
1779 w = (1.0 - specular_factor) / 2.0;
1780 //
1781 if (nbeams == 3) {
1782 surface_los(2, joker) = los2(0, joker);
1783 for (Index p = 0; p < rmatrix2.npages(); p++) {
1784 for (Index r = 0; r < rmatrix2.nrows(); r++) {
1785 surface_emission(p, r) += w * emission2(p, r);
1786 for (Index c = 0; c < rmatrix1.ncols(); c++) {
1787 surface_rmatrix(2, p, r, c) = w * rmatrix2(0, p, r, c);
1788 }
1789 }
1790 }
1791 }
1792
1793 // Do higher za and put in as index 0 (reusing variables for beam 2)
1794 //
1795 Vector los_new = rtp_los;
1796 los_new[0] += sign(rtp_los[0]) * dza; // Sign to also handle 2D negative za
1797 adjust_los(los_new, atmosphere_dim);
1799 skin_t_dummy,
1800 emission2,
1801 los2,
1802 rmatrix2,
1803 f_grid,
1804 rtp_pos,
1805 los_new,
1806 surface_rtprop_sub_agenda);
1807 //
1808 surface_los(0, joker) = los2(0, joker);
1809 for (Index p = 0; p < rmatrix2.npages(); p++) {
1810 for (Index r = 0; r < rmatrix2.nrows(); r++) {
1811 surface_emission(p, r) += w * emission2(p, r);
1812 for (Index c = 0; c < rmatrix1.ncols(); c++) {
1813 surface_rmatrix(0, p, r, c) = w * rmatrix2(0, p, r, c);
1814 }
1815 }
1816 }
1817}
1818
1819/* Workspace method: Doxygen documentation will be auto-generated */
1821 Tensor4& surface_rmatrix,
1822 const Index& atmosphere_dim,
1823 const Vector& rtp_pos,
1824 const Vector& rtp_los,
1825 const Numeric& specular_factor,
1826 const Numeric& dza,
1827 const Verbosity&) {
1828 chk_rte_pos(atmosphere_dim, rtp_pos);
1829 chk_rte_los(atmosphere_dim, rtp_los);
1830
1831 // Check that input surface data are of specular type
1832 if (surface_los.nrows() != 1)
1833 throw runtime_error(
1834 "Input surface data must be of specular type. That is, "
1835 "*surface_los* must contain a single direction.");
1836 if (surface_rmatrix.nbooks() != 1)
1837 throw runtime_error(
1838 "*surface_rmatrix* describes a different number of "
1839 "directions than *surface_los*.");
1840
1841 // Checks of GIN variables
1842 if (specular_factor > 1 || specular_factor < 1.0 / 3.0)
1843 throw runtime_error("The valid range for *specular_factor* is [1/3,1].");
1844 if (dza > 45 || dza <= 0)
1845 throw runtime_error("The valid range for *dza* is ]0,45].");
1846
1847 // Make copies of input data
1848 const Matrix los1 = surface_los;
1849 const Tensor4 rmatrix1 = surface_rmatrix;
1850
1851 // Use abs(za) in all expressions below, to also handle 2D
1852
1853 // Calculate highest possible za for downwelling radiation, with 1 degree
1854 // margin to the surface. This can be derived from za in rtp_los and surface_los.
1855 // (The directions in surface_los are not allowed to point into the surface)
1856 const Numeric za_max = 89 + (180 - abs(los1(0, 0)) - abs(rtp_los[0])) / 2.0;
1857
1858 // Number of downwelling beams
1859 Index nbeams = 3;
1860 if (abs(los1(0, 0)) > za_max) {
1861 nbeams = 2;
1862 }
1863
1864 // New los-s
1865 //
1866 surface_los.resize(nbeams, los1.ncols());
1867 //
1868 for (Index r = 0; r < nbeams; r++) {
1869 surface_los(r, 0) = ((Numeric)r - 1.0) * dza + abs(los1(0, 0));
1870 if (r == 2 && surface_los(r, 0) > za_max) {
1871 surface_los(r, 0) = za_max;
1872 }
1873 for (Index c = 1; c < los1.ncols(); c++) {
1874 surface_los(r, c) = los1(0, c);
1875 }
1876 }
1877
1878 // New rmatrix
1879 //
1880 surface_rmatrix.resize(
1881 nbeams, rmatrix1.npages(), rmatrix1.nrows(), rmatrix1.ncols());
1882 //
1883 for (Index b = 0; b < nbeams; b++) {
1884 Numeric w;
1885 if (b == 1 && nbeams == 3) // Specular direction with nbeams==3
1886 {
1887 w = specular_factor;
1888 } else if (b == 1) // Specular direction with nbeams==2
1889 {
1890 w = specular_factor + (1.0 - specular_factor) / 2.0;
1891 } else // Side directions
1892 {
1893 w = (1.0 - specular_factor) / 2.0;
1894 }
1895
1896 for (Index p = 0; p < rmatrix1.npages(); p++) {
1897 for (Index r = 0; r < rmatrix1.nrows(); r++) {
1898 for (Index c = 0; c < rmatrix1.ncols(); c++) {
1899 surface_rmatrix(b, p, r, c) = w * rmatrix1(0, p, r, c);
1900 }
1901 }
1902 }
1903 }
1904
1905 // Handle sign of za
1906 if (atmosphere_dim == 1) {
1907 // We only need to make sure that first direction has positive za
1908 surface_los(0, 0) = abs(surface_los(0, 0));
1909 } else if (atmosphere_dim == 2) {
1910 // Change sign if specular direction has za < 0
1911 if (los1(0, 0) < 0) {
1912 for (Index r = 0; r < rmatrix1.nrows(); r++) {
1913 surface_los(r, 0) = -surface_los(r, 0);
1914 }
1915 }
1916 } else if (atmosphere_dim == 1) {
1917 // We only need to make sure that first direction has positive za
1918 if (surface_los(0, 0) < 0) {
1919 surface_los(0, 0) = -surface_los(0, 0);
1920 surface_los(0, 1) += 180;
1921 if (surface_los(0, 1) > 180) {
1922 surface_los(0, 1) -= 360;
1923 }
1924 }
1925 }
1926}
1927
1928/* Workspace method: Doxygen documentation will be auto-generated */
1930 GriddedField3& surface_complex_refr_index,
1931 const Index& atmosphere_dim,
1932 const Vector& lat_grid,
1933 const Vector& lat_true,
1934 const Vector& lon_true,
1935 const Vector& rtp_pos,
1936 const GriddedField5& complex_n_field,
1937 const Verbosity&) {
1938 // Set expected order of grids
1939 Index gfield_fID = 0;
1940 Index gfield_tID = 1;
1941 Index gfield_compID = 2;
1942 Index gfield_latID = 3;
1943 Index gfield_lonID = 4;
1944
1945 // Basic checks and sizes
1946 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1947 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
1948 chk_rte_pos(atmosphere_dim, rtp_pos);
1949 complex_n_field.checksize_strict();
1950 //
1951 chk_griddedfield_gridname(complex_n_field, gfield_fID, "Frequency");
1952 chk_griddedfield_gridname(complex_n_field, gfield_tID, "Temperature");
1953 chk_griddedfield_gridname(complex_n_field, gfield_compID, "Complex");
1954 chk_griddedfield_gridname(complex_n_field, gfield_latID, "Latitude");
1955 chk_griddedfield_gridname(complex_n_field, gfield_lonID, "Longitude");
1956 //
1957 const Index nf = complex_n_field.data.nshelves();
1958 const Index nt = complex_n_field.data.nbooks();
1959 const Index nn = complex_n_field.data.npages();
1960 const Index nlat = complex_n_field.data.nrows();
1961 const Index nlon = complex_n_field.data.ncols();
1962 //
1963 if (nlat < 2 || nlon < 2) {
1964 ostringstream os;
1965 os << "The data in *complex_refr_index_field* must span a geographical "
1966 << "region. That is,\nthe latitude and longitude grids must have a "
1967 << "length >= 2.";
1968 throw runtime_error(os.str());
1969 }
1970 //
1971 if (nn != 2) {
1972 ostringstream os;
1973 os << "The data in *complex_refr_index_field* must have exactly two "
1974 << "pages. One page each\nfor the real and imaginary part of the "
1975 << "complex refractive index.";
1976 throw runtime_error(os.str());
1977 }
1978
1979 const Vector& GFlat = complex_n_field.get_numeric_grid(gfield_latID);
1980 const Vector& GFlon = complex_n_field.get_numeric_grid(gfield_lonID);
1981
1982 // Determine true geographical position
1983 Vector lat(1), lon(1);
1985 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
1986
1987 // Ensure correct coverage of lon grid
1988 Vector lon_shifted;
1989 lon_shiftgrid(lon_shifted, GFlon, lon[0]);
1990
1991 // Check if lat/lon we need are actually covered
1992 chk_if_in_range("rtp_pos.lat", lat[0], GFlat[0], GFlat[nlat - 1]);
1993 chk_if_in_range("rtp_pos.lon", lon[0], lon_shifted[0], lon_shifted[nlon - 1]);
1994
1995 // Size and fills grids of *surface_complex_refr_index*
1996 surface_complex_refr_index.resize(nf, nt, 2);
1997 surface_complex_refr_index.set_grid_name(0, "Frequency");
1998 surface_complex_refr_index.set_grid(
1999 0, complex_n_field.get_numeric_grid(gfield_fID));
2000 surface_complex_refr_index.set_grid_name(1, "Temperature");
2001 surface_complex_refr_index.set_grid(
2002 1, complex_n_field.get_numeric_grid(gfield_tID));
2003 surface_complex_refr_index.set_grid_name(2, "Complex");
2004 surface_complex_refr_index.set_grid(2, {"real", "imaginary"});
2005
2006 // Interpolate in lat and lon
2007 //
2008 GridPos gp_lat, gp_lon;
2009 gridpos(gp_lat, GFlat, lat[0]);
2010 gridpos(gp_lon, lon_shifted, lon[0]);
2011 Vector itw(4);
2012 interpweights(itw, gp_lat, gp_lon);
2013 //
2014 for (Index iv = 0; iv < nf; iv++) {
2015 for (Index it = 0; it < nt; it++) {
2016 surface_complex_refr_index.data(iv, it, 0) = interp(
2017 itw, complex_n_field.data(iv, it, 0, joker, joker), gp_lat, gp_lon);
2018 surface_complex_refr_index.data(iv, it, 1) = interp(
2019 itw, complex_n_field.data(iv, it, 1, joker, joker), gp_lat, gp_lon);
2020 }
2021 }
2022}
2023
2024/* Workspace method: Doxygen documentation will be auto-generated */
2026 const Index& stokes_dim,
2027 const Vector& f_grid,
2028 const Index& atmosphere_dim,
2029 const Vector& lat_grid,
2030 const Vector& lat_true,
2031 const Vector& lon_true,
2032 const Vector& rtp_pos,
2033 const Vector& rtp_los,
2034 const GriddedField6& r_field,
2035 const Verbosity&) {
2036 // Basic checks and sizes
2037 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2038 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2039 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
2040 chk_rte_pos(atmosphere_dim, rtp_pos);
2041 chk_rte_los(atmosphere_dim, rtp_los);
2042 r_field.checksize_strict();
2043 chk_griddedfield_gridname(r_field, 0, "Frequency");
2044 chk_griddedfield_gridname(r_field, 1, "Stokes element");
2045 chk_griddedfield_gridname(r_field, 2, "Stokes element");
2046 chk_griddedfield_gridname(r_field, 3, "Incidence angle");
2047 chk_griddedfield_gridname(r_field, 4, "Latitude");
2048 chk_griddedfield_gridname(r_field, 5, "Longitude");
2049 //
2050 const Index nf_in = r_field.data.nvitrines();
2051 const Index ns2 = r_field.data.nshelves();
2052 const Index ns1 = r_field.data.nbooks();
2053 const Index nza = r_field.data.npages();
2054 const Index nlat = r_field.data.nrows();
2055 const Index nlon = r_field.data.ncols();
2056 //
2057 if (nlat < 2 || nlon < 2) {
2058 ostringstream os;
2059 os << "The data in *r_field* must span a geographical region. That is,\n"
2060 << "the latitude and longitude grids must have a length >= 2.";
2061 throw runtime_error(os.str());
2062 }
2063 //
2064 if (nza < 2) {
2065 ostringstream os;
2066 os << "The data in *r_field* must span a range of zenith angles. That\n"
2067 << "is the zenith angle grid must have a length >= 2.";
2068 throw runtime_error(os.str());
2069 }
2070 if (ns1 < stokes_dim || ns2 < stokes_dim || ns1 > 4 || ns2 > 4) {
2071 ostringstream os;
2072 os << "The \"Stokes dimensions\" must have a size that is >= "
2073 << "*stokes_dim* (but not exceeding 4).";
2074 throw runtime_error(os.str());
2075 }
2076
2077 // Determine true geographical position
2078 Vector lat(1), lon(1);
2080 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
2081
2082 // Ensure correct coverage of lon grid
2083 Vector lon_shifted;
2084 lon_shiftgrid(lon_shifted, r_field.get_numeric_grid(5), lon[0]);
2085
2086 // Interpolate in lat and lon
2087 Tensor4 r_f_za(nf_in, stokes_dim, stokes_dim, nza);
2088 {
2090 "Latitude interpolation", r_field.get_numeric_grid(4), lat[0]);
2091 chk_interpolation_grids("Longitude interpolation", lon_shifted, lon[0]);
2092 GridPos gp_lat, gp_lon;
2093 gridpos(gp_lat, r_field.get_numeric_grid(4), lat[0]);
2094 gridpos(gp_lon, lon_shifted, lon[0]);
2095 Vector itw(4);
2096 interpweights(itw, gp_lat, gp_lon);
2097 for (Index iv = 0; iv < nf_in; iv++) {
2098 for (Index iz = 0; iz < nza; iz++) {
2099 for (Index is1 = 0; is1 < stokes_dim; is1++) {
2100 for (Index is2 = 0; is2 < stokes_dim; is2++) {
2101 r_f_za(iv, is1, is2, iz) =
2102 interp(itw,
2103 r_field.data(iv, is1, is2, iz, joker, joker),
2104 gp_lat,
2105 gp_lon);
2106 }
2107 }
2108 }
2109 }
2110 }
2111
2112 // Interpolate in incidence angle, cubic if possible
2113 Tensor3 r_f(nf_in, stokes_dim, stokes_dim);
2114 Index order = 3;
2115 if (nza < 4) {
2116 order = 1;
2117 }
2118 {
2120 "Incidence angle interpolation", r_field.get_numeric_grid(3), 180 - rtp_los[0]);
2121 const LagrangeInterpolation lag(0, 180 - rtp_los[0], r_field.get_numeric_grid(3), order);
2122 const auto itw = interpweights(lag);
2123 //
2124 for (Index i = 0; i < nf_in; i++) {
2125 for (Index is1 = 0; is1 < stokes_dim; is1++) {
2126 for (Index is2 = 0; is2 < stokes_dim; is2++) {
2127 r_f(i, is1, is2) = interp(r_f_za(i, is1, is2, joker), itw, lag);
2128 }
2129 }
2130 }
2131 }
2132
2133 // Extract or interpolate in frequency
2134 //
2135 if (nf_in == 1) {
2136 surface_reflectivity = r_f;
2137 } else {
2139 "Frequency interpolation", r_field.get_numeric_grid(0), f_grid);
2140 const Index nf_out = f_grid.nelem();
2141 surface_reflectivity.resize(nf_out, stokes_dim, stokes_dim);
2142 //
2143 ArrayOfGridPos gp(nf_out);
2144 Matrix itw(nf_out, 2);
2145 gridpos(gp, r_field.get_numeric_grid(0), f_grid);
2146 interpweights(itw, gp);
2147 for (Index is1 = 0; is1 < stokes_dim; is1++) {
2148 for (Index is2 = 0; is2 < stokes_dim; is2++) {
2149 interp(surface_reflectivity(joker, is1, is2),
2150 itw,
2151 r_f(joker, is1, is2),
2152 gp);
2153 }
2154 }
2155 }
2156}
2157
2158/* Workspace method: Doxygen documentation will be auto-generated */
2160 Vector& surface_scalar_reflectivity,
2161 const Index& stokes_dim,
2162 const Vector& f_grid,
2163 const Index& atmosphere_dim,
2164 const Vector& lat_grid,
2165 const Vector& lat_true,
2166 const Vector& lon_true,
2167 const Vector& rtp_pos,
2168 const Vector& rtp_los,
2169 const GriddedField4& r_field,
2170 const Verbosity&) {
2171 // Basic checks and sizes
2172 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2173 chk_if_in_range("stokes_dim", stokes_dim, 1, 1);
2174 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
2175 chk_rte_pos(atmosphere_dim, rtp_pos);
2176 chk_rte_los(atmosphere_dim, rtp_los);
2177 r_field.checksize_strict();
2178 chk_griddedfield_gridname(r_field, 0, "Frequency");
2179 chk_griddedfield_gridname(r_field, 1, "Incidence angle");
2180 chk_griddedfield_gridname(r_field, 2, "Latitude");
2181 chk_griddedfield_gridname(r_field, 3, "Longitude");
2182 //
2183 const Index nf_in = r_field.data.nbooks();
2184 const Index nza = r_field.data.npages();
2185 const Index nlat = r_field.data.nrows();
2186 const Index nlon = r_field.data.ncols();
2187 //
2188 if (nlat < 2 || nlon < 2) {
2189 ostringstream os;
2190 os << "The data in *r_field* must span a geographical region. That is,\n"
2191 << "the latitude and longitude grids must have a length >= 2.";
2192 throw runtime_error(os.str());
2193 }
2194 //
2195 if (nza < 2) {
2196 ostringstream os;
2197 os << "The data in *r_field* must span a range of zenith angles. That\n"
2198 << "is the zenith angle grid must have a length >= 2.";
2199 throw runtime_error(os.str());
2200 }
2201
2202 // Determine true geographical position
2203 Vector lat(1), lon(1);
2205 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
2206
2207 // Ensure correct coverage of lon grid
2208 Vector lon_shifted;
2209 lon_shiftgrid(lon_shifted, r_field.get_numeric_grid(3), lon[0]);
2210
2211 // Interpolate in lat and lon
2212 Matrix r_f_za(nf_in, nza);
2213 {
2215 "Latitude interpolation", r_field.get_numeric_grid(2), lat[0]);
2216 chk_interpolation_grids("Longitude interpolation", lon_shifted, lon[0]);
2217 GridPos gp_lat, gp_lon;
2218 gridpos(gp_lat, r_field.get_numeric_grid(2), lat[0]);
2219 gridpos(gp_lon, lon_shifted, lon[0]);
2220 Vector itw(4);
2221 interpweights(itw, gp_lat, gp_lon);
2222 for (Index iv = 0; iv < nf_in; iv++) {
2223 for (Index iz = 0; iz < nza; iz++) {
2224 r_f_za(iv, iz) =
2225 interp(itw, r_field.data(iv, iz, joker, joker), gp_lat, gp_lon);
2226 }
2227 }
2228 }
2229
2230 // Interpolate in incidence angle, cubic if possible
2231 Vector r_f(nf_in);
2232 Index order = 3;
2233 if (nza < 4) {
2234 order = 1;
2235 }
2236 {
2238 "Incidence angle interpolation", r_field.get_numeric_grid(1), 180 - rtp_los[0]);
2239 const LagrangeInterpolation lag(0, 180 - rtp_los[0], r_field.get_numeric_grid(1), order);
2240 const auto itw=interpweights(lag);
2241 for (Index i = 0; i < nf_in; i++) {
2242 r_f[i] = interp(r_f_za(i, joker), itw, lag);
2243 }
2244 }
2245
2246 // Extract or interpolate in frequency
2247 //
2248 if (nf_in == 1) {
2249 surface_scalar_reflectivity.resize(1);
2250 surface_scalar_reflectivity[0] = r_f[0];
2251 } else {
2253 "Frequency interpolation", r_field.get_numeric_grid(0), f_grid);
2254 const Index nf_out = f_grid.nelem();
2255 surface_scalar_reflectivity.resize(nf_out);
2256 //
2257 ArrayOfGridPos gp(nf_out);
2258 Matrix itw(nf_out, 2);
2259 gridpos(gp, r_field.get_numeric_grid(0), f_grid);
2260 interpweights(itw, gp);
2261 interp(surface_scalar_reflectivity, itw, r_f, gp);
2262 }
2263}
2264
2265/* Workspace method: Doxygen documentation will be auto-generated */
2267 Vector& surface_scalar_reflectivity,
2268 const Tensor4& surface_rmatrix,
2269 const Verbosity&) {
2270 const Index nf = surface_rmatrix.npages();
2271 const Index nlos = surface_rmatrix.nbooks();
2272
2273 surface_scalar_reflectivity.resize(nf);
2274 surface_scalar_reflectivity = 0;
2275
2276 for (Index i = 0; i < nf; i++) {
2277 for (Index l = 0; l < nlos; l++) {
2278 surface_scalar_reflectivity[i] += surface_rmatrix(l, i, 0, 0);
2279 }
2280 }
2281}
2282
2283/* Workspace method: Doxygen documentation will be auto-generated */
2285 Vector& surface_types_aux,
2286 Vector& surface_types_weights,
2287 const Index& atmosphere_dim,
2288 const Vector& lat_grid,
2289 const Vector& lat_true,
2290 const Vector& lon_true,
2291 const Vector& rtp_pos,
2292 const GriddedField2& surface_type_mask,
2293 const String& method,
2294 const Verbosity&) {
2295 // Set expected order of grids
2296 Index gfield_latID = 0;
2297 Index gfield_lonID = 1;
2298
2299 // Basic checks and sizes
2300 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2301 chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
2302 chk_rte_pos(atmosphere_dim, rtp_pos);
2303 surface_type_mask.checksize_strict();
2304 //
2305 chk_griddedfield_gridname(surface_type_mask, gfield_latID, "Latitude");
2306 chk_griddedfield_gridname(surface_type_mask, gfield_lonID, "Longitude");
2307 //
2308 const Index nlat = surface_type_mask.data.nrows();
2309 const Index nlon = surface_type_mask.data.ncols();
2310 //
2311 if (nlat < 2 || nlon < 2) {
2312 ostringstream os;
2313 os << "The data in *surface_type_mask* must span a geographical "
2314 << "region. Accordingly,\nthe latitude and longitude grids must "
2315 << "both have a length >= 2.";
2316 throw runtime_error(os.str());
2317 }
2318
2319 const Vector& GFlat = surface_type_mask.get_numeric_grid(gfield_latID);
2320 const Vector& GFlon = surface_type_mask.get_numeric_grid(gfield_lonID);
2321
2322 // Determine true geographical position
2323 Vector lat(1), lon(1);
2325 lat[0], lon[0], atmosphere_dim, lat_grid, lat_true, lon_true, rtp_pos);
2326
2327 // Ensure correct coverage of lon grid
2328 Vector lon_shifted;
2329 lon_shiftgrid(lon_shifted, GFlon, lon[0]);
2330
2331 // Check if lat/lon we need are actually covered
2332 chk_if_in_range("rtp_pos.lat", lat[0], GFlat[0], GFlat[nlat - 1]);
2333 chk_if_in_range("rtp_pos.lon", lon[0], lon_shifted[0], lon_shifted[nlon - 1]);
2334
2335 // Grid positions
2336 GridPos gp_lat, gp_lon;
2337 gridpos(gp_lat, GFlat, lat[0]);
2338 gridpos(gp_lon, lon_shifted, lon[0]);
2339
2340 if (method == "nearest" ) {
2341 // Extract closest point
2342 Index ilat, ilon;
2343 if (gp_lat.fd[0] < 0.5) {
2344 ilat = gp_lat.idx;
2345 } else {
2346 ilat = gp_lat.idx + 1;
2347 }
2348 if (gp_lon.fd[0] < 0.5) {
2349 ilon = gp_lon.idx;
2350 } else {
2351 ilon = gp_lon.idx + 1;
2352 }
2353 //
2354 surface_types.resize(1);
2355 surface_types_aux.resize(1);
2356 surface_types_weights.resize(1);
2357 surface_types[0] = (Index)floor(surface_type_mask.data(ilat, ilon));
2358 surface_types_aux[0] = surface_type_mask.data(ilat, ilon) -
2359 Numeric(surface_types[0]);
2360 surface_types_weights[0] = 1.0;
2361 }
2362
2363 else if (method == "linear" ) {
2364 // Determine types involved
2365 ArrayOfIndex types0(4), types(4);
2366 Index ntypes = 1;
2367 //
2368 types0[0] = (Index)floor(surface_type_mask.data(gp_lat.idx,gp_lon.idx));
2369 types0[1] = (Index)floor(surface_type_mask.data(gp_lat.idx,gp_lon.idx+1));
2370 types0[2] = (Index)floor(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx));
2371 types0[3] = (Index)floor(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx+1));
2372 //
2373 types[0] = types0[0];
2374 for (Index t=1; t<4; t++) {
2375 bool unique = true;
2376 for (Index n=0; n<ntypes && unique; n++)
2377 if (types0[t] == types[n] ) { unique = false; }
2378 if (unique) {
2379 types[ntypes] = types0[t];
2380 ntypes += 1;
2381 }
2382 }
2383 //
2384 surface_types.resize(ntypes);
2385 surface_types_aux.resize(ntypes);
2386 surface_types_weights.resize(ntypes);
2387
2388 // Interpolation weights
2389 Vector itw(4);
2390 interpweights(itw, gp_lat, gp_lon);
2391
2392 // Determine weight for each type, and make an interpolation of aux value
2393 for (Index t=0; t<ntypes; t++ ) {
2394 Numeric wtot = 0, auxtot = 0;
2395 Numeric ntype = (Numeric)types[t];
2396 if (types[t] == types0[0]) { wtot += itw[0];
2397 auxtot += itw[0]*(surface_type_mask.data(gp_lat.idx,gp_lon.idx)-ntype);
2398 };
2399 if (types[t] == types0[1]) { wtot += itw[1];
2400 auxtot += itw[1]*(surface_type_mask.data(gp_lat.idx,gp_lon.idx+1)-ntype);
2401 };
2402 if (types[t] == types0[2]) { wtot += itw[2];
2403 auxtot += itw[2]*(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx)-ntype);
2404 };
2405 if (types[t] == types0[3]) { wtot += itw[3];
2406 auxtot += itw[3]*(surface_type_mask.data(gp_lat.idx+1,gp_lon.idx+1)-ntype);
2407 };
2408 //
2409 surface_types[t] = types[t];
2410 surface_types_aux[t] = auxtot / wtot;
2411 surface_types_weights[t] = wtot;
2412 }
2413 }
2414
2415 else {
2416 ostringstream os;
2417 os << "The allowed options for *method are: \"nearest\" and \"linear\",\n"
2418 << "but you have selected: " << method;
2419 throw runtime_error(os.str());
2420 }
2421}
2422
2423/* Workspace method: Doxygen documentation will be auto-generated */
2425 Numeric& surface_skin_t,
2426 Matrix& surface_los,
2427 Tensor4& surface_rmatrix,
2428 Matrix& surface_emission,
2429 const Vector& f_grid,
2430 const Vector& rtp_pos,
2431 const Vector& rtp_los,
2432 const ArrayOfAgenda& surface_rtprop_agenda_array,
2433 const ArrayOfIndex& surface_types,
2434 const Vector& surface_types_aux,
2435 const Vector& surface_types_weights,
2436 const Verbosity&) {
2437 if (surface_types.nelem() != 1) {
2438 ostringstream os;
2439 os << "This method requires that *surface_types* have length 1.\n"
2440 << "If you are trying to apply a mixture model, you have to use\n"
2441 << "a setup based on *iySurfaceCallAgendaX* instead.";
2442 throw runtime_error(os.str());
2443 }
2444 if (surface_types[0] < 0)
2445 throw runtime_error("*surface_type* is not allowed to be negative.");
2446 if (surface_types[0] >= surface_rtprop_agenda_array.nelem()) {
2447 ostringstream os;
2448 os << "*surface_rtprop_agenda_array* has only "
2449 << surface_rtprop_agenda_array.nelem()
2450 << " elements,\n while you have selected element " << surface_types[0];
2451 throw runtime_error(os.str());
2452 }
2453 if (abs(surface_types_weights[0]-1)>1e-4)
2454 throw runtime_error("Sum of *surface_types_weights* deviates from 1.");
2455
2457 surface_skin_t,
2458 surface_emission,
2459 surface_los,
2460 surface_rmatrix,
2461 surface_types[0],
2462 f_grid,
2463 rtp_pos,
2464 rtp_los,
2465 surface_types_aux[0],
2466 surface_rtprop_agenda_array);
2467}
2468
2469/* Workspace method: Doxygen documentation will be auto-generated */
2470void SurfaceBlackbody(Matrix& surface_los,
2471 Tensor4& surface_rmatrix,
2472 ArrayOfTensor4& dsurface_rmatrix_dx,
2473 Matrix& surface_emission,
2474 ArrayOfMatrix& dsurface_emission_dx,
2475 const Index& stokes_dim,
2476 const Index& atmosphere_dim,
2477 const Vector& lat_grid,
2478 const Vector& lon_grid,
2479 const Vector& f_grid,
2480 const Vector& rtp_pos,
2481 const Vector& rtp_los,
2482 const Tensor3& surface_props_data,
2483 const ArrayOfString& surface_props_names,
2484 const ArrayOfString& dsurface_names,
2485 const Index& jacobian_do,
2486 const Verbosity& verbosity) {
2487 // Check surface_data
2488 surface_props_check(atmosphere_dim,
2489 lat_grid,
2490 lon_grid,
2491 surface_props_data,
2492 surface_props_names);
2493
2494 // Interplation grid positions and weights
2495 ArrayOfGridPos gp_lat(1), gp_lon(1);
2496 Matrix itw;
2498 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
2499 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
2500
2501 // Skin temperature
2502 Vector skin_t(1);
2503 surface_props_interp(skin_t,
2504 "Skin temperature",
2505 atmosphere_dim,
2506 gp_lat,
2507 gp_lon,
2508 itw,
2509 surface_props_data,
2510 surface_props_names);
2511
2512 surfaceBlackbody(surface_los,
2513 surface_rmatrix,
2514 surface_emission,
2515 atmosphere_dim,
2516 f_grid,
2517 stokes_dim,
2518 rtp_pos,
2519 rtp_los,
2520 skin_t[0],
2521 verbosity);
2522
2523 surface_rmatrix.resize(1, f_grid.nelem(), stokes_dim, stokes_dim);
2524 surface_rmatrix = 0.0;
2525
2526 // Jacobian part
2527 if (jacobian_do) {
2528 dsurface_check(surface_props_names,
2529 dsurface_names,
2530 dsurface_rmatrix_dx,
2531 dsurface_emission_dx);
2532
2533 Index irq;
2534
2535 // Skin temperature
2536 irq = find_first(dsurface_names, String("Skin temperature"));
2537 if (irq >= 0) {
2538 Matrix dbdt(f_grid.nelem(), 1);
2539 dplanck_dt(dbdt(joker, 0), f_grid, skin_t[0]);
2540 dsurface_emission_dx[irq] = dbdt;
2541 dsurface_rmatrix_dx[irq].resize(surface_rmatrix.nbooks(),
2542 surface_rmatrix.npages(),
2543 surface_rmatrix.nrows(),
2544 surface_rmatrix.ncols());
2545 dsurface_rmatrix_dx[irq] = 0;
2546 }
2547 }
2548}
2549
2550/* Workspace method: Doxygen documentation will be auto-generated */
2551void SurfaceDummy(ArrayOfTensor4& dsurface_rmatrix_dx,
2552 ArrayOfMatrix& dsurface_emission_dx,
2553 const Index& atmosphere_dim,
2554 const Vector& lat_grid,
2555 const Vector& lon_grid,
2556 const Tensor3& surface_props_data,
2557 const ArrayOfString& surface_props_names,
2558 const ArrayOfString& dsurface_names,
2559 const Index& jacobian_do,
2560 const Verbosity&) {
2561 if (surface_props_names.nelem())
2562 throw runtime_error(
2563 "When calling this method, *surface_props_names* should be empty.");
2564
2565 surface_props_check(atmosphere_dim,
2566 lat_grid,
2567 lon_grid,
2568 surface_props_data,
2569 surface_props_names);
2570
2571 if (jacobian_do) {
2572 dsurface_check(surface_props_names,
2573 dsurface_names,
2574 dsurface_rmatrix_dx,
2575 dsurface_emission_dx);
2576 }
2577}
2578
2579/* Workspace method: Doxygen documentation will be auto-generated */
2580void SurfaceFastem(Matrix& surface_los,
2581 Tensor4& surface_rmatrix,
2582 ArrayOfTensor4& dsurface_rmatrix_dx,
2583 Matrix& surface_emission,
2584 ArrayOfMatrix& dsurface_emission_dx,
2585 const Index& stokes_dim,
2586 const Index& atmosphere_dim,
2587 const Vector& lat_grid,
2588 const Vector& lon_grid,
2589 const Vector& f_grid,
2590 const Vector& rtp_pos,
2591 const Vector& rtp_los,
2592 const Tensor3& surface_props_data,
2593 const ArrayOfString& surface_props_names,
2594 const ArrayOfString& dsurface_names,
2595 const Index& jacobian_do,
2596 const Vector& transmittance,
2597 const Index& fastem_version,
2598 const Verbosity& verbosity) {
2599 // Check surface_data
2600 surface_props_check(atmosphere_dim,
2601 lat_grid,
2602 lon_grid,
2603 surface_props_data,
2604 surface_props_names);
2605
2606 // Interplation grid positions and weights
2607 ArrayOfGridPos gp_lat(1), gp_lon(1);
2608 Matrix itw;
2610 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
2611 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
2612
2613 // Skin temperature
2614 Vector skin_t(1);
2615 surface_props_interp(skin_t,
2616 "Water skin temperature",
2617 atmosphere_dim,
2618 gp_lat,
2619 gp_lon,
2620 itw,
2621 surface_props_data,
2622 surface_props_names);
2623
2624 // Wind speed
2625 Vector wind_speed(1);
2626 surface_props_interp(wind_speed,
2627 "Wind speed",
2628 atmosphere_dim,
2629 gp_lat,
2630 gp_lon,
2631 itw,
2632 surface_props_data,
2633 surface_props_names);
2634
2635 // Wind direction
2636 Vector wind_direction(1);
2637 surface_props_interp(wind_direction,
2638 "Wind direction",
2639 atmosphere_dim,
2640 gp_lat,
2641 gp_lon,
2642 itw,
2643 surface_props_data,
2644 surface_props_names);
2645
2646 // Salinity
2647 Vector salinity(1);
2648 surface_props_interp(salinity,
2649 "Salinity",
2650 atmosphere_dim,
2651 gp_lat,
2652 gp_lon,
2653 itw,
2654 surface_props_data,
2655 surface_props_names);
2656
2657 // Call FASTEM
2658 surfaceFastem(surface_los,
2659 surface_rmatrix,
2660 surface_emission,
2661 atmosphere_dim,
2662 stokes_dim,
2663 f_grid,
2664 rtp_pos,
2665 rtp_los,
2666 skin_t[0],
2667 salinity[0],
2668 wind_speed[0],
2669 wind_direction[0],
2670 transmittance,
2671 fastem_version,
2672 verbosity);
2673
2674 // Jacobian part
2675 if (jacobian_do) {
2676 dsurface_check(surface_props_names,
2677 dsurface_names,
2678 dsurface_rmatrix_dx,
2679 dsurface_emission_dx);
2680
2681 Index irq;
2682
2683 // Skin temperature
2684 irq = find_first(dsurface_names, String("Water skin temperature"));
2685 if (irq >= 0) {
2686 const Numeric dd = 0.1;
2687 Matrix surface_los2;
2688 surfaceFastem(surface_los2,
2689 dsurface_rmatrix_dx[irq],
2690 dsurface_emission_dx[irq],
2691 atmosphere_dim,
2692 stokes_dim,
2693 f_grid,
2694 rtp_pos,
2695 rtp_los,
2696 skin_t[0] + dd,
2697 salinity[0],
2698 wind_speed[0],
2699 wind_direction[0],
2700 transmittance,
2701 fastem_version,
2702 verbosity);
2703 //
2704 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
2705 dsurface_rmatrix_dx[irq] /= dd;
2706 dsurface_emission_dx[irq] -= surface_emission;
2707 dsurface_emission_dx[irq] /= dd;
2708 }
2709
2710 // Wind speed
2711 irq = find_first(dsurface_names, String("Wind speed"));
2712 if (irq >= 0) {
2713 const Numeric dd = 0.1;
2714 Matrix surface_los2;
2715 surfaceFastem(surface_los2,
2716 dsurface_rmatrix_dx[irq],
2717 dsurface_emission_dx[irq],
2718 atmosphere_dim,
2719 stokes_dim,
2720 f_grid,
2721 rtp_pos,
2722 rtp_los,
2723 skin_t[0],
2724 salinity[0],
2725 wind_speed[0] + dd,
2726 wind_direction[0],
2727 transmittance,
2728 fastem_version,
2729 verbosity);
2730 //
2731 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
2732 dsurface_rmatrix_dx[irq] /= dd;
2733 dsurface_emission_dx[irq] -= surface_emission;
2734 dsurface_emission_dx[irq] /= dd;
2735 }
2736
2737 // Wind direction
2738 irq = find_first(dsurface_names, String("Wind direction"));
2739 if (irq >= 0) {
2740 const Numeric dd = 1;
2741 Matrix surface_los2;
2742 surfaceFastem(surface_los2,
2743 dsurface_rmatrix_dx[irq],
2744 dsurface_emission_dx[irq],
2745 atmosphere_dim,
2746 stokes_dim,
2747 f_grid,
2748 rtp_pos,
2749 rtp_los,
2750 skin_t[0],
2751 salinity[0],
2752 wind_speed[0],
2753 wind_direction[0] + dd,
2754 transmittance,
2755 fastem_version,
2756 verbosity);
2757 //
2758 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
2759 dsurface_rmatrix_dx[irq] /= dd;
2760 dsurface_emission_dx[irq] -= surface_emission;
2761 dsurface_emission_dx[irq] /= dd;
2762 }
2763
2764 // Salinity
2765 irq = find_first(dsurface_names, String("Salinity"));
2766 if (irq >= 0) {
2767 const Numeric dd = 0.0005;
2768 Matrix surface_los2;
2769 surfaceFastem(surface_los2,
2770 dsurface_rmatrix_dx[irq],
2771 dsurface_emission_dx[irq],
2772 atmosphere_dim,
2773 stokes_dim,
2774 f_grid,
2775 rtp_pos,
2776 rtp_los,
2777 skin_t[0],
2778 salinity[0] + dd,
2779 wind_speed[0],
2780 wind_direction[0],
2781 transmittance,
2782 fastem_version,
2783 verbosity);
2784 //
2785 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
2786 dsurface_rmatrix_dx[irq] /= dd;
2787 dsurface_emission_dx[irq] -= surface_emission;
2788 dsurface_emission_dx[irq] /= dd;
2789 }
2790 }
2791}
2792
2793
2794/* Workspace method: Doxygen documentation will be auto-generated */
2796 Tensor4& surface_rmatrix,
2797 ArrayOfTensor4& dsurface_rmatrix_dx,
2798 Matrix& surface_emission,
2799 ArrayOfMatrix& dsurface_emission_dx,
2800 const Index& stokes_dim,
2801 const Index& atmosphere_dim,
2802 const Vector& lat_grid,
2803 const Vector& lon_grid,
2804 const Vector& f_grid,
2805 const Vector& rtp_pos,
2806 const Vector& rtp_los,
2807 const Vector& specular_los,
2808 const Tensor3& surface_props_data,
2809 const ArrayOfString& surface_props_names,
2810 const ArrayOfString& dsurface_names,
2811 const Index& jacobian_do,
2812 const Vector& f_reflectivities,
2813 const Verbosity& verbosity) {
2814 // Check surface_data
2815 surface_props_check(atmosphere_dim,
2816 lat_grid,
2817 lon_grid,
2818 surface_props_data,
2819 surface_props_names);
2820
2821 // Check of GINs
2822 chk_if_increasing("GIN *f_reflectivities*", f_reflectivities);
2823
2824 // Sizes
2825 const Index nr = f_reflectivities.nelem();
2826 const Index nf = f_grid.nelem();
2827
2828 // Interplation grid positions and weights
2829 ArrayOfGridPos gp_lat(1), gp_lon(1);
2830 Matrix itw;
2832 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
2833 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
2834
2835 // Skin temperature
2836 Vector surface_skin_t(1);
2837 surface_props_interp(surface_skin_t,
2838 "Skin temperature",
2839 atmosphere_dim,
2840 gp_lat,
2841 gp_lon,
2842 itw,
2843 surface_props_data,
2844 surface_props_names);
2845
2846 // Reflectivities
2847 Vector reflectivities(nr);
2848 //
2849 for (Index i=0; i<nr; i++) {
2850 Vector rv(1);
2851 ostringstream sstr;
2852 sstr << "Scalar reflectivity " << i;
2854 sstr.str(),
2855 atmosphere_dim,
2856 gp_lat,
2857 gp_lon,
2858 itw,
2859 surface_props_data,
2860 surface_props_names);
2861 reflectivities[i] = rv[0];
2862 }
2863
2864 // Create surface_scalar_reflectivity
2865 Vector surface_scalar_reflectivity(nf);
2866 ArrayOfGridPos gp(nf);
2867 if (nr==1) {
2868 gp4length1grid(gp);
2869 } else {
2870 gridpos(gp, f_reflectivities, f_grid, 1e99);
2872 }
2873 Matrix itw2(nf, 2);
2874 interpweights(itw2, gp);
2875 interp( surface_scalar_reflectivity, itw2, reflectivities, gp);
2876
2877 // Call non-Jacobian method
2879 surface_rmatrix,
2880 surface_emission,
2881 f_grid,
2882 stokes_dim,
2883 atmosphere_dim,
2884 rtp_pos,
2885 rtp_los,
2886 specular_los,
2887 surface_skin_t[0],
2888 surface_scalar_reflectivity,
2889 verbosity);
2890
2891 // Jacobian part
2892 if (jacobian_do) {
2893 dsurface_check(surface_props_names,
2894 dsurface_names,
2895 dsurface_rmatrix_dx,
2896 dsurface_emission_dx);
2897
2898 Index irq;
2899
2900 // Skin temperature
2901 irq = find_first(dsurface_names, String("Skin temperature"));
2902 if (irq >= 0) {
2903 const Numeric dd = 0.1;
2904 Matrix surface_los2;
2905 surfaceFlatScalarReflectivity(surface_los2,
2906 dsurface_rmatrix_dx[irq],
2907 dsurface_emission_dx[irq],
2908 f_grid,
2909 stokes_dim,
2910 atmosphere_dim,
2911 rtp_pos,
2912 rtp_los,
2913 specular_los,
2914 surface_skin_t[0]+dd,
2915 surface_scalar_reflectivity,
2916 verbosity);
2917 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
2918 dsurface_rmatrix_dx[irq] /= dd;
2919 dsurface_emission_dx[irq] -= surface_emission;
2920 dsurface_emission_dx[irq] /= dd;
2921 }
2922
2923 // Reflectivities
2924 for (Index i=0; i<nr; i++) {
2925 ostringstream sstr;
2926 sstr << "Scalar reflectivity " << i;
2927 irq = find_first(dsurface_names, String(sstr.str()));
2928 if (irq >= 0) {
2929 const Numeric dd = 1e-4;
2930 Vector rpert = reflectivities;
2931 rpert[i] += dd;
2932 interp( surface_scalar_reflectivity, itw2, rpert, gp);
2933 //
2934 Matrix surface_los2;
2935 surfaceFlatScalarReflectivity(surface_los2,
2936 dsurface_rmatrix_dx[irq],
2937 dsurface_emission_dx[irq],
2938 f_grid,
2939 stokes_dim,
2940 atmosphere_dim,
2941 rtp_pos,
2942 rtp_los,
2943 specular_los,
2944 surface_skin_t[0],
2945 surface_scalar_reflectivity,
2946 verbosity);
2947 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
2948 dsurface_rmatrix_dx[irq] /= dd;
2949 dsurface_emission_dx[irq] -= surface_emission;
2950 dsurface_emission_dx[irq] /= dd;
2951 }
2952 }
2953 }
2954}
2955
2956/* Workspace method: Doxygen documentation will be auto-generated */
2957void SurfaceTessem(Matrix& surface_los,
2958 Tensor4& surface_rmatrix,
2959 ArrayOfTensor4& dsurface_rmatrix_dx,
2960 Matrix& surface_emission,
2961 ArrayOfMatrix& dsurface_emission_dx,
2962 const Index& stokes_dim,
2963 const Index& atmosphere_dim,
2964 const Vector& lat_grid,
2965 const Vector& lon_grid,
2966 const Vector& f_grid,
2967 const Vector& rtp_pos,
2968 const Vector& rtp_los,
2969 const TessemNN& net_h,
2970 const TessemNN& net_v,
2971 const Tensor3& surface_props_data,
2972 const ArrayOfString& surface_props_names,
2973 const ArrayOfString& dsurface_names,
2974 const Index& jacobian_do,
2975 const Verbosity& verbosity) {
2976 // Check surface_data
2977 surface_props_check(atmosphere_dim,
2978 lat_grid,
2979 lon_grid,
2980 surface_props_data,
2981 surface_props_names);
2982
2983 // Interplation grid positions and weights
2984 ArrayOfGridPos gp_lat(1), gp_lon(1);
2985 Matrix itw;
2987 gp_lat[0], gp_lon[0], atmosphere_dim, lat_grid, lon_grid, rtp_pos);
2988 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
2989
2990 // Skin temperature
2991 Vector skin_t(1);
2992 surface_props_interp(skin_t,
2993 "Water skin temperature",
2994 atmosphere_dim,
2995 gp_lat,
2996 gp_lon,
2997 itw,
2998 surface_props_data,
2999 surface_props_names);
3000
3001 // Wind speed
3002 Vector wind_speed(1);
3003 surface_props_interp(wind_speed,
3004 "Wind speed",
3005 atmosphere_dim,
3006 gp_lat,
3007 gp_lon,
3008 itw,
3009 surface_props_data,
3010 surface_props_names);
3011
3012 // Salinity
3013 Vector salinity(1);
3014 surface_props_interp(salinity,
3015 "Salinity",
3016 atmosphere_dim,
3017 gp_lat,
3018 gp_lon,
3019 itw,
3020 surface_props_data,
3021 surface_props_names);
3022
3023 // Call TESSEM
3024 surfaceTessem(surface_los,
3025 surface_rmatrix,
3026 surface_emission,
3027 atmosphere_dim,
3028 stokes_dim,
3029 f_grid,
3030 rtp_pos,
3031 rtp_los,
3032 skin_t[0],
3033 net_h,
3034 net_v,
3035 salinity[0],
3036 wind_speed[0],
3037 verbosity);
3038
3039 // Jacobian part
3040 if (jacobian_do) {
3041 dsurface_check(surface_props_names,
3042 dsurface_names,
3043 dsurface_rmatrix_dx,
3044 dsurface_emission_dx);
3045
3046 Index irq;
3047
3048 // Skin temperature
3049 irq = find_first(dsurface_names, String("Water skin temperature"));
3050 if (irq >= 0) {
3051 const Numeric dd = 0.1;
3052 Matrix surface_los2;
3053 surfaceTessem(surface_los2,
3054 dsurface_rmatrix_dx[irq],
3055 dsurface_emission_dx[irq],
3056 atmosphere_dim,
3057 stokes_dim,
3058 f_grid,
3059 rtp_pos,
3060 rtp_los,
3061 skin_t[0] + dd,
3062 net_h,
3063 net_v,
3064 salinity[0],
3065 wind_speed[0],
3066 verbosity);
3067 //
3068 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3069 dsurface_rmatrix_dx[irq] /= dd;
3070 dsurface_emission_dx[irq] -= surface_emission;
3071 dsurface_emission_dx[irq] /= dd;
3072 }
3073
3074 // Wind speed
3075 irq = find_first(dsurface_names, String("Wind speed"));
3076 if (irq >= 0) {
3077 const Numeric dd = 0.1;
3078 Matrix surface_los2;
3079 surfaceTessem(surface_los2,
3080 dsurface_rmatrix_dx[irq],
3081 dsurface_emission_dx[irq],
3082 atmosphere_dim,
3083 stokes_dim,
3084 f_grid,
3085 rtp_pos,
3086 rtp_los,
3087 skin_t[0],
3088 net_h,
3089 net_v,
3090 salinity[0],
3091 wind_speed[0] + dd,
3092 verbosity);
3093 //
3094 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3095 dsurface_rmatrix_dx[irq] /= dd;
3096 dsurface_emission_dx[irq] -= surface_emission;
3097 dsurface_emission_dx[irq] /= dd;
3098 }
3099
3100 // Salinity
3101 irq = find_first(dsurface_names, String("Salinity"));
3102 if (irq >= 0) {
3103 const Numeric dd = 0.0005;
3104 Matrix surface_los2;
3105 surfaceTessem(surface_los2,
3106 dsurface_rmatrix_dx[irq],
3107 dsurface_emission_dx[irq],
3108 atmosphere_dim,
3109 stokes_dim,
3110 f_grid,
3111 rtp_pos,
3112 rtp_los,
3113 skin_t[0],
3114 net_h,
3115 net_v,
3116 salinity[0] + dd,
3117 wind_speed[0],
3118 verbosity);
3119 //
3120 dsurface_rmatrix_dx[irq] -= surface_rmatrix;
3121 dsurface_rmatrix_dx[irq] /= dd;
3122 dsurface_emission_dx[irq] -= surface_emission;
3123 dsurface_emission_dx[irq] /= dd;
3124 }
3125 }
3126}
3127
3128/* Workspace method: Doxygen documentation will be auto-generated */
3130 const ArrayOfString& iy_aux_vars,
3131 const ArrayOfMatrix& iy_aux,
3132 const Verbosity&) {
3133 Index ihit = -1;
3134
3135 for (Index i = 0; i < iy_aux_vars.nelem(); i++) {
3136 if (iy_aux_vars[i] == "Optical depth") {
3137 ihit = i;
3138 break;
3139 }
3140 }
3141
3142 if (ihit < 0)
3143 throw runtime_error("No element in *iy_aux* holds optical depths.");
3144
3145 const Index n = iy_aux[ihit].nrows();
3146
3147 transmittance.resize(n);
3148
3149 for (Index i = 0; i < n; i++) {
3150 transmittance[i] = exp(-iy_aux[ihit](i, 0));
3151 }
3152}
This file contains the definition of Array.
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:294
The global header file for ARTS.
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:25199
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:25266
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, 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:24276
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:24537
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:25148
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_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:51
String name() const
Agenda name.
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
Index nrows() const noexcept
Definition: matpackI.h:1061
Index ncols() const noexcept
Definition: matpackI.h:1062
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:541
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:1270
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1063
The range class.
Definition: matpackI.h:166
A telsem atlas.
Definition: telsem.h:59
Index calc_cellnum(Numeric lat, Numeric lon) const
Definition: telsem.cc:142
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:170
Vector get_emis_h(Index cellnum) const
Definition: telsem.h:159
std::pair< Numeric, Numeric > get_coordinates(Index cellnum) const
Definition: telsem.cc:225
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:287
Vector get_emis_v(Index i) const
Definition: telsem.h:137
The Tensor3 class.
Definition: matpackIII.h:344
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:427
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The VectorView class.
Definition: matpackI.h:658
The Vector class.
Definition: matpackI.h:908
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
Workspace class.
Definition: workspace_ng.h:40
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:107
This file contains functions that are adapted from FASTEM code which is used to calculate surface emi...
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1335
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
Definition: geodetic.cc:1408
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1304
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:474
void jacobian_type_extrapol(ArrayOfGridPos &gp)
Adopts grid positions to extrapolation used for jacobians.
Definition: jacobian.cc:835
#define abs(x)
#define q
#define min(a, b)
#define max(a, b)
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:1085
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:1820
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:2159
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:1929
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:200
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:1192
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:2424
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:1673
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:767
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:915
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:1024
Numeric EARTH_RADIUS
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:1271
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:2470
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:730
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:876
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 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:571
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:2580
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:1505
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:2284
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:2025
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:1341
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:1425
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:247
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 &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:443
void transmittanceFromIy_aux(Vector &transmittance, const ArrayOfString &iy_aux_vars, const ArrayOfMatrix &iy_aux, const Verbosity &)
WORKSPACE METHOD: transmittanceFromIy_aux.
Definition: m_surface.cc:3129
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:1571
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:326
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:2266
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:2957
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:2551
const Numeric DEG2RAD
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:141
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:66
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:2795
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:440
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1401
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:207
#define CREATE_OUT2
Definition: messages.h:206
Array< String > ArrayOfString
An array of Strings.
Definition: mystring.h:290
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:287
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:312
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:299
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:1067
void resolve_lon(Numeric &lon, const Numeric &lon5, const Numeric &lon6)
Resolves which longitude angle that shall be used.
Definition: ppath.cc:500
Numeric plevel_angletilt(const Numeric &r, const Numeric &c1)
Calculates the angular tilt of the surface or a pressure level.
Definition: ppath.cc:616
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:579
Propagation path structure and functions.
#define d
#define w
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:747
void mueller_stokes2modif(Matrix &Cm, const Index &stokes_dim)
mueller_stokes2modif
Definition: rte.cc:1941
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:1962
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:1803
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:83
void mueller_modif2stokes(Matrix &Cs, const Index &stokes_dim)
mueller_modif2stokes
Definition: rte.cc:1893
void mueller_rotation(Matrix &L, const Index &stokes_dim, const Numeric &rotangle)
mueller_rotation
Definition: rte.cc:1914
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.
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_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...