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