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