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