ARTS  2.4.0(git:4fb77825)
rt4.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016 Jana Mendrok <jana.mendrok@gmail.com>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA.
17 */
18 
28 /*===========================================================================
29  === External declarations
30  ===========================================================================*/
31 
32 #include "config.h"
33 
34 #ifdef ENABLE_RT4
35 
36 #include <complex.h>
37 #include <cfloat>
38 #include <stdexcept>
39 #include "disort.h"
40 #include "interpolation.h"
41 #include "m_xml.h"
42 #include "optproperties.h"
43 #include "physics_funcs.h"
44 #include "rt4.h"
45 #include "rte.h"
46 
47 using std::ostringstream;
48 using std::runtime_error;
49 
50 extern const Numeric PI;
51 extern const Numeric DEG2RAD;
52 extern const Numeric RAD2DEG;
53 extern const Numeric SPEED_OF_LIGHT;
54 extern const Numeric COSMIC_BG_TEMP;
55 
56 void check_rt4_input( // Output
57  Index& nhstreams,
58  Index& nhza,
59  Index& nummu,
60  // Input
61  const Index& cloudbox_on,
62  const Index& atmfields_checked,
63  const Index& atmgeom_checked,
64  const Index& cloudbox_checked,
65  const Index& scat_data_checked,
68  const Index& atmosphere_dim,
69  const Index& stokes_dim,
70  const Index& nstreams,
71  const String& quad_type,
72  const Index& add_straight_angles,
73  const Index& pnd_ncols) {
74  // Don't do anything if there's no cloudbox defined.
75  if (!cloudbox_on) {
76  throw runtime_error(
77  "Cloudbox is off, no scattering calculations to be"
78  " performed.");
79  }
80 
81  if (atmfields_checked != 1)
82  throw runtime_error(
83  "The atmospheric fields must be flagged to have"
84  " passed a consistency check (atmfields_checked=1).");
85  if (atmgeom_checked != 1)
86  throw runtime_error(
87  "The atmospheric geometry must be flagged to have"
88  " passed a consistency check (atmgeom_checked=1).");
89  if (cloudbox_checked != 1)
90  throw runtime_error(
91  "The cloudbox must be flagged to have"
92  " passed a consistency check (cloudbox_checked=1).");
93  if (scat_data_checked != 1)
94  throw runtime_error(
95  "The scat_data must be flagged to have "
96  "passed a consistency check (scat_data_checked=1).");
97 
98  if (atmosphere_dim != 1)
99  throw runtime_error(
100  "For running RT4, atmospheric dimensionality"
101  " must be 1.\n");
102 
103  if (stokes_dim < 0 || stokes_dim > 2)
104  throw runtime_error(
105  "For running RT4, the dimension of stokes vector"
106  " must be 1 or 2.\n");
107 
108  if (cloudbox_limits[0] != 0) {
109  ostringstream os;
110  os << "RT4 calculations currently only possible with"
111  << " lower cloudbox limit\n"
112  << "at 0th atmospheric level"
113  << " (assumes surface there, ignoring z_surface).\n";
114  throw runtime_error(os.str());
115  }
116 
117  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
118  throw runtime_error(
119  "*cloudbox_limits* is a vector which contains the"
120  " upper and lower limit of the cloud for all"
121  " atmospheric dimensions. So its dimension must"
122  " be 2 x *atmosphere_dim*");
123 
124  if (scat_data.empty())
125  throw runtime_error(
126  "No single scattering data present.\n"
127  "See documentation of WSV *scat_data* for options to"
128  " make single scattering data available.\n");
129 
130  if (pnd_ncols != 1)
131  throw runtime_error(
132  "*pnd_field* is not 1D! \n"
133  "RT4 can only be used for 1D!\n");
134 
135  if (quad_type.length() > 1) {
136  ostringstream os;
137  os << "Input parameter *quad_type* not allowed to contain more than a"
138  << " single character.\n"
139  << "Yours has " << quad_type.length() << ".\n";
140  throw runtime_error(os.str());
141  }
142 
143  if (quad_type == "D" || quad_type == "G") {
144  if (add_straight_angles)
145  nhza = 1;
146  else
147  nhza = 0;
148  } else if (quad_type == "L") {
149  nhza = 0;
150  } else {
151  ostringstream os;
152  os << "Unknown quadrature type: " << quad_type
153  << ".\nOnly D, G, and L allowed.\n";
154  throw runtime_error(os.str());
155  }
156 
157  // RT4 actually uses number of angles in single hemisphere. However, we don't
158  // want a bunch of different approaches used in the interface, so we apply the
159  // DISORT (and ARTS) way of total number of angles here. Hence, we have to
160  // ensure here that total number is even.
161  if (nstreams / 2 * 2 != nstreams) {
162  ostringstream os;
163  os << "RT4 requires an even number of streams, but yours is " << nstreams
164  << ".\n";
165  throw runtime_error(os.str());
166  }
167  nhstreams = nstreams / 2;
168  // nummu is the total number of angles in one hemisphere, including both
169  // the quadrature angles and the "extra" angles.
170  nummu = nhstreams + nhza;
171 
172  // RT4 can only completely or azimuthally randomly oriented particles.
173  bool no_arb_ori = true;
174  for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
175  for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++)
176  if (scat_data[i_ss][i_se].ptype != PTYPE_TOTAL_RND &&
177  scat_data[i_ss][i_se].ptype != PTYPE_AZIMUTH_RND)
178  no_arb_ori = false;
179  if (!no_arb_ori) {
180  ostringstream os;
181  os << "RT4 can only handle scattering elements of type " << PTYPE_TOTAL_RND
182  << " (" << PTypeToString(PTYPE_TOTAL_RND) << ") and\n"
184  << "),\n"
185  << "but at least one element of other type (" << PTYPE_GENERAL << "="
186  << PTypeToString(PTYPE_GENERAL) << ") is present.\n";
187  throw runtime_error(os.str());
188  }
189 
190  //------------- end of checks ---------------------------------------
191 }
192 
193 void get_quad_angles( // Output
194  VectorView mu_values,
195  VectorView quad_weights,
196  Vector& za_grid,
197  Vector& aa_grid,
198  // Input
199  const String& quad_type,
200  const Index& nhstreams,
201  const Index& nhza,
202  const Index& nummu) {
203  if (quad_type == "D") {
204  double_gauss_quadrature_(
205  nhstreams, mu_values.get_c_array(), quad_weights.get_c_array());
206  } else if (quad_type == "G") {
207  gauss_legendre_quadrature_(
208  nhstreams, mu_values.get_c_array(), quad_weights.get_c_array());
209  } else //if( quad_type=="L" )
210  {
211  lobatto_quadrature_(
212  nhstreams, mu_values.get_c_array(), quad_weights.get_c_array());
213  }
214 
215  // Set "extra" angle (at 0deg) if quad_type!="L" && !add_straight_angles
216  if (nhza > 0) mu_values[nhstreams] = 1.;
217 
218  // FIXME: we should be able to avoid setting za_grid here in one way,
219  // and resetting in another before leaving the WSM. This, however, requires
220  // rearranging the angle order and angle assignment in the RT4-SSP prep
221  // routines.
222  za_grid.resize(2 * nummu);
223  for (Index imu = 0; imu < nummu; imu++) {
224  za_grid[imu] = acos(mu_values[imu]) * RAD2DEG;
225  za_grid[nummu + imu] = 180. - za_grid[imu];
226  }
227  aa_grid.resize(1);
228  aa_grid[0] = 0.;
229 }
230 
231 void get_rt4surf_props( // Output
232  Vector& ground_albedo,
233  Tensor3& ground_reflec,
234  ComplexVector& ground_index,
235  // Input
237  const String& ground_type,
238  const Numeric& surface_skin_t,
242  const Index& stokes_dim) {
243  if (surface_skin_t < 0. || surface_skin_t > 1000.) {
244  ostringstream os;
245  os << "Surface temperature is set to " << surface_skin_t << " K,\n"
246  << "which is not considered a meaningful value.\n";
247  throw runtime_error(os.str());
248  }
249 
250  const Index nf = f_grid.nelem();
251 
252  if (ground_type == "L") // RT4's proprietary Lambertian
253  {
254  if (min(surface_scalar_reflectivity) < 0 ||
256  throw runtime_error(
257  "All values in *surface_scalar_reflectivity* must be inside [0,1].");
258  }
259 
260  // surface albedo
261  if (surface_scalar_reflectivity.nelem() == f_grid.nelem())
262  ground_albedo = surface_scalar_reflectivity;
263  else if (surface_scalar_reflectivity.nelem() == 1)
264  ground_albedo += surface_scalar_reflectivity[0];
265  else {
266  ostringstream os;
267  os << "For Lambertian surface reflection, the number of elements in\n"
268  << "*surface_scalar_reflectivity* needs to match the length of\n"
269  << "*f_grid* or be 1."
270  << "\n length of *f_grid* : " << f_grid.nelem()
271  << "\n length of *surface_scalar_reflectivity* : "
272  << surface_scalar_reflectivity.nelem() << "\n";
273  throw runtime_error(os.str());
274  }
275 
276  } else if (ground_type == "S") // RT4's 'proprietary' Specular
277  {
278  const Index ref_sto = surface_reflectivity.nrows();
279 
280  chk_if_in_range("surface_reflectivity's stokes_dim", ref_sto, 1, 4);
281  if (ref_sto != surface_reflectivity.ncols()) {
282  ostringstream os;
283  os << "The number of rows and columnss in *surface_reflectivity*\n"
284  << "must match each other.";
285  throw runtime_error(os.str());
286  }
287 
288  if (min(surface_reflectivity(joker, 0, 0)) < 0 ||
289  max(surface_reflectivity(joker, 0, 0)) > 1) {
290  throw runtime_error(
291  "All r11 values in *surface_reflectivity* must be inside [0,1].");
292  }
293 
294  // surface reflectivity
295  if (surface_reflectivity.npages() == f_grid.nelem())
296  if (ref_sto < stokes_dim)
297  ground_reflec(joker, Range(0, ref_sto), Range(0, ref_sto)) =
299  else
300  ground_reflec = surface_reflectivity(
302  else if (surface_reflectivity.npages() == 1)
303  if (ref_sto < stokes_dim)
304  for (Index f_index = 0; f_index < nf; f_index++)
305  ground_reflec(f_index, Range(0, ref_sto), Range(0, ref_sto)) +=
307  else
308  for (Index f_index = 0; f_index < nf; f_index++)
309  ground_reflec(f_index, joker, joker) += surface_reflectivity(
310  0, Range(0, stokes_dim), Range(0, stokes_dim));
311  else {
312  ostringstream os;
313  os << "For specular surface reflection, the number of elements in\n"
314  << "*surface_reflectivity* needs to match the length of\n"
315  << "*f_grid* or be 1."
316  << "\n length of *f_grid* : " << f_grid.nelem()
317  << "\n length of *surface_reflectivity* : "
318  << surface_reflectivity.npages() << "\n";
319  throw runtime_error(os.str());
320  }
321 
322  } else if (ground_type == "F") // RT4's proprietary Fresnel
323  {
324  //extract/interpolate from GriddedField
325  Matrix n_real(nf, 1), n_imag(nf, 1);
326  complex_n_interp(n_real,
327  n_imag,
329  "surface_complex_refr_index",
330  f_grid,
331  Vector(1, surface_skin_t));
332  //ground_index = Complex(n_real(joker,0),n_imag(joker,0));
333  for (Index f_index = 0; f_index < nf; f_index++) {
334  ground_index[f_index] = Complex(n_real(f_index, 0), n_imag(f_index, 0));
335  }
336  } else {
337  ostringstream os;
338  os << "Unknown surface type.\n";
339  throw runtime_error(os.str());
340  }
341 }
342 
343 void run_rt4(Workspace& ws,
344  // Output
346  Vector& za_grid,
347  // Input
350  ConstVectorView z_profile,
351  ConstVectorView t_profile,
352  ConstMatrixView vmr_profiles,
353  ConstMatrixView pnd_profiles,
357  const Index& stokes_dim,
358  const Index& nummu,
359  const Index& nhza,
360  const String& ground_type,
361  const Numeric& surface_skin_t,
362  ConstVectorView ground_albedo,
363  ConstTensor3View ground_reflec,
364  ConstComplexVectorView ground_index,
365  ConstTensor5View surf_refl_mat,
366  ConstTensor3View surf_emis_vec,
368  const Numeric& surf_altitude,
369  const String& quad_type,
370  Vector& mu_values,
371  ConstVectorView quad_weights,
372  const Index& auto_inc_nstreams,
373  const Index& robust,
374  const Index& za_interp_order,
375  const Index& cos_za_interp,
376  const String& pfct_method,
377  const Index& pfct_aa_grid_size,
378  const Numeric& pfct_threshold,
379  const Numeric& max_delta_tau,
380  const Verbosity& verbosity) {
381  // Create an atmosphere starting at z_surface
382  Vector p, z, t;
383  Matrix vmr, pnd;
384  ArrayOfIndex cboxlims;
385  Index ncboxremoved;
386  //
387  reduced_1datm(p,
388  z,
389  t,
390  vmr,
391  pnd,
392  cboxlims,
393  ncboxremoved,
394  p_grid,
395  z_profile,
396  surf_altitude,
397  t_profile,
398  vmr_profiles,
399  pnd_profiles,
401 
402  // Input variables for RT4
403  Index num_layers = p.nelem() - 1;
404 
405  // Top of the atmosphere temperature
406  // FIXME: so far hard-coded to cosmic background. However, change that to set
407  // according to/from space_agenda.
408  // to do so, we need to hand over sky_radiance instead of sky_temp as
409  // Tensor3(2,stokes_dim,nummu) per frequency. That is, for properly using
410  // iy_space_agenda, we need to recall the agenda over the stream angles (not
411  // sure what to do with the upwelling ones. according to the RT4-internal
412  // sizing, sky_radiance contains even those. but they might not be used
413  // (check!) and hence be set arbitrary.
414  const Numeric sky_temp = COSMIC_BG_TEMP;
415 
416  // Data fields
417  Vector height(num_layers + 1);
418  Vector temperatures(num_layers + 1);
419  for (Index i = 0; i < height.nelem(); i++) {
420  height[i] = z[num_layers - i];
421  temperatures[i] = t[num_layers - i];
422  }
423 
424  // this indexes all cloudbox layers as cloudy layers.
425  // optional FIXME: to use the power of RT4 (faster solving scheme for
426  // individual non-cloudy layers), one should consider non-cloudy layers within
427  // cloudbox. That requires some kind of recognition and index setting based on
428  // pnd_field or (derived) cloud layer extinction or scattering.
429  // we use something similar with iyHybrid. have a look there...
430  const Index num_scatlayers = pnd.ncols() - 1;
431  Vector scatlayers(num_layers, 0.);
432  Vector gas_extinct(num_layers, 0.);
433  Tensor6 scatter_matrix(
434  num_scatlayers, 4, nummu, stokes_dim, nummu, stokes_dim, 0.);
435  Tensor6 extinct_matrix(
436  1, num_scatlayers, 2, nummu, stokes_dim, stokes_dim, 0.);
437  Tensor5 emis_vector(1, num_scatlayers, 2, nummu, stokes_dim, 0.);
438 
439  // if there is no scatt particle at all, we don't need to calculate the
440  // scat properties (FIXME: that should rather be done by a proper setting
441  // of scat_layers).
442  Vector pnd_per_level(pnd.ncols());
443  for (Index clev = 0; clev < pnd.ncols(); clev++)
444  pnd_per_level[clev] = pnd(joker, clev).sum();
445  Numeric pndtot = pnd_per_level.sum();
446 
447  for (Index i = 0; i < cboxlims[1] - cboxlims[0]; i++) {
448  scatlayers[num_layers - 1 - cboxlims[0] - i] = float(i + 1);
449  }
450 
451  // Output variables
452  Tensor3 up_rad(num_layers + 1, nummu, stokes_dim, 0.);
453  Tensor3 down_rad(num_layers + 1, nummu, stokes_dim, 0.);
454 
455  Tensor6 extinct_matrix_allf;
456  Tensor5 emis_vector_allf;
457  if (!auto_inc_nstreams) {
458  extinct_matrix_allf.resize(
459  f_grid.nelem(), num_scatlayers, 2, nummu, stokes_dim, stokes_dim);
460  emis_vector_allf.resize(
461  f_grid.nelem(), num_scatlayers, 2, nummu, stokes_dim);
462  par_optpropCalc(emis_vector_allf,
463  extinct_matrix_allf,
464  //scatlayers,
465  scat_data,
466  za_grid,
467  -1,
468  pnd,
469  t,
470  cboxlims,
471  stokes_dim);
472  if (emis_vector_allf.nshelves() == 1) // scat_data had just a single freq
473  // point. copy into emis/ext here and
474  // don't touch anymore later on.
475  {
476  emis_vector = emis_vector_allf(Range(0, 1), joker, joker, joker, joker);
477  extinct_matrix =
478  extinct_matrix_allf(Range(0, 1), joker, joker, joker, joker, joker);
479  }
480  }
481 
482  // FIXME: once, all old optprop scheme incl. the applied agendas is removed,
483  // we can remove this as well.
484  Vector za_grid_orig;
485  if (auto_inc_nstreams) {
486  // For the WSV za_grid, we need to reset these grids instead of
487  // creating a new container. this because further down some agendas are
488  // used that access scat_za/aa_grid through the workspace.
489  // Later on, we need to reconstruct the original setting, hence backup
490  // that here.
491  za_grid_orig = za_grid;
492  }
493 
494  Index nummu_new = 0;
495  // Loop over frequencies
496  for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
497  // Wavelength [um]
498  Numeric wavelength;
499  wavelength = 1e6 * SPEED_OF_LIGHT / f_grid[f_index];
500 
501  Matrix groundreflec = ground_reflec(f_index, joker, joker);
502  Tensor4 surfreflmat = surf_refl_mat(f_index, joker, joker, joker, joker);
503  Matrix surfemisvec = surf_emis_vec(f_index, joker, joker);
504  //Vector muvalues=mu_values;
505 
506  // only update gas_extinct if there is any gas absorption at all (since
507  // vmr_field is not freq-dependent, gas_extinct will remain as above
508  // initialized (with 0) for all freqs, ie we can rely on that it wasn't
509  // changed).
510  if (vmr.ncols() > 0) {
511  gas_optpropCalc(ws,
512  gas_extinct,
514  t[Range(0, num_layers + 1)],
515  vmr(joker, Range(0, num_layers + 1)),
516  p[Range(0, num_layers + 1)],
517  f_grid[Range(f_index, 1)]);
518  }
519 
520  Index pfct_failed = 0;
521  if (pndtot != 0) {
522  if (nummu_new < nummu) {
523  if (!auto_inc_nstreams) // all freq calculated before. just copy
524  // here. but only if needed.
525  {
526  if (emis_vector_allf.nshelves() != 1) {
527  emis_vector =
528  emis_vector_allf(Range(f_index, 1), joker, joker, joker, joker);
529  extinct_matrix = extinct_matrix_allf(
531  }
532  } else {
533  par_optpropCalc(emis_vector,
534  extinct_matrix,
535  //scatlayers,
536  scat_data,
537  za_grid,
538  f_index,
539  pnd,
540  t[Range(0, num_layers + 1)],
541  cboxlims,
542  stokes_dim);
543  }
544  sca_optpropCalc(scatter_matrix,
545  pfct_failed,
546  emis_vector(0, joker, joker, joker, joker),
547  extinct_matrix(0, joker, joker, joker, joker, joker),
548  f_index,
549  scat_data,
550  pnd,
551  stokes_dim,
552  za_grid,
553  quad_weights,
554  pfct_method,
555  pfct_aa_grid_size,
556  pfct_threshold,
557  auto_inc_nstreams,
558  verbosity);
559  } else {
560  pfct_failed = 1;
561  }
562  }
563 
564  if (!pfct_failed) {
565 #pragma omp critical(fortran_rt4)
566  {
567  // Call RT4
568  radtrano_(stokes_dim,
569  nummu,
570  nhza,
571  max_delta_tau,
572  quad_type.c_str(),
574  ground_type.c_str(),
575  ground_albedo[f_index],
576  ground_index[f_index],
577  groundreflec.get_c_array(),
578  surfreflmat.get_c_array(),
579  surfemisvec.get_c_array(),
580  sky_temp,
581  wavelength,
582  num_layers,
583  height.get_c_array(),
584  temperatures.get_c_array(),
585  gas_extinct.get_c_array(),
586  num_scatlayers,
587  scatlayers.get_c_array(),
588  extinct_matrix.get_c_array(),
589  emis_vector.get_c_array(),
590  scatter_matrix.get_c_array(),
591  //noutlevels,
592  //outlevels.get_c_array(),
593  mu_values.get_c_array(),
594  up_rad.get_c_array(),
595  down_rad.get_c_array());
596  }
597 
598  } else { // if (auto_inc_nstreams)
599 
600  if (nummu_new < nummu) nummu_new = nummu + 1;
601 
602  Index nhstreams_new;
603  Vector mu_values_new, quad_weights_new, aa_grid_new;
604  Tensor6 scatter_matrix_new;
605  Tensor6 extinct_matrix_new;
606  Tensor5 emis_vector_new;
607  Tensor4 surfreflmat_new;
608  Matrix surfemisvec_new;
609 
610  while (pfct_failed && (2 * nummu_new) <= auto_inc_nstreams) {
611  // resize and recalc nstream-affected/determined variables:
612  // - mu_values, quad_weights (resize & recalc)
613  nhstreams_new = nummu_new - nhza;
614  mu_values_new.resize(nummu_new);
615  mu_values_new = 0.;
616  quad_weights_new.resize(nummu_new);
617  quad_weights_new = 0.;
618  get_quad_angles(mu_values_new,
619  quad_weights_new,
620  za_grid,
621  aa_grid_new,
622  quad_type,
623  nhstreams_new,
624  nhza,
625  nummu_new);
626 
627  // - resize & recalculate emis_vector, extinct_matrix (as input to scatter_matrix calc)
628  extinct_matrix_new.resize(
629  1, num_scatlayers, 2, nummu_new, stokes_dim, stokes_dim);
630  extinct_matrix_new = 0.;
631  emis_vector_new.resize(1, num_scatlayers, 2, nummu_new, stokes_dim);
632  emis_vector_new = 0.;
633  // FIXME: So far, outside-of-freq-loop calculated optprops will fall
634  // back to in-loop-calculated ones in case of auto-increasing stream
635  // numbers. There might be better options, but I (JM) couldn't come up
636  // with or decide for one so far (we could recalc over all freqs. but
637  // that would unnecessarily recalc lower-freq optprops, too, which are
638  // not needed anymore. which could likely take more time than we
639  // potentially safe through all-at-once temperature and direction
640  // interpolations.
641  par_optpropCalc(emis_vector_new,
642  extinct_matrix_new,
643  //scatlayers,
644  scat_data,
645  za_grid,
646  f_index,
647  pnd,
648  t[Range(0, num_layers + 1)],
649  cboxlims,
650  stokes_dim);
651 
652  // - resize & recalc scatter_matrix
653  scatter_matrix_new.resize(
654  num_scatlayers, 4, nummu_new, stokes_dim, nummu_new, stokes_dim);
655  scatter_matrix_new = 0.;
656  pfct_failed = 0;
657  sca_optpropCalc(
658  scatter_matrix_new,
659  pfct_failed,
660  emis_vector_new(0, joker, joker, joker, joker),
661  extinct_matrix_new(0, joker, joker, joker, joker, joker),
662  f_index,
663  scat_data,
664  pnd,
665  stokes_dim,
666  za_grid,
667  quad_weights_new,
668  pfct_method,
669  pfct_aa_grid_size,
670  pfct_threshold,
671  auto_inc_nstreams,
672  verbosity);
673 
674  if (pfct_failed) nummu_new = nummu_new + 1;
675  }
676 
677  if (pfct_failed) {
678  nummu_new = nummu_new - 1;
679  ostringstream os;
680  os << "Could not increase nstreams sufficiently (current: "
681  << 2 * nummu_new << ")\n"
682  << "to satisfy scattering matrix norm at f[" << f_index
683  << "]=" << f_grid[f_index] * 1e-9 << " GHz.\n";
684  if (!robust) {
685  // couldn't find a nstreams within the limits of auto_inc_nstremas
686  // (aka max. nstreams) that satisfies the scattering matrix norm.
687  // Hence fail completely.
688  os << "Try higher maximum number of allowed streams (ie. higher"
689  << " auto_inc_nstreams than " << auto_inc_nstreams << ").";
690  throw runtime_error(os.str());
691  } else {
692  CREATE_OUT1;
693  os << "Continuing with nstreams=" << 2 * nummu_new
694  << ". Output for this frequency might be erroneous.";
695  out1 << os.str();
696  pfct_failed = -1;
697  sca_optpropCalc(
698  scatter_matrix_new,
699  pfct_failed,
700  emis_vector_new(0, joker, joker, joker, joker),
701  extinct_matrix_new(0, joker, joker, joker, joker, joker),
702  f_index,
703  scat_data,
704  pnd,
705  stokes_dim,
706  za_grid,
707  quad_weights_new,
708  pfct_method,
709  pfct_aa_grid_size,
710  pfct_threshold,
711  0,
712  verbosity);
713  }
714  }
715 
716  // resize and calc remaining nstream-affected variables:
717  // - in case of surface_rtprop_agenda driven surface: surfreflmat, surfemisvec
718  if (ground_type == "A") // surface_rtprop_agenda driven surface
719  {
720  Tensor5 srm_new(1, nummu_new, stokes_dim, nummu_new, stokes_dim, 0.);
721  Tensor3 sev_new(1, nummu_new, stokes_dim, 0.);
722  surf_optpropCalc(ws,
723  srm_new,
724  sev_new,
726  f_grid[Range(f_index, 1)],
727  za_grid,
728  mu_values_new,
729  quad_weights_new,
730  stokes_dim,
731  surf_altitude);
732  surfreflmat_new = srm_new(0, joker, joker, joker, joker);
733  surfemisvec_new = sev_new(0, joker, joker);
734  }
735  // - up/down_rad (resize only)
736  Tensor3 up_rad_new(num_layers + 1, nummu_new, stokes_dim, 0.);
737  Tensor3 down_rad_new(num_layers + 1, nummu_new, stokes_dim, 0.);
738  //
739  // run radtrano_
740 #pragma omp critical(fortran_rt4)
741  {
742  // Call RT4
743  radtrano_(stokes_dim,
744  nummu_new,
745  nhza,
746  max_delta_tau,
747  quad_type.c_str(),
749  ground_type.c_str(),
750  ground_albedo[f_index],
751  ground_index[f_index],
752  groundreflec.get_c_array(),
753  surfreflmat_new.get_c_array(),
754  surfemisvec_new.get_c_array(),
755  sky_temp,
756  wavelength,
757  num_layers,
758  height.get_c_array(),
759  temperatures.get_c_array(),
760  gas_extinct.get_c_array(),
761  num_scatlayers,
762  scatlayers.get_c_array(),
763  extinct_matrix_new(0, joker, joker, joker, joker, joker)
764  .get_c_array(),
765  emis_vector_new(0, joker, joker, joker, joker).get_c_array(),
766  scatter_matrix_new.get_c_array(),
767  //noutlevels,
768  //outlevels.get_c_array(),
769  mu_values_new.get_c_array(),
770  up_rad_new.get_c_array(),
771  down_rad_new.get_c_array());
772  }
773  // back-interpolate nstream_new fields to nstreams
774  // (possible to use iyCloudboxInterp agenda? nja, not really a good
775  // idea. too much overhead there (checking, 3D+2ang interpol). rather
776  // use interp_order as additional user parameter.
777  // extrapol issues shouldn't occur as we go from finer to coarser
778  // angular grid)
779  // - loop over nummu:
780  // - determine weights per ummu ang (should be valid for both up and
781  // down)
782  // - loop over num_layers and stokes_dim:
783  // - apply weights
784  for (Index j = 0; j < nummu; j++) {
785  GridPosPoly gp_za;
786  if (cos_za_interp) {
787  gridpos_poly(
788  gp_za, mu_values_new, mu_values[j], za_interp_order, 0.5);
789  } else {
790  gridpos_poly(gp_za,
791  za_grid[Range(0, nummu_new)],
792  za_grid_orig[j],
793  za_interp_order,
794  0.5);
795  }
796  Vector itw(gp_za.idx.nelem());
797  interpweights(itw, gp_za);
798 
799  for (Index k = 0; k < num_layers + 1; k++)
800  for (Index ist = 0; ist < stokes_dim; ist++) {
801  up_rad(k, j, ist) = interp(itw, up_rad_new(k, joker, ist), gp_za);
802  down_rad(k, j, ist) =
803  interp(itw, down_rad_new(k, joker, ist), gp_za);
804  }
805  }
806 
807  // reconstruct za_grid
808  za_grid = za_grid_orig;
809  }
810 
811  // RT4 rad output is in wavelength units, nominally in W/(m2 sr um), where
812  // wavelength input is required in um.
813  // FIXME: When using wavelength input in m, output should be in W/(m2 sr
814  // m). However, check this. So, at first we use wavelength in um. Then
815  // change and compare.
816  //
817  // FIXME: if ever we allow the cloudbox to be not directly at the surface
818  // (at atm level #0, respectively), the assigning from up/down_rad to
819  // cloudbox_field needs to checked. there seems some offsetting going on
820  // (test example: TestDOIT.arts. if kept like below, cloudbox_field at
821  // top-of-cloudbox seems to actually be from somewhere within the
822  // cloud(box) indicated by downwelling being to high and downwelling
823  // exhibiting a non-zero polarisation signature (which it wouldn't with
824  // only scalar gas abs above).
825  //
826  Numeric rad_l2f = wavelength / f_grid[f_index];
827  // down/up_rad contain the radiances in order from slant (90deg) to steep
828  // (0 and 180deg, respectively) streams,then the possible extra angle(s).
829  // We need to resort them properly into cloudbox_field, such that order is
830  // from 0 to 180deg.
831  for (Index j = 0; j < nummu; j++) {
832  for (Index ist = 0; ist < stokes_dim; ist++) {
833  for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
834  cloudbox_field(f_index, k + ncboxremoved, 0, 0, nummu + j, 0, ist) =
835  up_rad(num_layers - k, j, ist) * rad_l2f;
837  f_index, k + ncboxremoved, 0, 0, nummu - 1 - j, 0, ist) =
838  down_rad(num_layers - k, j, ist) * rad_l2f;
839  }
840  // To avoid potential numerical problems at interpolation of the field,
841  // we copy the surface field to underground altitudes
842  for (Index k = ncboxremoved - 1; k >= 0; k--) {
843  cloudbox_field(f_index, k, 0, 0, nummu + j, 0, ist) =
844  cloudbox_field(f_index, k + 1, 0, 0, nummu + j, 0, ist);
845  cloudbox_field(f_index, k, 0, 0, nummu - 1 + j, 0, ist) =
846  cloudbox_field(f_index, k + 1, 0, 0, nummu - 1 + j, 0, ist);
847  }
848  }
849  }
850  }
851 }
852 
853 void za_grid_adjust( // Output
854  Vector& za_grid,
855  // Input
856  ConstVectorView mu_values,
857  const Index& nummu) {
858  for (Index j = 0; j < nummu; j++) {
859  za_grid[nummu - 1 - j] = acos(mu_values[j]) * RAD2DEG;
860  za_grid[nummu + j] = 180. - acos(mu_values[j]) * RAD2DEG;
861  //cout << "setting scat_za[" << nummu-1-j << "]=" << za_grid[nummu-1-j]
862  // << " and [" << nummu+j << "]=" << za_grid[nummu+j]
863  // << " from mu[" << j << "]=" << mu_values[j] << "\n";
864  }
865 }
866 
867 void gas_optpropCalc(Workspace& ws,
868  VectorView gas_extinct,
870  ConstVectorView t_profile,
871  ConstMatrixView vmr_profiles,
873  ConstVectorView f_mono) {
874  // Initialization
875  gas_extinct = 0.;
876 
877  const Index Np = p_grid.nelem();
878 
879  assert(gas_extinct.nelem() == Np - 1);
880 
881  // Local variables to be used in agendas
882  Numeric rtp_temperature_local;
883  Numeric rtp_pressure_local;
884  ArrayOfPropagationMatrix propmat_clearsky_local;
885  Vector rtp_vmr_local(vmr_profiles.nrows());
886 
887  const EnergyLevelMap rtp_nlte_local_dummy;
888 
889  // Calculate layer averaged gaseous extinction
890  for (Index i = 0; i < Np - 1; i++) {
891  rtp_pressure_local = 0.5 * (p_grid[i] + p_grid[i + 1]);
892  rtp_temperature_local = 0.5 * (t_profile[i] + t_profile[i + 1]);
893 
894  // Average vmrs
895  for (Index j = 0; j < vmr_profiles.nrows(); j++)
896  rtp_vmr_local[j] = 0.5 * (vmr_profiles(j, i) + vmr_profiles(j, i + 1));
897 
898  const Vector rtp_mag_dummy(3, 0);
899  const Vector ppath_los_dummy;
900 
901  //FIXME: do this right?
902  ArrayOfStokesVector nlte_dummy;
903  // This is right since there should be only clearsky partials
904  ArrayOfPropagationMatrix partial_dummy;
905  ArrayOfStokesVector partial_source_dummy, partial_nlte_dummy;
907  propmat_clearsky_local,
908  nlte_dummy,
909  partial_dummy,
910  partial_source_dummy,
911  partial_nlte_dummy,
913  f_mono, // monochromatic calculation
914  rtp_mag_dummy,
915  ppath_los_dummy,
916  rtp_pressure_local,
917  rtp_temperature_local,
918  rtp_nlte_local_dummy,
919  rtp_vmr_local,
921 
922  //Assuming non-polarized light and only one frequency
923  if (propmat_clearsky_local.nelem()) {
924  gas_extinct[Np - 2 - i] = propmat_clearsky_local[0].Kjj()[0];
925  for (Index j = 1; j < propmat_clearsky_local.nelem(); j++) {
926  gas_extinct[Np - 2 - i] += propmat_clearsky_local[j].Kjj()[0];
927  }
928  }
929  }
930 }
931 
932 void par_optpropCalc(Tensor5View emis_vector,
933  Tensor6View extinct_matrix,
934  //VectorView scatlayers,
936  const Vector& za_grid,
937  const Index& f_index,
938  ConstMatrixView pnd_profiles,
939  ConstVectorView t_profile,
941  const Index& stokes_dim) {
942  // Initialization
943  extinct_matrix = 0.;
944  emis_vector = 0.;
945 
946  const Index Np_cloud = pnd_profiles.ncols();
947  const Index nummu = za_grid.nelem() / 2;
948 
949  assert(emis_vector.nbooks() == Np_cloud - 1);
950  assert(extinct_matrix.nshelves() == Np_cloud - 1);
951 
952  // preparing input data
953  Vector T_array = t_profile[Range(cloudbox_limits[0], Np_cloud)];
954  Matrix dir_array(za_grid.nelem(), 2, 0.);
955  dir_array(joker, 0) = za_grid;
956 
957  // making output containers
958  ArrayOfArrayOfTensor5 ext_mat_Nse;
959  ArrayOfArrayOfTensor4 abs_vec_Nse;
960  ArrayOfArrayOfIndex ptypes_Nse;
961  Matrix t_ok;
962  ArrayOfTensor5 ext_mat_ssbulk;
963  ArrayOfTensor4 abs_vec_ssbulk;
964  ArrayOfIndex ptype_ssbulk;
965  Tensor5 ext_mat_bulk;
966  Tensor4 abs_vec_bulk;
967  Index ptype_bulk;
968 
969  opt_prop_NScatElems(ext_mat_Nse,
970  abs_vec_Nse,
971  ptypes_Nse,
972  t_ok,
973  scat_data,
974  stokes_dim,
975  T_array,
976  dir_array,
977  f_index);
978  opt_prop_ScatSpecBulk(ext_mat_ssbulk,
979  abs_vec_ssbulk,
980  ptype_ssbulk,
981  ext_mat_Nse,
982  abs_vec_Nse,
983  ptypes_Nse,
984  pnd_profiles(joker, joker),
985  t_ok);
986  opt_prop_Bulk(ext_mat_bulk,
987  abs_vec_bulk,
988  ptype_bulk,
989  ext_mat_ssbulk,
990  abs_vec_ssbulk,
991  ptype_ssbulk);
992 
993  // Calculate layer averaged extinction and absorption and sort into RT4-format
994  // data tensors
995  for (Index ipc = 0; ipc < Np_cloud - 1; ipc++) {
996  for (Index fi = 0; fi < abs_vec_bulk.nbooks(); fi++) {
997  for (Index imu = 0; imu < nummu; imu++) {
998  for (Index ist1 = 0; ist1 < stokes_dim; ist1++) {
999  for (Index ist2 = 0; ist2 < stokes_dim; ist2++) {
1000  extinct_matrix(fi, ipc, 0, imu, ist2, ist1) =
1001  .5 * (ext_mat_bulk(fi, ipc, imu, ist1, ist2) +
1002  ext_mat_bulk(fi, ipc + 1, imu, ist1, ist2));
1003  extinct_matrix(fi, ipc, 1, imu, ist2, ist1) =
1004  .5 * (ext_mat_bulk(fi, ipc, nummu + imu, ist1, ist2) +
1005  ext_mat_bulk(fi, ipc + 1, nummu + imu, ist1, ist2));
1006  }
1007  emis_vector(fi, ipc, 0, imu, ist1) =
1008  .5 * (abs_vec_bulk(fi, ipc, imu, ist1) +
1009  abs_vec_bulk(fi, ipc + 1, imu, ist1));
1010  emis_vector(fi, ipc, 1, imu, ist1) =
1011  .5 * (abs_vec_bulk(fi, ipc, nummu + imu, ist1) +
1012  abs_vec_bulk(fi, ipc + 1, nummu + imu, ist1));
1013  }
1014  }
1015  }
1016  }
1017 }
1018 
1019 void sca_optpropCalc( //Output
1020  Tensor6View scatter_matrix,
1021  Index& pfct_failed,
1022  //Input
1023  ConstTensor4View emis_vector,
1024  ConstTensor5View extinct_matrix,
1025  const Index& f_index,
1027  ConstMatrixView pnd_profiles,
1028  const Index& stokes_dim,
1029  const Vector& za_grid,
1030  ConstVectorView quad_weights,
1031  const String& pfct_method,
1032  const Index& pfct_aa_grid_size,
1033  const Numeric& pfct_threshold,
1034  const Index& auto_inc_nstreams,
1035  const Verbosity& verbosity) {
1036  // FIXME: this whole funtions needs revision/optimization regarding
1037  // - temperature dependence (using new-type scat_data)
1038  // - using redundancies in sca_mat data at least for totally random
1039  // orientation particles (are there some for az. random, too? like
1040  // upper/lower hemisphere equiv?) - we might at least have a flag
1041  // (transported down from calling function?) whether we deal exclusively
1042  // with totally random orient particles (and then take some shortcuts...)
1043 
1044  // FIXME: do we have numerical issues, too, here in case of tiny pnd? check
1045  // with Patrick's Disort-issue case.
1046 
1047  // Initialization
1048  scatter_matrix = 0.;
1049 
1050  const Index N_se = pnd_profiles.nrows();
1051  const Index Np_cloud = pnd_profiles.ncols();
1052  const Index nza_rt = za_grid.nelem();
1053 
1054  assert(scatter_matrix.nvitrines() == Np_cloud - 1);
1055 
1056  // Check that total number of scattering elements in scat_data and pnd_profiles
1057  // agree.
1058  // FIXME: this should be done earlier. in calling method. outside freq- and
1059  // other possible loops. rather assert than runtime error.
1060  if (TotalNumberOfElements(scat_data) != N_se) {
1061  ostringstream os;
1062  os << "Total number of scattering elements in scat_data ("
1063  << TotalNumberOfElements(scat_data) << ") and pnd_profiles (" << N_se
1064  << ") disagree.";
1065  throw runtime_error(os.str());
1066  }
1067 
1068  if (pfct_aa_grid_size < 2) {
1069  ostringstream os;
1070  os << "Azimuth grid size for scatt matrix extraction"
1071  << " (*pfct_aa_grid_size*) must be >1.\n"
1072  << "Yours is " << pfct_aa_grid_size << ".\n";
1073  throw runtime_error(os.str());
1074  }
1075  Vector aa_grid;
1076  nlinspace(aa_grid, 0, 180, pfct_aa_grid_size);
1077 
1078  Index i_se_flat = 0;
1079  Tensor5 sca_mat(N_se, nza_rt, nza_rt, stokes_dim, stokes_dim, 0.);
1080  Matrix ext_fixT_spt(N_se, nza_rt, 0.), abs_fixT_spt(N_se, nza_rt, 0.);
1081 
1082  // Precalculate azimuth integration weights for totally randomly oriented
1083  // (they are only determined by pfct_aa_grid_size)
1084  Numeric daa_totrand =
1085  1. / float(pfct_aa_grid_size - 1); // 2*180./360./(pfct_aa_grid_size-1)
1086 
1087  // first we extract Z at one T, integrate the azimuth data at each
1088  // za_inc/za_sca combi (to get the Fourier series 0.th mode), then interpolate
1089  // to the mu/mu' combis we need in RT.
1090  //
1091  // FIXME: are the stokes component assignments correct? for totally random?
1092  // and azimuthally random? (as they are done differently...)
1093  for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++) {
1094  for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
1095  SingleScatteringData ssd = scat_data[i_ss][i_se];
1096  Index this_f_index = (ssd.pha_mat_data.nlibraries() == 1 ? 0 : f_index);
1097  Index i_pfct;
1098  if (pfct_method == "low")
1099  i_pfct = 0;
1100  else if (pfct_method == "high")
1101  i_pfct = ssd.T_grid.nelem() - 1;
1102  else //if( pfct_method=="median" )
1103  i_pfct = ssd.T_grid.nelem() / 2;
1104 
1105  if (ssd.ptype == PTYPE_TOTAL_RND) {
1107  for (Index iza = 0; iza < nza_rt; iza++) {
1108  for (Index sza = 0; sza < nza_rt; sza++) {
1109  Matrix pha_mat_int(stokes_dim, stokes_dim, 0.);
1110  for (Index saa = 0; saa < pfct_aa_grid_size; saa++) {
1112  pha_mat(joker, joker),
1113  ssd.pha_mat_data(
1114  this_f_index, i_pfct, joker, joker, joker, joker, joker),
1115  ssd.za_grid,
1116  ssd.aa_grid,
1117  ssd.ptype,
1118  sza,
1119  saa,
1120  iza,
1121  0,
1122  za_grid,
1123  aa_grid,
1124  verbosity);
1125 
1126  if (saa == 0 || saa == pfct_aa_grid_size - 1)
1127  pha_mat *= (daa_totrand / 2.);
1128  else
1129  pha_mat *= daa_totrand;
1130  pha_mat_int += pha_mat;
1131  }
1132  sca_mat(i_se_flat, iza, sza, joker, joker) = pha_mat_int;
1133  }
1134  ext_fixT_spt(i_se_flat, iza) =
1135  ssd.ext_mat_data(this_f_index, i_pfct, 0, 0, 0);
1136  abs_fixT_spt(i_se_flat, iza) =
1137  ssd.abs_vec_data(this_f_index, i_pfct, 0, 0, 0);
1138  }
1139  } else if (ssd.ptype == PTYPE_AZIMUTH_RND) {
1140  Index nza_se = ssd.za_grid.nelem();
1141  Index naa_se = ssd.aa_grid.nelem();
1142  Tensor4 pha_mat_int(nza_se, nza_se, stokes_dim, stokes_dim, 0.);
1143  ConstVectorView za_datagrid = ssd.za_grid;
1144  ConstVectorView aa_datagrid = ssd.aa_grid;
1145  assert(aa_datagrid[0] == 0.);
1146  assert(aa_datagrid[naa_se - 1] == 180.);
1147  Vector daa(naa_se);
1148 
1149  // Precalculate azimuth integration weights for this azimuthally
1150  // randomly oriented scat element
1151  // (need to do this per scat element as ssd.aa_grid is scat
1152  // element specific (might change between elements) and need to do
1153  // this on actual grid instead of grid number since the grid can,
1154  // at least theoretically be non-equidistant)
1155  daa[0] = (aa_datagrid[1] - aa_datagrid[0]) / 360.;
1156  for (Index saa = 1; saa < naa_se - 1; saa++)
1157  daa[saa] = (aa_datagrid[saa + 1] - aa_datagrid[saa - 1]) / 360.;
1158  daa[naa_se - 1] =
1159  (aa_datagrid[naa_se - 1] - aa_datagrid[naa_se - 2]) / 360.;
1160 
1161  // first, extracting the phase matrix at the scatt elements own
1162  // polar angle grid, deriving their respective azimuthal (Fourier
1163  // series) 0-mode
1164  for (Index iza = 0; iza < nza_se; iza++)
1165  for (Index sza = 0; sza < nza_se; sza++) {
1166  for (Index saa = 0; saa < naa_se; saa++) {
1167  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1168  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1169  pha_mat_int(sza, iza, ist1, ist2) +=
1170  daa[saa] * ssd.pha_mat_data(this_f_index,
1171  i_pfct,
1172  sza,
1173  saa,
1174  iza,
1175  0,
1176  ist1 * 4 + ist2);
1177  }
1178  }
1179 
1180  // second, interpolating the extracted azimuthal mode to the RT4
1181  // solver polar angles
1182  for (Index iza = 0; iza < nza_rt; iza++) {
1183  for (Index sza = 0; sza < nza_rt; sza++) {
1184  GridPos za_sca_gp;
1185  GridPos za_inc_gp;
1186  Vector itw(4);
1187  Matrix pha_mat_lab(stokes_dim, stokes_dim, 0.);
1188  Numeric za_sca = za_grid[sza];
1189  Numeric za_inc = za_grid[iza];
1190 
1191  gridpos(za_inc_gp, za_datagrid, za_inc);
1192  gridpos(za_sca_gp, za_datagrid, za_sca);
1193  interpweights(itw, za_sca_gp, za_inc_gp);
1194 
1195  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1196  for (Index ist2 = 0; ist2 < stokes_dim; ist2++) {
1197  pha_mat_lab(ist1, ist2) =
1198  interp(itw,
1199  pha_mat_int(Range(joker), Range(joker), ist1, ist2),
1200  za_sca_gp,
1201  za_inc_gp);
1202  //if (ist1+ist2==1)
1203  // pha_mat_lab(ist1,ist2) *= -1.;
1204  }
1205 
1206  sca_mat(i_se_flat, iza, sza, joker, joker) = pha_mat_lab;
1207  }
1208  ext_fixT_spt(i_se_flat, iza) =
1209  ssd.ext_mat_data(this_f_index, i_pfct, iza, 0, 0);
1210  abs_fixT_spt(i_se_flat, iza) =
1211  ssd.abs_vec_data(this_f_index, i_pfct, iza, 0, 0);
1212  }
1213  } else {
1214  ostringstream os;
1215  os << "Unsuitable particle type encountered.";
1216  throw runtime_error(os.str());
1217  }
1218  i_se_flat++;
1219  }
1220  }
1221 
1222  assert(i_se_flat == N_se);
1223  // now we sum up the Z(mu,mu') over the scattering elements weighted by the
1224  // pnd_field data, deriving Z(z,mu,mu') and sorting this into
1225  // scatter_matrix
1226  // FIXME: it seems that at least for totally random, corresponding upward and
1227  // downward directions have exactly the same (0,0) elements. Can we use this
1228  // to reduce calc efforts (for sca and its normalization as well as for ext
1229  // and abs?)
1230  Index nummu = nza_rt / 2;
1231  for (Index scat_p_index_local = 0; scat_p_index_local < Np_cloud - 1;
1232  scat_p_index_local++) {
1233  Vector ext_fixT(nza_rt, 0.), abs_fixT(nza_rt, 0.);
1234 
1235  for (Index i_se = 0; i_se < N_se; i_se++) {
1236  Numeric pnd_mean = 0.5 * (pnd_profiles(i_se, scat_p_index_local + 1) +
1237  pnd_profiles(i_se, scat_p_index_local));
1238  if (pnd_mean != 0.)
1239  for (Index iza = 0; iza < nummu; iza++) {
1240  // JM171003: not clear to me anymore why this check. if pnd_mean
1241  // is non-zero, then extinction should also be non-zero by
1242  // default?
1243  //if ( (extinct_matrix(scat_p_index_local,0,iza,0,0)+
1244  // extinct_matrix(scat_p_index_local,1,iza,0,0)) > 0. )
1245  for (Index sza = 0; sza < nummu; sza++) {
1246  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1247  for (Index ist2 = 0; ist2 < stokes_dim; ist2++) {
1248  // we can't use stokes_dim jokers here since '*' doesn't
1249  // exist for Num*MatView. Also, order of stokes matrix
1250  // dimensions is inverted here (aka scat matrix is
1251  // transposed).
1252  scatter_matrix(scat_p_index_local, 0, iza, ist2, sza, ist1) +=
1253  pnd_mean * sca_mat(i_se, iza, sza, ist1, ist2);
1254  scatter_matrix(scat_p_index_local, 1, iza, ist2, sza, ist1) +=
1255  pnd_mean * sca_mat(i_se, nummu + iza, sza, ist1, ist2);
1256  scatter_matrix(scat_p_index_local, 2, iza, ist2, sza, ist1) +=
1257  pnd_mean * sca_mat(i_se, iza, nummu + sza, ist1, ist2);
1258  scatter_matrix(scat_p_index_local, 3, iza, ist2, sza, ist1) +=
1259  pnd_mean *
1260  sca_mat(i_se, nummu + iza, nummu + sza, ist1, ist2);
1261  }
1262  }
1263 
1264  ext_fixT[iza] += pnd_mean * ext_fixT_spt(i_se, iza);
1265  abs_fixT[iza] += pnd_mean * abs_fixT_spt(i_se, iza);
1266  }
1267  }
1268 
1269  for (Index iza = 0; iza < nummu; iza++) {
1270  for (Index ih = 0; ih < 2; ih++) {
1271  if (extinct_matrix(scat_p_index_local, ih, iza, 0, 0) > 0.) {
1272  Numeric sca_mat_integ = 0.;
1273 
1274  // We need to calculate the nominal values for the fixed T, we
1275  // used above in the pha_mat extraction. Only this tells us whether
1276  // angular resolution is sufficient.
1277  //
1278  //Numeric ext_nom = extinct_matrix(scat_p_index_local,ih,iza,0,0);
1279  //Numeric sca_nom = ext_nom-emis_vector(scat_p_index_local,ih,iza,0);
1280  Numeric ext_nom = ext_fixT[iza];
1281  Numeric sca_nom = ext_nom - abs_fixT[iza];
1282  Numeric w0_nom = sca_nom / ext_nom;
1283  assert(w0_nom >= 0.);
1284 
1285  for (Index sza = 0; sza < nummu; sza++) {
1286  // SUM2 = SUM2 + 2.0D0*PI*QUAD_WEIGHTS(J2)*
1287  // . ( SCATTER_MATRIX(1,J2,1,J1, L,TSL)
1288  // . + SCATTER_MATRIX(1,J2,1,J1, L+2,TSL) )
1289  sca_mat_integ +=
1290  quad_weights[sza] *
1291  (scatter_matrix(scat_p_index_local, ih, iza, 0, sza, 0) +
1292  scatter_matrix(scat_p_index_local, ih + 2, iza, 0, sza, 0));
1293  }
1294 
1295  // compare integrated scatt matrix with ext-abs for respective
1296  // incident polar angle - consistently with scat_dataCheck, we do
1297  // this in terms of albedo deviation (since PFCT deviations at
1298  // small albedos matter less than those at high albedos)
1299  // SUM1 = EMIS_VECTOR(1,J1,L,TSL)-EXTINCT_MATRIX(1,1,J1,L,TSL)
1300  Numeric w0_act = 2. * PI * sca_mat_integ / ext_nom;
1301  Numeric pfct_norm = 2. * PI * sca_mat_integ / sca_nom;
1302  /*
1303  cout << "sca_mat norm deviates " << 1e2*abs(1.-pfct_norm) << "%"
1304  << " (" << abs(w0_act-w0_nom) << " in albedo).\n";
1305  */
1306  Numeric sca_nom_paropt =
1307  extinct_matrix(scat_p_index_local, ih, iza, 0, 0) -
1308  emis_vector(scat_p_index_local, ih, iza, 0);
1309  //Numeric w0_nom_paropt = sca_nom_paropt /
1310  // extinct_matrix(scat_p_index_local,ih,iza,0,0);
1311  /*
1312  cout << "scat_p=" << scat_p_index_local
1313  << ", iza=" << iza << ", hem=" << ih << "\n";
1314  cout << " scaopt (paropt) w0_act= " << w0_act
1315  << ", w0_nom = " << w0_nom
1316  << " (" << w0_nom_paropt
1317  << "), diff=" << w0_act-w0_nom
1318  << " (" << w0_nom-w0_nom_paropt << ").\n";
1319  */
1320 
1321  if (abs(w0_act - w0_nom) > pfct_threshold) {
1322  if (pfct_failed >= 0) {
1323  if (auto_inc_nstreams) {
1324  pfct_failed = 1;
1325  return;
1326  } else {
1327  ostringstream os;
1328  os << "Bulk scattering matrix normalization deviates significantly\n"
1329  << "from expected value (" << 1e2 * abs(1. - pfct_norm)
1330  << "%,"
1331  << " resulting in albedo deviation of "
1332  << abs(w0_act - w0_nom) << ").\n"
1333  << "Something seems wrong with your scattering data "
1334  << " (did you run *scat_dataCheck*?)\n"
1335  << "or your RT4 setup (try increasing *nstreams* and in case"
1336  << " of randomly oriented particles possibly also"
1337  << " pfct_aa_grid_size).";
1338  throw runtime_error(os.str());
1339  }
1340  }
1341  } else if (abs(w0_act - w0_nom) > pfct_threshold * 0.1 ||
1342  abs(1. - pfct_norm) > 1e-2) {
1343  CREATE_OUT2;
1344  out2
1345  << "Warning: The bulk scattering matrix is not well normalized\n"
1346  << "Deviating from expected value by "
1347  << 1e2 * abs(1. - pfct_norm) << "% (and "
1348  << abs(w0_act - w0_nom) << " in terms of scattering albedo).\n";
1349  }
1350  // rescale scattering matrix to expected (0,0) value (and scale all
1351  // other elements accordingly)
1352  //
1353  // However, here we should not use the pfct_norm based on the
1354  // deviation from the fixed-temperature ext and abs. Instead, for
1355  // energy conservation reasons, this needs to be consistent with
1356  // extinct_matrix and emis_vector. Hence, recalculate pfct_norm
1357  // from them.
1358  //cout << " scaopt (paropt) pfct_norm dev = " << 1e2*pfct_norm-1e2;
1359  pfct_norm = 2. * PI * sca_mat_integ / sca_nom_paropt;
1360  //cout << " (" << 1e2*pfct_norm-1e2 << ")%.\n";
1361  //
1362  // FIXME: not fully clear whether applying the same rescaling
1363  // factor to all stokes elements is the correct way to do. check
1364  // out, e.g., Vasilieva (JQSRT 2006) for better approaches.
1365  // FIXME: rather rescale Z(0,0) only as we should have less issues
1366  // for other Z elements as they are less peaked/featured.
1367  scatter_matrix(scat_p_index_local, ih, iza, joker, joker, joker) /=
1368  pfct_norm;
1369  scatter_matrix(
1370  scat_p_index_local, ih + 2, iza, joker, joker, joker) /=
1371  pfct_norm;
1372  //if (scat_p_index_local==49)
1373  }
1374  }
1375  }
1376  }
1377 }
1378 
1379 void surf_optpropCalc(Workspace& ws,
1380  //Output
1381  Tensor5View surf_refl_mat,
1382  Tensor3View surf_emis_vec,
1383  //Input
1387  ConstVectorView mu_values,
1388  ConstVectorView quad_weights,
1389  const Index& stokes_dim,
1390  const Numeric& surf_alt) {
1391  // While proprietary RT4 - from the input/user control side - handles only
1392  // Lambertian and Fresnel, the Doubling&Adding solver core applies a surface
1393  // reflection matrix and a surface radiance term. The reflection matrix is
1394  // dependent on incident and reflection direction, ie forms a discrete
1395  // representation of the bidirectional reflection function; the radiance term
1396  // is dependent on outgoing polar angle. That is, the solver core can
1397  // basically handle ANY kind of surface reflection type.
1398  //
1399  // Here, we replace RT4's proprietary surface reflection and radiance term
1400  // calculation and use ARTS' surface_rtprop_agenda instead to set these
1401  // variables up.
1402  // That is, ARTS-RT4 is set up such that it can handle all surface reflection
1403  // types that ARTS itself can handle.
1404  // What is required here, is to derive the reflection matrix over all incident
1405  // and reflected polar angle directions. The surface_rtprop_agenda handles one
1406  // reflected direction at a time (rtp_los); it is sufficient here to simply
1407  // loop over all reflected angles as given by the stream directions. However,
1408  // the incident directions (surface_los) provided by surface_rtprop_agenda
1409  // depends on the specific, user-defined setup of the agenda. Out of what the
1410  // agenda call provides, we have to derive the reflection matrix entries for
1411  // the incident directions as defined by the RT4 stream directions. To be
1412  // completely general (handling everything from specular via semi-specular to
1413  // Lambertian) is a bit tricky. Here we decided to allow the ARTS-standard
1414  // 0.5-grid-spacing extrapolation and set everything outside that range to
1415  // zero.
1416  //
1417  // We do all frequencies here at once (assuming this is the faster variant as
1418  // the agenda anyway (can) provide output for full f_grid at once and as we
1419  // have to apply the same inter/extrapolation to all the frequencies).
1420  //
1421  // FIXME: Allow surface to be elsewhere than at lowest atm level (this
1422  // requires changes in the surface setting part and more extensive ones in the
1423  // atmospheric optical property prep part within the frequency loop further
1424  // below).
1425 
1426  chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
1427 
1428  const Index nf = f_grid.nelem();
1429  const Index nummu = za_grid.nelem() / 2;
1430  const String B_unit = "R";
1431 
1432  // Local input of surface_rtprop_agenda.
1433  Vector rtp_pos(1, surf_alt); //atmosphere_dim is 1
1434 
1435  for (Index rmu = 0; rmu < nummu; rmu++) {
1436  // Local output of surface_rtprop_agenda.
1441 
1442  // rtp_los is reflected direction, ie upwelling direction, which is >90deg
1443  // in ARTS. za_grid is sorted here as downwelling (90->0) in 1st
1444  // half, upwelling (90->180) in second half. that is, here we have to take
1445  // the second half grid or, alternatively, use 180deg-za[imu].
1446  Vector rtp_los(1, za_grid[nummu + rmu]);
1447 
1451  surface_los,
1453  f_grid,
1454  rtp_pos,
1455  rtp_los,
1457  Index nsl = surface_los.nrows();
1458  assert(surface_los.ncols() == 1 || nsl == 0);
1459 
1460  // ARTS' surface_emission is equivalent to RT4's gnd_radiance (here:
1461  // surf_emis_vec) except for ARTS using Planck in terms of frequency,
1462  // while RT4 uses it in terms of wavelengths.
1463  // To derive gnd_radiance, we can either use ARTS surface_emission and
1464  // rescale it by B_lambda(surface_skin_t)/B_freq(surface_skin_t). Or we
1465  // can create it from surface_rmatrix and B_lambda(surface_skin_t). The
1466  // latter is more direct regarding use of B(T), but is requires
1467  // (re-)deriving diffuse reflectivity stokes components (which should be
1468  // the sum of the incident polar angle dependent reflectivities. so, that
1469  // might ensure better consistency. but then we should calculate
1470  // gnd_radiance only after surface_los to za_grid conversion of the
1471  // reflection matrix). For now, we use the rescaling approach.
1472  for (Index f_index = 0; f_index < nf; f_index++) {
1473  Numeric freq = f_grid[f_index];
1474  Numeric B_freq = planck(freq, surface_skin_t);
1475  Numeric B_lambda;
1476  Numeric wave = 1e6 * SPEED_OF_LIGHT / freq;
1477  planck_function_(surface_skin_t, B_unit.c_str(), wave, B_lambda);
1478  Numeric B_ratio = B_lambda / B_freq;
1479  surf_emis_vec(f_index, rmu, joker) = surface_emission(f_index, joker);
1480  surf_emis_vec(f_index, rmu, joker) *= B_ratio;
1481  }
1482 
1483  // now we have to properly fill the RT4 stream directions in incident
1484  // angle dimension of the reflection matrix for RT4 with the agenda
1485  // output, which is given on/limited to surface_los. Only values inside
1486  // the (default) extrapolation range are filled; the rest is left as 0.
1487  // we do this for each incident stream separately as we need to check them
1488  // separately whether they are within allowed inter/extrapolation range
1489  // (ARTS can't do this for all elements of a vector at once. it will fail
1490  // for the whole vector if there's a single value outside.).
1491 
1492  // as we are rescaling surface_rmatrix within here, we need to keep its
1493  // original normalization for later renormalization. Create the container
1494  // here, fill in below.
1495  Vector R_arts(f_grid.nelem(), 0.);
1496 
1497  if (nsl > 1) // non-blackbody, non-specular reflection
1498  {
1499  for (Index f_index = 0; f_index < nf; f_index++)
1500  R_arts[f_index] = surface_rmatrix(joker, f_index, 0, 0).sum();
1501 
1502  // Determine angle range weights in surface_rmatrix and de-scale
1503  // surface_rmatrix with those.
1504  Vector surf_int_grid(nsl + 1);
1505  surf_int_grid[0] =
1506  surface_los(0, 0) - 0.5 * (surface_los(1, 0) - surface_los(0, 0));
1507  surf_int_grid[nsl] =
1508  surface_los(nsl - 1, 0) +
1509  0.5 * (surface_los(nsl - 1, 0) - surface_los(nsl - 2, 0));
1510  for (Index imu = 1; imu < nsl; imu++)
1511  surf_int_grid[imu] =
1512  0.5 * (surface_los(imu - 1, 0) + surface_los(imu, 0));
1513  surf_int_grid *= DEG2RAD;
1514  for (Index imu = 0; imu < nsl; imu++) {
1515  //Numeric coslow = cos(2.*surf_int_grid[imu]);
1516  //Numeric coshigh = cos(2.*surf_int_grid[imu+1]);
1517  //Numeric w = 0.5*(coslow-coshigh);
1518  Numeric w = 0.5 * (cos(2. * surf_int_grid[imu]) -
1519  cos(2. * surf_int_grid[imu + 1]));
1520  surface_rmatrix(imu, joker, joker, joker) /= w;
1521  }
1522 
1523  // Testing: interpolation in cos(za)
1524  //Vector mu_surf_los(nsl);
1525  //for (Index imu=0; imu<nsl; imu++)
1526  // mu_surf_los[imu] = cos(surface_los(imu,0)*DEG2RAD);
1527 
1528  for (Index imu = 0; imu < nummu; imu++) {
1529  try {
1530  GridPos gp_za;
1531  gridpos(gp_za, surface_los(joker, 0), za_grid[imu]);
1532  // Testing: interpolation in cos(za)
1533  //gridpos( gp_za, mu_surf_los, mu_values[imu] );
1534  Vector itw(2);
1535  interpweights(itw, gp_za);
1536 
1537  // now apply the interpolation weights. since we're not in
1538  // python, we can't do the whole tensor at once, but need to
1539  // loop over the dimensions.
1540  for (Index f_index = 0; f_index < nf; f_index++)
1541  for (Index sto1 = 0; sto1 < stokes_dim; sto1++)
1542  for (Index sto2 = 0; sto2 < stokes_dim; sto2++)
1543  surf_refl_mat(f_index, imu, sto2, rmu, sto1) = interp(
1544  itw, surface_rmatrix(joker, f_index, sto1, sto2), gp_za);
1545  // Apply new angle range weights - as this is for RT4, we apply
1546  // the actual RT4 angle (aka quadrature) weights:
1547  surf_refl_mat(joker, imu, joker, rmu, joker) *=
1548  (quad_weights[imu] * mu_values[imu]);
1549  } catch (const runtime_error& e) {
1550  // nothing to do here. we just leave the reflection matrix
1551  // entry at the 0.0 it was initialized with.
1552  }
1553  }
1554  } else if (nsl > 0) // specular reflection
1555  // no interpolation, no angle weight rescaling,
1556  // just setting diagonal elements of surf_refl_mat.
1557  {
1558  for (Index f_index = 0; f_index < nf; f_index++)
1559  R_arts[f_index] = surface_rmatrix(joker, f_index, 0, 0).sum();
1560 
1561  // surface_los angle should be identical to
1562  // 180. - (rtp_los=za_grid[nummu+rmu]=180.-za_grid[rmu])
1563  // = za_grid[rmu].
1564  // check that, and if so, sort values into refmat(rmu,rmu).
1565  assert(is_same_within_epsilon(surface_los(0, 0), za_grid[rmu], 1e-12));
1566  for (Index f_index = 0; f_index < nf; f_index++)
1567  for (Index sto1 = 0; sto1 < stokes_dim; sto1++)
1568  for (Index sto2 = 0; sto2 < stokes_dim; sto2++)
1569  surf_refl_mat(f_index, rmu, sto2, rmu, sto1) =
1570  surface_rmatrix(0, f_index, sto1, sto2);
1571  }
1572  //else {} // explicit blackbody
1573  // all surf_refl_mat elements to remain at 0., ie nothing to do.
1574 
1575  //eventually make sure the scaling of surf_refl_mat is correct
1576  for (Index f_index = 0; f_index < nf; f_index++) {
1577  Numeric R_scale = 1.;
1578  Numeric R_rt4 = surf_refl_mat(f_index, joker, 0, rmu, 0).sum();
1579  if (R_rt4 == 0.) {
1580  if (R_arts[f_index] != 0.) {
1581  ostringstream os;
1582  os << "Something went wrong.\n"
1583  << "At reflected stream #" << rmu
1584  << ", power reflection coefficient for RT4\n"
1585  << "became 0, although the one from surface_rtprop_agenda is "
1586  << R_arts << ".\n";
1587  throw runtime_error(os.str());
1588  }
1589  } else {
1590  R_scale = R_arts[f_index] / R_rt4;
1591  surf_refl_mat(f_index, joker, joker, rmu, joker) *= R_scale;
1592  }
1593  }
1594  }
1595 }
1596 
1597 void rt4_test(Tensor4& out_rad,
1598  const String& datapath,
1599  const Verbosity& verbosity) {
1600  //emissivity.resize(4);
1601  //reflectivity.resize(4);
1602 
1603  Index nstokes = 2;
1604  Index nummu = 8;
1605  Index nuummu = 0;
1606  Numeric max_delta_tau = 1.0E-6;
1607  String quad_type = "L";
1608  Numeric ground_temp = 300.;
1609  String ground_type = "L";
1610  Numeric ground_albedo = 0.05;
1611  Complex ground_index;
1612  Numeric sky_temp = 0.;
1613  Numeric wavelength = 880.;
1614  //Index noutlevels=1;
1615  //ArrayOfIndex outlevels(1);
1616  //outlevels[0]=1;
1617 
1618  Vector height, temperatures, gas_extinct;
1619  Tensor5 sca_data;
1620  Tensor4 ext_data;
1621  Tensor3 abs_data;
1622  ReadXML(height, "height", datapath + "z.xml", "", verbosity);
1623  ReadXML(temperatures, "temperatures", datapath + "T.xml", "", verbosity);
1624  ReadXML(gas_extinct, "gas_extinct", datapath + "abs_gas.xml", "", verbosity);
1625  ReadXML(abs_data, "abs_data", datapath + "abs_par.xml", "", verbosity);
1626  ReadXML(ext_data, "ext_data", datapath + "ext_par.xml", "", verbosity);
1627  ReadXML(sca_data, "sca_data", datapath + "sca_par.xml", "", verbosity);
1628  Index num_layers = height.nelem() - 1;
1629  Index num_scatlayers = 3;
1630  Vector scatlayers(num_layers, 0.);
1631  scatlayers[3] = 1.;
1632  scatlayers[4] = 2.;
1633  scatlayers[5] = 3.;
1634 
1635  // the read in sca/ext/abs_data is the complete set (and it's in the wrong
1636  // order for passing it directly to radtrano). before handing over to
1637  // fortran, we need to reduce it to the number of stokes elements to be
1638  // used. we can't use views here as all data needs to be continuous in
1639  // memory; that is, we have to explicitly copy the data we need.
1640  Tensor6 scatter_matrix(num_scatlayers, 4, nummu, nstokes, nummu, nstokes);
1641  for (Index ii = 0; ii < 4; ii++)
1642  for (Index ij = 0; ij < nummu; ij++)
1643  for (Index ik = 0; ik < nstokes; ik++)
1644  for (Index il = 0; il < nummu; il++)
1645  for (Index im = 0; im < nstokes; im++)
1646  for (Index in = 0; in < num_scatlayers; in++)
1647  scatter_matrix(in, ii, ij, ik, il, im) =
1648  sca_data(im, il, ik, ij, ii);
1649  Tensor5 extinct_matrix(num_scatlayers, 2, nummu, nstokes, nstokes);
1650  for (Index ii = 0; ii < 2; ii++)
1651  for (Index ij = 0; ij < nummu; ij++)
1652  for (Index ik = 0; ik < nstokes; ik++)
1653  for (Index il = 0; il < nstokes; il++)
1654  for (Index im = 0; im < num_scatlayers; im++)
1655  extinct_matrix(im, ii, ij, ik, il) = ext_data(il, ik, ij, ii);
1656  Tensor4 emis_vector(num_scatlayers, 2, nummu, nstokes);
1657  for (Index ii = 0; ii < 2; ii++)
1658  for (Index ij = 0; ij < nummu; ij++)
1659  for (Index ik = 0; ik < nstokes; ik++)
1660  for (Index il = 0; il < num_scatlayers; il++)
1661  emis_vector(il, ii, ij, ik) = abs_data(ik, ij, ii);
1662 
1663  // dummy parameters necessary due to modified, flexible surface handling
1664  Tensor4 surf_refl_mat(nummu, nstokes, nummu, nstokes, 0.);
1665  Matrix surf_emis_vec(nummu, nstokes, 0.);
1666  Matrix ground_reflec(nstokes, nstokes, 0.);
1667 
1668  // Output variables
1669  Vector mu_values(nummu);
1670  Tensor3 up_rad(num_layers + 1, nummu, nstokes, 0.);
1671  Tensor3 down_rad(num_layers + 1, nummu, nstokes, 0.);
1672 
1673  radtrano_(nstokes,
1674  nummu,
1675  nuummu,
1676  max_delta_tau,
1677  quad_type.c_str(),
1678  ground_temp,
1679  ground_type.c_str(),
1680  ground_albedo,
1681  ground_index,
1682  ground_reflec.get_c_array(),
1683  surf_refl_mat.get_c_array(),
1684  surf_emis_vec.get_c_array(),
1685  sky_temp,
1686  wavelength,
1687  num_layers,
1688  height.get_c_array(),
1689  temperatures.get_c_array(),
1690  gas_extinct.get_c_array(),
1691  num_scatlayers,
1692  scatlayers.get_c_array(),
1693  extinct_matrix.get_c_array(),
1694  emis_vector.get_c_array(),
1695  scatter_matrix.get_c_array(),
1696  //noutlevels,
1697  //outlevels.get_c_array(),
1698  mu_values.get_c_array(),
1699  up_rad.get_c_array(),
1700  down_rad.get_c_array());
1701 
1702  //so far, output is in
1703  // units W/m^2 um sr
1704  // dimensions [numlayers+1,nummu,nstokes]
1705  // sorted from high to low (altitudes) and 0 to |1| (mu)
1706  //WriteXML( "ascii", up_rad, "up_rad.xml", 0, "up_rad", "", "", verbosity );
1707  //WriteXML( "ascii", down_rad, "down_rad.xml", 0, "down_rad", "", "", verbosity );
1708 
1709  //to be able to compare with RT4 reference results, reshape output into
1710  //RT4-output type table (specifically, resort up_rad such that it runs from
1711  //zenith welling to horizontal, thus forms a continuous angle grid with
1712  //down_rad. if later changing up_rad/down_rad sorting such that it is in
1713  //line with cloudbox_field, then this has to be adapted as well...
1714  out_rad.resize(num_layers + 1, 2, nummu, nstokes);
1715  for (Index ii = 0; ii < nummu; ii++)
1716  out_rad(joker, 0, ii, joker) = up_rad(joker, nummu - 1 - ii, joker);
1717  //out_rad(joker,0,joker,joker) = up_rad;
1718  out_rad(joker, 1, joker, joker) = down_rad;
1719  //WriteXML( "ascii", out_rad, "out_rad.xml", 0, "out_rad", "", "", verbosity );
1720 }
1721 
1722 #endif /* ENABLE_RT4 */
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
is_same_within_epsilon
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:351
Tensor4::resize
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1064
ConstTensor6View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackVI.cc:45
ARTS::Var::rtp_los
Vector rtp_los(Workspace &ws) noexcept
Definition: autoarts.h:5715
ConstTensor5View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackV.cc:47
SingleScatteringData::za_grid
Vector za_grid
Definition: optproperties.h:85
PTYPE_GENERAL
@ PTYPE_GENERAL
Definition: optproperties.h:54
Tensor6::resize
void resize(Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVI.cc:2175
ARTS::Var::atmgeom_checked
Index atmgeom_checked(Workspace &ws) noexcept
Definition: autoarts.h:2489
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
complex_n_interp
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.
Definition: special_interp.cc:926
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
RAD2DEG
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
joker
const Joker joker
ARTS::Var::surface_rmatrix
Tensor4 surface_rmatrix(Workspace &ws) noexcept
Definition: autoarts.h:6800
interpolation.h
Header file for interpolation.cc.
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ComplexVector
The ComplexVector class.
Definition: complex.h:573
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
chk_not_empty
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:694
ARTS::Var::f_index
Index f_index(Workspace &ws) noexcept
Definition: autoarts.h:3464
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
SingleScatteringData::aa_grid
Vector aa_grid
Definition: optproperties.h:86
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
COSMIC_BG_TEMP
const Numeric COSMIC_BG_TEMP
Global constant, Planck temperature for cosmic background radiation [K].
opt_prop_ScatSpecBulk
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
Definition: optproperties.cc:148
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
ArrayOfRetrievalQuantity
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:402
DEG2RAD
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
reduced_1datm
void reduced_1datm(Vector &p, Vector &z, Vector &t, Matrix &vmr, Matrix &pnd, ArrayOfIndex &cboxlims, Index &ncboxremoved, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfIndex &cloudbox_limits)
reduced_1datm
Definition: disort.cc:693
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
ARTS::Var::surface_rtprop_agenda
Agenda surface_rtprop_agenda(Workspace &ws) noexcept
Definition: autoarts.h:6807
CREATE_OUT1
#define CREATE_OUT1
Definition: messages.h:205
Agenda
The Agenda class.
Definition: agenda_class.h:44
GriddedField3
Definition: gridded_fields.h:284
SingleScatteringData
Definition: optproperties.h:80
ARTS::Var::pha_mat
Tensor4 pha_mat(Workspace &ws) noexcept
Definition: autoarts.h:4881
ConstTensor7View::nlibraries
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:42
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:133
opt_prop_NScatElems
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
Definition: optproperties.cc:248
complex.h
A class implementing complex numbers for ARTS.
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
Array
This can be used to make arrays out of anything.
Definition: array.h:108
ReadXML
void ReadXML(T &v, const String &v_name, const String &f, const String &f_name, const Verbosity &verbosity)
WORKSPACE METHOD: ReadXML.
Definition: m_xml.h:41
PTypeToString
String PTypeToString(const PType &ptype)
Convert particle type enum value to String.
Definition: optproperties.cc:2512
ConstTensor6View::nvitrines
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVI.cc:42
min
#define min(a, b)
Definition: legacy_continua.cc:20628
SingleScatteringData::ptype
PType ptype
Definition: optproperties.h:81
Tensor5::resize
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1743
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
GridPosPoly
Structure to store a grid position for higher order interpolation.
Definition: interpolation_poly.h:64
my_basic_string< char >
Tensor5View::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackV.cc:675
SingleScatteringData::ext_mat_data
Tensor5 ext_mat_data
Definition: optproperties.h:88
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
surface_rtprop_agendaExecute
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:25033
physics_funcs.h
This file contains declerations of functions of physical character.
ARTS::Var::atmfields_checked
Index atmfields_checked(Workspace &ws) noexcept
Definition: autoarts.h:2473
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
optproperties.h
Scattering database structure and functions.
GridPosPoly::idx
ArrayOfIndex idx
Definition: interpolation_poly.h:67
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
Tensor5View
The Tensor5View class.
Definition: matpackV.h:333
propmat_clearsky_agendaExecute
void propmat_clearsky_agendaExecute(Workspace &ws, ArrayOfPropagationMatrix &propmat_clearsky, ArrayOfStokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_dx_source, ArrayOfStokesVector &nlte_dsource_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:23513
ARTS::Var::surface_scalar_reflectivity
Vector surface_scalar_reflectivity(Workspace &ws) noexcept
Definition: autoarts.h:6864
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
Tensor6View::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackVI.cc:1734
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
EnergyLevelMap
Definition: energylevelmap.h:60
max
#define max(a, b)
Definition: legacy_continua.cc:20629
ConstComplexVectorView
A constant view of a ComplexVector.
Definition: complex.h:268
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
m_xml.h
Workspace methods and template functions for supergeneric XML IO.
ARTS::Var::za_grid
Vector za_grid(Workspace &ws) noexcept
Definition: autoarts.h:7771
ARTS::Var::surface_emission
Matrix surface_emission(Workspace &ws) noexcept
Definition: autoarts.h:6700
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
PTYPE_AZIMUTH_RND
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:55
pha_matTransform
void pha_matTransform(MatrixView pha_mat_lab, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const PType &ptype, const Index &za_sca_idx, const Index &aa_sca_idx, const Index &za_inc_idx, const Index &aa_inc_idx, ConstVectorView za_grid, ConstVectorView aa_grid, const Verbosity &verbosity)
Transformation of phase matrix.
Definition: optproperties.cc:1692
ARTS::Var::surface_complex_refr_index
GriddedField3 surface_complex_refr_index(Workspace &ws) noexcept
Definition: autoarts.h:6688
Tensor5
The Tensor5 class.
Definition: matpackV.h:506
ARTS::Var::scat_data_checked
Index scat_data_checked(Workspace &ws) noexcept
Definition: autoarts.h:5850
planck
Numeric planck(const Numeric &f, const Numeric &t)
planck
Definition: physics_funcs.cc:269
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
Range
The range class.
Definition: matpackI.h:160
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
opt_prop_Bulk
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
Definition: optproperties.cc:102
ConstTensor5View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackV.cc:44
disort.h
Functions for disort interface.
rte.h
Declaration of functions in rte.cc.
Workspace
Workspace class.
Definition: workspace_ng.h:40
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
SingleScatteringData::abs_vec_data
Tensor5 abs_vec_data
Definition: optproperties.h:89
config.h
ARTS::Var::cloudbox_field
Tensor7 cloudbox_field(Workspace &ws) noexcept
Definition: autoarts.h:2676
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
ARTS::Var::propmat_clearsky_agenda
Agenda propmat_clearsky_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5405
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
chk_if_in_range
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
ARTS::Var::surface_skin_t
Numeric surface_skin_t(Workspace &ws) noexcept
Definition: autoarts.h:6877
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
gridpos_poly
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.
Definition: interpolation_poly.cc:120
MatrixView::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackI.cc:735
ARTS::Var::surface_los
Matrix surface_los(Workspace &ws) noexcept
Definition: autoarts.h:6714
Tensor3View::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackIII.cc:321
Vector
The Vector class.
Definition: matpackI.h:860
SingleScatteringData::pha_mat_data
Tensor7 pha_mat_data
Definition: optproperties.h:87
ARTS::Var::aa_grid
Vector aa_grid(Workspace &ws) noexcept
Definition: autoarts.h:1717
PTYPE_TOTAL_RND
@ PTYPE_TOTAL_RND
Definition: optproperties.h:56
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
SingleScatteringData::T_grid
Vector T_grid
Definition: optproperties.h:84
Tensor4View::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackIV.cc:358
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
ARTS::Var::cloudbox_checked
Index cloudbox_checked(Workspace &ws) noexcept
Definition: autoarts.h:2650
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
ARTS::Var::surface_reflectivity
Tensor3 surface_reflectivity(Workspace &ws) noexcept
Definition: autoarts.h:6778
ConstTensor5View
A constant view of a Tensor5.
Definition: matpackV.h:143
VectorView::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array, const-version.
Definition: matpackI.cc:272
PI
const Numeric PI
Global constant, pi.
Tensor7
The Tensor7 class.
Definition: matpackVII.h:2382
ARTS::Var::rtp_pos
Vector rtp_pos(Workspace &ws) noexcept
Definition: autoarts.h:5774
Tensor6View
The Tensor6View class.
Definition: matpackVI.h:621
rt4.h
Contains functions related to application of scattering solver RT4.
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782