ARTS  2.3.1285(git:92a29ea9-dirty)
m_doit.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012
2  Claudia Emde <claudia.emde@dlr.de>
3  Sreerekha T.R. <rekha@uni-bremen.de>
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 
34 /*===========================================================================
35  === External declarations
36  ===========================================================================*/
37 
38 #include <cmath>
39 #include <cstdlib>
40 #include <iostream>
41 #include <stdexcept>
42 #include "agenda_class.h"
43 #include "array.h"
44 #include "arts.h"
45 #include "auto_md.h"
46 #include "check_input.h"
47 #include "doit.h"
48 #include "geodetic.h"
49 #include "lin_alg.h"
50 #include "logic.h"
51 #include "m_general.h"
52 #include "math_funcs.h"
53 #include "matpackVII.h"
54 #include "messages.h"
55 #include "physics_funcs.h"
56 #include "ppath.h"
57 #include "rte.h"
58 #include "special_interp.h"
59 #include "wsv_aux.h"
60 #include "xml_io.h"
61 
62 extern const Numeric PI;
63 extern const Numeric RAD2DEG;
64 
65 /*===========================================================================
66  === The functions (in alphabetical order)
67  ===========================================================================*/
68 
69 /* Workspace method: Doxygen documentation will be auto-generated */
70 void DOAngularGridsSet( // WS Output:
71  Index& doit_za_grid_size,
72  Vector& aa_grid,
73  Vector& za_grid,
74  // Keywords:
75  const Index& N_za_grid,
76  const Index& N_aa_grid,
77  const String& za_grid_opt_file,
78  const Verbosity& verbosity) {
79  // Azimuth angle grid (the same is used for the scattering integral and
80  // for the radiative transfer.
81  if (N_aa_grid > 1)
82  nlinspace(aa_grid, 0, 360, N_aa_grid);
83  else if (N_aa_grid < 1) {
84  ostringstream os;
85  os << "N_aa_grid must be > 0 (even for 1D / DISORT cases).";
86  throw runtime_error(os.str());
87  } else {
88  aa_grid.resize(1);
89  aa_grid[0] = 0.;
90  }
91 
92  // Zenith angle grid:
93  // Number of zenith angle grid points (only for scattering integral):
94  if (N_za_grid < 0) {
95  ostringstream os;
96  os << "N_za_grid must be >= 0.";
97  throw runtime_error(os.str());
98  } else
99  doit_za_grid_size = N_za_grid;
100 
101  if (za_grid_opt_file == "")
102  if (N_za_grid == 0)
103  za_grid.resize(0);
104  else if (N_za_grid == 1) {
105  ostringstream os;
106  os << "N_za_grid must be >1 or =0 (the latter only allowed for RT4).";
107  throw runtime_error(os.str());
108  } else
109  nlinspace(za_grid, 0, 180, N_za_grid);
110  else
111  xml_read_from_file(za_grid_opt_file, za_grid, verbosity);
112 }
113 
114 /* Workspace method: Doxygen documentation will be auto-generated */
115 void doit_conv_flagAbs( //WS Input and Output:
116  Index& doit_conv_flag,
117  Index& doit_iteration_counter,
118  Tensor6& cloudbox_field_mono,
119  // WS Input:
120  const Tensor6& cloudbox_field_mono_old,
121  // Keyword:
122  const Vector& epsilon,
123  const Index& max_iterations,
124  const Index& throw_nonconv_error,
125  const Verbosity& verbosity) {
126  CREATE_OUT1;
127  CREATE_OUT2;
128 
129  //------------Check the input---------------------------------------
130  if (doit_conv_flag != 0)
131  throw runtime_error(
132  "Convergence flag is non-zero, which means that this\n"
133  "WSM is not used correctly. *doit_conv_flagAbs* should\n"
134  "be used only in *doit_conv_test_agenda*\n");
135 
136  const Index N_p = cloudbox_field_mono.nvitrines();
137  const Index N_lat = cloudbox_field_mono.nshelves();
138  const Index N_lon = cloudbox_field_mono.nbooks();
139  const Index N_za = cloudbox_field_mono.npages();
140  const Index N_aa = cloudbox_field_mono.nrows();
141  const Index stokes_dim = cloudbox_field_mono.ncols();
142 
143  // Check keyword "epsilon":
144  if (epsilon.nelem() != stokes_dim)
145  throw runtime_error(
146  "You have to specify limiting values for the "
147  "convergence test for each Stokes component "
148  "separately. That means that *epsilon* must "
149  "have *stokes_dim* elements!");
150 
151  // Check if cloudbox_field and cloudbox_field_old have the same dimensions:
152  if (!is_size(
153  cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim))
154  throw runtime_error(
155  "The fields (Tensor6) *cloudbox_field* and \n"
156  "*cloudbox_field_old* which are compared in the \n"
157  "convergence test do not have the same size.\n");
158 
159  //-----------End of checks-------------------------------------------------
160 
161  doit_iteration_counter += 1;
162  out2 << " Number of DOIT iteration: " << doit_iteration_counter << "\n";
163 
164  if (doit_iteration_counter > max_iterations) {
165  ostringstream out;
166  out << "Method does not converge (number of iterations \n"
167  << "is > " << max_iterations << "). Either the cloud "
168  << "particle number density \n"
169  << "is too large or the numerical setup for the DOIT \n"
170  << "calculation is not correct. In case of limb \n"
171  << "simulations please make sure that you use an \n"
172  << "optimized zenith angle grid. \n"
173  << "*cloudbox_field* might be wrong.\n";
174  if (throw_nonconv_error != 0) {
175  // FIXME: OLE: Remove this later
176  // ostringstream os;
177  // os << "Error in DOIT calculation:\n"
178  // << out.str();
179  // throw runtime_error( os.str() );
180  out1 << "Warning in DOIT calculation (output set to NaN):\n" << out.str();
181  cloudbox_field_mono = NAN;
182  doit_conv_flag = 1;
183  } else {
184  out1 << "Warning in DOIT calculation (output equals current status):\n"
185  << out.str();
186  doit_conv_flag = 1;
187  }
188  } else {
189  for (Index p_index = 0; p_index < N_p; p_index++) {
190  for (Index lat_index = 0; lat_index < N_lat; lat_index++) {
191  for (Index lon_index = 0; lon_index < N_lon; lon_index++) {
192  for (Index za_index = 0; za_index < N_za; za_index++) {
193  for (Index aa_index = 0; aa_index < N_aa; aa_index++) {
194  for (Index stokes_index = 0; stokes_index < stokes_dim;
195  stokes_index++) {
196  Numeric diff = (cloudbox_field_mono(p_index,
197  lat_index,
198  lon_index,
199  za_index,
200  aa_index,
201  stokes_index) -
202  cloudbox_field_mono_old(p_index,
203  lat_index,
204  lon_index,
205  za_index,
206  aa_index,
207  stokes_index));
208 
209  // If the absolute difference of the components
210  // is larger than the pre-defined values, return
211  // to *cloudbox_fieldIterarte* and do next iteration
212 
213  if (abs(diff) > epsilon[stokes_index]) {
214  out1 << "difference: " << diff << "\n";
215  return;
216  }
217 
218  } // End loop stokes_dom.
219  } // End loop aa_grid.
220  } // End loop za_grid.
221  } // End loop lon_grid.
222  } // End loop lat_grid.
223  } // End p_grid.
224 
225  // Convergence test has been successful, doit_conv_flag can be set to 1.
226  doit_conv_flag = 1;
227  }
228 }
229 
230 /* Workspace method: Doxygen documentation will be auto-generated */
231 void doit_conv_flagAbsBT( //WS Input and Output:
232  Index& doit_conv_flag,
233  Index& doit_iteration_counter,
234  Tensor6& cloudbox_field_mono,
235  // WS Input:
236  const Tensor6& cloudbox_field_mono_old,
237  const Vector& f_grid,
238  const Index& f_index,
239  // Keyword:
240  const Vector& epsilon,
241  const Index& max_iterations,
242  const Index& throw_nonconv_error,
243  const Verbosity& verbosity) {
244  CREATE_OUT1;
245  CREATE_OUT2;
246 
247  //------------Check the input---------------------------------------
248 
249  if (doit_conv_flag != 0)
250  throw runtime_error(
251  "Convergence flag is non-zero, which means that this \n"
252  "WSM is not used correctly. *doit_conv_flagAbs* should\n"
253  "be used only in *doit_conv_test_agenda*\n");
254 
255  const Index N_p = cloudbox_field_mono.nvitrines();
256  const Index N_lat = cloudbox_field_mono.nshelves();
257  const Index N_lon = cloudbox_field_mono.nbooks();
258  const Index N_za = cloudbox_field_mono.npages();
259  const Index N_aa = cloudbox_field_mono.nrows();
260  const Index stokes_dim = cloudbox_field_mono.ncols();
261 
262  // Check keyword "epsilon":
263  if (epsilon.nelem() != stokes_dim)
264  throw runtime_error(
265  "You have to specify limiting values for the "
266  "convergence test for each Stokes component "
267  "separately. That means that *epsilon* must "
268  "have *stokes_dim* elements!");
269 
270  // Check if cloudbox_field and cloudbox_field_old have the same dimensions:
271  if (!is_size(
272  cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim))
273  throw runtime_error(
274  "The fields (Tensor6) *cloudbox_field* and \n"
275  "*cloudbox_field_old* which are compared in the \n"
276  "convergence test do not have the same size.\n");
277 
278  // Frequency grid
279  //
280  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
281  chk_if_increasing("f_grid", f_grid);
282 
283  // Is the frequency index valid?
284  if (f_index >= f_grid.nelem())
285  throw runtime_error(
286  "*f_index* is greater than number of elements in the\n"
287  "frequency grid.\n");
288 
289  //-----------End of checks--------------------------------
290 
291  doit_iteration_counter += 1;
292  out2 << " Number of DOIT iteration: " << doit_iteration_counter << "\n";
293 
294  //Numeric max_diff_bt = 0.;
295  if (doit_iteration_counter > max_iterations) {
296  ostringstream out;
297  out << "At frequency " << f_grid[f_index] << " GHz \n"
298  << "method does not converge (number of iterations \n"
299  << "is > " << max_iterations << "). Either the particle"
300  << " number density \n"
301  << "is too large or the numerical setup for the DOIT \n"
302  << "calculation is not correct. In case of limb \n"
303  << "simulations please make sure that you use an \n"
304  << "optimized zenith angle grid. \n"
305  << "*cloudbox_field* might be wrong.\n";
306  if (throw_nonconv_error != 0) {
307  // FIXME: OLE: Remove this later
308  // ostringstream os;
309  // os << "Error in DOIT calculation:\n"
310  // << out.str();
311  // throw runtime_error( os.str() );
312  out1 << "Warning in DOIT calculation (output set to NaN):\n" << out.str();
313  cloudbox_field_mono = NAN;
314  doit_conv_flag = 1;
315  } else {
316  out1 << "Warning in DOIT calculation (output equals current status):\n"
317  << out.str();
318  doit_conv_flag = 1;
319  }
320  } else {
321  for (Index p_index = 0; p_index < N_p; p_index++) {
322  for (Index lat_index = 0; lat_index < N_lat; lat_index++) {
323  for (Index lon_index = 0; lon_index < N_lon; lon_index++) {
324  for (Index za_index = 0; za_index < N_za; za_index++) {
325  for (Index aa_index = 0; aa_index < N_aa; aa_index++) {
326  for (Index stokes_index = 0; stokes_index < stokes_dim;
327  stokes_index++) {
328  Numeric diff = cloudbox_field_mono(p_index,
329  lat_index,
330  lon_index,
331  za_index,
332  aa_index,
333  stokes_index) -
334  cloudbox_field_mono_old(p_index,
335  lat_index,
336  lon_index,
337  za_index,
338  aa_index,
339  stokes_index);
340 
341  // If the absolute difference of the components
342  // is larger than the pre-defined values, return
343  // to *cloudbox_fieldIterate* and do next iteration
344  Numeric diff_bt = invrayjean(diff, f_grid[f_index]);
345  // if( abs(diff_bt) > max_diff_bt )
346  // max_diff_bt = abs(diff_bt);
347  if (abs(diff_bt) > epsilon[stokes_index]) {
348  out1 << "BT difference: " << diff_bt << " in stokes dim "
349  << stokes_index << "\n";
350  // cout << "max BT difference in iteration #"
351  // << doit_iteration_counter << ": "
352  // << max_diff_bt << "\n";
353  return;
354  }
355  } // End loop stokes_dom.
356  } // End loop aa_grid.
357  } // End loop za_grid.
358  } // End loop lon_grid.
359  } // End loop lat_grid.
360  } // End p_grid.
361 
362  //cout << "max BT difference in iteration #" << doit_iteration_counter
363  // << ": " << max_diff_bt << "\n";
364  // Convergence test has been successful, doit_conv_flag can be set to 1.
365  doit_conv_flag = 1;
366  }
367 }
368 
369 /* Workspace method: Doxygen documentation will be auto-generated */
370 void doit_conv_flagLsq( //WS Output:
371  Index& doit_conv_flag,
372  Index& doit_iteration_counter,
373  Tensor6& cloudbox_field_mono,
374  // WS Input:
375  const Tensor6& cloudbox_field_mono_old,
376  const Vector& f_grid,
377  const Index& f_index,
378  // Keyword:
379  const Vector& epsilon,
380  const Index& max_iterations,
381  const Index& throw_nonconv_error,
382  const Verbosity& verbosity) {
383  CREATE_OUT1;
384  CREATE_OUT2;
385 
386  //------------Check the input---------------------------------------
387 
388  if (doit_conv_flag != 0)
389  throw runtime_error(
390  "Convergence flag is non-zero, which means that this \n"
391  "WSM is not used correctly. *doit_conv_flagAbs* should\n"
392  "be used only in *doit_conv_test_agenda*\n");
393 
394  const Index N_p = cloudbox_field_mono.nvitrines();
395  const Index N_lat = cloudbox_field_mono.nshelves();
396  const Index N_lon = cloudbox_field_mono.nbooks();
397  const Index N_za = cloudbox_field_mono.npages();
398  const Index N_aa = cloudbox_field_mono.nrows();
399  const Index stokes_dim = cloudbox_field_mono.ncols();
400 
401  // Check keyword "epsilon":
402  if (epsilon.nelem() != stokes_dim)
403  throw runtime_error(
404  "You have to specify limiting values for the "
405  "convergence test for each Stokes component "
406  "separately. That means that *epsilon* must "
407  "have *stokes_dim* elements!");
408 
409  // Check if cloudbox_field and cloudbox_field_old have the same dimensions:
410  if (!is_size(
411  cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim))
412  throw runtime_error(
413  "The fields (Tensor6) *cloudbox_field* and \n"
414  "*cloudbox_field_old* which are compared in the \n"
415  "convergence test do not have the same size.\n");
416 
417  // Frequency grid
418  //
419  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
420  chk_if_increasing("f_grid", f_grid);
421 
422  // Is the frequency index valid?
423  if (f_index >= f_grid.nelem())
424  throw runtime_error(
425  "*f_index* is greater than number of elements in the\n"
426  "frequency grid.\n");
427 
428  //-----------End of checks--------------------------------
429 
430  doit_iteration_counter += 1;
431  out2 << " Number of DOIT iteration: " << doit_iteration_counter << "\n";
432 
433  if (doit_iteration_counter > max_iterations) {
434  ostringstream out;
435  out << "Method does not converge (number of iterations \n"
436  << "is > " << max_iterations << "). Either the"
437  << " particle number density \n"
438  << "is too large or the numerical setup for the DOIT \n"
439  << "calculation is not correct. In case of limb \n"
440  << "simulations please make sure that you use an \n"
441  << "optimized zenith angle grid. \n";
442  if (throw_nonconv_error != 0) {
443  // FIXME: OLE: Remove this later
444  // ostringstream os;
445  // os << "Error in DOIT calculation:\n"
446  // << out.str();
447  // throw runtime_error( os.str() );
448  out1 << "Warning in DOIT calculation (output set to NaN):\n" << out.str();
449  cloudbox_field_mono = NAN;
450  doit_conv_flag = 1;
451  } else {
452  out1 << "Warning in DOIT calculation (output equals current status):\n"
453  << out.str();
454  doit_conv_flag = 1;
455  }
456  } else {
457  Vector lqs(4, 0.);
458 
459  // Will be set to zero if convergence not fullfilled
460  doit_conv_flag = 1;
461  for (Index i = 0; i < epsilon.nelem(); i++) {
462  for (Index p_index = 0; p_index < N_p; p_index++) {
463  for (Index lat_index = 0; lat_index < N_lat; lat_index++) {
464  for (Index lon_index = 0; lon_index < N_lon; lon_index++) {
465  for (Index za_index = 0; za_index < N_za; za_index++) {
466  for (Index aa_index = 0; aa_index < N_aa; aa_index++) {
467  lqs[i] += pow(
468  cloudbox_field_mono(
469  p_index, lat_index, lon_index, za_index, aa_index, i) -
470  cloudbox_field_mono_old(p_index,
471  lat_index,
472  lon_index,
473  za_index,
474  aa_index,
475  i),
476  2);
477  } // End loop aa_grid.
478  } // End loop za_grid.
479  } // End loop lon_grid.
480  } // End loop lat_grid.
481  } // End p_grid.
482 
483  lqs[i] = sqrt(lqs[i]);
484  lqs[i] /= (Numeric)(N_p * N_lat * N_lon * N_za * N_aa);
485 
486  // Convert difference to Rayleigh Jeans BT
487  lqs[i] = invrayjean(lqs[i], f_grid[f_index]);
488 
489  if (lqs[i] >= epsilon[i]) doit_conv_flag = 0;
490  }
491  // end loop stokes_index
492  out1 << "lqs [I]: " << lqs[0] << "\n";
493  }
494 }
495 
496 /* Workspace method: Doxygen documentation will be auto-generated */
497 void cloudbox_field_monoIterate(Workspace& ws,
498  // WS Input and Output:
499  Tensor6& cloudbox_field_mono,
500 
501  // WS Input:
502  const Agenda& doit_scat_field_agenda,
503  const Agenda& doit_rte_agenda,
504  const Agenda& doit_conv_test_agenda,
505  const Index& accelerated,
506  const Verbosity& verbosity)
507 
508 {
509  CREATE_OUT2;
510 
511  //---------------Check input---------------------------------
512  chk_not_empty("doit_scat_field_agenda", doit_scat_field_agenda);
513  chk_not_empty("doit_rte_agenda", doit_rte_agenda);
514  chk_not_empty("doit_conv_test_agenda", doit_conv_test_agenda);
515 
516  for (Index v = 0; v < cloudbox_field_mono.nvitrines(); v++)
517  for (Index s = 0; s < cloudbox_field_mono.nshelves(); s++)
518  for (Index b = 0; b < cloudbox_field_mono.nbooks(); b++)
519  for (Index p = 0; p < cloudbox_field_mono.npages(); p++)
520  for (Index r = 0; r < cloudbox_field_mono.nrows(); r++)
521  for (Index c = 0; c < cloudbox_field_mono.ncols(); c++)
522  if (std::isnan(cloudbox_field_mono(v, s, b, p, r, c)))
523  throw std::runtime_error(
524  "*cloudbox_field_mono* contains at least one NaN value.\n"
525  "This indicates an improper initialization of *cloudbox_field*.");
526 
527  //cloudbox_field_mono can not be further checked here, because there is no way
528  //to find out the size without including a lot more interface
529  //variables
530  //-----------End of checks--------------------------------------
531 
532  Tensor6 cloudbox_field_mono_old_local;
533  Index doit_conv_flag_local;
534  Index doit_iteration_counter_local;
535 
536  // Resize and initialize doit_scat_field,
537  // which has the same dimensions as cloudbox_field
538  Tensor6 doit_scat_field_local(cloudbox_field_mono.nvitrines(),
539  cloudbox_field_mono.nshelves(),
540  cloudbox_field_mono.nbooks(),
541  cloudbox_field_mono.npages(),
542  cloudbox_field_mono.nrows(),
543  cloudbox_field_mono.ncols(),
544  0.);
545 
546  doit_conv_flag_local = 0;
547  doit_iteration_counter_local = 0;
548  // Array to save the last iteration steps
549  ArrayOfTensor6 acceleration_input;
550  if (accelerated) {
551  acceleration_input.resize(4);
552  }
553  while (doit_conv_flag_local == 0) {
554  // 1. Copy cloudbox_field to cloudbox_field_old.
555  cloudbox_field_mono_old_local = cloudbox_field_mono;
556 
557  // 2.Calculate scattered field vector for all points in the cloudbox.
558 
559  // Calculate the scattered field.
560  out2 << " Execute doit_scat_field_agenda. \n";
561  doit_scat_field_agendaExecute(
562  ws, doit_scat_field_local, cloudbox_field_mono, doit_scat_field_agenda);
563 
564  // Update cloudbox_field.
565  out2 << " Execute doit_rte_agenda. \n";
566  doit_rte_agendaExecute(
567  ws, cloudbox_field_mono, doit_scat_field_local, doit_rte_agenda);
568 
569  //Convergence test.
570  doit_conv_test_agendaExecute(ws,
571  doit_conv_flag_local,
572  doit_iteration_counter_local,
573  cloudbox_field_mono,
574  cloudbox_field_mono_old_local,
575  doit_conv_test_agenda);
576 
577  // Convergence Acceleration, if wished.
578  if (accelerated > 0 && doit_conv_flag_local == 0) {
579  acceleration_input[(doit_iteration_counter_local - 1) % 4] =
580  cloudbox_field_mono;
581  // NG - Acceleration
582  if (doit_iteration_counter_local % 4 == 0) {
583  cloudbox_field_ngAcceleration(
584  cloudbox_field_mono, acceleration_input, accelerated, verbosity);
585  }
586  }
587  } //end of while loop, convergence is reached.
588 }
589 
590 /* Workspace method: Doxygen documentation will be auto-generated */
591 void cloudbox_fieldUpdate1D(
592  Workspace& ws,
593  // WS Input and Output:
594  Tensor6& cloudbox_field_mono,
595  // WS Input:
596  const Tensor6& doit_scat_field,
597  const ArrayOfIndex& cloudbox_limits,
598  // Calculate scalar gas absorption:
599  const Agenda& propmat_clearsky_agenda,
600  const Tensor4& vmr_field,
601  // Optical properties for individual scattering elements:
602  const Agenda& spt_calc_agenda,
603  const Vector& za_grid,
604  const Tensor4& pnd_field,
605  // Propagation path calculation:
606  const Agenda& ppath_step_agenda,
607  const Numeric& ppath_lmax,
608  const Numeric& ppath_lraytrace,
609  const Vector& p_grid,
610  const Tensor3& z_field,
611  const Vector& refellipsoid,
612  // Calculate thermal emission:
613  const Tensor3& t_field,
614  const Vector& f_grid,
615  const Index& f_index,
616  const Agenda& surface_rtprop_agenda,
617  const Index& doit_za_interp,
618  const Verbosity& verbosity) {
619  CREATE_OUT2;
620  CREATE_OUT3;
621 
622  out2
623  << " cloudbox_fieldUpdate1D: Radiative transfer calculation in cloudbox\n";
624  out2 << " ------------------------------------------------------------- \n";
625 
626  // ---------- Check the input ----------------------------------------
627 
628  // Agendas
629  chk_not_empty("spt_calc_agenda", spt_calc_agenda);
630  chk_not_empty("ppath_step_agenda", ppath_step_agenda);
631 
632  if (cloudbox_limits.nelem() != 2)
633  throw runtime_error(
634  "The cloudbox dimension is not 1D! \n"
635  "Do you really want to do a 1D calculation? \n"
636  "If not, use *cloudbox_fieldUpdateSeq3D*.\n");
637 
638  // Number of zenith angles.
639  const Index N_scat_za = za_grid.nelem();
640 
641  if (za_grid[0] != 0. || za_grid[N_scat_za - 1] != 180.)
642  throw runtime_error("The range of *za_grid* must [0 180].");
643 
644  if (p_grid.nelem() < 2)
645  throw runtime_error("The length of *p_grid* must be >= 2.");
646  chk_if_decreasing("p_grid", p_grid);
647 
648  chk_size("z_field", z_field, p_grid.nelem(), 1, 1);
649  chk_size("t_field", t_field, p_grid.nelem(), 1, 1);
650 
651  // Frequency grid
652  //
653  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
654  chk_if_increasing("f_grid", f_grid);
655 
656  // Is the frequency index valid?
657  if (f_index >= f_grid.nelem())
658  throw runtime_error(
659  "*f_index* is greater than number of elements in the\n"
660  "frequency grid.\n");
661 
662  if (!(doit_za_interp == 0 || doit_za_interp == 1))
663  throw runtime_error(
664  "Interpolation method is not defined. Use \n"
665  "*doit_za_interpSet*.\n");
666 
667  const Index stokes_dim = doit_scat_field.ncols();
668  assert(stokes_dim > 0 || stokes_dim < 5);
669 
670  // These variables are calculated internally, so assertions should be o.k.
671  assert(is_size(cloudbox_field_mono,
672  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
673  1,
674  1,
675  N_scat_za,
676  1,
677  stokes_dim));
678 
679  assert(is_size(doit_scat_field,
680  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
681  1,
682  1,
683  N_scat_za,
684  1,
685  stokes_dim));
686 
687  // FIXME: Check *vmr_field*
688 
689  // -------------- End of checks --------------------------------------
690 
691  //=======================================================================
692  // Calculate scattering coefficients for all positions in the cloudbox
693  //=======================================================================
694  out3 << "Calculate optical properties of individual scattering elements\n";
695 
696  // At this place only the particle properties are calculated. Gaseous
697  // absorption is calculated inside the radiative transfer part. Inter-
698  // polating absorption coefficients for gaseous species gives very bad
699  // results, so they are calulated for interpolated VMRs,
700  // temperature and pressure.
701 
702  // To use special interpolation functions for atmospheric fields we
703  // use ext_mat_field and abs_vec_field:
704  Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
705  1,
706  1,
707  stokes_dim,
708  stokes_dim,
709  0.);
710  Tensor4 abs_vec_field(
711  cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
712 
713  Tensor6 cloudbox_field_old(cloudbox_field_mono);
714 
715  //Only dummy variable:
716  Index aa_index_local = 0;
717 
718  //Loop over all directions, defined by za_grid
719  for (Index za_index_local = 0; za_index_local < N_scat_za; za_index_local++) {
720  // This function has to be called inside the angular loop, as
721  // spt_calc_agenda takes *za_index_local* and *aa_index*
722  // from the workspace.
723  cloud_fieldsCalc(ws,
724  ext_mat_field,
725  abs_vec_field,
726  spt_calc_agenda,
727  za_index_local,
728  aa_index_local,
729  cloudbox_limits,
730  t_field,
731  pnd_field,
732  verbosity);
733 
734  //======================================================================
735  // Radiative transfer inside the cloudbox
736  //=====================================================================
737 
738  for (Index p_index = cloudbox_limits[0]; p_index <= cloudbox_limits[1];
739  p_index++) {
740  if ((p_index != 0) || (za_grid[za_index_local] <= 90.)) {
741  cloud_ppath_update1D_noseq(ws,
742  cloudbox_field_mono,
743  p_index,
744  za_index_local,
745  za_grid,
746  cloudbox_limits,
747  cloudbox_field_old,
748  doit_scat_field,
749  propmat_clearsky_agenda,
750  vmr_field,
751  ppath_step_agenda,
752  ppath_lmax,
753  ppath_lraytrace,
754  p_grid,
755  z_field,
756  refellipsoid,
757  t_field,
758  f_grid,
759  f_index,
760  ext_mat_field,
761  abs_vec_field,
762  surface_rtprop_agenda,
763  doit_za_interp,
764  verbosity);
765  }
766  }
767  } // Closes loop over za_grid.
768 }
769 
770 /* Workspace method: Doxygen documentation will be auto-generated */
771 void cloudbox_fieldUpdateSeq1D(
772  Workspace& ws,
773  // WS Input and Output:
774  Tensor6& cloudbox_field_mono,
775  Tensor6& doit_scat_field,
776  // WS Input:
777  const ArrayOfIndex& cloudbox_limits,
778  // Calculate scalar gas absorption:
779  const Agenda& propmat_clearsky_agenda,
780  const Tensor4& vmr_field,
781  // Optical properties for individual scattering elements:
782  const Agenda& spt_calc_agenda,
783  const Vector& za_grid,
784  const Vector& aa_grid,
785  const Tensor4& pnd_field,
786  // Propagation path calculation:
787  const Agenda& ppath_step_agenda,
788  const Numeric& ppath_lmax,
789  const Numeric& ppath_lraytrace,
790  const Vector& p_grid,
791  const Tensor3& z_field,
792  const Vector& refellipsoid,
793  // Calculate thermal emission:
794  const Tensor3& t_field,
795  const Vector& f_grid,
796  const Index& f_index,
797  const Agenda& surface_rtprop_agenda, //STR
798  const Index& doit_za_interp,
799  const Index& normalize,
800  const Numeric& norm_error_threshold,
801  const Index& norm_debug,
802  const Verbosity& verbosity) {
803  CREATE_OUT2;
804  CREATE_OUT3;
805 
806  out2
807  << " cloudbox_fieldUpdateSeq1D: Radiative transfer calculation in cloudbox\n";
808  out2 << " ------------------------------------------------------------- \n";
809 
810  // ---------- Check the input ----------------------------------------
811 
812  // Agendas
813  chk_not_empty("propmat_clearsky_agenda", propmat_clearsky_agenda);
814  chk_not_empty("spt_calc_agenda", spt_calc_agenda);
815  chk_not_empty("ppath_step_agenda", ppath_step_agenda);
816 
817  if (cloudbox_limits.nelem() != 2)
818  throw runtime_error(
819  "The cloudbox dimension is not 1D! \n"
820  "Do you really want to do a 1D calculation? \n"
821  "For 3D use *cloudbox_fieldUpdateSeq3D*.\n");
822 
823  // Number of zenith angles.
824  const Index N_scat_za = za_grid.nelem();
825 
826  if (za_grid[0] != 0. || za_grid[N_scat_za - 1] != 180.)
827  throw runtime_error("The range of *za_grid* must [0 180].");
828 
829  if (p_grid.nelem() < 2)
830  throw runtime_error("The length of *p_grid* must be >= 2.");
831  chk_if_decreasing("p_grid", p_grid);
832 
833  chk_size("z_field", z_field, p_grid.nelem(), 1, 1);
834  chk_size("t_field", t_field, p_grid.nelem(), 1, 1);
835 
836  // Frequency grid
837  //
838  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
839  chk_if_increasing("f_grid", f_grid);
840 
841  // Is the frequency index valid?
842  if (f_index >= f_grid.nelem())
843  throw runtime_error(
844  "*f_index* is greater than number of elements in the\n"
845  "frequency grid.\n");
846 
847  if (!(doit_za_interp == 0 || doit_za_interp == 1))
848  throw runtime_error(
849  "Interpolation method is not defined. Use \n"
850  "*doit_za_interpSet*.\n");
851 
852  const Index stokes_dim = doit_scat_field.ncols();
853  assert(stokes_dim > 0 || stokes_dim < 5);
854 
855  // These variables are calculated internally, so assertions should be o.k.
856  assert(is_size(cloudbox_field_mono,
857  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
858  1,
859  1,
860  N_scat_za,
861  1,
862  stokes_dim));
863 
864  assert(is_size(doit_scat_field,
865  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
866  1,
867  1,
868  N_scat_za,
869  1,
870  stokes_dim));
871 
872  // FIXME: Check *vmr_field*
873 
874  // -------------- End of checks --------------------------------------
875 
876  //=======================================================================
877  // Calculate scattering coefficients for all positions in the cloudbox
878  //=======================================================================
879  out3 << "Calculate optical properties of individual scattering elements\n";
880 
881  // At this place only the particle properties are calculated. Gaseous
882  // absorption is calculated inside the radiative transfer part. Inter-
883  // polating absorption coefficients for gaseous species gives very bad
884  // results, so they are calulated for interpolated VMRs,
885  // temperature and pressure.
886 
887  // To use special interpolation functions for atmospheric fields we
888  // use ext_mat_field and abs_vec_field:
889  Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
890  1,
891  1,
892  stokes_dim,
893  stokes_dim,
894  0.);
895  Tensor4 abs_vec_field(
896  cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
897 
898  // If theta is between 90° and the limiting value, the intersection point
899  // is exactly at the same level as the starting point (cp. AUG)
900  Numeric theta_lim =
901  180. - asin((refellipsoid[0] + z_field(cloudbox_limits[0], 0, 0)) /
902  (refellipsoid[0] + z_field(cloudbox_limits[1], 0, 0))) *
903  RAD2DEG;
904 
905  // Epsilon for additional limb iterations
906  Vector epsilon(4);
907  epsilon[0] = 0.1;
908  epsilon[1] = 0.01;
909  epsilon[2] = 0.01;
910  epsilon[3] = 0.01;
911 
912  Matrix cloudbox_field_limb;
913 
914  //Only dummy variables:
915  Index aa_index_local = 0;
916 
917  if (normalize) {
918  Tensor4 si, sei, si_corr;
919  doit_scat_fieldNormalize(ws,
920  doit_scat_field,
921  cloudbox_field_mono,
922  cloudbox_limits,
923  spt_calc_agenda,
924  1,
925  za_grid,
926  aa_grid,
927  pnd_field,
928  t_field,
929  norm_error_threshold,
930  norm_debug,
931  verbosity);
932  }
933 
934  //Loop over all directions, defined by za_grid
935  for (Index za_index_local = 0; za_index_local < N_scat_za; za_index_local++) {
936  // This function has to be called inside the angular loop, as
937  // spt_calc_agenda takes *za_index* and *aa_index*
938  // from the workspace.
939  cloud_fieldsCalc(ws,
940  ext_mat_field,
941  abs_vec_field,
942  spt_calc_agenda,
943  za_index_local,
944  aa_index_local,
945  cloudbox_limits,
946  t_field,
947  pnd_field,
948  verbosity);
949 
950  //======================================================================
951  // Radiative transfer inside the cloudbox
952  //=====================================================================
953 
954  // Sequential update for uplooking angles
955  if (za_grid[za_index_local] <= 90.) {
956  // Loop over all positions inside the cloud box defined by the
957  // cloudbox_limits excluding the upper boundary. For uplooking
958  // directions, we start from cloudbox_limits[1]-1 and go down
959  // to cloudbox_limits[0] to do a sequential update of the
960  // radiation field
961  for (Index p_index = cloudbox_limits[1] - 1;
962  p_index >= cloudbox_limits[0];
963  p_index--) {
964  cloud_ppath_update1D(ws,
965  cloudbox_field_mono,
966  p_index,
967  za_index_local,
968  za_grid,
969  cloudbox_limits,
970  doit_scat_field,
971  propmat_clearsky_agenda,
972  vmr_field,
973  ppath_step_agenda,
974  ppath_lmax,
975  ppath_lraytrace,
976  p_grid,
977  z_field,
978  refellipsoid,
979  t_field,
980  f_grid,
981  f_index,
982  ext_mat_field,
983  abs_vec_field,
984  surface_rtprop_agenda,
985  doit_za_interp,
986  verbosity);
987  }
988  } else if (za_grid[za_index_local] >= theta_lim) {
989  //
990  // Sequential updating for downlooking angles
991  //
992  for (Index p_index = cloudbox_limits[0] + 1;
993  p_index <= cloudbox_limits[1];
994  p_index++) {
995  cloud_ppath_update1D(ws,
996  cloudbox_field_mono,
997  p_index,
998  za_index_local,
999  za_grid,
1000  cloudbox_limits,
1001  doit_scat_field,
1002  propmat_clearsky_agenda,
1003  vmr_field,
1004  ppath_step_agenda,
1005  ppath_lmax,
1006  ppath_lraytrace,
1007  p_grid,
1008  z_field,
1009  refellipsoid,
1010  t_field,
1011  f_grid,
1012  f_index,
1013  ext_mat_field,
1014  abs_vec_field,
1015  surface_rtprop_agenda,
1016  doit_za_interp,
1017  verbosity);
1018  } // Close loop over p_grid (inside cloudbox).
1019  } // end if downlooking.
1020 
1021  //
1022  // Limb looking:
1023  // We have to include a special case here, as we may miss the endpoints
1024  // when the intersection point is at the same level as the aactual point.
1025  // To be save we loop over the full cloudbox. Inside the function
1026  // cloud_ppath_update1D it is checked whether the intersection point is
1027  // inside the cloudbox or not.
1028  else {
1029  bool conv_flag = false;
1030  Index limb_it = 0;
1031  while (!conv_flag && limb_it < 10) {
1032  limb_it++;
1033  cloudbox_field_limb =
1034  cloudbox_field_mono(joker, 0, 0, za_index_local, 0, joker);
1035  for (Index p_index = cloudbox_limits[0]; p_index <= cloudbox_limits[1];
1036  p_index++) {
1037  // For this case the cloudbox goes down to the surface and we
1038  // look downwards. These cases are outside the cloudbox and
1039  // not needed. Switch is included here, as ppath_step_agenda
1040  // gives an error for such cases.
1041  if (p_index != 0) {
1042  cloud_ppath_update1D(ws,
1043  cloudbox_field_mono,
1044  p_index,
1045  za_index_local,
1046  za_grid,
1047  cloudbox_limits,
1048  doit_scat_field,
1049  propmat_clearsky_agenda,
1050  vmr_field,
1051  ppath_step_agenda,
1052  ppath_lmax,
1053  ppath_lraytrace,
1054  p_grid,
1055  z_field,
1056  refellipsoid,
1057  t_field,
1058  f_grid,
1059  f_index,
1060  ext_mat_field,
1061  abs_vec_field,
1062  surface_rtprop_agenda,
1063  doit_za_interp,
1064  verbosity);
1065  }
1066  }
1067 
1068  conv_flag = true;
1069  for (Index p_index = 0;
1070  conv_flag && p_index < cloudbox_field_mono.nvitrines();
1071  p_index++) {
1072  for (Index stokes_index = 0; conv_flag && stokes_index < stokes_dim;
1073  stokes_index++) {
1074  Numeric diff = cloudbox_field_mono(
1075  p_index, 0, 0, za_index_local, 0, stokes_index) -
1076  cloudbox_field_limb(p_index, stokes_index);
1077 
1078  // If the absolute difference of the components
1079  // is larger than the pre-defined values, continue with
1080  // another iteration
1081  Numeric diff_bt = invrayjean(diff, f_grid[f_index]);
1082  if (abs(diff_bt) > epsilon[stokes_index]) {
1083  out2 << "Limb BT difference: " << diff_bt << " in stokes dim "
1084  << stokes_index << "\n";
1085  conv_flag = false;
1086  }
1087  }
1088  }
1089  }
1090  out2 << "Limb iterations: " << limb_it << "\n";
1091  }
1092  } // Closes loop over za_grid.
1093 } // End of the function.
1094 
1095 /* Workspace method: Doxygen documentation will be auto-generated */
1096 void cloudbox_fieldUpdateSeq3D(
1097  Workspace& ws,
1098  // WS Output and Input:
1099  Tensor6& cloudbox_field_mono,
1100  // WS Input:
1101  const Tensor6& doit_scat_field,
1102  const ArrayOfIndex& cloudbox_limits,
1103  // Calculate scalar gas absorption:
1104  const Agenda& propmat_clearsky_agenda,
1105  const Tensor4& vmr_field,
1106  // Optical properties for individual scattering elements:
1107  const Agenda& spt_calc_agenda,
1108  const Vector& za_grid,
1109  const Vector& aa_grid,
1110  const Tensor4& pnd_field,
1111  // Propagation path calculation:
1112  const Agenda& ppath_step_agenda,
1113  const Numeric& ppath_lmax,
1114  const Numeric& ppath_lraytrace,
1115  const Vector& p_grid,
1116  const Vector& lat_grid,
1117  const Vector& lon_grid,
1118  const Tensor3& z_field,
1119  const Vector& refellipsoid,
1120  // Calculate thermal emission:
1121  const Tensor3& t_field,
1122  const Vector& f_grid,
1123  const Index& f_index,
1124  const Index& doit_za_interp,
1125  const Verbosity& verbosity) {
1126  CREATE_OUT2;
1127  CREATE_OUT3;
1128 
1129  out2
1130  << " cloudbox_fieldUpdateSeq3D: Radiative transfer calculatiuon in cloudbox.\n";
1131  out2 << " ------------------------------------------------------------- \n";
1132 
1133  // ---------- Check the input ----------------------------------------
1134 
1135  // Agendas
1136  chk_not_empty("propmat_clearsky_agenda", propmat_clearsky_agenda);
1137  chk_not_empty("spt_calc_agenda", spt_calc_agenda);
1138  chk_not_empty("ppath_step_agenda", ppath_step_agenda);
1139 
1140  if (cloudbox_limits.nelem() != 6)
1141  throw runtime_error(
1142  "The cloudbox dimension is not 3D! \n"
1143  "Do you really want to do a 3D calculation? \n"
1144  "For 1D use *cloudbox_fieldUpdateSeq1D*.\n");
1145 
1146  // Number of zenith angles.
1147  const Index N_scat_za = za_grid.nelem();
1148 
1149  if (za_grid[0] != 0. || za_grid[N_scat_za - 1] != 180.)
1150  throw runtime_error("The range of *za_grid* must [0 180].");
1151 
1152  // Number of azimuth angles.
1153  const Index N_scat_aa = aa_grid.nelem();
1154 
1155  if (aa_grid[0] != 0. || aa_grid[N_scat_aa - 1] != 360.)
1156  throw runtime_error("The range of *za_grid* must [0 360].");
1157 
1158  // Check atmospheric grids
1159  chk_atm_grids(3, p_grid, lat_grid, lon_grid);
1160 
1161  // Check atmospheric fields
1162  chk_size(
1163  "z_field", z_field, p_grid.nelem(), lat_grid.nelem(), lon_grid.nelem());
1164  chk_size(
1165  "t_field", t_field, p_grid.nelem(), lat_grid.nelem(), lon_grid.nelem());
1166 
1167  // Frequency grid
1168  //
1169  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
1170  chk_if_increasing("f_grid", f_grid);
1171 
1172  // Is the frequency index valid?
1173  if (f_index >= f_grid.nelem())
1174  throw runtime_error(
1175  "*f_index* is greater than number of elements in the\n"
1176  "frequency grid.\n");
1177 
1178  if (!(doit_za_interp == 0 || doit_za_interp == 1))
1179  throw runtime_error(
1180  "Interpolation method is not defined. Use \n"
1181  "*doit_za_interpSet*.\n");
1182 
1183  const Index stokes_dim = doit_scat_field.ncols();
1184  assert(stokes_dim > 0 || stokes_dim < 5);
1185 
1186  // These variables are calculated internally, so assertions should be o.k.
1187  assert(is_size(cloudbox_field_mono,
1188  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1189  (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1190  (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1191  N_scat_za,
1192  N_scat_aa,
1193  stokes_dim));
1194 
1195  assert(is_size(doit_scat_field,
1196  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1197  (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1198  (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1199  N_scat_za,
1200  N_scat_aa,
1201  stokes_dim));
1202 
1203  // FIXME: Check *vmr_field*
1204 
1205  // ---------- End of checks ------------------------------------------
1206 
1207  //=======================================================================
1208  // Calculate coefficients for all positions in the cloudbox
1209  //=======================================================================
1210  out3 << "Calculate optical properties of individual scattering elements\n";
1211 
1212  // At this place only the particle properties are calculated. Gaseous
1213  // absorption is calculated inside the radiative transfer part. Inter-
1214  // polating absorption coefficients for gaseous species gives very bad
1215  // results, so they are
1216  // calulated for interpolated VMRs, temperature and pressure.
1217 
1218  // Define shorter names for cloudbox_limits.
1219 
1220  const Index p_low = cloudbox_limits[0];
1221  const Index p_up = cloudbox_limits[1];
1222  const Index lat_low = cloudbox_limits[2];
1223  const Index lat_up = cloudbox_limits[3];
1224  const Index lon_low = cloudbox_limits[4];
1225  const Index lon_up = cloudbox_limits[5];
1226 
1227  // To use special interpolation functions for atmospheric fields we
1228  // use ext_mat_field and abs_vec_field:
1229  Tensor5 ext_mat_field(p_up - p_low + 1,
1230  lat_up - lat_low + 1,
1231  lon_up - lon_low + 1,
1232  stokes_dim,
1233  stokes_dim,
1234  0.);
1235  Tensor4 abs_vec_field(p_up - p_low + 1,
1236  lat_up - lat_low + 1,
1237  lon_up - lon_low + 1,
1238  stokes_dim,
1239  0.);
1240 
1241  //Loop over all directions, defined by za_grid
1242  for (Index za_index = 0; za_index < N_scat_za; za_index++) {
1243  //Loop over azimuth directions (aa_grid). First and last point in
1244  // azimuth angle grid are euqal. Start with second element.
1245  for (Index aa_index = 1; aa_index < N_scat_aa; aa_index++) {
1246  //==================================================================
1247  // Radiative transfer inside the cloudbox
1248  //==================================================================
1249 
1250  // This function has to be called inside the angular loop, as
1251  // it spt_calc_agenda takes *za_index* and *aa_index*
1252  // from the workspace.
1253  cloud_fieldsCalc(ws,
1254  ext_mat_field,
1255  abs_vec_field,
1256  spt_calc_agenda,
1257  za_index,
1258  aa_index,
1259  cloudbox_limits,
1260  t_field,
1261  pnd_field,
1262  verbosity);
1263 
1264  Vector stokes_vec(stokes_dim, 0.);
1265 
1266  Numeric theta_lim = 180. - asin((refellipsoid[0] + z_field(p_low, 0, 0)) /
1267  (refellipsoid[0] + z_field(p_up, 0, 0))) *
1268  RAD2DEG;
1269 
1270  // Sequential update for uplooking angles
1271  if (za_grid[za_index] <= 90.) {
1272  // Loop over all positions inside the cloud box defined by the
1273  // cloudbox_limits exculding the upper boundary. For uplooking
1274  // directions, we start from cloudbox_limits[1]-1 and go down
1275  // to cloudbox_limits[0] to do a sequential update of the
1276  // aradiation field
1277  for (Index p_index = p_up - 1; p_index >= p_low; p_index--) {
1278  for (Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1279  for (Index lon_index = lon_low; lon_index <= lon_up; lon_index++) {
1280  cloud_ppath_update3D(ws,
1281  cloudbox_field_mono,
1282  p_index,
1283  lat_index,
1284  lon_index,
1285  za_index,
1286  aa_index,
1287  za_grid,
1288  aa_grid,
1289  cloudbox_limits,
1290  doit_scat_field,
1291  propmat_clearsky_agenda,
1292  vmr_field,
1293  ppath_step_agenda,
1294  ppath_lmax,
1295  ppath_lraytrace,
1296  p_grid,
1297  lat_grid,
1298  lon_grid,
1299  z_field,
1300  refellipsoid,
1301  t_field,
1302  f_grid,
1303  f_index,
1304  ext_mat_field,
1305  abs_vec_field,
1306  doit_za_interp,
1307  verbosity);
1308  }
1309  }
1310  }
1311  } // close up-looking case
1312  else if (za_grid[za_index] > theta_lim) {
1313  //
1314  // Sequential updating for downlooking angles
1315  //
1316  for (Index p_index = p_low + 1; p_index <= p_up; p_index++) {
1317  for (Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1318  for (Index lon_index = lon_low; lon_index <= lon_up; lon_index++) {
1319  cloud_ppath_update3D(ws,
1320  cloudbox_field_mono,
1321  p_index,
1322  lat_index,
1323  lon_index,
1324  za_index,
1325  aa_index,
1326  za_grid,
1327  aa_grid,
1328  cloudbox_limits,
1329  doit_scat_field,
1330  propmat_clearsky_agenda,
1331  vmr_field,
1332  ppath_step_agenda,
1333  ppath_lmax,
1334  ppath_lraytrace,
1335  p_grid,
1336  lat_grid,
1337  lon_grid,
1338  z_field,
1339  refellipsoid,
1340  t_field,
1341  f_grid,
1342  f_index,
1343  ext_mat_field,
1344  abs_vec_field,
1345  doit_za_interp,
1346  verbosity);
1347  }
1348  }
1349  }
1350  } // end if downlooking.
1351 
1352  //
1353  // Limb looking:
1354  // We have to include a special case here, as we may miss the endpoints
1355  // when the intersection point is at the same level as the actual point.
1356  // To be save we loop over the full cloudbox. Inside the function
1357  // cloud_ppath_update3D it is checked whether the intersection point is
1358  // inside the cloudbox or not.
1359  else if (za_grid[za_index] > 90. && za_grid[za_index] < theta_lim) {
1360  for (Index p_index = p_low; p_index <= p_up; p_index++) {
1361  // For this case the cloudbox goes down to the surface an we
1362  // look downwards. These cases are outside the cloudbox and
1363  // not needed. Switch is included here, as ppath_step_agenda
1364  // gives an error for such cases.
1365  if (!(p_index == 0 && za_grid[za_index] > 90.)) {
1366  for (Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1367  for (Index lon_index = lon_low; lon_index <= lon_up;
1368  lon_index++) {
1369  cloud_ppath_update3D(ws,
1370  cloudbox_field_mono,
1371  p_index,
1372  lat_index,
1373  lon_index,
1374  za_index,
1375  aa_index,
1376  za_grid,
1377  aa_grid,
1378  cloudbox_limits,
1379  doit_scat_field,
1380  propmat_clearsky_agenda,
1381  vmr_field,
1382  ppath_step_agenda,
1383  ppath_lmax,
1384  ppath_lraytrace,
1385  p_grid,
1386  lat_grid,
1387  lon_grid,
1388  z_field,
1389  refellipsoid,
1390  t_field,
1391  f_grid,
1392  f_index,
1393  ext_mat_field,
1394  abs_vec_field,
1395  doit_za_interp,
1396  verbosity);
1397  }
1398  }
1399  }
1400  }
1401  }
1402  } // Closes loop over aa_grid.
1403  } // Closes loop over za_grid.
1404 
1405  cloudbox_field_mono(joker, joker, joker, joker, 0, joker) =
1406  cloudbox_field_mono(joker, joker, joker, joker, N_scat_aa - 1, joker);
1407 }
1408 
1409 /* Workspace method: Doxygen documentation will be auto-generated */
1410 void cloudbox_fieldUpdateSeq1DPP(
1411  Workspace& ws,
1412  // WS Output:
1413  Tensor6& cloudbox_field_mono,
1414  // spt_calc_agenda:
1415  Index& za_index,
1416  // WS Input:
1417  const Tensor6& doit_scat_field,
1418  const ArrayOfIndex& cloudbox_limits,
1419  // Calculate scalar gas absorption:
1420  const Agenda& propmat_clearsky_agenda,
1421  const Tensor4& vmr_field,
1422  // Optical properties for individual scattering elements:
1423  const Agenda& spt_calc_agenda,
1424  const Vector& za_grid,
1425  const Tensor4& pnd_field,
1426  // Propagation path calculation:
1427  const Vector& p_grid,
1428  const Tensor3& z_field,
1429  // Calculate thermal emission:
1430  const Tensor3& t_field,
1431  const Vector& f_grid,
1432  const Index& f_index,
1433  const Verbosity& verbosity) {
1434  CREATE_OUT2;
1435  CREATE_OUT3;
1436 
1437  out2
1438  << " cloudbox_fieldUpdateSeq1DPP: Radiative transfer calculation in cloudbox.\n";
1439  out2
1440  << " --------------------------------------------------------------------- \n";
1441 
1442  const Index stokes_dim = doit_scat_field.ncols();
1443  // const Index atmosphere_dim = 1;
1444 
1445  //Check the input
1446 
1447  if (stokes_dim < 0 || stokes_dim > 4)
1448  throw runtime_error(
1449  "The dimension of stokes vector must be"
1450  "1,2,3, or 4");
1451 
1452  assert(is_size(cloudbox_field_mono,
1453  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1454  1,
1455  1,
1456  za_grid.nelem(),
1457  1,
1458  stokes_dim));
1459 
1460  assert(is_size(doit_scat_field,
1461  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1462  1,
1463  1,
1464  za_grid.nelem(),
1465  1,
1466  stokes_dim));
1467 
1468  // Is the frequency index valid?
1469  assert(f_index <= f_grid.nelem());
1470 
1471  // End of checks
1472 
1473  // Number of zenith angles.
1474  const Index N_scat_za = za_grid.nelem();
1475 
1476  //=======================================================================
1477  // Calculate scattering coefficients for all positions in the cloudbox
1478  //=======================================================================
1479  out3 << "Calculate optical properties of individual scattering elements\n";
1480 
1481  // At this place only the particle properties are calculated. Gaseous
1482  // absorption is calculated inside the radiative transfer part. Inter-
1483  // polating absorption coefficients for gaseous species gives very bad
1484  // results, so they are
1485  // calulated for interpolated VMRs, temperature and pressure.
1486 
1487  // To use special interpolation functions for atmospheric fields we
1488  // use ext_mat_field and abs_vec_field:
1489 
1490  Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1491  1,
1492  1,
1493  stokes_dim,
1494  stokes_dim,
1495  0.);
1496  Tensor4 abs_vec_field(
1497  cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
1498 
1499  //Loop over all directions, defined by za_grid
1500  for (za_index = 0; za_index < N_scat_za; za_index++) {
1501  //Only dummy variables:
1502  Index aa_index = 0;
1503 
1504  cloud_fieldsCalc(ws,
1505  ext_mat_field,
1506  abs_vec_field,
1507  spt_calc_agenda,
1508  za_index,
1509  aa_index,
1510  cloudbox_limits,
1511  t_field,
1512  pnd_field,
1513  verbosity);
1514 
1515  //======================================================================
1516  // Radiative transfer inside the cloudbox
1517  //=====================================================================
1518 
1519  Vector stokes_vec(stokes_dim, 0.);
1520  // Sequential update for uplooking angles
1521  if (za_grid[za_index] <= 90) {
1522  // Loop over all positions inside the cloud box defined by the
1523  // cloudbox_limits exculding the upper boundary. For uplooking
1524  // directions, we start from cloudbox_limits[1]-1 and go down
1525  // to cloudbox_limits[0] to do a sequential update of the
1526  // aradiation field
1527 
1528  // Loop over all positions inside the cloudbox defined by the
1529  // cloudbox_limits.
1530  for (Index p_index = cloudbox_limits[1] - 1;
1531  p_index >= cloudbox_limits[0];
1532  p_index--) {
1533  cloud_ppath_update1D_planeparallel(ws,
1534  cloudbox_field_mono,
1535  p_index,
1536  za_index,
1537  za_grid,
1538  cloudbox_limits,
1539  doit_scat_field,
1540  propmat_clearsky_agenda,
1541  vmr_field,
1542  p_grid,
1543  z_field,
1544  t_field,
1545  f_grid,
1546  f_index,
1547  ext_mat_field,
1548  abs_vec_field,
1549  verbosity);
1550  }
1551  } else if (za_grid[za_index] > 90) {
1552  //
1553  // Sequential updating for downlooking angles
1554  //
1555  for (Index p_index = cloudbox_limits[0] + 1;
1556  p_index <= cloudbox_limits[1];
1557  p_index++) {
1558  cloud_ppath_update1D_planeparallel(ws,
1559  cloudbox_field_mono,
1560  p_index,
1561  za_index,
1562  za_grid,
1563  cloudbox_limits,
1564  doit_scat_field,
1565  propmat_clearsky_agenda,
1566  vmr_field,
1567  p_grid,
1568  z_field,
1569  t_field,
1570  f_grid,
1571  f_index,
1572  ext_mat_field,
1573  abs_vec_field,
1574  verbosity);
1575  } // Close loop over p_grid (inside cloudbox).
1576  } // end if downlooking.
1577 
1578  } // Closes loop over za_grid.
1579 }
1580 
1581 /* Workspace method: Doxygen documentation will be auto-generated */
1582 void DoitInit( //WS Output
1583  Tensor6& doit_scat_field,
1584  Tensor7& cloudbox_field,
1585  Index& doit_is_initialized,
1586  // WS Input
1587  const Index& stokes_dim,
1588  const Index& atmosphere_dim,
1589  const Vector& f_grid,
1590  const Vector& za_grid,
1591  const Vector& aa_grid,
1592  const Index& doit_za_grid_size,
1593  const Index& cloudbox_on,
1594  const ArrayOfIndex& cloudbox_limits,
1595  const Verbosity& verbosity) {
1596  if (!cloudbox_on) {
1597  CREATE_OUT0;
1598  doit_is_initialized = 0;
1599  out0 << " Cloudbox is off, DOIT calculation will be skipped.\n";
1600  return;
1601  //throw runtime_error( "Cloudbox is off, no scattering calculations to be"
1602  // "performed." );
1603  }
1604 
1605  // -------------- Check the input ------------------------------
1606 
1607  if (stokes_dim < 0 || stokes_dim > 4)
1608  throw runtime_error(
1609  "The dimension of stokes vector must be"
1610  "1,2,3, or 4");
1611 
1612  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1613 
1614  // Number of zenith angles.
1615  const Index N_scat_za = za_grid.nelem();
1616 
1617  // The recommended values were found by testing the accuracy and the speed of
1618  // 1D DOIT calculations for different grid sizes. For 3D calculations it can
1619  // be necessary to use more grid points.
1620  if (N_scat_za < 16)
1621  throw runtime_error(
1622  "For accurate results, za_grid must have "
1623  "more than 15 elements.");
1624  else if (N_scat_za > 100) {
1625  CREATE_OUT1;
1626  out1 << "Warning: za_grid is very large, which means that the \n"
1627  << "calculation will be very slow.\n";
1628  }
1629 
1630  if (za_grid[0] != 0. || za_grid[N_scat_za - 1] != 180.)
1631  throw runtime_error("The range of *za_grid* must [0 180].");
1632 
1633  if (!is_increasing(za_grid))
1634  throw runtime_error("*za_grid* must be increasing.");
1635 
1636  // Number of azimuth angles.
1637  const Index N_scat_aa = aa_grid.nelem();
1638 
1639  if (N_scat_aa < 6)
1640  throw runtime_error(
1641  "For accurate results, aa_grid must have "
1642  "more than 5 elements.");
1643  else if (N_scat_aa > 100) {
1644  CREATE_OUT1;
1645  out1 << "Warning: aa_grid is very large which means that the \n"
1646  << "calculation will be very slow.\n";
1647  }
1648 
1649  if (aa_grid[0] != 0. || aa_grid[N_scat_aa - 1] != 360.)
1650  throw runtime_error("The range of *aa_grid* must [0 360].");
1651 
1652  if (doit_za_grid_size < 16)
1653  throw runtime_error(
1654  "*doit_za_grid_size* must be greater than 15 for accurate results");
1655  else if (doit_za_grid_size > 100) {
1656  CREATE_OUT1;
1657  out1 << "Warning: doit_za_grid_size is very large which means that the \n"
1658  << "calculation will be very slow.\n";
1659  }
1660 
1661  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
1662  throw runtime_error(
1663  "*cloudbox_limits* is a vector which contains the"
1664  "upper and lower limit of the cloud for all "
1665  "atmospheric dimensions. So its dimension must"
1666  "be 2 x *atmosphere_dim*");
1667 
1668  //------------- end of checks ---------------------------------------
1669 
1670  const Index Nf = f_grid.nelem();
1671  const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1672  const Index Nza = za_grid.nelem();
1673  const Index Ns = stokes_dim;
1674 
1675  // Resize and initialize radiation field in the cloudbox
1676  if (atmosphere_dim == 1) {
1677  cloudbox_field.resize(Nf, Np_cloud, 1, 1, Nza, 1, Ns);
1678  doit_scat_field.resize(Np_cloud, 1, 1, Nza, 1, Ns);
1679  } else if (atmosphere_dim == 3) {
1680  const Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1681  const Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
1682  const Index Naa = aa_grid.nelem();
1683 
1684  cloudbox_field.resize(Nf, Np_cloud, Nlat_cloud, Nlon_cloud, Nza, Naa, Ns);
1685  doit_scat_field.resize(Np_cloud, Nlat_cloud, Nlon_cloud, Nza, Naa, Ns);
1686  } else {
1687  throw runtime_error(
1688  "Scattering calculations are not possible for a 2D"
1689  "atmosphere. If you want to do scattering calculations"
1690  "*atmosphere_dim* has to be either 1 or 3");
1691  }
1692 
1693  cloudbox_field = NAN;
1694  doit_scat_field = NAN;
1695  doit_is_initialized = 1;
1696 }
1697 
1698 /* Workspace method: Doxygen documentation will be auto-generated */
1699 void cloudbox_field_monoOptimizeReverse( //WS input
1700  Tensor6& cloudbox_field_mono,
1701  const Vector& p_grid_orig,
1702  const Vector& p_grid,
1703  const ArrayOfIndex& cloudbox_limits,
1704  const Verbosity& /* verbosity */) {
1705  Tensor6 cloudbox_field_mono_opt(
1706  p_grid_orig.nelem(), 1, 1, cloudbox_field_mono.npages(), 1, 1);
1707  ArrayOfGridPos Z_gp(p_grid_orig.nelem());
1708  Matrix itw_z(Z_gp.nelem(), 2);
1709  // We only need the p_grid inside the cloudbox as cloudbox_field_mono is only defined in the
1710  // cloudbox
1711  Vector p_grid_cloudbox = p_grid[Range(
1712  cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1713  ostringstream os;
1714  os << "There is a problem with the pressure grid interpolation";
1715  chk_interpolation_grids(os.str(), p_grid, p_grid_orig);
1716 
1717  // Gridpositions:
1718  gridpos(Z_gp, p_grid_cloudbox, p_grid_orig);
1719  // Interpolation weights:
1720  interpweights(itw_z, Z_gp);
1721  // Interpolate cloudbox_field_mono
1722  for (Index i = 0; i < cloudbox_field_mono.npages(); i++) {
1723  interp(cloudbox_field_mono_opt(joker, 0, 0, i, 0, 0),
1724  itw_z,
1725  cloudbox_field_mono(joker, 0, 0, i, 0, 0),
1726  Z_gp);
1727  }
1728  cloudbox_field_mono = cloudbox_field_mono_opt;
1729 }
1730 
1731 /* Workspace method: Doxygen documentation will be auto-generated */
1732 void OptimizeDoitPressureGrid(
1733  Workspace& ws,
1734  //WS input
1735  Vector& p_grid,
1736  Tensor4& pnd_field,
1737  Tensor3& t_field,
1738  ArrayOfArrayOfSingleScatteringData& scat_data_mono,
1739  Tensor3& z_field,
1740  ArrayOfIndex& cloudbox_limits,
1741  Tensor6& cloudbox_field_mono,
1742  Tensor7& pha_mat_doit,
1743  Tensor4& vmr_field,
1744  Vector& p_grid_orig,
1745  const Vector& f_grid,
1746  const Index& f_index,
1747  const Agenda& propmat_clearsky_agenda,
1748  const Numeric& tau_scat_max,
1749  const Numeric& sgl_alb_max,
1750  const Index& cloudbox_size_max,
1751  const Verbosity& verbosity) {
1752  CREATE_OUT2;
1753  CREATE_OUT3;
1754  // Make sure that stokes_dim = 1 and that ScatSpeciesMerge has been applied:
1755  if (cloudbox_field_mono.ncols() != 1)
1756  throw runtime_error(
1757  " This method currently only works for unpolarized radiation "
1758  "( stokes_dim = 1)");
1759  // If ScatSpeciesMerged has been applied, then scat_data_mono should have only one element, and
1760  // pnd_field should have the number of pressure grid points in the cloudbox as first dimension
1761  if (scat_data_mono.nelem() > 1 ||
1762  pnd_field.nbooks() != cloudbox_limits[1] - cloudbox_limits[0] + 1)
1763  throw runtime_error(
1764  " ScatSpeciesMerge has to be applied in order to use this method");
1765 
1766  bool was_too_much = false;
1767  Numeric tau_scat_max_internal = tau_scat_max;
1768  ArrayOfSingleScatteringData& scat_data_local = scat_data_mono[0];
1769  p_grid_orig = p_grid;
1770  vector<Numeric> z_grid_new;
1771  Vector ext_mat(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1772  Vector abs_vec(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1773  Vector scat_vec(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1774  ArrayOfIndex cloudbox_limits_opt(2);
1775  z_grid_new.reserve(1000);
1776  for (Index i = 0; i < cloudbox_limits[0]; i++)
1777  z_grid_new.push_back(z_field(i, 0, 0));
1778  //-----------------------------------------------
1779  // Calculate optical thicknesses of the layers
1780  //------------------------------------------------
1781 
1782  // Fields for scalar gas absorption
1783  const Vector rtp_mag_dummy(3, 0);
1784  const Vector ppath_los_dummy;
1785  ArrayOfStokesVector nlte_dummy;
1786  ArrayOfPropagationMatrix partial_dummy;
1787  ArrayOfStokesVector partial_source_dummy, partial_nlte_dummy;
1788  EnergyLevelMap rtp_nlte_dummy;
1789  ArrayOfPropagationMatrix cur_propmat_clearsky;
1790 
1791  Index scat_data_insert_offset = 0;
1792  Vector single_scattering_albedo(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1793  0.);
1794  for (Index k = cloudbox_limits[0]; k < cloudbox_limits[1] + 1; ++k) {
1795  Index cloudbox_index = k - cloudbox_limits[0];
1796  Numeric abs_coeff = 0;
1797  // Scattering particles
1798  ext_mat[cloudbox_index] =
1799  scat_data_local[cloudbox_index].ext_mat_data(0, 0, 0, 0, 0);
1800  abs_vec[cloudbox_index] =
1801  scat_data_local[cloudbox_index].abs_vec_data(0, 0, 0, 0, 0);
1802  scat_vec[cloudbox_index] =
1803  ext_mat[cloudbox_index] - abs_vec[cloudbox_index];
1804  // Calculate scalar gas absorption
1805  propmat_clearsky_agendaExecute(ws,
1806  cur_propmat_clearsky,
1807  nlte_dummy,
1808  partial_dummy,
1809  partial_source_dummy,
1810  partial_nlte_dummy,
1811  ArrayOfRetrievalQuantity(0),
1812  f_grid[Range(f_index, 1)],
1813  rtp_mag_dummy,
1814  ppath_los_dummy,
1815  p_grid[k],
1816  t_field(k, 0, 0),
1817  rtp_nlte_dummy,
1818  vmr_field(joker, k, 0, 0),
1819  propmat_clearsky_agenda);
1820  for (auto& pm : cur_propmat_clearsky) {
1821  abs_coeff += pm.Kjj()[0];
1822  }
1823  abs_coeff /= (Numeric)cur_propmat_clearsky.nelem();
1824  single_scattering_albedo[cloudbox_index] =
1825  scat_vec[cloudbox_index] / (ext_mat[cloudbox_index] + abs_coeff);
1826  }
1827  //
1828  // Try out the current settings of tau_scat_max and sgl_alb_max
1829  //
1830  do {
1831  scat_data_insert_offset = 0;
1832  for (Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1833  Index cloudbox_index = k - cloudbox_limits[0];
1834  const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1835  single_scattering_albedo[cloudbox_index + 1]) /
1836  2;
1837  const Numeric scat_opt_thk =
1838  (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1839  ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1840  if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1841  Index factor = (Index)ceil(scat_opt_thk / tau_scat_max_internal);
1842  for (Index j = 1; j < factor; j++) {
1843  scat_data_insert_offset++;
1844  }
1845  }
1846  }
1847  // If enhancement is too large, change tau_scat_max to a higher value:
1848  if (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] + 1 >
1849  cloudbox_size_max) {
1850  tau_scat_max_internal += 0.01;
1851  was_too_much = true;
1852  }
1853  } while (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] +
1854  1 >
1855  cloudbox_size_max);
1856  scat_data_insert_offset = 0;
1857  // Give warning if enhancement was too much and threshold had to be changed:
1858  if (was_too_much) {
1859  out2
1860  << "Warning: Pressure grid optimization with the thresholds tau_scat_max = "
1861  << tau_scat_max << " and sgl_slb_max = " << sgl_alb_max
1862  << " has lead to an enhancement larger than the value of"
1863  << " cloudbox_size_max = " << cloudbox_size_max
1864  << ". This is why I changed tau_scat_max to " << tau_scat_max_internal
1865  << ". This may lead to errors of a too coarse grid! \n";
1866  }
1867 
1868  //--------------------------------------
1869  //Optimize the altitude grid z_grid
1870  //---------------------------------------
1871  for (Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1872  Index cloudbox_index = k - cloudbox_limits[0];
1873  const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1874  single_scattering_albedo[cloudbox_index + 1]) /
1875  2;
1876  const Numeric scat_opt_thk =
1877  (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1878  ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1879  z_grid_new.push_back(z_field(k, 0, 0));
1880 
1881  if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1882  Index factor = (Index)ceil(scat_opt_thk / tau_scat_max_internal);
1883  Numeric step =
1884  (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) / (Numeric)factor;
1885  const SingleScatteringData nextLayer =
1886  scat_data_local[cloudbox_index + scat_data_insert_offset + 1];
1887  const SingleScatteringData currentLayer =
1888  scat_data_local[cloudbox_index + scat_data_insert_offset];
1889 
1890  for (Index j = 1; j < factor; j++) {
1891  z_grid_new.push_back(z_field(k, 0, 0) + (Numeric)j * step);
1892  // Perform manual interpolation of scat_data
1893  const Numeric weight = (Numeric)j / (Numeric)factor;
1894  SingleScatteringData newLayer = currentLayer;
1895  Tensor7 weightednextLayerPhamat = nextLayer.pha_mat_data;
1896  Tensor5 weightednextLayerExtmat = nextLayer.ext_mat_data;
1897  Tensor5 weightednextLayerAbsvec = nextLayer.abs_vec_data;
1898 
1899  weightednextLayerPhamat *= weight;
1900  weightednextLayerExtmat *= weight;
1901  weightednextLayerAbsvec *= weight;
1902 
1903  newLayer.pha_mat_data *= 1. - weight;
1904  newLayer.ext_mat_data *= 1. - weight;
1905  newLayer.abs_vec_data *= 1. - weight;
1906 
1907  newLayer.pha_mat_data += weightednextLayerPhamat;
1908  newLayer.ext_mat_data += weightednextLayerExtmat;
1909  newLayer.abs_vec_data += weightednextLayerAbsvec;
1910 
1911  // Optimize scat_data
1912  scat_data_local.insert(scat_data_local.begin() + cloudbox_index +
1913  scat_data_insert_offset + 1,
1914  std::move(newLayer));
1915 
1916  scat_data_insert_offset++;
1917  }
1918  }
1919  }
1920  // New cloudbox limits
1921  cloudbox_limits_opt[0] = cloudbox_limits[0];
1922  cloudbox_limits_opt[1] = scat_data_insert_offset + cloudbox_limits[1];
1923  const Index cloudbox_opt_size =
1924  cloudbox_limits_opt[1] - cloudbox_limits_opt[0] + 1;
1925  out3 << "Frequency: " << f_grid[f_index]
1926  << " old: " << cloudbox_limits[1] - cloudbox_limits[0] + 1
1927  << " new: " << cloudbox_opt_size << " Factor: "
1928  << (Numeric)cloudbox_opt_size /
1929  (Numeric)(cloudbox_limits[1] - cloudbox_limits[0] + 1)
1930  << "\n";
1931 
1932  for (Index i = cloudbox_limits[1]; i < z_field.npages(); i++)
1933  z_grid_new.push_back(z_field(i, 0, 0));
1934 
1935  Vector z_grid(z_grid_new.size());
1936  for (Index i = 0; i < z_grid.nelem(); i++) z_grid[i] = z_grid_new[i];
1937  p_grid_orig = p_grid[Range(cloudbox_limits[0],
1938  cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1939  // ---------------------------------------
1940  // Interpolate fields to new z_grid
1941  // ----------------------------------------
1942  ArrayOfArrayOfSingleScatteringData scat_data_new;
1943  Tensor3 t_field_new(z_grid.nelem(), 1, 1);
1944  Vector p_grid_opt(z_grid.nelem());
1945  Tensor6 cloudbox_field_mono_opt(
1946  cloudbox_opt_size, 1, 1, cloudbox_field_mono.npages(), 1, 1);
1947  Tensor7 pha_mat_doit_opt(cloudbox_opt_size,
1948  pha_mat_doit.nvitrines(),
1949  1,
1950  pha_mat_doit.nbooks(),
1951  pha_mat_doit.npages(),
1952  1,
1953  1);
1954  ArrayOfGridPos Z_gp(z_grid.nelem());
1955  Matrix itw_z(Z_gp.nelem(), 2);
1956  ostringstream os;
1957  os << "At the current frequency " << f_grid[f_index]
1958  << " there was an error while interpolating the fields to the new z_field";
1959  chk_interpolation_grids(os.str(), z_field(joker, 0, 0), z_grid);
1960 
1961  // Gridpositions of interpolation:
1962  gridpos(Z_gp, z_field(joker, 0, 0), z_grid);
1963  // Interpolation weights:
1964  interpweights(itw_z, Z_gp);
1965 
1966  // Interpolate Temperature
1967  interp(t_field_new(joker, 0, 0), itw_z, t_field(joker, 0, 0), Z_gp);
1968  // Write new Temperature to scat_data
1969  for (Index k = cloudbox_limits_opt[0]; k < cloudbox_limits_opt[1]; k++) {
1970  Index i = k - cloudbox_limits[0];
1971  scat_data_local[i].T_grid = t_field_new(i, 0, 0);
1972  }
1973 
1974  // Interpolate p_grid
1975  interp(p_grid_opt, itw_z, p_grid, Z_gp);
1976 
1977  // Interpolate vmr_field
1978  Tensor4 vmr_field_opt(vmr_field.nbooks(), p_grid_opt.nelem(), 1, 1);
1979  for (Index i = 0; i < vmr_field.nbooks(); i++)
1980  interp(
1981  vmr_field_opt(i, joker, 0, 0), itw_z, vmr_field(i, joker, 0, 0), Z_gp);
1982 
1983  // Interpolate cloudbox_field_mono and pha_mat_doit
1984  ArrayOfGridPos Z_gp_2(cloudbox_opt_size);
1985  Matrix itw_z_2(Z_gp_2.nelem(), 2);
1986  Range r1 =
1987  Range(cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1);
1988  Range r2 = Range(cloudbox_limits_opt[0], cloudbox_opt_size);
1989  chk_interpolation_grids(os.str(), z_field(r1, 0, 0), z_grid[r2]);
1990  gridpos(Z_gp_2,
1991  z_field(Range(cloudbox_limits[0],
1992  cloudbox_limits[1] - cloudbox_limits[0] + 1),
1993  0,
1994  0),
1995  z_grid[Range(cloudbox_limits_opt[0], cloudbox_opt_size)]);
1996  interpweights(itw_z_2, Z_gp_2);
1997 
1998  for (Index i = 0; i < cloudbox_field_mono.npages(); i++) {
1999  interp(cloudbox_field_mono_opt(joker, 0, 0, i, 0, 0),
2000  itw_z_2,
2001  cloudbox_field_mono(joker, 0, 0, i, 0, 0),
2002  Z_gp_2);
2003  }
2004  for (Index i = 0; i < pha_mat_doit.nvitrines(); i++) {
2005  for (Index j = 0; j < pha_mat_doit.nbooks(); j++) {
2006  for (Index k = 0; k < pha_mat_doit.npages(); k++) {
2007  interp(pha_mat_doit_opt(joker, i, 0, j, k, 0, 0),
2008  itw_z_2,
2009  pha_mat_doit(joker, i, 0, j, k, 0, 0),
2010  Z_gp_2);
2011  }
2012  }
2013  }
2014 
2015  // Interpolate pnd-field
2016  pnd_field.resize(cloudbox_opt_size, cloudbox_opt_size, 1, 1);
2017  pnd_field = 0.;
2018  for (Index i = 0; i < cloudbox_opt_size; i++) pnd_field(i, i, 0, 0) = 1.;
2019 
2020  //Write new fields
2021  p_grid = p_grid_opt;
2022  t_field = t_field_new;
2023  cloudbox_limits = cloudbox_limits_opt;
2024  cloudbox_field_mono = cloudbox_field_mono_opt;
2025  pha_mat_doit = pha_mat_doit_opt;
2026  z_field.resize(z_grid.nelem(), 1, 1);
2027  z_field(joker, 0, 0) = z_grid;
2028  vmr_field = vmr_field_opt;
2029 }
2030 
2031 /* Workspace method: Doxygen documentation will be auto-generated */
2032 void DoitWriteIterationFields( //WS input
2033  const Index& doit_iteration_counter,
2034  const Tensor6& cloudbox_field_mono,
2035  const Index& f_index,
2036  //Keyword:
2037  const ArrayOfIndex& iterations,
2038  const ArrayOfIndex& frequencies,
2039  const Verbosity& verbosity) {
2040  if (!frequencies.nelem() || !iterations.nelem()) return;
2041 
2042  // If the number of iterations is less than a number specified in the
2043  // keyword *iterations*, this number will be ignored.
2044 
2045  ostringstream os;
2046  os << "doit_iteration_f" << f_index << "_i" << doit_iteration_counter;
2047 
2048  // All iterations for all frequencies are written to files
2049  if (frequencies[0] == -1 && iterations[0] == -1) {
2050  xml_write_to_file(
2051  os.str() + ".xml", cloudbox_field_mono, FILE_TYPE_ASCII, 0, verbosity);
2052  }
2053 
2054  for (Index j = 0; j < frequencies.nelem(); j++) {
2055  if (f_index == frequencies[j] || (!j && frequencies[j] == -1)) {
2056  // All iterations are written to files
2057  if (iterations[0] == -1) {
2058  xml_write_to_file(os.str() + ".xml",
2059  cloudbox_field_mono,
2060  FILE_TYPE_ASCII,
2061  0,
2062  verbosity);
2063  }
2064 
2065  // Only the iterations given by the keyword are written to a file
2066  else {
2067  for (Index i = 0; i < iterations.nelem(); i++) {
2068  if (doit_iteration_counter == iterations[i])
2069  xml_write_to_file(os.str() + ".xml",
2070  cloudbox_field_mono,
2071  FILE_TYPE_ASCII,
2072  0,
2073  verbosity);
2074  }
2075  }
2076  }
2077  }
2078 }
2079 
2080 /* Workspace method: Doxygen documentation will be auto-generated */
2081 void doit_scat_fieldCalc(Workspace& ws,
2082  // WS Output and Input
2083  Tensor6& doit_scat_field,
2084  //WS Input:
2085  const Agenda& pha_mat_spt_agenda,
2086  const Tensor6& cloudbox_field_mono,
2087  const Tensor4& pnd_field,
2088  const Tensor3& t_field,
2089  const Index& atmosphere_dim,
2090  const ArrayOfIndex& cloudbox_limits,
2091  const Vector& za_grid,
2092  const Vector& aa_grid,
2093  const Index& doit_za_grid_size,
2094  const Tensor7& pha_mat_doit,
2095  const Verbosity& verbosity)
2096 
2097 {
2098  CREATE_OUT2;
2099  CREATE_OUT3;
2100 
2101  // ------------ Check the input -------------------------------
2102 
2103  // Agenda for calculation of phase matrix
2104  chk_not_empty("pha_mat_spt_agenda", pha_mat_spt_agenda);
2105 
2106  // Number of zenith angles.
2107  const Index Nza = za_grid.nelem();
2108 
2109  if (za_grid[0] != 0. || za_grid[Nza - 1] != 180.)
2110  throw runtime_error("The range of *za_grid* must [0 180].");
2111 
2112  // Number of azimuth angles.
2113  const Index Naa = aa_grid.nelem();
2114 
2115  if (Naa > 1 && (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.))
2116  throw runtime_error("The range of *aa_grid* must [0 360].");
2117 
2118  // Get stokes dimension from *doit_scat_field*:
2119  const Index stokes_dim = doit_scat_field.ncols();
2120  assert(stokes_dim > 0 || stokes_dim < 5);
2121 
2122  // Size of particle number density field can not be checked here,
2123  // because the function does not use the atmospheric grids.
2124  // Check will be included after fixing the size of pnd field, so
2125  // that it is defined only inside the cloudbox.
2126 
2127  // Check atmospheric dimension and dimensions of
2128  // radiation field (*cloudbox_field*) and scattering integral field
2129  // (*doit_scat_field*)
2130  if (atmosphere_dim == 1) {
2131  assert(is_size(cloudbox_field_mono,
2132  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2133  1,
2134  1,
2135  Nza,
2136  1,
2137  stokes_dim));
2138  assert(is_size(doit_scat_field,
2139  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2140  1,
2141  1,
2142  za_grid.nelem(),
2143  1,
2144  stokes_dim));
2145  } else if (atmosphere_dim == 3) {
2146  assert(is_size(cloudbox_field_mono,
2147  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2148  (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2149  (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2150  Nza,
2151  Naa,
2152  stokes_dim));
2153  assert(is_size(doit_scat_field,
2154  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2155  (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2156  (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2157  Nza,
2158  Naa,
2159  stokes_dim));
2160  } else {
2161  ostringstream os;
2162  os << "The atmospheric dimension must be 1D or 3D \n"
2163  << "for scattering calculations using the DOIT\n"
2164  << "module, but it is not. The value of *atmosphere_dim*\n"
2165  << "is " << atmosphere_dim << ".";
2166  throw runtime_error(os.str());
2167  }
2168 
2169  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
2170  throw runtime_error(
2171  "*cloudbox_limits* is a vector which contains the"
2172  "upper and lower limit of the cloud for all "
2173  "atmospheric dimensions. So its dimension must"
2174  "be 2 x *atmosphere_dim*");
2175 
2176  // This function should only be used for down-looking cases where no
2177  // optimized zenith angle grid is required.
2178  if (doit_za_grid_size != Nza)
2179  throw runtime_error(
2180  "The zenith angle grids for the computation of\n"
2181  "the scattering integral and the RT part must \n"
2182  "be equal. Check definitions in \n"
2183  "*DOAngularGridsSet*. The keyword \n"
2184  "'za_grid_opt_file' should be empty. \n");
2185 
2186  // ------ end of checks -----------------------------------------------
2187 
2188  // Initialize variables *pha_mat* and *pha_mat_spt*
2189  Tensor4 pha_mat_local(
2190  doit_za_grid_size, aa_grid.nelem(), stokes_dim, stokes_dim, 0.);
2191 
2192  Tensor5 pha_mat_spt_local(pnd_field.nbooks(),
2193  doit_za_grid_size,
2194  aa_grid.nelem(),
2195  stokes_dim,
2196  stokes_dim,
2197  0.);
2198 
2199  // Equidistant step size for integration
2200  Vector grid_stepsize(2);
2201  grid_stepsize[0] = 180. / (Numeric)(doit_za_grid_size - 1);
2202  if (Naa > 1) {
2203  grid_stepsize[1] = 360. / (Numeric)(Naa - 1);
2204  }
2205 
2206  Tensor3 product_field(Nza, Naa, stokes_dim, 0);
2207 
2208  out2 << " Calculate the scattered field\n";
2209 
2210  if (atmosphere_dim == 1) {
2211  // Get pha_mat at the grid positions
2212  // Since atmosphere_dim = 1, there is no loop over lat and lon grids
2213  for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2214  p_index++) {
2215  //There is only loop over zenith angle grid ; no azimuth angle grid.
2216  for (Index za_index_local = 0; za_index_local < Nza; za_index_local++) {
2217  // Calculate the phase matrix of individual scattering elements
2218  out3 << "Multiplication of phase matrix with incoming"
2219  << " intensities \n";
2220 
2221  product_field = 0;
2222 
2223  // za_in and aa_in are for incoming zenith and azimuth
2224  //angle direction for which pha_mat is calculated
2225  for (Index za_in = 0; za_in < Nza; ++za_in) {
2226  for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2227  // Multiplication of phase matrix with incoming
2228  // intensity field.
2229 
2230  for (Index i = 0; i < stokes_dim; i++) {
2231  for (Index j = 0; j < stokes_dim; j++) {
2232  product_field(za_in, aa_in, i) +=
2233  pha_mat_doit(
2234  p_index, za_index_local, 0, za_in, aa_in, i, j) *
2235  cloudbox_field_mono(p_index, 0, 0, za_in, 0, j);
2236  }
2237  }
2238 
2239  } //end aa_in loop
2240  } //end za_in loop
2241  //integration of the product of ifield_in and pha
2242  // over zenith angle and azimuth angle grid. It calls
2243  if (Naa == 1) {
2244  for (Index i = 0; i < stokes_dim; i++) {
2245  doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2246  AngIntegrate_trapezoid(product_field(joker, 0, i), za_grid) /
2247  2 / PI;
2248  } //end i loop
2249  } else {
2250  for (Index i = 0; i < stokes_dim; i++) {
2251  doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2252  AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2253  za_grid,
2254  aa_grid,
2255  grid_stepsize);
2256  }
2257  }
2258  } //end za_prop loop
2259  } //end p_index loop
2260  } //end atmosphere_dim = 1
2261 
2262  //atmosphere_dim = 3
2263  else if (atmosphere_dim == 3) {
2264  /*there is a loop over pressure, latitude and longitudeindex
2265  when we calculate the pha_mat from pha_mat_spt and pnd_field
2266  using the method pha_matCalc. */
2267 
2268  for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2269  p_index++) {
2270  for (Index lat_index = 0;
2271  lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2272  lat_index++) {
2273  for (Index lon_index = 0;
2274  lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2275  lon_index++) {
2276  Numeric rtp_temperature_local =
2277  t_field(p_index + cloudbox_limits[0],
2278  lat_index + cloudbox_limits[2],
2279  lon_index + cloudbox_limits[4]);
2280 
2281  for (Index aa_index_local = 1; aa_index_local < Naa;
2282  aa_index_local++) {
2283  for (Index za_index_local = 0; za_index_local < Nza;
2284  za_index_local++) {
2285  out3 << "Calculate phase matrix \n";
2286  pha_mat_spt_agendaExecute(ws,
2287  pha_mat_spt_local,
2288  za_index_local,
2289  lat_index,
2290  lon_index,
2291  p_index,
2292  aa_index_local,
2293  rtp_temperature_local,
2294  pha_mat_spt_agenda);
2295 
2296  pha_matCalc(pha_mat_local,
2297  pha_mat_spt_local,
2298  pnd_field,
2299  atmosphere_dim,
2300  p_index,
2301  lat_index,
2302  lon_index,
2303  verbosity);
2304 
2305  product_field = 0;
2306 
2307  //za_in and aa_in are the incoming directions
2308  //for which pha_mat_spt is calculated
2309  for (Index za_in = 0; za_in < Nza; ++za_in) {
2310  for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2311  // Multiplication of phase matrix
2312  // with incloming intensity field.
2313  for (Index i = 0; i < stokes_dim; i++) {
2314  for (Index j = 0; j < stokes_dim; j++) {
2315  product_field(za_in, aa_in, i) +=
2316  pha_mat_local(za_in, aa_in, i, j) *
2317  cloudbox_field_mono(p_index,
2318  lat_index,
2319  lon_index,
2320  za_index_local,
2321  aa_index_local,
2322  j);
2323  }
2324  }
2325  } //end aa_in loop
2326  } //end za_in loop
2327  //integration of the product of ifield_in and pha
2328  //over zenith angle and azimuth angle grid. It
2329  //calls here the integration routine
2330  //AngIntegrate_trapezoid_opti
2331  for (Index i = 0; i < stokes_dim; i++) {
2332  doit_scat_field(p_index,
2333  lat_index,
2334  lon_index,
2335  za_index_local,
2336  aa_index_local,
2337  i) =
2338  AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2339  za_grid,
2340  aa_grid,
2341  grid_stepsize);
2342  } //end i loop
2343  } //end aa_prop loop
2344  } //end za_prop loop
2345  } //end lon loop
2346  } // end lat loop
2347  } // end p loop
2348  // aa = 0 is the same as aa = 180:
2349  doit_scat_field(joker, joker, joker, joker, 0, joker) =
2350  doit_scat_field(joker, joker, joker, joker, Naa - 1, joker);
2351  } // end atmosphere_dim = 3
2352 }
2353 
2354 /* Workspace method: Doxygen documentation will be auto-generated */
2355 void doit_scat_fieldCalcLimb(Workspace& ws,
2356  // WS Output and Input
2357  Tensor6& doit_scat_field,
2358  //WS Input:
2359  const Agenda& pha_mat_spt_agenda,
2360  const Tensor6& cloudbox_field_mono,
2361  const Tensor4& pnd_field,
2362  const Tensor3& t_field,
2363  const Index& atmosphere_dim,
2364  const ArrayOfIndex& cloudbox_limits,
2365  const Vector& za_grid,
2366  const Vector& aa_grid,
2367  const Index& doit_za_grid_size,
2368  const Index& doit_za_interp,
2369  const Tensor7& pha_mat_doit,
2370  const Verbosity& verbosity) {
2371  CREATE_OUT2;
2372  CREATE_OUT3;
2373 
2374  // ------------ Check the input -------------------------------
2375 
2376  // Agenda for calculation of phase matrix
2377  chk_not_empty("pha_mat_spt_agenda", pha_mat_spt_agenda);
2378 
2379  // Number of zenith angles.
2380  const Index Nza = za_grid.nelem();
2381 
2382  if (za_grid[0] != 0. || za_grid[Nza - 1] != 180.)
2383  throw runtime_error("The range of *za_grid* must [0 180].");
2384 
2385  // Number of azimuth angles.
2386  const Index Naa = aa_grid.nelem();
2387 
2388  if (Naa > 1 && (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.))
2389  throw runtime_error("The range of *aa_grid* must [0 360].");
2390 
2391  // Get stokes dimension from *doit_scat_field*:
2392  const Index stokes_dim = doit_scat_field.ncols();
2393  assert(stokes_dim > 0 || stokes_dim < 5);
2394 
2395  // Size of particle number density field can not be checked here,
2396  // because the function does not use the atmospheric grids.
2397  // Check will be included after fixing the size of pnd field, so
2398  // that it is defined only inside the cloudbox.
2399 
2400  // Check atmospheric dimension and dimensions of
2401  // radiation field (*cloudbox_field*) and scattering integral field
2402  // (*doit_scat_field*)
2403  if (atmosphere_dim == 1) {
2404  assert(is_size(cloudbox_field_mono,
2405  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2406  1,
2407  1,
2408  Nza,
2409  1,
2410  stokes_dim));
2411  assert(is_size(doit_scat_field,
2412  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2413  1,
2414  1,
2415  za_grid.nelem(),
2416  1,
2417  stokes_dim));
2418  } else if (atmosphere_dim == 3) {
2419  assert(is_size(cloudbox_field_mono,
2420  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2421  (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2422  (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2423  Nza,
2424  Naa,
2425  stokes_dim));
2426  assert(is_size(doit_scat_field,
2427  (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2428  (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2429  (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2430  Nza,
2431  Naa,
2432  stokes_dim));
2433  } else {
2434  ostringstream os;
2435  os << "The atmospheric dimension must be 1D or 3D \n"
2436  << "for scattering calculations using the DOIT\n"
2437  << "module, but it is not. The value of *atmosphere_dim*\n"
2438  << "is " << atmosphere_dim << ".";
2439  throw runtime_error(os.str());
2440  }
2441 
2442  if (!(doit_za_interp == 0 || doit_za_interp == 1))
2443  throw runtime_error(
2444  "Interpolation method is not defined. Use \n"
2445  "*doit_za_interpSet*.\n");
2446 
2447  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
2448  throw runtime_error(
2449  "*cloudbox_limits* is a vector which contains the"
2450  "upper and lower limit of the cloud for all "
2451  "atmospheric dimensions. So its dimension must"
2452  "be 2 x *atmosphere_dim*");
2453 
2454  if (doit_za_grid_size < 16)
2455  throw runtime_error(
2456  "*doit_za_grid_size* must be greater than 15 for"
2457  "accurate results");
2458  else if (doit_za_grid_size > 100) {
2459  CREATE_OUT1;
2460  out1 << "Warning: doit_za_grid_size is very large which means that the \n"
2461  << "calculation will be very slow. The recommended value is 19.\n";
2462  }
2463 
2464  // ------ end of checks -----------------------------------------------
2465 
2466  // Initialize variables *pha_mat* and *pha_mat_spt*
2467  Tensor4 pha_mat_local(
2468  doit_za_grid_size, aa_grid.nelem(), stokes_dim, stokes_dim, 0.);
2469 
2470  Tensor5 pha_mat_spt_local(pnd_field.nbooks(),
2471  doit_za_grid_size,
2472  aa_grid.nelem(),
2473  stokes_dim,
2474  stokes_dim,
2475  0.);
2476 
2477  // Create the grids for the calculation of the scattering integral.
2478  Vector za_g;
2479  nlinspace(za_g, 0, 180, doit_za_grid_size);
2480 
2481  // Two interpolations are required. First we have to interpolate the
2482  // intensity field on the equidistant grid:
2483  ArrayOfGridPos gp_za_i(doit_za_grid_size);
2484  gridpos(gp_za_i, za_grid, za_g);
2485 
2486  Matrix itw_za_i(doit_za_grid_size, 2);
2487  interpweights(itw_za_i, gp_za_i);
2488 
2489  // Intensity field interpolated on equidistant grid.
2490  Matrix cloudbox_field_int(doit_za_grid_size, stokes_dim, 0);
2491 
2492  // Second, we have to interpolate the scattering integral on the RT
2493  // zenith angle grid.
2494  ArrayOfGridPos gp_za(Nza);
2495  gridpos(gp_za, za_g, za_grid);
2496 
2497  Matrix itw_za(Nza, 2);
2498  interpweights(itw_za, gp_za);
2499 
2500  // Original scattered field, on equidistant zenith angle grid.
2501  Matrix doit_scat_field_org(doit_za_grid_size, stokes_dim, 0);
2502 
2503  // Grid stepsize of zenith and azimuth angle grid, these are needed for the
2504  // integration function.
2505  Vector grid_stepsize(2);
2506  grid_stepsize[0] = 180. / (Numeric)(doit_za_grid_size - 1);
2507  if (Naa > 1) {
2508  grid_stepsize[1] = 360. / (Numeric)(Naa - 1);
2509  }
2510 
2511  Tensor3 product_field(doit_za_grid_size, Naa, stokes_dim, 0);
2512 
2513  if (atmosphere_dim == 1) {
2514  // Get pha_mat at the grid positions
2515  // Since atmosphere_dim = 1, there is no loop over lat and lon grids
2516  for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2517  p_index++) {
2518  // Interpolate intensity field:
2519  for (Index i = 0; i < stokes_dim; i++) {
2520  if (doit_za_interp == 0) {
2521  interp(cloudbox_field_int(joker, i),
2522  itw_za_i,
2523  cloudbox_field_mono(p_index, 0, 0, joker, 0, i),
2524  gp_za_i);
2525  } else if (doit_za_interp == 1) {
2526  // Polynomial
2527  for (Index za = 0; za < za_g.nelem(); za++) {
2528  cloudbox_field_int(za, i) =
2529  interp_poly(za_grid,
2530  cloudbox_field_mono(p_index, 0, 0, joker, 0, i),
2531  za_g[za],
2532  gp_za_i[za]);
2533  }
2534  }
2535  // doit_za_interp must be 0 or 1 (linear or polynomial)!!!
2536  else
2537  assert(false);
2538  }
2539 
2540  //There is only loop over zenith angle grid; no azimuth angle grid.
2541  for (Index za_index_local = 0; za_index_local < doit_za_grid_size;
2542  za_index_local++) {
2543  out3 << "Multiplication of phase matrix with incoming"
2544  << " intensities \n";
2545 
2546  product_field = 0;
2547 
2548  // za_in and aa_in are for incoming zenith and azimuth
2549  // angle direction for which pha_mat is calculated
2550  for (Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2551  for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2552  // Multiplication of phase matrix with incoming
2553  // intensity field.
2554 
2555  for (Index i = 0; i < stokes_dim; i++) {
2556  for (Index j = 0; j < stokes_dim; j++) {
2557  product_field(za_in, aa_in, i) +=
2558  pha_mat_doit(
2559  p_index, za_index_local, 0, za_in, aa_in, i, j) *
2560  cloudbox_field_int(za_in, j);
2561  }
2562  }
2563 
2564  } //end aa_in loop
2565  } //end za_in loop
2566 
2567  out3 << "Compute integral. \n";
2568  if (Naa == 1) {
2569  for (Index i = 0; i < stokes_dim; i++) {
2570  doit_scat_field_org(za_index_local, i) =
2571  AngIntegrate_trapezoid(product_field(joker, 0, i), za_g) / 2 /
2572  PI;
2573  } //end i loop
2574  } else {
2575  for (Index i = 0; i < stokes_dim; i++) {
2576  doit_scat_field_org(za_index_local, i) =
2577  AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2578  za_g,
2579  aa_grid,
2580  grid_stepsize);
2581  }
2582  }
2583 
2584  } //end za_prop loop
2585 
2586  // Interpolation on za_grid, which is used in
2587  // radiative transfer part.
2588  for (Index i = 0; i < stokes_dim; i++) {
2589  if (doit_za_interp == 0) // linear interpolation
2590  {
2591  interp(doit_scat_field(p_index, 0, 0, joker, 0, i),
2592  itw_za,
2593  doit_scat_field_org(joker, i),
2594  gp_za);
2595  } else // polynomial interpolation
2596  {
2597  for (Index za = 0; za < za_grid.nelem(); za++) {
2598  doit_scat_field(p_index, 0, 0, za, 0, i) = interp_poly(
2599  za_g, doit_scat_field_org(joker, i), za_grid[za], gp_za[za]);
2600  }
2601  }
2602  }
2603  } //end p_index loop
2604 
2605  } //end atmosphere_dim = 1
2606 
2607  else if (atmosphere_dim == 3) {
2608  // Loop over all positions
2609  for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2610  p_index++) {
2611  for (Index lat_index = 0;
2612  lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2613  lat_index++) {
2614  for (Index lon_index = 0;
2615  lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2616  lon_index++) {
2617  Numeric rtp_temperature_local =
2618  t_field(p_index + cloudbox_limits[0],
2619  lat_index + cloudbox_limits[2],
2620  lon_index + cloudbox_limits[4]);
2621 
2622  // Loop over scattered directions
2623  for (Index aa_index_local = 1; aa_index_local < Naa;
2624  aa_index_local++) {
2625  // Interpolate intensity field:
2626  for (Index i = 0; i < stokes_dim; i++) {
2627  interp(
2628  cloudbox_field_int(joker, i),
2629  itw_za_i,
2630  cloudbox_field_mono(
2631  p_index, lat_index, lon_index, joker, aa_index_local, i),
2632  gp_za_i);
2633  }
2634 
2635  for (Index za_index_local = 0; za_index_local < doit_za_grid_size;
2636  za_index_local++) {
2637  out3 << "Calculate phase matrix \n";
2638  pha_mat_spt_agendaExecute(ws,
2639  pha_mat_spt_local,
2640  za_index_local,
2641  lat_index,
2642  lon_index,
2643  p_index,
2644  aa_index_local,
2645  rtp_temperature_local,
2646  pha_mat_spt_agenda);
2647 
2648  pha_matCalc(pha_mat_local,
2649  pha_mat_spt_local,
2650  pnd_field,
2651  atmosphere_dim,
2652  p_index,
2653  lat_index,
2654  lon_index,
2655  verbosity);
2656 
2657  product_field = 0;
2658 
2659  //za_in and aa_in are the incoming directions
2660  //for which pha_mat_spt is calculated
2661  out3 << "Multiplication of phase matrix with"
2662  << "incoming intensity \n";
2663 
2664  for (Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2665  for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2666  // Multiplication of phase matrix
2667  // with incloming intensity field.
2668  for (Index i = 0; i < stokes_dim; i++) {
2669  for (Index j = 0; j < stokes_dim; j++) {
2670  product_field(za_in, aa_in, i) +=
2671  pha_mat_local(za_in, aa_in, i, j) *
2672  cloudbox_field_int(za_in, j);
2673  }
2674  }
2675  } //end aa_in loop
2676  } //end za_in loop
2677 
2678  out3 << "Compute the integral \n";
2679 
2680  for (Index i = 0; i < stokes_dim; i++) {
2681  doit_scat_field_org(za_index_local, i) =
2682  AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2683  za_grid,
2684  aa_grid,
2685  grid_stepsize);
2686  } //end stokes_dim loop
2687 
2688  } //end za_prop loop
2689  //Interpolate on original za_grid.
2690  for (Index i = 0; i < stokes_dim; i++) {
2691  interp(
2692  doit_scat_field(
2693  p_index, lat_index, lon_index, joker, aa_index_local, i),
2694  itw_za,
2695  doit_scat_field_org(joker, i),
2696  gp_za);
2697  }
2698  } // end aa_prop loop
2699  } //end lon loop
2700  } //end lat loop
2701  } // end p loop
2702  doit_scat_field(joker, joker, joker, joker, 0, joker) =
2703  doit_scat_field(joker, joker, joker, joker, Naa - 1, joker);
2704  } // end atm_dim=3
2705  out2 << " Finished scattered field.\n";
2706 }
2707 
2708 /* Workspace method: Doxygen documentation will be auto-generated */
2709 void doit_za_grid_optCalc( //WS Output
2710  Vector& doit_za_grid_opt,
2711  // WS Input:
2712  const Tensor6& cloudbox_field_mono,
2713  const Vector& za_grid,
2714  const Index& doit_za_interp,
2715  //Keywords:
2716  const Numeric& acc,
2717  const Verbosity& verbosity) {
2718  CREATE_OUT1;
2719  //-------- Check the input ---------------------------------
2720 
2721  // Here it is checked whether cloudbox_field is 1D and whether it is
2722  // consistent with za_grid. The number of pressure levels and the
2723  // number of stokes components does not matter.
2724  chk_size("cloudbox_field",
2725  cloudbox_field_mono,
2726  cloudbox_field_mono.nvitrines(),
2727  1,
2728  1,
2729  za_grid.nelem(),
2730  1,
2731  cloudbox_field_mono.ncols());
2732 
2733  if (cloudbox_field_mono.ncols() < 1 || cloudbox_field_mono.ncols() > 4)
2734  throw runtime_error(
2735  "The last dimension of *cloudbox_field* corresponds\n"
2736  "to the Stokes dimension, therefore the number of\n"
2737  "columns in *cloudbox_field* must be a number between\n"
2738  "1 and 4, but it is not!");
2739 
2740  if (!(doit_za_interp == 0 || doit_za_interp == 1))
2741  throw runtime_error(
2742  "Interpolation method is not defined. Use \n"
2743  "*doit_za_interpSet*.\n");
2744 
2745  if (za_grid.nelem() < 500) {
2746  out1 << "Warning: The fine grid (*za_grid*) has less than\n"
2747  << "500 grid points which is likely not sufficient for\n"
2748  << "grid_optimization\n";
2749  /* throw runtime_error("The fine grid (*za_grid*) has less than \n"
2750  "500 grid points which is not sufficient for \n"
2751  "grid_optimization");
2752 */
2753  }
2754  // ------------- end of checks -------------------------------------
2755 
2756  // Here only used as dummy variable.
2757  Matrix cloudbox_field_opt_mat;
2758  cloudbox_field_opt_mat = 0.;
2759 
2760  // Optimize zenith angle grid.
2761  za_gridOpt(doit_za_grid_opt,
2762  cloudbox_field_opt_mat,
2763  za_grid,
2764  cloudbox_field_mono,
2765  acc,
2766  doit_za_interp);
2767 }
2768 
2769 /* Workspace method: Doxygen documentation will be auto-generated */
2770 void doit_za_interpSet(Index& doit_za_interp,
2771  const Index& atmosphere_dim,
2772  //Keyword
2773  const String& method,
2774  const Verbosity&) {
2775  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2776 
2777  if (atmosphere_dim != 1 && method == "polynomial")
2778  throw runtime_error(
2779  "Polynomial interpolation is only implemented for\n"
2780  "1D DOIT calculations as \n"
2781  "in 3D there can be numerical problems.\n"
2782  "Please use 'linear' interpolation method.");
2783 
2784  if (method == "linear")
2785  doit_za_interp = 0;
2786  else if (method == "polynomial")
2787  doit_za_interp = 1;
2788  else
2789  throw runtime_error(
2790  "Possible interpolation methods are 'linear' "
2791  "and 'polynomial'.\n");
2792 }
2793 
2794 /* Workspace method: Doxygen documentation will be auto-generated */
2795 void DoitCalc(Workspace& ws,
2796  Tensor7& cloudbox_field,
2797  const Index& atmfields_checked,
2798  const Index& atmgeom_checked,
2799  const Index& cloudbox_checked,
2800  const Index& scat_data_checked,
2801  const Index& cloudbox_on,
2802  const Vector& f_grid,
2803  const Agenda& doit_mono_agenda,
2804  const Index& doit_is_initialized,
2805  const Verbosity& verbosity)
2806 
2807 {
2808  CREATE_OUT2;
2809 
2810  if (!cloudbox_on) {
2811  CREATE_OUT0;
2812  out0 << " Cloudbox is off, DOIT calculation will be skipped.\n";
2813  return;
2814  //throw runtime_error( "Cloudbox is off, no scattering calculations to be"
2815  // "performed." );
2816  }
2817 
2818  //-------- Check input -------------------------------------------
2819 
2820  if (atmfields_checked != 1)
2821  throw runtime_error(
2822  "The atmospheric fields must be flagged to have "
2823  "passed a consistency check (atmfields_checked=1).");
2824  if (atmgeom_checked != 1)
2825  throw runtime_error(
2826  "The atmospheric geometry must be flagged to have "
2827  "passed a consistency check (atmgeom_checked=1).");
2828  if (cloudbox_checked != 1)
2829  throw runtime_error(
2830  "The cloudbox must be flagged to have "
2831  "passed a consistency check (cloudbox_checked=1).");
2832 
2833  // Don't do anything if there's no cloudbox defined.
2834  if (!cloudbox_on) return;
2835 
2836  if (scat_data_checked != 1)
2837  throw runtime_error(
2838  "The scattering data must be flagged to have "
2839  "passed a consistency check (scat_data_checked=1).");
2840 
2841  chk_not_empty("doit_mono_agenda", doit_mono_agenda);
2842 
2843  // Frequency grid
2844  //
2845  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
2846  chk_if_increasing("f_grid", f_grid);
2847 
2848  // Check whether DoitInit was executed
2849  if (!doit_is_initialized) {
2850  ostringstream os;
2851  os << "Initialization method *DoitInit* has to be called before "
2852  << "*DoitCalc*";
2853  throw runtime_error(os.str());
2854  }
2855 
2856  //-------- end of checks ----------------------------------------
2857 
2858  // We have to make a local copy of the Workspace and the agendas because
2859  // only non-reference types can be declared firstprivate in OpenMP
2860  Workspace l_ws(ws);
2861  Agenda l_doit_mono_agenda(doit_mono_agenda);
2862 
2863  // OMP likes simple loop end conditions, so we make a local copy here:
2864  const Index nf = f_grid.nelem();
2865 
2866  if (nf) {
2867  String fail_msg;
2868  bool failed = false;
2869 
2870 #pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) \
2871  firstprivate(l_ws, l_doit_mono_agenda)
2872  for (Index f_index = 0; f_index < nf; f_index++) {
2873  if (failed) {
2874  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) = NAN;
2875  continue;
2876  }
2877 
2878  try {
2879  ostringstream os;
2880  os << "Frequency: " << f_grid[f_index] / 1e9 << " GHz \n";
2881  out2 << os.str();
2882 
2883  Tensor6 cloudbox_field_mono_local =
2884  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker);
2885  doit_mono_agendaExecute(l_ws,
2886  cloudbox_field_mono_local,
2887  f_grid,
2888  f_index,
2889  l_doit_mono_agenda);
2890  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) =
2891  cloudbox_field_mono_local;
2892  } catch (const std::exception& e) {
2893  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) = NAN;
2894  ostringstream os;
2895  os << "Error for f_index = " << f_index << " (" << f_grid[f_index]
2896  << " Hz)" << endl
2897  << e.what();
2898 #pragma omp critical(DoitCalc_fail)
2899  {
2900  failed = true;
2901  fail_msg = os.str();
2902  }
2903  continue;
2904  }
2905  }
2906 
2907  if (failed) throw runtime_error(fail_msg);
2908  }
2909 }
2910 
2911 /* Workspace method: Doxygen documentation will be auto-generated */
2912 void DoitGetIncoming(Workspace& ws,
2913  Tensor7& cloudbox_field,
2914  const Index& atmfields_checked,
2915  const Index& atmgeom_checked,
2916  const Index& cloudbox_checked,
2917  const Index& doit_is_initialized,
2918  const Agenda& iy_main_agenda,
2919  const Index& atmosphere_dim,
2920  const Vector& lat_grid,
2921  const Vector& lon_grid,
2922  const Tensor3& z_field,
2923  const EnergyLevelMap& nlte_field,
2924  const Index& cloudbox_on,
2925  const ArrayOfIndex& cloudbox_limits,
2926  const Vector& f_grid,
2927  const Index& stokes_dim,
2928  const Vector& za_grid,
2929  const Vector& aa_grid,
2930  const Index& rigorous,
2931  const Numeric& maxratio,
2932  const Verbosity&) {
2933  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2934  if (atmfields_checked != 1)
2935  throw runtime_error(
2936  "The atmospheric fields must be flagged to have "
2937  "passed a consistency check (atmfields_checked=1).");
2938  if (atmgeom_checked != 1)
2939  throw runtime_error(
2940  "The atmospheric geometry must be flagged to have "
2941  "passed a consistency check (atmgeom_checked=1).");
2942  if (cloudbox_checked != 1)
2943  throw runtime_error(
2944  "The cloudbox must be flagged to have "
2945  "passed a consistency check (cloudbox_checked=1).");
2946 
2947  // Don't do anything if there's no cloudbox defined.
2948  if (!cloudbox_on) return;
2949 
2950  // Check whether DoitInit was executed
2951  if (!doit_is_initialized) {
2952  ostringstream os;
2953  os << "Initialization method *DoitInit* has to be called before "
2954  << "*DoitGetIncoming*.";
2955  throw runtime_error(os.str());
2956  }
2957 
2958  // iy_unit hard.coded to "1" here
2959  const String iy_unit = "1";
2960 
2961  Index Nf = f_grid.nelem();
2962  Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2963  Index Nza = za_grid.nelem();
2964  Matrix iy;
2965  Ppath ppath;
2966 
2967  //--- Check input ----------------------------------------------------------
2968  if (!(atmosphere_dim == 1 || atmosphere_dim == 3))
2969  throw runtime_error("The atmospheric dimensionality must be 1 or 3.");
2970  if (za_grid[0] != 0. || za_grid[Nza - 1] != 180.)
2971  throw runtime_error(
2972  "*za_grid* must include 0 and 180 degrees as endpoints.");
2973  //--------------------------------------------------------------------------
2974 
2975  if (atmosphere_dim == 1) {
2976  //Define the variables for position and direction.
2977  Vector los(1), pos(1);
2978 
2979  //--- Get cloudbox_field at lower and upper boundary
2980  // (boundary=0: lower, boundary=1: upper)
2981  for (Index boundary = 0; boundary <= 1; boundary++) {
2982  const Index boundary_index =
2983  boundary ? cloudbox_field.nvitrines() - 1 : 0;
2984  pos[0] = z_field(cloudbox_limits[boundary], 0, 0);
2985 
2986  // doing the first angle separately for allowing dy between 2 angles
2987  // in the loop
2988  los[0] = za_grid[0];
2989  get_iy(ws,
2990  iy,
2991  0,
2992  f_grid,
2993  nlte_field,
2994  pos,
2995  los,
2996  Vector(0),
2997  iy_unit,
2998  iy_main_agenda);
2999  cloudbox_field(joker, boundary_index, 0, 0, 0, 0, joker) = iy;
3000 
3001  for (Index za_index = 1; za_index < Nza; za_index++) {
3002  los[0] = za_grid[za_index];
3003 
3004  get_iy(ws,
3005  iy,
3006  0,
3007  f_grid,
3008  nlte_field,
3009  pos,
3010  los,
3011  Vector(0),
3012  iy_unit,
3013  iy_main_agenda);
3014 
3015  cloudbox_field(joker, boundary_index, 0, 0, za_index, 0, joker) = iy;
3016 
3017  if (rigorous) {
3018  for (Index fi = 0; fi < Nf; fi++) {
3019  if (cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0) /
3020  cloudbox_field(
3021  fi, boundary_index, 0, 0, za_index, 0, 0) >
3022  maxratio ||
3023  cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0) /
3024  cloudbox_field(
3025  fi, boundary_index, 0, 0, za_index, 0, 0) <
3026  1 / maxratio) {
3027  ostringstream os;
3028  os << "ERROR: Radiance difference between "
3029  << "interpolation points is too large (factor " << maxratio
3030  << ")\n"
3031  << "to safely interpolate. This might be due to "
3032  << "za_grid being too coarse or the radiance field "
3033  << "being a step-like function.\n"
3034  << "Happens at boundary " << boundary_index
3035  << " between zenith angles " << za_grid[za_index - 1]
3036  << " and " << za_grid[za_index] << "deg\n"
3037  << "for frequency #" << fi << ", where radiances are "
3038  << cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0)
3039  << " and "
3040  << cloudbox_field(fi, boundary_index, 0, 0, za_index, 0, 0)
3041  << " W/(sr m2 Hz).";
3042  throw runtime_error(os.str());
3043  }
3044  }
3045  }
3046  }
3047  }
3048  }
3049 
3050  //--- atmosphere_dim = 3: --------------------------------------------------
3051  else {
3052  Index Naa = aa_grid.nelem();
3053 
3054  if (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.)
3055  throw runtime_error(
3056  "*aa_grid* must include 0 and 360 degrees as endpoints.");
3057 
3058  Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3059  Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3060 
3061  // Convert aa_grid to "sensor coordinates"
3062  // (-180° < azimuth angle < 180°)
3063  //
3064  Vector aa_g(Naa);
3065  for (Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3066 
3067  // Define the variables for position and direction.
3068  Vector los(2), pos(3);
3069 
3070  //--- Get cloudbox_field at lower and upper boundary
3071  // (boundary=0: lower, boundary=1: upper)
3072  for (Index boundary = 0; boundary <= 1; boundary++) {
3073  const Index boundary_index =
3074  boundary ? cloudbox_field.nvitrines() - 1 : 0;
3075  for (Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3076  for (Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3077  pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3078  pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3079  pos[0] = z_field(cloudbox_limits[boundary],
3080  lat_index + cloudbox_limits[2],
3081  lon_index + cloudbox_limits[4]);
3082 
3083  for (Index za_index = 0; za_index < Nza; za_index++) {
3084  for (Index aa_index = 0; aa_index < Naa; aa_index++) {
3085  los[0] = za_grid[za_index];
3086  los[1] = aa_g[aa_index];
3087 
3088  // For end points of za_index (0 & 180deg), we
3089  // only need to perform calculations for one scat_aa
3090  // and set the others to same value
3091  if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3092  get_iy(ws,
3093  iy,
3094  0,
3095  f_grid,
3096  nlte_field,
3097  pos,
3098  los,
3099  Vector(0),
3100  iy_unit,
3101  iy_main_agenda);
3102  }
3103 
3104  cloudbox_field(joker,
3105  boundary_index,
3106  lat_index,
3107  lon_index,
3108  za_index,
3109  aa_index,
3110  joker) = iy;
3111  }
3112  }
3113  }
3114  }
3115  }
3116 
3117  //--- Get scat_i_lat (2nd and 3rd boundary)
3118  for (Index boundary = 0; boundary <= 1; boundary++) {
3119  const Index boundary_index = boundary ? cloudbox_field.nshelves() - 1 : 0;
3120  for (Index p_index = 0; p_index < Np_cloud; p_index++) {
3121  for (Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3122  pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3123  pos[1] = lat_grid[cloudbox_limits[boundary + 2]];
3124  pos[0] = z_field(p_index + cloudbox_limits[0],
3125  cloudbox_limits[boundary + 2],
3126  lon_index + cloudbox_limits[4]);
3127 
3128  for (Index za_index = 0; za_index < Nza; za_index++) {
3129  for (Index aa_index = 0; aa_index < Naa; aa_index++) {
3130  los[0] = za_grid[za_index];
3131  los[1] = aa_g[aa_index];
3132 
3133  // For end points of za_index, we need only to
3134  // perform calculations for first scat_aa
3135  if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3136  get_iy(ws,
3137  iy,
3138  0,
3139  f_grid,
3140  nlte_field,
3141  pos,
3142  los,
3143  Vector(0),
3144  iy_unit,
3145  iy_main_agenda);
3146  }
3147 
3148  cloudbox_field(joker,
3149  p_index,
3150  boundary_index,
3151  lon_index,
3152  za_index,
3153  aa_index,
3154  joker) = iy;
3155  }
3156  }
3157  }
3158  }
3159  }
3160 
3161  //--- Get scat_i_lon (1st and 2nd boundary):
3162  for (Index boundary = 0; boundary <= 1; boundary++) {
3163  const Index boundary_index = boundary ? cloudbox_field.nbooks() - 1 : 0;
3164  for (Index p_index = 0; p_index < Np_cloud; p_index++) {
3165  for (Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3166  pos[2] = lon_grid[cloudbox_limits[boundary + 4]];
3167  pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3168  pos[0] = z_field(p_index + cloudbox_limits[0],
3169  lat_index + cloudbox_limits[2],
3170  cloudbox_limits[boundary + 4]);
3171 
3172  for (Index za_index = 0; za_index < Nza; za_index++) {
3173  for (Index aa_index = 0; aa_index < Naa; aa_index++) {
3174  los[0] = za_grid[za_index];
3175  los[1] = aa_g[aa_index];
3176 
3177  // For end points of za_index, we need only to
3178  // perform calculations for first scat_aa
3179  if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3180  get_iy(ws,
3181  iy,
3182  0,
3183  f_grid,
3184  nlte_field,
3185  pos,
3186  los,
3187  Vector(0),
3188  iy_unit,
3189  iy_main_agenda);
3190  }
3191 
3192  cloudbox_field(joker,
3193  p_index,
3194  lat_index,
3195  boundary_index,
3196  za_index,
3197  aa_index,
3198  joker) = iy;
3199  }
3200  }
3201  }
3202  }
3203  }
3204  } // End atmosphere_dim = 3.
3205 }
3206 
3207 /* Workspace method: Doxygen documentation will be auto-generated */
3208 void DoitGetIncoming1DAtm(Workspace& ws,
3209  Tensor7& cloudbox_field,
3210  Index& cloudbox_on,
3211  const Index& atmfields_checked,
3212  const Index& atmgeom_checked,
3213  const Index& cloudbox_checked,
3214  const Index& doit_is_initialized,
3215  const Agenda& iy_main_agenda,
3216  const Index& atmosphere_dim,
3217  const Vector& lat_grid,
3218  const Vector& lon_grid,
3219  const Tensor3& z_field,
3220  const EnergyLevelMap& nlte_field,
3221  const ArrayOfIndex& cloudbox_limits,
3222  const Vector& f_grid,
3223  const Index& stokes_dim,
3224  const Vector& za_grid,
3225  const Vector& aa_grid,
3226  const Verbosity&) {
3227  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
3228  if (atmfields_checked != 1)
3229  throw runtime_error(
3230  "The atmospheric fields must be flagged to have "
3231  "passed a consistency check (atmfields_checked=1).");
3232  if (atmgeom_checked != 1)
3233  throw runtime_error(
3234  "The atmospheric geometry must be flagged to have "
3235  "passed a consistency check (atmgeom_checked=1).");
3236  if (cloudbox_checked != 1)
3237  throw runtime_error(
3238  "The cloudbox must be flagged to have "
3239  "passed a consistency check (cloudbox_checked=1).");
3240 
3241  // Don't do anything if there's no cloudbox defined.
3242  if (!cloudbox_on) return;
3243 
3244  // Check whether DoitInit was executed
3245  if (!doit_is_initialized) {
3246  ostringstream os;
3247  os << "Initialization method *DoitInit* has to be called before "
3248  << "*DoitGetIncoming1DAtm*.";
3249  throw runtime_error(os.str());
3250  }
3251 
3252  // iy_unit hard.coded to "1" here
3253  const String iy_unit = "1";
3254 
3255  Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3256  Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3257  Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3258  Index Nza = za_grid.nelem();
3259  Index Naa = aa_grid.nelem();
3260  Matrix iy;
3261  Ppath ppath;
3262 
3263  //--- Check input ----------------------------------------------------------
3264  if (atmosphere_dim != 3)
3265  throw runtime_error("The atmospheric dimensionality must be 3.");
3266  if (za_grid[0] != 0. || za_grid[Nza - 1] != 180.)
3267  throw runtime_error(
3268  "*za_grid* must include 0 and 180 degrees as endpoints.");
3269  if (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.)
3270  throw runtime_error(
3271  "*aa_grid* must include 0 and 360 degrees as endpoints.");
3272  //--------------------------------------------------------------------------
3273 
3274  // Dummy variable for flag cloudbox_on. It has to be 0 here not to get
3275  // stuck in an infinite loop (if some propagation path hits the cloud
3276  // box at some other position.
3277  cloudbox_on = 0;
3278 
3279  // Convert za_grid to "sensor coordinates"
3280  //(-180° < azimuth angle < 180°)
3281  //
3282  Vector aa_g(Naa);
3283  for (Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3284 
3285  // As the atmosphere is spherically symmetric we only have to calculate
3286  // one azimuth angle.
3287  Index aa_index = 0;
3288 
3289  // Define the variables for position and direction.
3290  Vector los(2), pos(3);
3291 
3292  // These variables are constant for all calculations:
3293  pos[1] = lat_grid[cloudbox_limits[2]];
3294  pos[2] = lon_grid[cloudbox_limits[4]];
3295  los[1] = aa_g[aa_index];
3296 
3297  // Calculate cloudbox_field on boundary
3298  for (Index p_index = 0; p_index < Np_cloud; p_index++) {
3299  pos[0] = z_field(
3300  cloudbox_limits[0] + p_index, cloudbox_limits[2], cloudbox_limits[4]);
3301 
3302  for (Index za_index = 0; za_index < Nza; za_index++) {
3303  los[0] = za_grid[za_index];
3304 
3305  get_iy(ws,
3306  iy,
3307  0,
3308  f_grid,
3309  nlte_field,
3310  pos,
3311  los,
3312  Vector(0),
3313  iy_unit,
3314  iy_main_agenda);
3315 
3316  for (Index aa = 0; aa < Naa; aa++) {
3317  // cloudbox_field lower boundary
3318  if (p_index == 0) {
3319  for (Index lat = 0; lat < Nlat_cloud; lat++) {
3320  for (Index lon = 0; lon < Nlon_cloud; lon++) {
3321  cloudbox_field(joker, 0, lat, lon, za_index, aa, joker) = iy;
3322  }
3323  }
3324  }
3325  //cloudbox_field at upper boundary
3326  else if (p_index == Np_cloud - 1)
3327  for (Index lat = 0; lat < Nlat_cloud; lat++) {
3328  for (Index lon = 0; lon < Nlon_cloud; lon++) {
3329  cloudbox_field(joker,
3330  cloudbox_field.nvitrines(),
3331  lat,
3332  lon,
3333  za_index,
3334  aa,
3335  joker) = iy;
3336  }
3337  }
3338 
3339  // scat_i_lat (both boundaries)
3340  for (Index lat = 0; lat < 2; lat++) {
3341  for (Index lon = 0; lon < Nlon_cloud; lon++) {
3342  const Index boundary_index =
3343  lat ? cloudbox_field.nshelves() - 1 : 0;
3344  cloudbox_field(
3345  joker, p_index, boundary_index, lon, za_index, aa, joker) = iy;
3346  }
3347  }
3348 
3349  // scat_i_lon (both boundaries)
3350  for (Index lat = 0; lat < Nlat_cloud; lat++) {
3351  for (Index lon = 0; lon < 2; lon++) {
3352  const Index boundary_index = lon ? cloudbox_field.nbooks() - 1 : 0;
3353  cloudbox_field(
3354  joker, p_index, lat, boundary_index, za_index, aa, joker) = iy;
3355  }
3356  }
3357  }
3358  }
3359  }
3360  cloudbox_on = 1;
3361 }
3362 
3363 /* Workspace method: Doxygen documentation will be auto-generated */
3364 void cloudbox_fieldSetFromPrecalc(Tensor7& cloudbox_field,
3365  const Vector& za_grid,
3366  const Vector& f_grid,
3367  const Index& atmosphere_dim,
3368  const Index& stokes_dim,
3369  const ArrayOfIndex& cloudbox_limits,
3370  const Index& doit_is_initialized,
3371  const Tensor7& cloudbox_field_precalc,
3372  const Verbosity&) //verbosity)
3373 {
3374  // this is only for 1D atmo!
3375  if (atmosphere_dim != 1) {
3376  ostringstream os;
3377  os << "This method is currently only implemented for 1D atmospheres!\n";
3378  throw runtime_error(os.str());
3379  }
3380 
3381  // Check whether DoitInit was executed
3382  if (!doit_is_initialized) {
3383  ostringstream os;
3384  os << "Initialization method *DoitInit* has to be called before "
3385  << "*cloudbox_fieldSetFromPrecalc*";
3386  throw runtime_error(os.str());
3387  }
3388 
3389  // check dimensions of cloudbox_field_precalc
3390  Index nf = f_grid.nelem();
3391  Index nza = za_grid.nelem();
3392  Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3393 
3394  if (nf != cloudbox_field_precalc.nlibraries()) {
3395  ostringstream os;
3396  os << "cloudbox_field_precalc has wrong size in frequency dimension.\n"
3397  << nf << " frequency points are expected, but cloudbox_field_precalc "
3398  << "contains " << cloudbox_field_precalc.nlibraries()
3399  << "frequency points.\n";
3400  throw runtime_error(os.str());
3401  }
3402  if (np != cloudbox_field_precalc.nvitrines()) {
3403  ostringstream os;
3404  os << "cloudbox_field_precalc has wrong size in pressure level dimension.\n"
3405  << np << " pressure levels expected, but cloudbox_field_precalc "
3406  << "contains " << cloudbox_field_precalc.nvitrines()
3407  << "pressure levels.\n";
3408  throw runtime_error(os.str());
3409  }
3410  if (nza != cloudbox_field_precalc.npages()) {
3411  ostringstream os;
3412  os << "cloudbox_field_precalc has wrong size in polar angle dimension.\n"
3413  << nza << " angles expected, but cloudbox_field_precalc "
3414  << "contains " << cloudbox_field_precalc.npages() << "angles.\n";
3415  throw runtime_error(os.str());
3416  }
3417  if (stokes_dim != cloudbox_field_precalc.ncols()) {
3418  ostringstream os;
3419  os << "cloudbox_field_precalc has wrong stokes dimension.\n"
3420  << "Dimension " << stokes_dim
3421  << " expected, but cloudbox_field_precalc is dimesnion "
3422  << cloudbox_field_precalc.ncols() << ".\n";
3423  throw runtime_error(os.str());
3424  }
3425 
3426  // We have to update cloudbox incoming (!) field - this is because
3427  // compared to our first guess initialization, the clearsky field around might
3428  // have changed.
3429 
3430  // Find which za_grid entries describe upwelling, which downwelling
3431  // radiation.
3432  Index first_upwell = 0;
3433  assert(za_grid[0] < 90.);
3434  assert(za_grid[za_grid.nelem() - 1] > 90.);
3435  while (za_grid[first_upwell] < 90.) first_upwell++;
3436 
3437  Range downwell(0, first_upwell);
3438  Range upwell(first_upwell, za_grid.nelem() - first_upwell);
3439 
3440  // Copy everything inside the field
3441  cloudbox_field(joker, Range(1, np - 2), 0, 0, joker, 0, joker) =
3442  cloudbox_field_precalc(joker, Range(1, np - 2), 0, 0, joker, 0, joker);
3443 
3444  // At boundaries we need to be a bit careful. We shouldn't overwrite the
3445  // boundary conditions.
3446 
3447  // Copy only upwelling radiation at upper boundary from precalc field
3448  // (Downwelling is "boundary condition" and has been set by DoitGetIncoming)
3449  cloudbox_field(joker, np - 1, 0, 0, upwell, 0, joker) =
3450  cloudbox_field_precalc(joker, np - 1, 0, 0, upwell, 0, joker);
3451 
3452  if (cloudbox_limits[0] != 0)
3453  // Copy only downwelling radiation at lower boundary from precalc field
3454  // (Upwelling is "boundary condition" and has been set by DoitGetIncoming)
3455  cloudbox_field(joker, 0, 0, 0, downwell, 0, joker) =
3456  cloudbox_field_precalc(joker, 0, 0, 0, downwell, 0, joker);
3457  else
3458  // Copy all directions at lower boundary from precalc field
3459  // (when lower boundary at surface, the upwelling field is fixed "boundary
3460  // condition", but part of the field getting iteratively updated according
3461  // to the cloudbox-dependent surface reflection contribution)
3462  cloudbox_field(joker, 0, 0, 0, joker, 0, joker) =
3463  cloudbox_field_precalc(joker, 0, 0, 0, joker, 0, joker);
3464 }
3465 
3466 /* Workspace method: Doxygen documentation will be auto-generated */
3467 void cloudbox_fieldSetClearsky(Tensor7& cloudbox_field,
3468  const Vector& p_grid,
3469  const Vector& lat_grid,
3470  const Vector& lon_grid,
3471  const ArrayOfIndex& cloudbox_limits,
3472  const Index& atmosphere_dim,
3473  const Index& cloudbox_on,
3474  const Index& doit_is_initialized,
3475  const Index& all_frequencies,
3476  const Verbosity& verbosity) {
3477  CREATE_OUT2;
3478 
3479  // Don't do anything if there's no cloudbox defined.
3480  if (!cloudbox_on) return;
3481 
3482  // Check whether DoitInit was executed
3483  if (!doit_is_initialized) {
3484  ostringstream os;
3485  os << "Initialization method *DoitInit* has to be called before "
3486  << "*cloudbox_fieldSetClearsky*.";
3487  throw runtime_error(os.str());
3488  }
3489 
3490  out2
3491  << " Interpolate boundary clearsky field to obtain the initial field.\n";
3492 
3493  // Initial field only needs to be calculated from clearsky field for the
3494  // first frequency. For the next frequencies the solution field from the
3495  // previous frequencies is used.
3496  if (atmosphere_dim == 1) {
3497  const Index nf = all_frequencies ? cloudbox_field.nlibraries() : 1;
3498 
3499  for (Index f_index = 0; f_index < nf; f_index++) {
3500  Index N_za = cloudbox_field.npages();
3501  Index N_aa = cloudbox_field.nrows();
3502  Index N_i = cloudbox_field.ncols();
3503 
3504  //1. interpolation - pressure grid
3505 
3506  /*the old grid is having only two elements, corresponding to the
3507  cloudbox_limits and the new grid have elements corresponding to
3508  all grid points inside the cloudbox plus the cloud_box_limits*/
3509 
3510  ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3511 
3512  p2gridpos(p_gp,
3513  p_grid[Range(cloudbox_limits[0],
3514  2,
3515  (cloudbox_limits[1] - cloudbox_limits[0]))],
3516  p_grid[Range(cloudbox_limits[0],
3517  (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3518 
3519  Matrix itw((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3520  interpweights(itw, p_gp);
3521 
3522  Tensor6 scat_i_p(2, 1, 1, N_za, 1, N_i);
3523  scat_i_p(0, joker, joker, joker, joker, joker) =
3524  cloudbox_field(f_index, 0, joker, joker, joker, joker, joker);
3525  scat_i_p(1, joker, joker, joker, joker, joker) =
3526  cloudbox_field(f_index,
3527  cloudbox_field.nvitrines() - 1,
3528  joker,
3529  joker,
3530  joker,
3531  joker,
3532  joker);
3533 
3534  for (Index za_index = 0; za_index < N_za; ++za_index) {
3535  for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3536  for (Index i = 0; i < N_i; ++i) {
3537  VectorView target_field = cloudbox_field(
3538  f_index, Range(joker), 0, 0, za_index, aa_index, i);
3539 
3540  ConstVectorView source_field =
3541  scat_i_p(Range(joker), 0, 0, za_index, aa_index, i);
3542 
3543  interp(target_field, itw, source_field, p_gp);
3544  }
3545  }
3546  }
3547  }
3548  } else if (atmosphere_dim == 3) {
3549  if (all_frequencies == false)
3550  throw runtime_error(
3551  "Error in cloudbox_fieldSetClearsky: For 3D "
3552  "all_frequencies option is not implemented \n");
3553 
3554  for (Index f_index = 0; f_index < cloudbox_field.nvitrines(); f_index++) {
3555  Index N_p = cloudbox_field.nvitrines();
3556  Index N_lat = cloudbox_field.nshelves();
3557  Index N_lon = cloudbox_field.nbooks();
3558  Index N_za = cloudbox_field.npages();
3559  Index N_aa = cloudbox_field.nrows();
3560  Index N_i = cloudbox_field.ncols();
3561 
3562  Tensor6 scat_i_p(2, N_lat, N_lon, N_za, N_aa, N_i);
3563  scat_i_p(0, joker, joker, joker, joker, joker) =
3564  cloudbox_field(f_index, 0, joker, joker, joker, joker, joker);
3565  scat_i_p(1, joker, joker, joker, joker, joker) =
3566  cloudbox_field(f_index, N_p - 1, joker, joker, joker, joker, joker);
3567 
3568  Tensor6 scat_i_lat(N_p, 2, N_lon, N_za, N_aa, N_i);
3569  scat_i_lat(joker, 0, joker, joker, joker, joker) =
3570  cloudbox_field(f_index, joker, 0, joker, joker, joker, joker);
3571  scat_i_lat(joker, 1, joker, joker, joker, joker) =
3572  cloudbox_field(f_index, joker, N_lat - 1, joker, joker, joker, joker);
3573 
3574  Tensor6 scat_i_lon(N_p, N_lat, 2, N_za, N_aa, N_i);
3575  scat_i_lon(joker, joker, 0, joker, joker, joker) =
3576  cloudbox_field(f_index, joker, joker, 0, joker, joker, joker);
3577  scat_i_lon(joker, joker, 1, joker, joker, joker) =
3578  cloudbox_field(f_index, joker, joker, N_lon - 1, joker, joker, joker);
3579 
3580  //1. interpolation - pressure grid, latitude grid and longitude grid
3581 
3582  ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3583  ArrayOfGridPos lat_gp((cloudbox_limits[3] - cloudbox_limits[2]) + 1);
3584  ArrayOfGridPos lon_gp((cloudbox_limits[5] - cloudbox_limits[4]) + 1);
3585 
3586  /*the old grid is having only two elements, corresponding to the
3587  cloudbox_limits and the new grid have elements corresponding to
3588  all grid points inside the cloudbox plus the cloud_box_limits*/
3589 
3590  p2gridpos(p_gp,
3591  p_grid[Range(cloudbox_limits[0],
3592  2,
3593  (cloudbox_limits[1] - cloudbox_limits[0]))],
3594  p_grid[Range(cloudbox_limits[0],
3595  (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3596  gridpos(lat_gp,
3597  lat_grid[Range(cloudbox_limits[2],
3598  2,
3599  (cloudbox_limits[3] - cloudbox_limits[2]))],
3600  lat_grid[Range(cloudbox_limits[2],
3601  (cloudbox_limits[3] - cloudbox_limits[2]) + 1)]);
3602  gridpos(lon_gp,
3603  lon_grid[Range(cloudbox_limits[4],
3604  2,
3605  (cloudbox_limits[5] - cloudbox_limits[4]))],
3606  lon_grid[Range(cloudbox_limits[4],
3607  (cloudbox_limits[5] - cloudbox_limits[4]) + 1)]);
3608 
3609  //interpolation weights corresponding to pressure, latitude and
3610  //longitude grids.
3611 
3612  Matrix itw_p((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3613  Matrix itw_lat((cloudbox_limits[3] - cloudbox_limits[2]) + 1, 2);
3614  Matrix itw_lon((cloudbox_limits[5] - cloudbox_limits[4]) + 1, 2);
3615 
3616  interpweights(itw_p, p_gp);
3617  interpweights(itw_lat, lat_gp);
3618  interpweights(itw_lon, lon_gp);
3619 
3620  // interpolation - pressure grid
3621  for (Index lat_index = 0;
3622  lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3623  ++lat_index) {
3624  for (Index lon_index = 0;
3625  lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3626  ++lon_index) {
3627  for (Index za_index = 0; za_index < N_za; ++za_index) {
3628  for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3629  for (Index i = 0; i < N_i; ++i) {
3630  VectorView target_field = cloudbox_field(f_index,
3631  Range(joker),
3632  lat_index,
3633  lon_index,
3634  za_index,
3635  aa_index,
3636  i);
3637 
3638  ConstVectorView source_field = scat_i_p(
3639  Range(joker), lat_index, lon_index, za_index, aa_index, i);
3640 
3641  interp(target_field, itw_p, source_field, p_gp);
3642  }
3643  }
3644  }
3645  }
3646  }
3647  //interpolation latitude
3648  for (Index p_index = 0;
3649  p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3650  ++p_index) {
3651  for (Index lon_index = 0;
3652  lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3653  ++lon_index) {
3654  for (Index za_index = 0; za_index < N_za; ++za_index) {
3655  for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3656  for (Index i = 0; i < N_i; ++i) {
3657  VectorView target_field = cloudbox_field(f_index,
3658  p_index,
3659  Range(joker),
3660  lon_index,
3661  za_index,
3662  aa_index,
3663  i);
3664 
3665  ConstVectorView source_field = scat_i_lat(
3666  p_index, Range(joker), lon_index, za_index, aa_index, i);
3667 
3668  interp(target_field, itw_lat, source_field, lat_gp);
3669  }
3670  }
3671  }
3672  }
3673  }
3674  //interpolation -longitude
3675  for (Index p_index = 0;
3676  p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3677  ++p_index) {
3678  for (Index lat_index = 0;
3679  lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3680  ++lat_index) {
3681  for (Index za_index = 0; za_index < N_za; ++za_index) {
3682  for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3683  for (Index i = 0; i < N_i; ++i) {
3684  VectorView target_field = cloudbox_field(f_index,
3685  p_index,
3686  lat_index,
3687  Range(joker),
3688  za_index,
3689  aa_index,
3690  i);
3691 
3692  ConstVectorView source_field = scat_i_lon(
3693  p_index, lat_index, Range(joker), za_index, aa_index, i);
3694 
3695  interp(target_field, itw_lon, source_field, lon_gp);
3696  }
3697  }
3698  }
3699  }
3700  } //end of interpolation
3701  } // end of frequency loop
3702  } //ends atmosphere_dim = 3
3703 }
3704 
3705 /* Workspace method: Doxygen documentation will be auto-generated */
3706 void cloudbox_fieldSetConst( //WS Output:
3707  Tensor7& cloudbox_field,
3708  //WS Input:
3709  const Vector& p_grid,
3710  const Vector& lat_grid,
3711  const Vector& lon_grid,
3712  const ArrayOfIndex& cloudbox_limits,
3713  const Index& atmosphere_dim,
3714  const Index& stokes_dim,
3715  // Keyword
3716  const Vector& cloudbox_field_values,
3717  const Verbosity& verbosity) {
3718  CREATE_OUT2;
3719 
3720  Tensor6 cloudbox_field_mono(cloudbox_field.nvitrines(),
3721  cloudbox_field.nshelves(),
3722  cloudbox_field.nbooks(),
3723  cloudbox_field.npages(),
3724  cloudbox_field.nrows(),
3725  cloudbox_field.ncols());
3726 
3727  for (Index f_index = 0; f_index < cloudbox_field.nlibraries(); f_index++) {
3728  cloudbox_field_mono =
3729  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker);
3730 
3731  cloudbox_field_monoSetConst(cloudbox_field_mono,
3732  p_grid,
3733  lat_grid,
3734  lon_grid,
3735  cloudbox_limits,
3736  atmosphere_dim,
3737  stokes_dim,
3738  cloudbox_field_values,
3739  verbosity);
3740 
3741  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) =
3742  cloudbox_field_mono;
3743  }
3744 }
3745 
3746 /* Workspace method: Doxygen documentation will be auto-generated */
3747 void cloudbox_fieldSetConstPerFreq( //WS Output:
3748  Tensor7& cloudbox_field,
3749  //WS Input:
3750  const Vector& p_grid,
3751  const Vector& lat_grid,
3752  const Vector& lon_grid,
3753  const ArrayOfIndex& cloudbox_limits,
3754  const Index& atmosphere_dim,
3755  const Index& stokes_dim,
3756  // Keyword
3757  const Matrix& cloudbox_field_values,
3758  const Verbosity& verbosity) {
3759  CREATE_OUT2;
3760 
3761  if (cloudbox_field.nlibraries() != cloudbox_field_values.nrows())
3762  throw runtime_error(
3763  "Number of rows in *cloudbox_field_values* has to be equal"
3764  " the frequency dimension of *cloudbox_field*.");
3765 
3766  Tensor6 cloudbox_field_mono(cloudbox_field.nvitrines(),
3767  cloudbox_field.nshelves(),
3768  cloudbox_field.nbooks(),
3769  cloudbox_field.npages(),
3770  cloudbox_field.nrows(),
3771  cloudbox_field.ncols());
3772 
3773  for (Index f_index = 0; f_index < cloudbox_field.nlibraries(); f_index++) {
3774  cloudbox_field_mono =
3775  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker);
3776 
3777  cloudbox_field_monoSetConst(cloudbox_field_mono,
3778  p_grid,
3779  lat_grid,
3780  lon_grid,
3781  cloudbox_limits,
3782  atmosphere_dim,
3783  stokes_dim,
3784  cloudbox_field_values(f_index, joker),
3785  verbosity);
3786 
3787  cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) =
3788  cloudbox_field_mono;
3789  }
3790 }
3791 
3792 /* Workspace method: Doxygen documentation will be auto-generated */
3793 void cloudbox_field_monoSetConst( //WS Output:
3794  Tensor6& cloudbox_field_mono,
3795  //WS Input:
3796  const Vector& p_grid,
3797  const Vector& lat_grid,
3798  const Vector& lon_grid,
3799  const ArrayOfIndex& cloudbox_limits,
3800  const Index& atmosphere_dim,
3801  const Index& stokes_dim,
3802  // Keyword
3803  const Vector& cloudbox_field_values,
3804  const Verbosity& verbosity) {
3805  CREATE_OUT2;
3806  CREATE_OUT3;
3807 
3808  out2 << " Set initial field to constant values: " << cloudbox_field_values
3809  << "\n";
3810 
3811  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3812 
3813  // Grids have to be adapted to atmosphere_dim.
3814  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3815 
3816  // Check the input:
3817  if (stokes_dim < 0 || stokes_dim > 4)
3818  throw runtime_error(
3819  "The dimension of stokes vector must be"
3820  "1,2,3, or 4.");
3821 
3822  if (stokes_dim != cloudbox_field_values.nelem())
3823  throw runtime_error(
3824  "Length of *cloudbox_field_values* has to be equal"
3825  " *stokes_dim*.");
3826  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
3827  throw runtime_error(
3828  "*cloudbox_limits* is a vector which contains the"
3829  "upper and lower limit of the cloud for all "
3830  "atmospheric dimensions. So its dimension must"
3831  "be 2 x *atmosphere_dim*.");
3832 
3833  for (Index i = 0; i < stokes_dim; i++) {
3834  cloudbox_field_mono(joker, joker, joker, joker, joker, i) =
3835  cloudbox_field_values[i];
3836  }
3837 }
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
m_general.h
Template functions for general supergeneric ws methods.
cloudbox_field_monoSetConst
void cloudbox_field_monoSetConst(Tensor6 &cloudbox_field_mono, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &cloudbox_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_field_monoSetConst.
Definition: m_doit.cc:3793
interp_poly
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
Definition: interpolation.cc:2723
VectorView
The VectorView class.
Definition: matpackI.h:610
DOAngularGridsSet
void DOAngularGridsSet(Index &doit_za_grid_size, Vector &aa_grid, Vector &za_grid, const Index &N_za_grid, const Index &N_aa_grid, const String &za_grid_opt_file, const Verbosity &verbosity)
WORKSPACE METHOD: DOAngularGridsSet.
Definition: m_doit.cc:70
doit_mono_agendaExecute
void doit_mono_agendaExecute(Workspace &ws, Tensor6 &cloudbox_field_mono, const Vector &f_grid, const Index f_index, const Agenda &input_agenda)
Definition: auto_md.cc:23717
doit_scat_field_agendaExecute
void doit_scat_field_agendaExecute(Workspace &ws, Tensor6 &doit_scat_field, const Tensor6 &cloudbox_field_mono, const Agenda &input_agenda)
Definition: auto_md.cc:23756
Agenda
The Agenda class.
Definition: agenda_class.h:44
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
chk_if_decreasing
void chk_if_decreasing(const String &x_name, ConstVectorView x)
chk_if_decreasing
Definition: check_input.cc:358
messages.h
Declarations having to do with the four output streams.
doit_conv_flagLsq
void doit_conv_flagLsq(Index &doit_conv_flag, Index &doit_iteration_counter, Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Vector &f_grid, const Index &f_index, const Vector &epsilon, const Index &max_iterations, const Index &throw_nonconv_error, const Verbosity &verbosity)
WORKSPACE METHOD: doit_conv_flagLsq.
Definition: m_doit.cc:370
invrayjean
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
Definition: physics_funcs.cc:192
chk_if_increasing
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:117
AngIntegrate_trapezoid_opti
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:340
Vector
The Vector class.
Definition: matpackI.h:860
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ConstTensor6View::nvitrines
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVI.cc:42
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
RAD2DEG
const Numeric RAD2DEG
DoitGetIncoming1DAtm
void DoitGetIncoming1DAtm(Workspace &ws, Tensor7 &cloudbox_field, Index &cloudbox_on, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &doit_is_initialized, const Agenda &iy_main_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const EnergyLevelMap &nlte_field, const ArrayOfIndex &cloudbox_limits, const Vector &f_grid, const Index &stokes_dim, const Vector &za_grid, const Vector &aa_grid, const Verbosity &)
WORKSPACE METHOD: DoitGetIncoming1DAtm.
Definition: m_doit.cc:3208
is_increasing
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
Range
The range class.
Definition: matpackI.h:160
FILE_TYPE_ASCII
Definition: xml_io.h:38
ArrayOfRetrievalQuantity
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:402
ConstTensor7View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackVII.cc:57
lin_alg.h
Linear algebra functions.
DoitWriteIterationFields
void DoitWriteIterationFields(const Index &doit_iteration_counter, const Tensor6 &cloudbox_field_mono, const Index &f_index, const ArrayOfIndex &iterations, const ArrayOfIndex &frequencies, const Verbosity &verbosity)
WORKSPACE METHOD: DoitWriteIterationFields.
Definition: m_doit.cc:2032
ConstVectorView::empty
bool empty() const
Returns true if variable size is zero.
Definition: matpackI.cc:49
cloudbox_field_ngAcceleration
void cloudbox_field_ngAcceleration(Tensor6 &cloudbox_field_mono, const ArrayOfTensor6 &acceleration_input, const Index &accelerated, const Verbosity &)
Convergence acceleration.
Definition: doit.cc:1620
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
doit_za_interpSet
void doit_za_interpSet(Index &doit_za_interp, const Index &atmosphere_dim, const String &method, const Verbosity &)
WORKSPACE METHOD: doit_za_interpSet.
Definition: m_doit.cc:2770
xml_io.h
This file contains basic functions to handle XML data files.
p2gridpos
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
Definition: special_interp.cc:759
Tensor7
The Tensor7 class.
Definition: matpackVII.h:2382
basic_ostringstream
Definition: sstream.h:171
cloud_ppath_update3D
void cloud_ppath_update3D(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &lat_index, const Index &lon_index, const Index &za_index, const Index &aa_index, ConstVectorView za_grid, ConstVectorView aa_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Index &, const Verbosity &verbosity)
Radiative transfer calculation along a path inside the cloudbox (3D).
Definition: doit.cc:1108
cloudbox_field_monoOptimizeReverse
void cloudbox_field_monoOptimizeReverse(Tensor6 &cloudbox_field_mono, const Vector &p_grid_orig, const Vector &p_grid, const ArrayOfIndex &cloudbox_limits, const Verbosity &)
WORKSPACE METHOD: cloudbox_field_monoOptimizeReverse.
Definition: m_doit.cc:1699
cloudbox_fieldUpdate1D
void cloudbox_fieldUpdate1D(Workspace &ws, Tensor6 &cloudbox_field_mono, const Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Tensor4 &pnd_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Vector &p_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Agenda &surface_rtprop_agenda, const Index &doit_za_interp, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdate1D.
Definition: m_doit.cc:591
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
chk_interpolation_grids
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:989
DoitGetIncoming
void DoitGetIncoming(Workspace &ws, Tensor7 &cloudbox_field, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &doit_is_initialized, const Agenda &iy_main_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &f_grid, const Index &stokes_dim, const Vector &za_grid, const Vector &aa_grid, const Index &rigorous, const Numeric &maxratio, const Verbosity &)
WORKSPACE METHOD: DoitGetIncoming.
Definition: m_doit.cc:2912
array.h
This file contains the definition of Array.
chk_not_empty
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:694
cloudbox_fieldSetClearsky
void cloudbox_fieldSetClearsky(Tensor7 &cloudbox_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &cloudbox_on, const Index &doit_is_initialized, const Index &all_frequencies, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldSetClearsky.
Definition: m_doit.cc:3467
pha_mat_spt_agendaExecute
void pha_mat_spt_agendaExecute(Workspace &ws, Tensor5 &pha_mat_spt, const Index za_index, const Index scat_lat_index, const Index scat_lon_index, const Index scat_p_index, const Index aa_index, const Numeric rtp_temperature, const Agenda &input_agenda)
Definition: auto_md.cc:24645
DoitCalc
void DoitCalc(Workspace &ws, Tensor7 &cloudbox_field, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &cloudbox_on, const Vector &f_grid, const Agenda &doit_mono_agenda, const Index &doit_is_initialized, const Verbosity &verbosity)
WORKSPACE METHOD: DoitCalc.
Definition: m_doit.cc:2795
cloudbox_fieldUpdateSeq3D
void cloudbox_fieldUpdateSeq3D(Workspace &ws, Tensor6 &cloudbox_field_mono, const Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Vector &aa_grid, const Tensor4 &pnd_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Index &doit_za_interp, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdateSeq3D.
Definition: m_doit.cc:1096
my_basic_string< char >
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
arts.h
The global header file for ARTS.
cloudbox_fieldUpdateSeq1D
void cloudbox_fieldUpdateSeq1D(Workspace &ws, Tensor6 &cloudbox_field_mono, Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Vector &aa_grid, const Tensor4 &pnd_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Vector &p_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Agenda &surface_rtprop_agenda, const Index &doit_za_interp, const Index &normalize, const Numeric &norm_error_threshold, const Index &norm_debug, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdateSeq1D.
Definition: m_doit.cc:771
doit_scat_fieldCalc
void doit_scat_fieldCalc(Workspace &ws, Tensor6 &doit_scat_field, const Agenda &pha_mat_spt_agenda, const Tensor6 &cloudbox_field_mono, const Tensor4 &pnd_field, const Tensor3 &t_field, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Vector &za_grid, const Vector &aa_grid, const Index &doit_za_grid_size, const Tensor7 &pha_mat_doit, const Verbosity &verbosity)
WORKSPACE METHOD: doit_scat_fieldCalc.
Definition: m_doit.cc:2081
_CS_basic_sstream_base::str
_CS_string_type str() const
Definition: sstream.h:491
DoitInit
void DoitInit(Tensor6 &doit_scat_field, Tensor7 &cloudbox_field, Index &doit_is_initialized, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &f_grid, const Vector &za_grid, const Vector &aa_grid, const Index &doit_za_grid_size, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Verbosity &verbosity)
WORKSPACE METHOD: DoitInit.
Definition: m_doit.cc:1582
agenda_class.h
Declarations for agendas.
doit_conv_flagAbsBT
void doit_conv_flagAbsBT(Index &doit_conv_flag, Index &doit_iteration_counter, Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Vector &f_grid, const Index &f_index, const Vector &epsilon, const Index &max_iterations, const Index &throw_nonconv_error, const Verbosity &verbosity)
WORKSPACE METHOD: doit_conv_flagAbsBT.
Definition: m_doit.cc:231
cloudbox_fieldSetConstPerFreq
void cloudbox_fieldSetConstPerFreq(Tensor7 &cloudbox_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Matrix &cloudbox_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldSetConstPerFreq.
Definition: m_doit.cc:3747
math_funcs.h
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
doit_conv_test_agendaExecute
void doit_conv_test_agendaExecute(Workspace &ws, Index &doit_conv_flag, Index &doit_iteration_counter, const Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Agenda &input_agenda)
Definition: auto_md.cc:23675
cloudbox_field_monoIterate
void cloudbox_field_monoIterate(Workspace &ws, Tensor6 &cloudbox_field_mono, const Agenda &doit_scat_field_agenda, const Agenda &doit_rte_agenda, const Agenda &doit_conv_test_agenda, const Index &accelerated, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_field_monoIterate.
Definition: m_doit.cc:497
get_iy
void get_iy(Workspace &ws, Matrix &iy, const Index &cloudbox_on, ConstVectorView f_grid, const EnergyLevelMap &nlte_field, ConstVectorView rte_pos, ConstVectorView rte_los, ConstVectorView rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda)
Basic call of iy_main_agenda.
Definition: rte.cc:877
CREATE_OUT1
#define CREATE_OUT1
Definition: messages.h:205
cloud_fieldsCalc
void cloud_fieldsCalc(Workspace &ws, Tensor5View ext_mat_field, Tensor4View abs_vec_field, const Agenda &spt_calc_agenda, const Index &za_index, const Index &aa_index, const ArrayOfIndex &cloudbox_limits, ConstTensor3View t_field, ConstTensor4View pnd_field, const Verbosity &verbosity)
Calculate ext_mat, abs_vec for all points inside the cloudbox for one.
Definition: doit.cc:163
ConstTensor6View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackVI.cc:54
doit_scat_fieldCalcLimb
void doit_scat_fieldCalcLimb(Workspace &ws, Tensor6 &doit_scat_field, const Agenda &pha_mat_spt_agenda, const Tensor6 &cloudbox_field_mono, const Tensor4 &pnd_field, const Tensor3 &t_field, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Vector &za_grid, const Vector &aa_grid, const Index &doit_za_grid_size, const Index &doit_za_interp, const Tensor7 &pha_mat_doit, const Verbosity &verbosity)
WORKSPACE METHOD: doit_scat_fieldCalcLimb.
Definition: m_doit.cc:2355
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
xml_read_from_file
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:901
joker
const Joker joker
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Matrix
The Matrix class.
Definition: matpackI.h:1193
cloudbox_fieldSetConst
void cloudbox_fieldSetConst(Tensor7 &cloudbox_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &cloudbox_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldSetConst.
Definition: m_doit.cc:3706
ConstTensor7View::nlibraries
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:42
auto_md.h
special_interp.h
Header file for special_interp.cc.
ppath.h
Propagation path structure and functions.
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
Tensor3::resize
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:664
logic.h
Header file for logic.cc.
doit.h
Radiative transfer in cloudbox.
doit_za_grid_optCalc
void doit_za_grid_optCalc(Vector &doit_za_grid_opt, const Tensor6 &cloudbox_field_mono, const Vector &za_grid, const Index &doit_za_interp, const Numeric &acc, const Verbosity &verbosity)
WORKSPACE METHOD: doit_za_grid_optCalc.
Definition: m_doit.cc:2709
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Array
This can be used to make arrays out of anything.
Definition: array.h:40
ConstTensor6View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackVI.cc:45
ConstTensor6View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackVI.cc:57
OptimizeDoitPressureGrid
void OptimizeDoitPressureGrid(Workspace &ws, Vector &p_grid, Tensor4 &pnd_field, Tensor3 &t_field, ArrayOfArrayOfSingleScatteringData &scat_data_mono, Tensor3 &z_field, ArrayOfIndex &cloudbox_limits, Tensor6 &cloudbox_field_mono, Tensor7 &pha_mat_doit, Tensor4 &vmr_field, Vector &p_grid_orig, const Vector &f_grid, const Index &f_index, const Agenda &propmat_clearsky_agenda, const Numeric &tau_scat_max, const Numeric &sgl_alb_max, const Index &cloudbox_size_max, const Verbosity &verbosity)
WORKSPACE METHOD: OptimizeDoitPressureGrid.
Definition: m_doit.cc:1732
chk_size
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:466
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
doit_conv_flagAbs
void doit_conv_flagAbs(Index &doit_conv_flag, Index &doit_iteration_counter, Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Vector &epsilon, const Index &max_iterations, const Index &throw_nonconv_error, const Verbosity &verbosity)
WORKSPACE METHOD: doit_conv_flagAbs.
Definition: m_doit.cc:115
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
pha_matCalc
void pha_matCalc(Tensor4 &pha_mat, const Tensor5 &pha_mat_spt, const Tensor4 &pnd_field, const Index &atmosphere_dim, const Index &scat_p_index, const Index &scat_lat_index, const Index &scat_lon_index, const Verbosity &)
WORKSPACE METHOD: pha_matCalc.
Definition: m_optproperties.cc:1015
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
doit_scat_fieldNormalize
void doit_scat_fieldNormalize(Workspace &ws, Tensor6 &doit_scat_field, const Tensor6 &cloudbox_field_mono, const ArrayOfIndex &cloudbox_limits, const Agenda &spt_calc_agenda, const Index &atmosphere_dim, const Vector &za_grid, const Vector &aa_grid, const Tensor4 &pnd_field, const Tensor3 &t_field, const Numeric &norm_error_threshold, const Index &norm_debug, const Verbosity &verbosity)
Normalization of scattered field.
Definition: doit.cc:1959
ConstTensor6View::npages
Index npages() const
Returns the number of pages.
Definition: matpackVI.cc:51
cloudbox_fieldUpdateSeq1DPP
void cloudbox_fieldUpdateSeq1DPP(Workspace &ws, Tensor6 &cloudbox_field_mono, Index &za_index, const Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Tensor4 &pnd_field, const Vector &p_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdateSeq1DPP.
Definition: m_doit.cc:1410
ConstTensor7View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:48
cloud_ppath_update1D_noseq
void cloud_ppath_update1D_noseq(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &za_index, ConstVectorView za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View cloudbox_field_mono_old, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, ConstVectorView p_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Agenda &surface_rtprop_agenda, const Index &scat_za_interp, const Verbosity &verbosity)
Calculation of radiation field along a propagation path step for specified zenith direction and press...
Definition: doit.cc:457
ConstTensor7View::npages
Index npages() const
Returns the number of pages.
Definition: matpackVII.cc:54
Workspace
Workspace class.
Definition: workspace_ng.h:40
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
Tensor6::resize
void resize(Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVI.cc:2175
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:204
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:207
doit_rte_agendaExecute
void doit_rte_agendaExecute(Workspace &ws, Tensor6 &cloudbox_field_mono, const Tensor6 &doit_scat_field, const Agenda &input_agenda)
Definition: auto_md.cc:23792
cloudbox_fieldSetFromPrecalc
void cloudbox_fieldSetFromPrecalc(Tensor7 &cloudbox_field, const Vector &za_grid, const Vector &f_grid, const Index &atmosphere_dim, const Index &stokes_dim, const ArrayOfIndex &cloudbox_limits, const Index &doit_is_initialized, const Tensor7 &cloudbox_field_precalc, const Verbosity &)
WORKSPACE METHOD: cloudbox_fieldSetFromPrecalc.
Definition: m_doit.cc:3364
is_size
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
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:23495
ConstTensor7View::nvitrines
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVII.cc:45
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
ConstTensor7View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:60
cloud_ppath_update1D_planeparallel
void cloud_ppath_update1D_planeparallel(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &za_index, ConstVectorView za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, ConstVectorView p_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Verbosity &verbosity)
Radiative transfer calculation inside cloudbox for planeparallel case.
Definition: doit.cc:615
PI
const Numeric PI
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
Tensor7::resize
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5484
chk_atm_grids
void chk_atm_grids(const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_grids
Definition: check_input.cc:1201
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
xml_write_to_file
void xml_write_to_file(const String &filename, const T &type, const FileType ftype, const Index no_clobber, const Verbosity &verbosity)
Write data to XML file.
Definition: xml_io.cc:972
cloud_ppath_update1D
void cloud_ppath_update1D(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &za_index, ConstVectorView za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, ConstVectorView p_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Agenda &surface_rtprop_agenda, const Index &scat_za_interp, const Verbosity &verbosity)
Calculates radiation field along a propagation path step for specified zenith direction and pressure ...
Definition: doit.cc:301
physics_funcs.h
This file contains declerations of functions of physical character.
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
ConstTensor7View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:51
Tensor4::resize
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1064
AngIntegrate_trapezoid
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:296
Tensor5
The Tensor5 class.
Definition: matpackV.h:506
rte.h
Declaration of functions in rte.cc.
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
za_gridOpt
void za_gridOpt(Vector &za_grid_opt, Matrix &cloudbox_field_opt, ConstVectorView za_grid_fine, ConstTensor6View cloudbox_field_mono, const Numeric &acc, const Index &scat_za_interp)
Optimize the zenith angle grid.
Definition: doit.cc:1853
wsv_aux.h
Auxiliary header stuff related to workspace variable groups.
check_input.h
ConstTensor6View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackVI.cc:48