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