ARTS 2.5.0 (git: 9ee3ac6c)
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 f_grid[Range(f_index, 1)],
1777 rtp_mag_dummy,
1778 ppath_los_dummy,
1779 p_grid[k],
1780 t_field(k, 0, 0),
1781 rtp_nlte_dummy,
1782 vmr_field(joker, k, 0, 0),
1783 propmat_clearsky_agenda);
1784 abs_coeff += cur_propmat_clearsky.Kjj()[0];
1785 abs_coeff /= (Numeric)vmr_field.nbooks(); // FIXME: Is this really as intended???
1786 single_scattering_albedo[cloudbox_index] =
1787 scat_vec[cloudbox_index] / (ext_mat[cloudbox_index] + abs_coeff);
1788 }
1789 //
1790 // Try out the current settings of tau_scat_max and sgl_alb_max
1791 //
1792 do {
1793 scat_data_insert_offset = 0;
1794 for (Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1795 Index cloudbox_index = k - cloudbox_limits[0];
1796 const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1797 single_scattering_albedo[cloudbox_index + 1]) /
1798 2;
1799 const Numeric scat_opt_thk =
1800 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1801 ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1802 if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1803 Index factor = (Index)ceil(scat_opt_thk / tau_scat_max_internal);
1804 for (Index j = 1; j < factor; j++) {
1805 scat_data_insert_offset++;
1806 }
1807 }
1808 }
1809 // If enhancement is too large, change tau_scat_max to a higher value:
1810 if (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] + 1 >
1811 cloudbox_size_max) {
1812 tau_scat_max_internal += 0.01;
1813 was_too_much = true;
1814 }
1815 } while (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] +
1816 1 >
1817 cloudbox_size_max);
1818 scat_data_insert_offset = 0;
1819 // Give warning if enhancement was too much and threshold had to be changed:
1820 if (was_too_much) {
1821 out2
1822 << "Warning: Pressure grid optimization with the thresholds tau_scat_max = "
1823 << tau_scat_max << " and sgl_slb_max = " << sgl_alb_max
1824 << " has lead to an enhancement larger than the value of"
1825 << " cloudbox_size_max = " << cloudbox_size_max
1826 << ". This is why I changed tau_scat_max to " << tau_scat_max_internal
1827 << ". This may lead to errors of a too coarse grid! \n";
1828 }
1829
1830 //--------------------------------------
1831 //Optimize the altitude grid z_grid
1832 //---------------------------------------
1833 for (Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1834 Index cloudbox_index = k - cloudbox_limits[0];
1835 const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1836 single_scattering_albedo[cloudbox_index + 1]) /
1837 2;
1838 const Numeric scat_opt_thk =
1839 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1840 ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1841 z_grid_new.push_back(z_field(k, 0, 0));
1842
1843 if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1844 Index factor = (Index)ceil(scat_opt_thk / tau_scat_max_internal);
1845 Numeric step =
1846 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) / (Numeric)factor;
1847 const SingleScatteringData nextLayer =
1848 scat_data_local[cloudbox_index + scat_data_insert_offset + 1];
1849 const SingleScatteringData currentLayer =
1850 scat_data_local[cloudbox_index + scat_data_insert_offset];
1851
1852 for (Index j = 1; j < factor; j++) {
1853 z_grid_new.push_back(z_field(k, 0, 0) + (Numeric)j * step);
1854 // Perform manual interpolation of scat_data
1855 const Numeric weight = (Numeric)j / (Numeric)factor;
1856 SingleScatteringData newLayer = currentLayer;
1857 Tensor7 weightednextLayerPhamat = nextLayer.pha_mat_data;
1858 Tensor5 weightednextLayerExtmat = nextLayer.ext_mat_data;
1859 Tensor5 weightednextLayerAbsvec = nextLayer.abs_vec_data;
1860
1861 weightednextLayerPhamat *= weight;
1862 weightednextLayerExtmat *= weight;
1863 weightednextLayerAbsvec *= weight;
1864
1865 newLayer.pha_mat_data *= 1. - weight;
1866 newLayer.ext_mat_data *= 1. - weight;
1867 newLayer.abs_vec_data *= 1. - weight;
1868
1869 newLayer.pha_mat_data += weightednextLayerPhamat;
1870 newLayer.ext_mat_data += weightednextLayerExtmat;
1871 newLayer.abs_vec_data += weightednextLayerAbsvec;
1872
1873 // Optimize scat_data
1874 scat_data_local.insert(scat_data_local.begin() + cloudbox_index +
1875 scat_data_insert_offset + 1,
1876 std::move(newLayer));
1877
1878 scat_data_insert_offset++;
1879 }
1880 }
1881 }
1882 // New cloudbox limits
1883 cloudbox_limits_opt[0] = cloudbox_limits[0];
1884 cloudbox_limits_opt[1] = scat_data_insert_offset + cloudbox_limits[1];
1885 const Index cloudbox_opt_size =
1886 cloudbox_limits_opt[1] - cloudbox_limits_opt[0] + 1;
1887 out3 << "Frequency: " << f_grid[f_index]
1888 << " old: " << cloudbox_limits[1] - cloudbox_limits[0] + 1
1889 << " new: " << cloudbox_opt_size << " Factor: "
1890 << (Numeric)cloudbox_opt_size /
1891 (Numeric)(cloudbox_limits[1] - cloudbox_limits[0] + 1)
1892 << "\n";
1893
1894 for (Index i = cloudbox_limits[1]; i < z_field.npages(); i++)
1895 z_grid_new.push_back(z_field(i, 0, 0));
1896
1897 Vector z_grid(z_grid_new.size());
1898 for (Index i = 0; i < z_grid.nelem(); i++) z_grid[i] = z_grid_new[i];
1899 p_grid_orig = p_grid[Range(cloudbox_limits[0],
1900 cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1901 // ---------------------------------------
1902 // Interpolate fields to new z_grid
1903 // ----------------------------------------
1905 Tensor3 t_field_new(z_grid.nelem(), 1, 1);
1906 Vector p_grid_opt(z_grid.nelem());
1907 Tensor6 cloudbox_field_mono_opt(
1908 cloudbox_opt_size, 1, 1, cloudbox_field_mono.npages(), 1, 1);
1909 Tensor7 pha_mat_doit_opt(cloudbox_opt_size,
1910 pha_mat_doit.nvitrines(),
1911 1,
1912 pha_mat_doit.nbooks(),
1913 pha_mat_doit.npages(),
1914 1,
1915 1);
1916 ArrayOfGridPos Z_gp(z_grid.nelem());
1917 Matrix itw_z(Z_gp.nelem(), 2);
1918 ostringstream os;
1919 os << "At the current frequency " << f_grid[f_index]
1920 << " there was an error while interpolating the fields to the new z_field";
1921 chk_interpolation_grids(os.str(), z_field(joker, 0, 0), z_grid);
1922
1923 // Gridpositions of interpolation:
1924 gridpos(Z_gp, z_field(joker, 0, 0), z_grid);
1925 // Interpolation weights:
1926 interpweights(itw_z, Z_gp);
1927
1928 // Interpolate Temperature
1929 interp(t_field_new(joker, 0, 0), itw_z, t_field(joker, 0, 0), Z_gp);
1930 // Write new Temperature to scat_data
1931 for (Index k = cloudbox_limits_opt[0]; k < cloudbox_limits_opt[1]; k++) {
1932 Index i = k - cloudbox_limits[0];
1933 scat_data_local[i].T_grid = t_field_new(i, 0, 0);
1934 }
1935
1936 // Interpolate p_grid
1937 interp(p_grid_opt, itw_z, p_grid, Z_gp);
1938
1939 // Interpolate vmr_field
1940 Tensor4 vmr_field_opt(vmr_field.nbooks(), p_grid_opt.nelem(), 1, 1);
1941 for (Index i = 0; i < vmr_field.nbooks(); i++)
1942 interp(
1943 vmr_field_opt(i, joker, 0, 0), itw_z, vmr_field(i, joker, 0, 0), Z_gp);
1944
1945 // Interpolate cloudbox_field_mono and pha_mat_doit
1946 ArrayOfGridPos Z_gp_2(cloudbox_opt_size);
1947 Matrix itw_z_2(Z_gp_2.nelem(), 2);
1948 Range r1 =
1949 Range(cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1);
1950 Range r2 = Range(cloudbox_limits_opt[0], cloudbox_opt_size);
1951 chk_interpolation_grids(os.str(), z_field(r1, 0, 0), z_grid[r2]);
1952 gridpos(Z_gp_2,
1953 z_field(Range(cloudbox_limits[0],
1954 cloudbox_limits[1] - cloudbox_limits[0] + 1),
1955 0,
1956 0),
1957 z_grid[Range(cloudbox_limits_opt[0], cloudbox_opt_size)]);
1958 interpweights(itw_z_2, Z_gp_2);
1959
1960 for (Index i = 0; i < cloudbox_field_mono.npages(); i++) {
1961 interp(cloudbox_field_mono_opt(joker, 0, 0, i, 0, 0),
1962 itw_z_2,
1963 cloudbox_field_mono(joker, 0, 0, i, 0, 0),
1964 Z_gp_2);
1965 }
1966 for (Index i = 0; i < pha_mat_doit.nvitrines(); i++) {
1967 for (Index j = 0; j < pha_mat_doit.nbooks(); j++) {
1968 for (Index k = 0; k < pha_mat_doit.npages(); k++) {
1969 interp(pha_mat_doit_opt(joker, i, 0, j, k, 0, 0),
1970 itw_z_2,
1971 pha_mat_doit(joker, i, 0, j, k, 0, 0),
1972 Z_gp_2);
1973 }
1974 }
1975 }
1976
1977 // Interpolate pnd-field
1978 pnd_field.resize(cloudbox_opt_size, cloudbox_opt_size, 1, 1);
1979 pnd_field = 0.;
1980 for (Index i = 0; i < cloudbox_opt_size; i++) pnd_field(i, i, 0, 0) = 1.;
1981
1982 //Write new fields
1983 p_grid = p_grid_opt;
1984 t_field = t_field_new;
1985 cloudbox_limits = cloudbox_limits_opt;
1986 cloudbox_field_mono = cloudbox_field_mono_opt;
1987 pha_mat_doit = pha_mat_doit_opt;
1988 z_field.resize(z_grid.nelem(), 1, 1);
1989 z_field(joker, 0, 0) = z_grid;
1990 vmr_field = vmr_field_opt;
1991}
1992
1993/* Workspace method: Doxygen documentation will be auto-generated */
1995 const Index& doit_iteration_counter,
1996 const Tensor6& cloudbox_field_mono,
1997 const Index& f_index,
1998 //Keyword:
1999 const ArrayOfIndex& iterations,
2000 const ArrayOfIndex& frequencies,
2001 const Verbosity& verbosity) {
2002 if (!frequencies.nelem() || !iterations.nelem()) return;
2003
2004 // If the number of iterations is less than a number specified in the
2005 // keyword *iterations*, this number will be ignored.
2006
2007 ostringstream os;
2008 os << "doit_iteration_f" << f_index << "_i" << doit_iteration_counter;
2009
2010 // All iterations for all frequencies are written to files
2011 if (frequencies[0] == -1 && iterations[0] == -1) {
2013 os.str() + ".xml", cloudbox_field_mono, FILE_TYPE_ASCII, 0, verbosity);
2014 }
2015
2016 for (Index j = 0; j < frequencies.nelem(); j++) {
2017 if (f_index == frequencies[j] || (!j && frequencies[j] == -1)) {
2018 // All iterations are written to files
2019 if (iterations[0] == -1) {
2020 xml_write_to_file(os.str() + ".xml",
2021 cloudbox_field_mono,
2023 0,
2024 verbosity);
2025 }
2026
2027 // Only the iterations given by the keyword are written to a file
2028 else {
2029 for (Index i = 0; i < iterations.nelem(); i++) {
2030 if (doit_iteration_counter == iterations[i])
2031 xml_write_to_file(os.str() + ".xml",
2032 cloudbox_field_mono,
2034 0,
2035 verbosity);
2036 }
2037 }
2038 }
2039 }
2040}
2041
2042/* Workspace method: Doxygen documentation will be auto-generated */
2044 // WS Output and Input
2045 Tensor6& doit_scat_field,
2046 //WS Input:
2047 const Agenda& pha_mat_spt_agenda,
2048 const Tensor6& cloudbox_field_mono,
2049 const Tensor4& pnd_field,
2050 const Tensor3& t_field,
2051 const Index& atmosphere_dim,
2052 const ArrayOfIndex& cloudbox_limits,
2053 const Vector& za_grid,
2054 const Vector& aa_grid,
2055 const Index& doit_za_grid_size,
2056 const Tensor7& pha_mat_doit,
2057 const Verbosity& verbosity)
2058
2059{
2062
2063 // ------------ Check the input -------------------------------
2064
2065 // Agenda for calculation of phase matrix
2066 chk_not_empty("pha_mat_spt_agenda", pha_mat_spt_agenda);
2067
2068 // Number of zenith angles.
2069 const Index Nza = za_grid.nelem();
2070
2071 ARTS_USER_ERROR_IF (za_grid[0] != 0. || za_grid[Nza - 1] != 180.,
2072 "The range of *za_grid* must [0 180].");
2073
2074 // Number of azimuth angles.
2075 const Index Naa = aa_grid.nelem();
2076
2077 ARTS_USER_ERROR_IF (Naa > 1 && (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.),
2078 "The range of *aa_grid* must [0 360].");
2079
2080 // Get stokes dimension from *doit_scat_field*:
2081 const Index stokes_dim = doit_scat_field.ncols();
2082 ARTS_ASSERT(stokes_dim > 0 || stokes_dim < 5);
2083
2084 // Size of particle number density field can not be checked here,
2085 // because the function does not use the atmospheric grids.
2086 // Check will be included after fixing the size of pnd field, so
2087 // that it is defined only inside the cloudbox.
2088
2089 // Check atmospheric dimension and dimensions of
2090 // radiation field (*cloudbox_field*) and scattering integral field
2091 // (*doit_scat_field*)
2092 if (atmosphere_dim == 1) {
2093 ARTS_ASSERT(is_size(cloudbox_field_mono,
2094 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2095 1,
2096 1,
2097 Nza,
2098 1,
2099 stokes_dim));
2100 ARTS_ASSERT(is_size(doit_scat_field,
2101 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2102 1,
2103 1,
2104 za_grid.nelem(),
2105 1,
2106 stokes_dim));
2107 } else if (atmosphere_dim == 3) {
2108 ARTS_ASSERT(is_size(cloudbox_field_mono,
2109 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2110 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2111 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2112 Nza,
2113 Naa,
2114 stokes_dim));
2115 ARTS_ASSERT(is_size(doit_scat_field,
2116 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2117 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2118 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2119 Nza,
2120 Naa,
2121 stokes_dim));
2122 } else {
2124 "The atmospheric dimension must be 1D or 3D \n"
2125 "for scattering calculations using the DOIT\n"
2126 "module, but it is not. The value of *atmosphere_dim*\n"
2127 "is ", atmosphere_dim, ".")
2128 }
2129
2130 ARTS_USER_ERROR_IF (cloudbox_limits.nelem() != 2 * atmosphere_dim,
2131 "*cloudbox_limits* is a vector which contains the"
2132 "upper and lower limit of the cloud for all "
2133 "atmospheric dimensions. So its dimension must"
2134 "be 2 x *atmosphere_dim*");
2135
2136 // This function should only be used for down-looking cases where no
2137 // optimized zenith angle grid is required.
2138 ARTS_USER_ERROR_IF (doit_za_grid_size != Nza,
2139 "The zenith angle grids for the computation of\n"
2140 "the scattering integral and the RT part must \n"
2141 "be equal. Check definitions in \n"
2142 "*DOAngularGridsSet*. The keyword \n"
2143 "'za_grid_opt_file' should be empty. \n");
2144
2145 // ------ end of checks -----------------------------------------------
2146
2147 // Initialize variables *pha_mat* and *pha_mat_spt*
2148 Tensor4 pha_mat_local(
2149 doit_za_grid_size, aa_grid.nelem(), stokes_dim, stokes_dim, 0.);
2150
2151 Tensor5 pha_mat_spt_local(pnd_field.nbooks(),
2152 doit_za_grid_size,
2153 aa_grid.nelem(),
2154 stokes_dim,
2155 stokes_dim,
2156 0.);
2157
2158 // Equidistant step size for integration
2159 Vector grid_stepsize(2);
2160 grid_stepsize[0] = 180. / (Numeric)(doit_za_grid_size - 1);
2161 if (Naa > 1) {
2162 grid_stepsize[1] = 360. / (Numeric)(Naa - 1);
2163 }
2164
2165 Tensor3 product_field(Nza, Naa, stokes_dim, 0);
2166
2167 out2 << " Calculate the scattered field\n";
2168
2169 if (atmosphere_dim == 1) {
2170 // Get pha_mat at the grid positions
2171 // Since atmosphere_dim = 1, there is no loop over lat and lon grids
2172 for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2173 p_index++) {
2174 //There is only loop over zenith angle grid ; no azimuth angle grid.
2175 for (Index za_index_local = 0; za_index_local < Nza; za_index_local++) {
2176 // Calculate the phase matrix of individual scattering elements
2177 out3 << "Multiplication of phase matrix with incoming"
2178 << " intensities \n";
2179
2180 product_field = 0;
2181
2182 // za_in and aa_in are for incoming zenith and azimuth
2183 //angle direction for which pha_mat is calculated
2184 for (Index za_in = 0; za_in < Nza; ++za_in) {
2185 for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2186 // Multiplication of phase matrix with incoming
2187 // intensity field.
2188
2189 for (Index i = 0; i < stokes_dim; i++) {
2190 for (Index j = 0; j < stokes_dim; j++) {
2191 product_field(za_in, aa_in, i) +=
2192 pha_mat_doit(
2193 p_index, za_index_local, 0, za_in, aa_in, i, j) *
2194 cloudbox_field_mono(p_index, 0, 0, za_in, 0, j);
2195 }
2196 }
2197
2198 } //end aa_in loop
2199 } //end za_in loop
2200 //integration of the product of ifield_in and pha
2201 // over zenith angle and azimuth angle grid. It calls
2202 if (Naa == 1) {
2203 for (Index i = 0; i < stokes_dim; i++) {
2204 doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2205 AngIntegrate_trapezoid(product_field(joker, 0, i), za_grid) /
2206 2 / PI;
2207 } //end i loop
2208 } else {
2209 for (Index i = 0; i < stokes_dim; i++) {
2210 doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2211 AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2212 za_grid,
2213 aa_grid,
2214 grid_stepsize);
2215 }
2216 }
2217 } //end za_prop loop
2218 } //end p_index loop
2219 } //end atmosphere_dim = 1
2220
2221 //atmosphere_dim = 3
2222 else if (atmosphere_dim == 3) {
2223 /*there is a loop over pressure, latitude and longitudeindex
2224 when we calculate the pha_mat from pha_mat_spt and pnd_field
2225 using the method pha_matCalc. */
2226
2227 for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2228 p_index++) {
2229 for (Index lat_index = 0;
2230 lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2231 lat_index++) {
2232 for (Index lon_index = 0;
2233 lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2234 lon_index++) {
2235 Numeric rtp_temperature_local =
2236 t_field(p_index + cloudbox_limits[0],
2237 lat_index + cloudbox_limits[2],
2238 lon_index + cloudbox_limits[4]);
2239
2240 for (Index aa_index_local = 1; aa_index_local < Naa;
2241 aa_index_local++) {
2242 for (Index za_index_local = 0; za_index_local < Nza;
2243 za_index_local++) {
2244 out3 << "Calculate phase matrix \n";
2246 pha_mat_spt_local,
2247 za_index_local,
2248 lat_index,
2249 lon_index,
2250 p_index,
2251 aa_index_local,
2252 rtp_temperature_local,
2253 pha_mat_spt_agenda);
2254
2255 pha_matCalc(pha_mat_local,
2256 pha_mat_spt_local,
2257 pnd_field,
2258 atmosphere_dim,
2259 p_index,
2260 lat_index,
2261 lon_index,
2262 verbosity);
2263
2264 product_field = 0;
2265
2266 //za_in and aa_in are the incoming directions
2267 //for which pha_mat_spt is calculated
2268 for (Index za_in = 0; za_in < Nza; ++za_in) {
2269 for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2270 // Multiplication of phase matrix
2271 // with incloming intensity field.
2272 for (Index i = 0; i < stokes_dim; i++) {
2273 for (Index j = 0; j < stokes_dim; j++) {
2274 product_field(za_in, aa_in, i) +=
2275 pha_mat_local(za_in, aa_in, i, j) *
2276 cloudbox_field_mono(p_index,
2277 lat_index,
2278 lon_index,
2279 za_index_local,
2280 aa_index_local,
2281 j);
2282 }
2283 }
2284 } //end aa_in loop
2285 } //end za_in loop
2286 //integration of the product of ifield_in and pha
2287 //over zenith angle and azimuth angle grid. It
2288 //calls here the integration routine
2289 //AngIntegrate_trapezoid_opti
2290 for (Index i = 0; i < stokes_dim; i++) {
2291 doit_scat_field(p_index,
2292 lat_index,
2293 lon_index,
2294 za_index_local,
2295 aa_index_local,
2296 i) =
2297 AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2298 za_grid,
2299 aa_grid,
2300 grid_stepsize);
2301 } //end i loop
2302 } //end aa_prop loop
2303 } //end za_prop loop
2304 } //end lon loop
2305 } // end lat loop
2306 } // end p loop
2307 // aa = 0 is the same as aa = 180:
2308 doit_scat_field(joker, joker, joker, joker, 0, joker) =
2309 doit_scat_field(joker, joker, joker, joker, Naa - 1, joker);
2310 } // end atmosphere_dim = 3
2311}
2312
2313/* Workspace method: Doxygen documentation will be auto-generated */
2315 // WS Output and Input
2316 Tensor6& doit_scat_field,
2317 //WS Input:
2318 const Agenda& pha_mat_spt_agenda,
2319 const Tensor6& cloudbox_field_mono,
2320 const Tensor4& pnd_field,
2321 const Tensor3& t_field,
2322 const Index& atmosphere_dim,
2323 const ArrayOfIndex& cloudbox_limits,
2324 const Vector& za_grid,
2325 const Vector& aa_grid,
2326 const Index& doit_za_grid_size,
2327 const Index& doit_za_interp,
2328 const Tensor7& pha_mat_doit,
2329 const Verbosity& verbosity) {
2332
2333 // ------------ Check the input -------------------------------
2334
2335 // Agenda for calculation of phase matrix
2336 chk_not_empty("pha_mat_spt_agenda", pha_mat_spt_agenda);
2337
2338 // Number of zenith angles.
2339 const Index Nza = za_grid.nelem();
2340
2341 ARTS_USER_ERROR_IF (za_grid[0] != 0. || za_grid[Nza - 1] != 180.,
2342 "The range of *za_grid* must [0 180].");
2343
2344 // Number of azimuth angles.
2345 const Index Naa = aa_grid.nelem();
2346
2347 ARTS_USER_ERROR_IF (Naa > 1 && (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.),
2348 "The range of *aa_grid* must [0 360].");
2349
2350 // Get stokes dimension from *doit_scat_field*:
2351 const Index stokes_dim = doit_scat_field.ncols();
2352 ARTS_ASSERT(stokes_dim > 0 || stokes_dim < 5);
2353
2354 // Size of particle number density field can not be checked here,
2355 // because the function does not use the atmospheric grids.
2356 // Check will be included after fixing the size of pnd field, so
2357 // that it is defined only inside the cloudbox.
2358
2359 // Check atmospheric dimension and dimensions of
2360 // radiation field (*cloudbox_field*) and scattering integral field
2361 // (*doit_scat_field*)
2362 if (atmosphere_dim == 1) {
2363 ARTS_ASSERT(is_size(cloudbox_field_mono,
2364 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2365 1,
2366 1,
2367 Nza,
2368 1,
2369 stokes_dim));
2370 ARTS_ASSERT(is_size(doit_scat_field,
2371 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2372 1,
2373 1,
2374 za_grid.nelem(),
2375 1,
2376 stokes_dim));
2377 } else if (atmosphere_dim == 3) {
2378 ARTS_ASSERT(is_size(cloudbox_field_mono,
2379 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2380 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2381 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2382 Nza,
2383 Naa,
2384 stokes_dim));
2385 ARTS_ASSERT(is_size(doit_scat_field,
2386 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2387 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2388 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2389 Nza,
2390 Naa,
2391 stokes_dim));
2392 } else {
2394 "The atmospheric dimension must be 1D or 3D \n"
2395 "for scattering calculations using the DOIT\n"
2396 "module, but it is not. The value of *atmosphere_dim*\n"
2397 "is ", atmosphere_dim, ".")
2398 }
2399
2400 ARTS_USER_ERROR_IF (!(doit_za_interp == 0 || doit_za_interp == 1),
2401 "Interpolation method is not defined. Use \n"
2402 "*doit_za_interpSet*.\n");
2403
2404 ARTS_USER_ERROR_IF (cloudbox_limits.nelem() != 2 * atmosphere_dim,
2405 "*cloudbox_limits* is a vector which contains the"
2406 "upper and lower limit of the cloud for all "
2407 "atmospheric dimensions. So its dimension must"
2408 "be 2 x *atmosphere_dim*");
2409
2410 ARTS_USER_ERROR_IF (doit_za_grid_size < 16,
2411 "*doit_za_grid_size* must be greater than 15 for"
2412 "accurate results");
2413 if (doit_za_grid_size > 100) {
2415 out1 << "Warning: doit_za_grid_size is very large which means that the \n"
2416 << "calculation will be very slow. The recommended value is 19.\n";
2417 }
2418
2419 // ------ end of checks -----------------------------------------------
2420
2421 // Initialize variables *pha_mat* and *pha_mat_spt*
2422 Tensor4 pha_mat_local(
2423 doit_za_grid_size, aa_grid.nelem(), stokes_dim, stokes_dim, 0.);
2424
2425 Tensor5 pha_mat_spt_local(pnd_field.nbooks(),
2426 doit_za_grid_size,
2427 aa_grid.nelem(),
2428 stokes_dim,
2429 stokes_dim,
2430 0.);
2431
2432 // Create the grids for the calculation of the scattering integral.
2433 Vector za_g;
2434 nlinspace(za_g, 0, 180, doit_za_grid_size);
2435
2436 // Two interpolations are required. First we have to interpolate the
2437 // intensity field on the equidistant grid:
2438 ArrayOfGridPos gp_za_i(doit_za_grid_size);
2439 gridpos(gp_za_i, za_grid, za_g);
2440
2441 Matrix itw_za_i(doit_za_grid_size, 2);
2442 interpweights(itw_za_i, gp_za_i);
2443
2444 // Intensity field interpolated on equidistant grid.
2445 Matrix cloudbox_field_int(doit_za_grid_size, stokes_dim, 0);
2446
2447 // Second, we have to interpolate the scattering integral on the RT
2448 // zenith angle grid.
2449 ArrayOfGridPos gp_za(Nza);
2450 gridpos(gp_za, za_g, za_grid);
2451
2452 Matrix itw_za(Nza, 2);
2453 interpweights(itw_za, gp_za);
2454
2455 // Original scattered field, on equidistant zenith angle grid.
2456 Matrix doit_scat_field_org(doit_za_grid_size, stokes_dim, 0);
2457
2458 // Grid stepsize of zenith and azimuth angle grid, these are needed for the
2459 // integration function.
2460 Vector grid_stepsize(2);
2461 grid_stepsize[0] = 180. / (Numeric)(doit_za_grid_size - 1);
2462 if (Naa > 1) {
2463 grid_stepsize[1] = 360. / (Numeric)(Naa - 1);
2464 }
2465
2466 Tensor3 product_field(doit_za_grid_size, Naa, stokes_dim, 0);
2467
2468 if (atmosphere_dim == 1) {
2469 // Get pha_mat at the grid positions
2470 // Since atmosphere_dim = 1, there is no loop over lat and lon grids
2471 for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2472 p_index++) {
2473 // Interpolate intensity field:
2474 for (Index i = 0; i < stokes_dim; i++) {
2475 if (doit_za_interp == 0) {
2476 interp(cloudbox_field_int(joker, i),
2477 itw_za_i,
2478 cloudbox_field_mono(p_index, 0, 0, joker, 0, i),
2479 gp_za_i);
2480 } else if (doit_za_interp == 1) {
2481 // Polynomial
2482 for (Index za = 0; za < za_g.nelem(); za++) {
2483 cloudbox_field_int(za, i) =
2484 interp_poly(za_grid,
2485 cloudbox_field_mono(p_index, 0, 0, joker, 0, i),
2486 za_g[za],
2487 gp_za_i[za]);
2488 }
2489 }
2490 // doit_za_interp must be 0 or 1 (linear or polynomial)!!!
2491 else
2492 ARTS_ASSERT(false);
2493 }
2494
2495 //There is only loop over zenith angle grid; no azimuth angle grid.
2496 for (Index za_index_local = 0; za_index_local < doit_za_grid_size;
2497 za_index_local++) {
2498 out3 << "Multiplication of phase matrix with incoming"
2499 << " intensities \n";
2500
2501 product_field = 0;
2502
2503 // za_in and aa_in are for incoming zenith and azimuth
2504 // angle direction for which pha_mat is calculated
2505 for (Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2506 for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2507 // Multiplication of phase matrix with incoming
2508 // intensity field.
2509
2510 for (Index i = 0; i < stokes_dim; i++) {
2511 for (Index j = 0; j < stokes_dim; j++) {
2512 product_field(za_in, aa_in, i) +=
2513 pha_mat_doit(
2514 p_index, za_index_local, 0, za_in, aa_in, i, j) *
2515 cloudbox_field_int(za_in, j);
2516 }
2517 }
2518
2519 } //end aa_in loop
2520 } //end za_in loop
2521
2522 out3 << "Compute integral. \n";
2523 if (Naa == 1) {
2524 for (Index i = 0; i < stokes_dim; i++) {
2525 doit_scat_field_org(za_index_local, i) =
2526 AngIntegrate_trapezoid(product_field(joker, 0, i), za_g) / 2 /
2527 PI;
2528 } //end i loop
2529 } else {
2530 for (Index i = 0; i < stokes_dim; i++) {
2531 doit_scat_field_org(za_index_local, i) =
2532 AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2533 za_g,
2534 aa_grid,
2535 grid_stepsize);
2536 }
2537 }
2538
2539 } //end za_prop loop
2540
2541 // Interpolation on za_grid, which is used in
2542 // radiative transfer part.
2543 for (Index i = 0; i < stokes_dim; i++) {
2544 if (doit_za_interp == 0) // linear interpolation
2545 {
2546 interp(doit_scat_field(p_index, 0, 0, joker, 0, i),
2547 itw_za,
2548 doit_scat_field_org(joker, i),
2549 gp_za);
2550 } else // polynomial interpolation
2551 {
2552 for (Index za = 0; za < za_grid.nelem(); za++) {
2553 doit_scat_field(p_index, 0, 0, za, 0, i) = interp_poly(
2554 za_g, doit_scat_field_org(joker, i), za_grid[za], gp_za[za]);
2555 }
2556 }
2557 }
2558 } //end p_index loop
2559
2560 } //end atmosphere_dim = 1
2561
2562 else if (atmosphere_dim == 3) {
2563 // Loop over all positions
2564 for (Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2565 p_index++) {
2566 for (Index lat_index = 0;
2567 lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2568 lat_index++) {
2569 for (Index lon_index = 0;
2570 lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2571 lon_index++) {
2572 Numeric rtp_temperature_local =
2573 t_field(p_index + cloudbox_limits[0],
2574 lat_index + cloudbox_limits[2],
2575 lon_index + cloudbox_limits[4]);
2576
2577 // Loop over scattered directions
2578 for (Index aa_index_local = 1; aa_index_local < Naa;
2579 aa_index_local++) {
2580 // Interpolate intensity field:
2581 for (Index i = 0; i < stokes_dim; i++) {
2582 interp(
2583 cloudbox_field_int(joker, i),
2584 itw_za_i,
2585 cloudbox_field_mono(
2586 p_index, lat_index, lon_index, joker, aa_index_local, i),
2587 gp_za_i);
2588 }
2589
2590 for (Index za_index_local = 0; za_index_local < doit_za_grid_size;
2591 za_index_local++) {
2592 out3 << "Calculate phase matrix \n";
2594 pha_mat_spt_local,
2595 za_index_local,
2596 lat_index,
2597 lon_index,
2598 p_index,
2599 aa_index_local,
2600 rtp_temperature_local,
2601 pha_mat_spt_agenda);
2602
2603 pha_matCalc(pha_mat_local,
2604 pha_mat_spt_local,
2605 pnd_field,
2606 atmosphere_dim,
2607 p_index,
2608 lat_index,
2609 lon_index,
2610 verbosity);
2611
2612 product_field = 0;
2613
2614 //za_in and aa_in are the incoming directions
2615 //for which pha_mat_spt is calculated
2616 out3 << "Multiplication of phase matrix with"
2617 << "incoming intensity \n";
2618
2619 for (Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2620 for (Index aa_in = 0; aa_in < Naa; ++aa_in) {
2621 // Multiplication of phase matrix
2622 // with incloming intensity field.
2623 for (Index i = 0; i < stokes_dim; i++) {
2624 for (Index j = 0; j < stokes_dim; j++) {
2625 product_field(za_in, aa_in, i) +=
2626 pha_mat_local(za_in, aa_in, i, j) *
2627 cloudbox_field_int(za_in, j);
2628 }
2629 }
2630 } //end aa_in loop
2631 } //end za_in loop
2632
2633 out3 << "Compute the integral \n";
2634
2635 for (Index i = 0; i < stokes_dim; i++) {
2636 doit_scat_field_org(za_index_local, i) =
2637 AngIntegrate_trapezoid_opti(product_field(joker, joker, i),
2638 za_grid,
2639 aa_grid,
2640 grid_stepsize);
2641 } //end stokes_dim loop
2642
2643 } //end za_prop loop
2644 //Interpolate on original za_grid.
2645 for (Index i = 0; i < stokes_dim; i++) {
2646 interp(
2647 doit_scat_field(
2648 p_index, lat_index, lon_index, joker, aa_index_local, i),
2649 itw_za,
2650 doit_scat_field_org(joker, i),
2651 gp_za);
2652 }
2653 } // end aa_prop loop
2654 } //end lon loop
2655 } //end lat loop
2656 } // end p loop
2657 doit_scat_field(joker, joker, joker, joker, 0, joker) =
2658 doit_scat_field(joker, joker, joker, joker, Naa - 1, joker);
2659 } // end atm_dim=3
2660 out2 << " Finished scattered field.\n";
2661}
2662
2663/* Workspace method: Doxygen documentation will be auto-generated */
2664void doit_za_grid_optCalc( //WS Output
2665 Vector& doit_za_grid_opt,
2666 // WS Input:
2667 const Tensor6& cloudbox_field_mono,
2668 const Vector& za_grid,
2669 const Index& doit_za_interp,
2670 //Keywords:
2671 const Numeric& acc,
2672 const Verbosity& verbosity) {
2674 //-------- Check the input ---------------------------------
2675
2676 // Here it is checked whether cloudbox_field is 1D and whether it is
2677 // consistent with za_grid. The number of pressure levels and the
2678 // number of stokes components does not matter.
2679 chk_size("cloudbox_field",
2680 cloudbox_field_mono,
2681 cloudbox_field_mono.nvitrines(),
2682 1,
2683 1,
2684 za_grid.nelem(),
2685 1,
2686 cloudbox_field_mono.ncols());
2687
2688 ARTS_USER_ERROR_IF (cloudbox_field_mono.ncols() < 1 || cloudbox_field_mono.ncols() > 4,
2689 "The last dimension of *cloudbox_field* corresponds\n"
2690 "to the Stokes dimension, therefore the number of\n"
2691 "columns in *cloudbox_field* must be a number between\n"
2692 "1 and 4, but it is not!");
2693
2694 ARTS_USER_ERROR_IF (!(doit_za_interp == 0 || doit_za_interp == 1),
2695 "Interpolation method is not defined. Use \n"
2696 "*doit_za_interpSet*.\n");
2697
2698 if (za_grid.nelem() < 500) {
2699 out1 << "Warning: The fine grid (*za_grid*) has less than\n"
2700 << "500 grid points which is likely not sufficient for\n"
2701 << "grid_optimization\n";
2702 /* throw runtime_error("The fine grid (*za_grid*) has less than \n"
2703 "500 grid points which is not sufficient for \n"
2704 "grid_optimization");
2705*/
2706 }
2707 // ------------- end of checks -------------------------------------
2708
2709 // Here only used as dummy variable.
2710 Matrix cloudbox_field_opt_mat;
2711 cloudbox_field_opt_mat = 0.;
2712
2713 // Optimize zenith angle grid.
2714 za_gridOpt(doit_za_grid_opt,
2715 cloudbox_field_opt_mat,
2716 za_grid,
2717 cloudbox_field_mono,
2718 acc,
2719 doit_za_interp);
2720}
2721
2722/* Workspace method: Doxygen documentation will be auto-generated */
2723void doit_za_interpSet(Index& doit_za_interp,
2724 const Index& atmosphere_dim,
2725 //Keyword
2726 const String& method,
2727 const Verbosity&) {
2728 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2729
2730 ARTS_USER_ERROR_IF (atmosphere_dim != 1 && method == "polynomial",
2731 "Polynomial interpolation is only implemented for\n"
2732 "1D DOIT calculations as \n"
2733 "in 3D there can be numerical problems.\n"
2734 "Please use 'linear' interpolation method.");
2735
2736 if (method == "linear")
2737 doit_za_interp = 0;
2738 else if (method == "polynomial")
2739 doit_za_interp = 1;
2740 else {
2742 "Possible interpolation methods are 'linear' "
2743 "and 'polynomial'.\n");
2744 }
2745}
2746
2747/* Workspace method: Doxygen documentation will be auto-generated */
2749 Tensor7& cloudbox_field,
2750 const Index& atmfields_checked,
2751 const Index& atmgeom_checked,
2752 const Index& cloudbox_checked,
2753 const Index& scat_data_checked,
2754 const Index& cloudbox_on,
2755 const Vector& f_grid,
2756 const Agenda& doit_mono_agenda,
2757 const Index& doit_is_initialized,
2758 const Verbosity& verbosity)
2759
2760{
2762
2763 if (!cloudbox_on) {
2765 out0 << " Cloudbox is off, DOIT calculation will be skipped.\n";
2766 return;
2767 //throw runtime_error( "Cloudbox is off, no scattering calculations to be"
2768 // "performed." );
2769 }
2770
2771 //-------- Check input -------------------------------------------
2772
2773 ARTS_USER_ERROR_IF (atmfields_checked != 1,
2774 "The atmospheric fields must be flagged to have "
2775 "passed a consistency check (atmfields_checked=1).");
2776 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
2777 "The atmospheric geometry must be flagged to have "
2778 "passed a consistency check (atmgeom_checked=1).");
2779 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
2780 "The cloudbox must be flagged to have "
2781 "passed a consistency check (cloudbox_checked=1).");
2782
2783 // Don't do anything if there's no cloudbox defined.
2784 if (!cloudbox_on) return;
2785
2786 ARTS_USER_ERROR_IF (scat_data_checked != 1,
2787 "The scattering data must be flagged to have "
2788 "passed a consistency check (scat_data_checked=1).");
2789
2790 chk_not_empty("doit_mono_agenda", doit_mono_agenda);
2791
2792 // Frequency grid
2793 //
2794 ARTS_USER_ERROR_IF (f_grid.empty(), "The frequency grid is empty.");
2795 chk_if_increasing("f_grid", f_grid);
2796
2797 // Check whether DoitInit was executed
2798 ARTS_USER_ERROR_IF (!doit_is_initialized,
2799 "Initialization method *DoitInit* has to be called before "
2800 "*DoitCalc*")
2801
2802 //-------- end of checks ----------------------------------------
2803
2804 // We have to make a local copy of the Workspace and the agendas because
2805 // only non-reference types can be declared firstprivate in OpenMP
2806 Workspace l_ws(ws);
2807 Agenda l_doit_mono_agenda(doit_mono_agenda);
2808
2809 // OMP likes simple loop end conditions, so we make a local copy here:
2810 const Index nf = f_grid.nelem();
2811
2812 if (nf) {
2813 String fail_msg;
2814 bool failed = false;
2815
2816#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) \
2817 firstprivate(l_ws, l_doit_mono_agenda)
2818 for (Index f_index = 0; f_index < nf; f_index++) {
2819 if (failed) {
2820 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) = NAN;
2821 continue;
2822 }
2823
2824 try {
2825 ostringstream os;
2826 os << "Frequency: " << f_grid[f_index] / 1e9 << " GHz \n";
2827 out2 << os.str();
2828
2829 Tensor6 cloudbox_field_mono_local =
2830 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker);
2832 cloudbox_field_mono_local,
2833 f_grid,
2834 f_index,
2835 l_doit_mono_agenda);
2836 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) =
2837 cloudbox_field_mono_local;
2838 } catch (const std::exception& e) {
2839 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) = NAN;
2840 ostringstream os;
2841 os << "Error for f_index = " << f_index << " (" << f_grid[f_index]
2842 << " Hz)" << endl
2843 << e.what();
2844#pragma omp critical(DoitCalc_fail)
2845 {
2846 failed = true;
2847 fail_msg = os.str();
2848 }
2849 continue;
2850 }
2851 }
2852
2853 ARTS_USER_ERROR_IF (failed, fail_msg);
2854 }
2855}
2856
2857/* Workspace method: Doxygen documentation will be auto-generated */
2859 Tensor7& cloudbox_field,
2860 const Index& atmfields_checked,
2861 const Index& atmgeom_checked,
2862 const Index& cloudbox_checked,
2863 const Index& doit_is_initialized,
2864 const Agenda& iy_main_agenda,
2865 const Index& atmosphere_dim,
2866 const Vector& lat_grid,
2867 const Vector& lon_grid,
2868 const Tensor3& z_field,
2869 const EnergyLevelMap& nlte_field,
2870 const Index& cloudbox_on,
2871 const ArrayOfIndex& cloudbox_limits,
2872 const Vector& f_grid,
2873 const Index& stokes_dim,
2874 const Vector& za_grid,
2875 const Vector& aa_grid,
2876 const Index& rigorous,
2877 const Numeric& maxratio,
2878 const Verbosity&) {
2879 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2880 ARTS_USER_ERROR_IF (atmfields_checked != 1,
2881 "The atmospheric fields must be flagged to have "
2882 "passed a consistency check (atmfields_checked=1).");
2883 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
2884 "The atmospheric geometry must be flagged to have "
2885 "passed a consistency check (atmgeom_checked=1).");
2886 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
2887 "The cloudbox must be flagged to have "
2888 "passed a consistency check (cloudbox_checked=1).");
2889
2890 // Don't do anything if there's no cloudbox defined.
2891 if (!cloudbox_on) return;
2892
2893 // Check whether DoitInit was executed
2894 ARTS_USER_ERROR_IF (!doit_is_initialized,
2895 "Initialization method *DoitInit* has to be called before "
2896 "*DoitGetIncoming*.")
2897
2898 // iy_unit hard.coded to "1" here
2899 const String iy_unit = "1";
2900
2901 Index Nf = f_grid.nelem();
2902 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2903 Index Nza = za_grid.nelem();
2904
2905 Matrix iy;
2906
2907 //dummy variables needed for the call of iy_main_agendaExecute
2908 ArrayOfMatrix iy_aux;
2909 ArrayOfTensor3 diy_dx;
2910 Ppath ppath;
2911 Tensor3 iy_transmittance(0, 0, 0);
2912 const ArrayOfString iy_aux_vars(0);
2913
2914
2915
2916 //--- Check input ----------------------------------------------------------
2917 ARTS_USER_ERROR_IF (!(atmosphere_dim == 1 || atmosphere_dim == 3),
2918 "The atmospheric dimensionality must be 1 or 3.");
2919 ARTS_USER_ERROR_IF (za_grid[0] != 0. || za_grid[Nza - 1] != 180.,
2920 "*za_grid* must include 0 and 180 degrees as endpoints.");
2921 //--------------------------------------------------------------------------
2922
2923 if (atmosphere_dim == 1) {
2924 //Define the variables for position and direction.
2925 Vector los(1), pos(1);
2926
2927 //--- Get cloudbox_field at lower and upper boundary
2928 // (boundary=0: lower, boundary=1: upper)
2929 for (Index boundary = 0; boundary <= 1; boundary++) {
2930 const Index boundary_index =
2931 boundary ? cloudbox_field.nvitrines() - 1 : 0;
2932 pos[0] = z_field(cloudbox_limits[boundary], 0, 0);
2933
2934 // doing the first angle separately for allowing dy between 2 angles
2935 // in the loop
2936 los[0] = za_grid[0];
2937
2939 iy,
2940 iy_aux,
2941 ppath,
2942 diy_dx,
2943 1,
2944 iy_transmittance,
2945 iy_aux_vars,
2946 0,
2947 iy_unit,
2948 0,
2949 0,
2950 f_grid,
2951 nlte_field,
2952 pos,
2953 los,
2954 Vector(0),
2955 iy_main_agenda);
2956
2957 cloudbox_field(joker, boundary_index, 0, 0, 0, 0, joker) = iy;
2958
2959 for (Index za_index = 1; za_index < Nza; za_index++) {
2960 los[0] = za_grid[za_index];
2961
2963 iy,
2964 iy_aux,
2965 ppath,
2966 diy_dx,
2967 1,
2968 iy_transmittance,
2969 iy_aux_vars,
2970 0,
2971 iy_unit,
2972 0,
2973 0,
2974 f_grid,
2975 nlte_field,
2976 pos,
2977 los,
2978 Vector(0),
2979 iy_main_agenda);
2980
2981 cloudbox_field(joker, boundary_index, 0, 0, za_index, 0, joker) = iy;
2982
2983 if (rigorous) {
2984 for (Index fi = 0; fi < Nf; fi++) {
2985 ARTS_USER_ERROR_IF (cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0) /
2986 cloudbox_field(
2987 fi, boundary_index, 0, 0, za_index, 0, 0) >
2988 maxratio ||
2989 cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0) /
2990 cloudbox_field(
2991 fi, boundary_index, 0, 0, za_index, 0, 0) <
2992 1 / maxratio,
2993 "ERROR: Radiance difference between "
2994 "interpolation points is too large (factor ", maxratio,
2995 ")\n"
2996 "to safely interpolate. This might be due to "
2997 "za_grid being too coarse or the radiance field "
2998 "being a step-like function.\n"
2999 "Happens at boundary ", boundary_index,
3000 " between zenith angles ", za_grid[za_index - 1],
3001 " and ", za_grid[za_index], "deg\n"
3002 "for frequency #", fi, ", where radiances are ",
3003 cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0),
3004 " and ",
3005 cloudbox_field(fi, boundary_index, 0, 0, za_index, 0, 0),
3006 " W/(sr m2 Hz).")
3007 }
3008 }
3009 }
3010 }
3011 }
3012
3013 //--- atmosphere_dim = 3: --------------------------------------------------
3014 else {
3015 Index Naa = aa_grid.nelem();
3016
3017 ARTS_USER_ERROR_IF (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.,
3018 "*aa_grid* must include 0 and 360 degrees as endpoints.");
3019
3020 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3021 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3022
3023 // Convert aa_grid to "sensor coordinates"
3024 // (-180° < azimuth angle < 180°)
3025 //
3026 Vector aa_g(Naa);
3027 for (Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3028
3029 // Define the variables for position and direction.
3030 Vector los(2), pos(3);
3031
3032 //--- Get cloudbox_field at lower and upper boundary
3033 // (boundary=0: lower, boundary=1: upper)
3034 for (Index boundary = 0; boundary <= 1; boundary++) {
3035 const Index boundary_index =
3036 boundary ? cloudbox_field.nvitrines() - 1 : 0;
3037 for (Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3038 for (Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3039 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3040 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3041 pos[0] = z_field(cloudbox_limits[boundary],
3042 lat_index + cloudbox_limits[2],
3043 lon_index + cloudbox_limits[4]);
3044
3045 for (Index za_index = 0; za_index < Nza; za_index++) {
3046 for (Index aa_index = 0; aa_index < Naa; aa_index++) {
3047 los[0] = za_grid[za_index];
3048 los[1] = aa_g[aa_index];
3049
3050 // For end points of za_index (0 & 180deg), we
3051 // only need to perform calculations for one scat_aa
3052 // and set the others to same value
3053 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3055 iy,
3056 iy_aux,
3057 ppath,
3058 diy_dx,
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 1,
3112 iy_transmittance,
3113 iy_aux_vars,
3114 0,
3115 iy_unit,
3116 0,
3117 0,
3118 f_grid,
3119 nlte_field,
3120 pos,
3121 los,
3122 Vector(0),
3123 iy_main_agenda);
3124 }
3125
3126 cloudbox_field(joker,
3127 p_index,
3128 boundary_index,
3129 lon_index,
3130 za_index,
3131 aa_index,
3132 joker) = iy;
3133 }
3134 }
3135 }
3136 }
3137 }
3138
3139 //--- Get scat_i_lon (1st and 2nd boundary):
3140 for (Index boundary = 0; boundary <= 1; boundary++) {
3141 const Index boundary_index = boundary ? cloudbox_field.nbooks() - 1 : 0;
3142 for (Index p_index = 0; p_index < Np_cloud; p_index++) {
3143 for (Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3144 pos[2] = lon_grid[cloudbox_limits[boundary + 4]];
3145 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3146 pos[0] = z_field(p_index + cloudbox_limits[0],
3147 lat_index + cloudbox_limits[2],
3148 cloudbox_limits[boundary + 4]);
3149
3150 for (Index za_index = 0; za_index < Nza; za_index++) {
3151 for (Index aa_index = 0; aa_index < Naa; aa_index++) {
3152 los[0] = za_grid[za_index];
3153 los[1] = aa_g[aa_index];
3154
3155 // For end points of za_index, we need only to
3156 // perform calculations for first scat_aa
3157 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3159 iy,
3160 iy_aux,
3161 ppath,
3162 diy_dx,
3163 1,
3164 iy_transmittance,
3165 iy_aux_vars,
3166 0,
3167 iy_unit,
3168 0,
3169 0,
3170 f_grid,
3171 nlte_field,
3172 pos,
3173 los,
3174 Vector(0),
3175 iy_main_agenda);
3176 }
3177
3178 cloudbox_field(joker,
3179 p_index,
3180 lat_index,
3181 boundary_index,
3182 za_index,
3183 aa_index,
3184 joker) = iy;
3185 }
3186 }
3187 }
3188 }
3189 }
3190 } // End atmosphere_dim = 3.
3191}
3192
3193/* Workspace method: Doxygen documentation will be auto-generated */
3195 Tensor7& cloudbox_field,
3196 Index& cloudbox_on,
3197 const Index& atmfields_checked,
3198 const Index& atmgeom_checked,
3199 const Index& cloudbox_checked,
3200 const Index& doit_is_initialized,
3201 const Agenda& iy_main_agenda,
3202 const Index& atmosphere_dim,
3203 const Vector& lat_grid,
3204 const Vector& lon_grid,
3205 const Tensor3& z_field,
3206 const EnergyLevelMap& nlte_field,
3207 const ArrayOfIndex& cloudbox_limits,
3208 const Vector& f_grid,
3209 const Index& stokes_dim,
3210 const Vector& za_grid,
3211 const Vector& aa_grid,
3212 const Verbosity&) {
3213 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
3214 ARTS_USER_ERROR_IF (atmfields_checked != 1,
3215 "The atmospheric fields must be flagged to have "
3216 "passed a consistency check (atmfields_checked=1).");
3217 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
3218 "The atmospheric geometry must be flagged to have "
3219 "passed a consistency check (atmgeom_checked=1).");
3220 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
3221 "The cloudbox must be flagged to have "
3222 "passed a consistency check (cloudbox_checked=1).");
3223
3224 // Don't do anything if there's no cloudbox defined.
3225 if (!cloudbox_on) return;
3226
3227 // Check whether DoitInit was executed
3228 ARTS_USER_ERROR_IF (!doit_is_initialized,
3229 "Initialization method *DoitInit* has to be called before "
3230 "*DoitGetIncoming1DAtm*.")
3231
3232 // iy_unit hard.coded to "1" here
3233 const String iy_unit = "1";
3234
3235 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3236 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3237 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3238 Index Nza = za_grid.nelem();
3239 Index Naa = aa_grid.nelem();
3240
3241 Matrix iy;
3242
3243 //dummy variables needed for the call of iy_main_agendaExecute
3244 ArrayOfMatrix iy_aux;
3245 ArrayOfTensor3 diy_dx;
3246 Ppath ppath;
3247 Tensor3 iy_transmittance(0, 0, 0);
3248 const ArrayOfString iy_aux_vars(0);
3249
3250 //--- Check input ----------------------------------------------------------
3251 ARTS_USER_ERROR_IF (atmosphere_dim != 3,
3252 "The atmospheric dimensionality must be 3.");
3253 ARTS_USER_ERROR_IF (za_grid[0] != 0. || za_grid[Nza - 1] != 180.,
3254 "*za_grid* must include 0 and 180 degrees as endpoints.");
3255 ARTS_USER_ERROR_IF (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.,
3256 "*aa_grid* must include 0 and 360 degrees as endpoints.");
3257 //--------------------------------------------------------------------------
3258
3259 // Dummy variable for flag cloudbox_on. It has to be 0 here not to get
3260 // stuck in an infinite loop (if some propagation path hits the cloud
3261 // box at some other position.
3262 cloudbox_on = 0;
3263
3264 // Convert za_grid to "sensor coordinates"
3265 //(-180° < azimuth angle < 180°)
3266 //
3267 Vector aa_g(Naa);
3268 for (Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3269
3270 // As the atmosphere is spherically symmetric we only have to calculate
3271 // one azimuth angle.
3272 Index aa_index = 0;
3273
3274 // Define the variables for position and direction.
3275 Vector los(2), pos(3);
3276
3277 // These variables are constant for all calculations:
3278 pos[1] = lat_grid[cloudbox_limits[2]];
3279 pos[2] = lon_grid[cloudbox_limits[4]];
3280 los[1] = aa_g[aa_index];
3281
3282 // Calculate cloudbox_field on boundary
3283 for (Index p_index = 0; p_index < Np_cloud; p_index++) {
3284 pos[0] = z_field(
3285 cloudbox_limits[0] + p_index, cloudbox_limits[2], cloudbox_limits[4]);
3286
3287 for (Index za_index = 0; za_index < Nza; za_index++) {
3288 los[0] = za_grid[za_index];
3289
3291 iy,
3292 iy_aux,
3293 ppath,
3294 diy_dx,
3295 1,
3296 iy_transmittance,
3297 iy_aux_vars,
3298 0,
3299 iy_unit,
3300 0,
3301 0,
3302 f_grid,
3303 nlte_field,
3304 pos,
3305 los,
3306 Vector(0),
3307 iy_main_agenda);
3308
3309 for (Index aa = 0; aa < Naa; aa++) {
3310 // cloudbox_field lower boundary
3311 if (p_index == 0) {
3312 for (Index lat = 0; lat < Nlat_cloud; lat++) {
3313 for (Index lon = 0; lon < Nlon_cloud; lon++) {
3314 cloudbox_field(joker, 0, lat, lon, za_index, aa, joker) = iy;
3315 }
3316 }
3317 }
3318 //cloudbox_field at upper boundary
3319 else if (p_index == Np_cloud - 1)
3320 for (Index lat = 0; lat < Nlat_cloud; lat++) {
3321 for (Index lon = 0; lon < Nlon_cloud; lon++) {
3322 cloudbox_field(joker,
3323 cloudbox_field.nvitrines(),
3324 lat,
3325 lon,
3326 za_index,
3327 aa,
3328 joker) = iy;
3329 }
3330 }
3331
3332 // scat_i_lat (both boundaries)
3333 for (Index lat = 0; lat < 2; lat++) {
3334 for (Index lon = 0; lon < Nlon_cloud; lon++) {
3335 const Index boundary_index =
3336 lat ? cloudbox_field.nshelves() - 1 : 0;
3337 cloudbox_field(
3338 joker, p_index, boundary_index, lon, za_index, aa, joker) = iy;
3339 }
3340 }
3341
3342 // scat_i_lon (both boundaries)
3343 for (Index lat = 0; lat < Nlat_cloud; lat++) {
3344 for (Index lon = 0; lon < 2; lon++) {
3345 const Index boundary_index = lon ? cloudbox_field.nbooks() - 1 : 0;
3346 cloudbox_field(
3347 joker, p_index, lat, boundary_index, za_index, aa, joker) = iy;
3348 }
3349 }
3350 }
3351 }
3352 }
3353 cloudbox_on = 1;
3354}
3355
3356/* Workspace method: Doxygen documentation will be auto-generated */
3358 const Vector& za_grid,
3359 const Vector& f_grid,
3360 const Index& atmosphere_dim,
3361 const Index& stokes_dim,
3362 const ArrayOfIndex& cloudbox_limits,
3363 const Index& doit_is_initialized,
3364 const Tensor7& cloudbox_field_precalc,
3365 const Verbosity&) //verbosity)
3366{
3367 // this is only for 1D atmo!
3368 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
3369 "This method is currently only implemented for 1D atmospheres!\n")
3370
3371 // Check whether DoitInit was executed
3372 ARTS_USER_ERROR_IF (!doit_is_initialized,
3373 "Initialization method *DoitInit* has to be called before "
3374 "*cloudbox_fieldSetFromPrecalc*")
3375
3376 // check dimensions of cloudbox_field_precalc
3377 Index nf = f_grid.nelem();
3378 Index nza = za_grid.nelem();
3379 Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3380
3381 ARTS_USER_ERROR_IF (nf != cloudbox_field_precalc.nlibraries(),
3382 "cloudbox_field_precalc has wrong size in frequency dimension.\n",
3383 nf, " frequency points are expected, but cloudbox_field_precalc "
3384 "contains ", cloudbox_field_precalc.nlibraries(),
3385 "frequency points.\n")
3386 ARTS_USER_ERROR_IF (np != cloudbox_field_precalc.nvitrines(),
3387 "cloudbox_field_precalc has wrong size in pressure level dimension.\n",
3388 np, " pressure levels expected, but cloudbox_field_precalc "
3389 "contains ", cloudbox_field_precalc.nvitrines(),
3390 "pressure levels.\n")
3391 ARTS_USER_ERROR_IF (nza != cloudbox_field_precalc.npages(),
3392 "cloudbox_field_precalc has wrong size in polar angle dimension.\n",
3393 nza, " angles expected, but cloudbox_field_precalc "
3394 "contains ", cloudbox_field_precalc.npages(), "angles.\n")
3395 ARTS_USER_ERROR_IF (stokes_dim != cloudbox_field_precalc.ncols(),
3396 "cloudbox_field_precalc has wrong stokes dimension.\n"
3397 "Dimension ", stokes_dim,
3398 " expected, but cloudbox_field_precalc is dimesnion ",
3399 cloudbox_field_precalc.ncols(), ".\n")
3400
3401 // We have to update cloudbox incoming (!) field - this is because
3402 // compared to our first guess initialization, the clearsky field around might
3403 // have changed.
3404
3405 // Find which za_grid entries describe upwelling, which downwelling
3406 // radiation.
3407 Index first_upwell = 0;
3408 ARTS_ASSERT(za_grid[0] < 90.);
3409 ARTS_ASSERT(za_grid[za_grid.nelem() - 1] > 90.);
3410 while (za_grid[first_upwell] < 90.) first_upwell++;
3411
3412 Range downwell(0, first_upwell);
3413 Range upwell(first_upwell, za_grid.nelem() - first_upwell);
3414
3415 // Copy everything inside the field
3416 cloudbox_field(joker, Range(1, np - 2), 0, 0, joker, 0, joker) =
3417 cloudbox_field_precalc(joker, Range(1, np - 2), 0, 0, joker, 0, joker);
3418
3419 // At boundaries we need to be a bit careful. We shouldn't overwrite the
3420 // boundary conditions.
3421
3422 // Copy only upwelling radiation at upper boundary from precalc field
3423 // (Downwelling is "boundary condition" and has been set by DoitGetIncoming)
3424 cloudbox_field(joker, np - 1, 0, 0, upwell, 0, joker) =
3425 cloudbox_field_precalc(joker, np - 1, 0, 0, upwell, 0, joker);
3426
3427 if (cloudbox_limits[0] != 0)
3428 // Copy only downwelling radiation at lower boundary from precalc field
3429 // (Upwelling is "boundary condition" and has been set by DoitGetIncoming)
3430 cloudbox_field(joker, 0, 0, 0, downwell, 0, joker) =
3431 cloudbox_field_precalc(joker, 0, 0, 0, downwell, 0, joker);
3432 else
3433 // Copy all directions at lower boundary from precalc field
3434 // (when lower boundary at surface, the upwelling field is fixed "boundary
3435 // condition", but part of the field getting iteratively updated according
3436 // to the cloudbox-dependent surface reflection contribution)
3437 cloudbox_field(joker, 0, 0, 0, joker, 0, joker) =
3438 cloudbox_field_precalc(joker, 0, 0, 0, joker, 0, joker);
3439}
3440
3441/* Workspace method: Doxygen documentation will be auto-generated */
3443 const Vector& p_grid,
3444 const Vector& lat_grid,
3445 const Vector& lon_grid,
3446 const ArrayOfIndex& cloudbox_limits,
3447 const Index& atmosphere_dim,
3448 const Index& cloudbox_on,
3449 const Index& doit_is_initialized,
3450 const Index& all_frequencies,
3451 const Verbosity& verbosity) {
3453
3454 // Don't do anything if there's no cloudbox defined.
3455 if (!cloudbox_on) return;
3456
3457 // Check whether DoitInit was executed
3458 ARTS_USER_ERROR_IF (!doit_is_initialized,
3459 "Initialization method *DoitInit* has to be called before "
3460 "*cloudbox_fieldSetClearsky*.")
3461
3462 out2
3463 << " Interpolate boundary clearsky field to obtain the initial field.\n";
3464
3465 // Initial field only needs to be calculated from clearsky field for the
3466 // first frequency. For the next frequencies the solution field from the
3467 // previous frequencies is used.
3468 if (atmosphere_dim == 1) {
3469 const Index nf = all_frequencies ? cloudbox_field.nlibraries() : 1;
3470
3471 for (Index f_index = 0; f_index < nf; f_index++) {
3472 Index N_za = cloudbox_field.npages();
3473 Index N_aa = cloudbox_field.nrows();
3474 Index N_i = cloudbox_field.ncols();
3475
3476 //1. interpolation - pressure grid
3477
3478 /*the old grid is having only two elements, corresponding to the
3479 cloudbox_limits and the new grid have elements corresponding to
3480 all grid points inside the cloudbox plus the cloud_box_limits*/
3481
3482 ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3483
3484 p2gridpos(p_gp,
3485 p_grid[Range(cloudbox_limits[0],
3486 2,
3487 (cloudbox_limits[1] - cloudbox_limits[0]))],
3488 p_grid[Range(cloudbox_limits[0],
3489 (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3490
3491 Matrix itw((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3492 interpweights(itw, p_gp);
3493
3494 Tensor6 scat_i_p(2, 1, 1, N_za, 1, N_i);
3495 scat_i_p(0, joker, joker, joker, joker, joker) =
3496 cloudbox_field(f_index, 0, joker, joker, joker, joker, joker);
3497 scat_i_p(1, joker, joker, joker, joker, joker) =
3498 cloudbox_field(f_index,
3499 cloudbox_field.nvitrines() - 1,
3500 joker,
3501 joker,
3502 joker,
3503 joker,
3504 joker);
3505
3506 for (Index za_index = 0; za_index < N_za; ++za_index) {
3507 for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3508 for (Index i = 0; i < N_i; ++i) {
3509 VectorView target_field = cloudbox_field(
3510 f_index, Range(joker), 0, 0, za_index, aa_index, i);
3511
3512 ConstVectorView source_field =
3513 scat_i_p(Range(joker), 0, 0, za_index, aa_index, i);
3514
3515 interp(target_field, itw, source_field, p_gp);
3516 }
3517 }
3518 }
3519 }
3520 } else if (atmosphere_dim == 3) {
3521 ARTS_USER_ERROR_IF (all_frequencies == false,
3522 "Error in cloudbox_fieldSetClearsky: For 3D "
3523 "all_frequencies option is not implemented \n");
3524
3525 for (Index f_index = 0; f_index < cloudbox_field.nvitrines(); f_index++) {
3526 Index N_p = cloudbox_field.nvitrines();
3527 Index N_lat = cloudbox_field.nshelves();
3528 Index N_lon = cloudbox_field.nbooks();
3529 Index N_za = cloudbox_field.npages();
3530 Index N_aa = cloudbox_field.nrows();
3531 Index N_i = cloudbox_field.ncols();
3532
3533 Tensor6 scat_i_p(2, N_lat, N_lon, N_za, N_aa, N_i);
3534 scat_i_p(0, joker, joker, joker, joker, joker) =
3535 cloudbox_field(f_index, 0, joker, joker, joker, joker, joker);
3536 scat_i_p(1, joker, joker, joker, joker, joker) =
3537 cloudbox_field(f_index, N_p - 1, joker, joker, joker, joker, joker);
3538
3539 Tensor6 scat_i_lat(N_p, 2, N_lon, N_za, N_aa, N_i);
3540 scat_i_lat(joker, 0, joker, joker, joker, joker) =
3541 cloudbox_field(f_index, joker, 0, joker, joker, joker, joker);
3542 scat_i_lat(joker, 1, joker, joker, joker, joker) =
3543 cloudbox_field(f_index, joker, N_lat - 1, joker, joker, joker, joker);
3544
3545 Tensor6 scat_i_lon(N_p, N_lat, 2, N_za, N_aa, N_i);
3546 scat_i_lon(joker, joker, 0, joker, joker, joker) =
3547 cloudbox_field(f_index, joker, joker, 0, joker, joker, joker);
3548 scat_i_lon(joker, joker, 1, joker, joker, joker) =
3549 cloudbox_field(f_index, joker, joker, N_lon - 1, joker, joker, joker);
3550
3551 //1. interpolation - pressure grid, latitude grid and longitude grid
3552
3553 ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3554 ArrayOfGridPos lat_gp((cloudbox_limits[3] - cloudbox_limits[2]) + 1);
3555 ArrayOfGridPos lon_gp((cloudbox_limits[5] - cloudbox_limits[4]) + 1);
3556
3557 /*the old grid is having only two elements, corresponding to the
3558 cloudbox_limits and the new grid have elements corresponding to
3559 all grid points inside the cloudbox plus the cloud_box_limits*/
3560
3561 p2gridpos(p_gp,
3562 p_grid[Range(cloudbox_limits[0],
3563 2,
3564 (cloudbox_limits[1] - cloudbox_limits[0]))],
3565 p_grid[Range(cloudbox_limits[0],
3566 (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3567 gridpos(lat_gp,
3568 lat_grid[Range(cloudbox_limits[2],
3569 2,
3570 (cloudbox_limits[3] - cloudbox_limits[2]))],
3571 lat_grid[Range(cloudbox_limits[2],
3572 (cloudbox_limits[3] - cloudbox_limits[2]) + 1)]);
3573 gridpos(lon_gp,
3574 lon_grid[Range(cloudbox_limits[4],
3575 2,
3576 (cloudbox_limits[5] - cloudbox_limits[4]))],
3577 lon_grid[Range(cloudbox_limits[4],
3578 (cloudbox_limits[5] - cloudbox_limits[4]) + 1)]);
3579
3580 //interpolation weights corresponding to pressure, latitude and
3581 //longitude grids.
3582
3583 Matrix itw_p((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3584 Matrix itw_lat((cloudbox_limits[3] - cloudbox_limits[2]) + 1, 2);
3585 Matrix itw_lon((cloudbox_limits[5] - cloudbox_limits[4]) + 1, 2);
3586
3587 interpweights(itw_p, p_gp);
3588 interpweights(itw_lat, lat_gp);
3589 interpweights(itw_lon, lon_gp);
3590
3591 // interpolation - pressure grid
3592 for (Index lat_index = 0;
3593 lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3594 ++lat_index) {
3595 for (Index lon_index = 0;
3596 lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3597 ++lon_index) {
3598 for (Index za_index = 0; za_index < N_za; ++za_index) {
3599 for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3600 for (Index i = 0; i < N_i; ++i) {
3601 VectorView target_field = cloudbox_field(f_index,
3602 Range(joker),
3603 lat_index,
3604 lon_index,
3605 za_index,
3606 aa_index,
3607 i);
3608
3609 ConstVectorView source_field = scat_i_p(
3610 Range(joker), lat_index, lon_index, za_index, aa_index, i);
3611
3612 interp(target_field, itw_p, source_field, p_gp);
3613 }
3614 }
3615 }
3616 }
3617 }
3618 //interpolation latitude
3619 for (Index p_index = 0;
3620 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3621 ++p_index) {
3622 for (Index lon_index = 0;
3623 lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3624 ++lon_index) {
3625 for (Index za_index = 0; za_index < N_za; ++za_index) {
3626 for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3627 for (Index i = 0; i < N_i; ++i) {
3628 VectorView target_field = cloudbox_field(f_index,
3629 p_index,
3630 Range(joker),
3631 lon_index,
3632 za_index,
3633 aa_index,
3634 i);
3635
3636 ConstVectorView source_field = scat_i_lat(
3637 p_index, Range(joker), lon_index, za_index, aa_index, i);
3638
3639 interp(target_field, itw_lat, source_field, lat_gp);
3640 }
3641 }
3642 }
3643 }
3644 }
3645 //interpolation -longitude
3646 for (Index p_index = 0;
3647 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3648 ++p_index) {
3649 for (Index lat_index = 0;
3650 lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3651 ++lat_index) {
3652 for (Index za_index = 0; za_index < N_za; ++za_index) {
3653 for (Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3654 for (Index i = 0; i < N_i; ++i) {
3655 VectorView target_field = cloudbox_field(f_index,
3656 p_index,
3657 lat_index,
3658 Range(joker),
3659 za_index,
3660 aa_index,
3661 i);
3662
3663 ConstVectorView source_field = scat_i_lon(
3664 p_index, lat_index, Range(joker), za_index, aa_index, i);
3665
3666 interp(target_field, itw_lon, source_field, lon_gp);
3667 }
3668 }
3669 }
3670 }
3671 } //end of interpolation
3672 } // end of frequency loop
3673 } //ends atmosphere_dim = 3
3674}
3675
3676/* Workspace method: Doxygen documentation will be auto-generated */
3677void cloudbox_fieldSetConst( //WS Output:
3678 Tensor7& cloudbox_field,
3679 //WS Input:
3680 const Vector& p_grid,
3681 const Vector& lat_grid,
3682 const Vector& lon_grid,
3683 const ArrayOfIndex& cloudbox_limits,
3684 const Index& atmosphere_dim,
3685 const Index& stokes_dim,
3686 // Keyword
3687 const Vector& cloudbox_field_values,
3688 const Verbosity& verbosity) {
3690
3691 Tensor6 cloudbox_field_mono(cloudbox_field.nvitrines(),
3692 cloudbox_field.nshelves(),
3693 cloudbox_field.nbooks(),
3694 cloudbox_field.npages(),
3695 cloudbox_field.nrows(),
3696 cloudbox_field.ncols());
3697
3698 for (Index f_index = 0; f_index < cloudbox_field.nlibraries(); f_index++) {
3699 cloudbox_field_mono =
3700 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker);
3701
3702 cloudbox_field_monoSetConst(cloudbox_field_mono,
3703 p_grid,
3704 lat_grid,
3705 lon_grid,
3706 cloudbox_limits,
3707 atmosphere_dim,
3708 stokes_dim,
3709 cloudbox_field_values,
3710 verbosity);
3711
3712 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) =
3713 cloudbox_field_mono;
3714 }
3715}
3716
3717/* Workspace method: Doxygen documentation will be auto-generated */
3719 Tensor7& cloudbox_field,
3720 //WS Input:
3721 const Vector& p_grid,
3722 const Vector& lat_grid,
3723 const Vector& lon_grid,
3724 const ArrayOfIndex& cloudbox_limits,
3725 const Index& atmosphere_dim,
3726 const Index& stokes_dim,
3727 // Keyword
3728 const Matrix& cloudbox_field_values,
3729 const Verbosity& verbosity) {
3731
3732 ARTS_USER_ERROR_IF (cloudbox_field.nlibraries() != cloudbox_field_values.nrows(),
3733 "Number of rows in *cloudbox_field_values* has to be equal"
3734 " the frequency dimension of *cloudbox_field*.");
3735
3736 Tensor6 cloudbox_field_mono(cloudbox_field.nvitrines(),
3737 cloudbox_field.nshelves(),
3738 cloudbox_field.nbooks(),
3739 cloudbox_field.npages(),
3740 cloudbox_field.nrows(),
3741 cloudbox_field.ncols());
3742
3743 for (Index f_index = 0; f_index < cloudbox_field.nlibraries(); f_index++) {
3744 cloudbox_field_mono =
3745 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker);
3746
3747 cloudbox_field_monoSetConst(cloudbox_field_mono,
3748 p_grid,
3749 lat_grid,
3750 lon_grid,
3751 cloudbox_limits,
3752 atmosphere_dim,
3753 stokes_dim,
3754 cloudbox_field_values(f_index, joker),
3755 verbosity);
3756
3757 cloudbox_field(f_index, joker, joker, joker, joker, joker, joker) =
3758 cloudbox_field_mono;
3759 }
3760}
3761
3762/* Workspace method: Doxygen documentation will be auto-generated */
3764 Tensor6& cloudbox_field_mono,
3765 //WS Input:
3766 const Vector& p_grid,
3767 const Vector& lat_grid,
3768 const Vector& lon_grid,
3769 const ArrayOfIndex& cloudbox_limits,
3770 const Index& atmosphere_dim,
3771 const Index& stokes_dim,
3772 // Keyword
3773 const Vector& cloudbox_field_values,
3774 const Verbosity& verbosity) {
3777
3778 out2 << " Set initial field to constant values: " << cloudbox_field_values
3779 << "\n";
3780
3781 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3782
3783 // Grids have to be adapted to atmosphere_dim.
3784 chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3785
3786 // Check the input:
3787 ARTS_USER_ERROR_IF (stokes_dim < 0 || stokes_dim > 4,
3788 "The dimension of stokes vector must be"
3789 "1,2,3, or 4.");
3790
3791 ARTS_USER_ERROR_IF (stokes_dim != cloudbox_field_values.nelem(),
3792 "Length of *cloudbox_field_values* has to be equal"
3793 " *stokes_dim*.");
3794 ARTS_USER_ERROR_IF (cloudbox_limits.nelem() != 2 * atmosphere_dim,
3795 "*cloudbox_limits* is a vector which contains the"
3796 "upper and lower limit of the cloud for all "
3797 "atmospheric dimensions. So its dimension must"
3798 "be 2 x *atmosphere_dim*.");
3799
3800 for (Index i = 0; i < stokes_dim; i++) {
3801 cloudbox_field_mono(joker, joker, joker, joker, joker, i) =
3802 cloudbox_field_values[i];
3803 }
3804}
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:23109
void doit_scat_field_agendaExecute(Workspace &ws, Tensor6 &doit_scat_field, const Tensor6 &cloudbox_field_mono, const Agenda &input_agenda)
Definition: auto_md.cc:23073
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:23034
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 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:22824
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, 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:23520
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:22992
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:24022
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 &verbosity)
WORKSPACE METHOD: pha_matCalc.
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:44
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
Index ncols() const
Returns the number of columns.
Definition: matpackVI.cc:57
Index npages() const
Returns the number of pages.
Definition: matpackVI.cc:51
Index nrows() const
Returns the number of rows.
Definition: matpackVI.cc:54
Index nshelves() const
Returns the number of shelves.
Definition: matpackVI.cc:45
Index nbooks() const
Returns the number of books.
Definition: matpackVI.cc:48
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVI.cc:42
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:48
Index nrows() const
Returns the number of rows.
Definition: matpackVII.cc:57
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVII.cc:45
Index npages() const
Returns the number of pages.
Definition: matpackVII.cc:54
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:60
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:42
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:51
A constant view of a Vector.
Definition: matpackI.h:489
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:46
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The Matrix class.
Definition: matpackI.h:1225
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The range class.
Definition: matpackI.h:165
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:339
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:658
The Tensor4 class.
Definition: matpackIV.h:421
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1058
The Tensor5 class.
Definition: matpackV.h:506
The Tensor6 class.
Definition: matpackVI.h:1088
void resize(Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVI.cc:2171
The Tensor7 class.
Definition: matpackVII.h:2382
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5482
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
#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:1842
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:1948
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:1609
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:1102
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:560
#define abs(x)
Linear algebra functions.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
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:2748
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:2664
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:3194
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:1994
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:2723
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:2314
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:3718
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:3357
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:2858
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:3763
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:3677
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:3442
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:2043
Template functions for general supergeneric ws methods.
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:290
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:334
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition: messages.h:205
#define CREATE_OUT3
Definition: messages.h:207
#define CREATE_OUT2
Definition: messages.h:206
#define CREATE_OUT0
Definition: messages.h:204
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:39
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:747
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
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.h:48
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:181
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.h:160
@ FILE_TYPE_ASCII
Definition: xml_io_base.h:41