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