ARTS 2.5.11 (git: 725533f0)
doit.cc
Go to the documentation of this file.
1
11/*===========================================================================
12 === External declarations
13 ===========================================================================*/
14
15#include "doit.h"
16#include <cmath>
17#include <cstdlib>
18#include <iostream>
19#include <stdexcept>
20#include "agenda_class.h"
21#include "array.h"
22#include "arts_constants.h"
23#include "arts_conversions.h"
24#include "auto_md.h"
25#include "check_input.h"
26#include "cloudbox.h"
27#include "geodetic.h"
28#include "lin_alg.h"
29#include "logic.h"
30#include "math_funcs.h"
31#include "matpack_data.h"
32#include "messages.h"
33#include "physics_funcs.h"
34#include "ppath.h"
35#include "propagationmatrix.h"
36#include "rte.h"
37#include "sorting.h"
38#include "special_interp.h"
39#include "xml_io.h"
40
41inline constexpr Numeric PI=Constant::pi;
42inline constexpr Numeric RAD2DEG=Conversion::rad2deg(1);
43
44//FIXME function name of 'rte_step_doit_replacement' should be replaced by
45//proper name
46void rte_step_doit_replacement( //Output and Input:
47 VectorView stokes_vec,
48 MatrixView trans_mat,
49 //Input
50 const PropagationMatrix& ext_mat_av,
51 const StokesVector& abs_vec_av,
52 ConstVectorView sca_vec_av,
53 const Numeric& lstep,
54 const Numeric& rtp_planck_value,
55 const bool& trans_is_precalc) {
56 //Stokes dimension:
57 Index stokes_dim = stokes_vec.nelem();
58
59 //Test sizes
60 ARTS_ASSERT(ext_mat_av.NumberOfFrequencies() == 1 and
61 abs_vec_av.NumberOfFrequencies() == 1);
62
63 //Check inputs:
64 ARTS_ASSERT(is_size(trans_mat, stokes_dim, stokes_dim));
65 ARTS_ASSERT(stokes_dim == ext_mat_av.StokesDimensions() and
66 stokes_dim == abs_vec_av.StokesDimensions());
67 ARTS_ASSERT(is_size(sca_vec_av, stokes_dim));
68 ARTS_ASSERT(rtp_planck_value >= 0);
69 ARTS_ASSERT(lstep >= 0);
70 //ARTS_ASSERT (not ext_mat_av.AnySingular()); This is ARTS_ASSERTed at a later time in this version...
71
72 // Check, if only the first component of abs_vec is non-zero:
73// const bool unpol_abs_vec = abs_vec_av.IsUnpolarized(0);
74
75// bool unpol_sca_vec = true;
76
77// for (Index i = 1; i < stokes_dim; i++)
78// if (sca_vec_av[i] != 0) unpol_sca_vec = false;
79
80 // Calculate transmission by general function, if not precalculated
81// Index extmat_case = 0;
82 if (!trans_is_precalc) {
83 compute_transmission_matrix_from_averaged_matrix_at_frequency(
84 trans_mat, lstep, ext_mat_av, 0);
85 }
86
87 //--- Scalar case: ---------------------------------------------------------
88 if (stokes_dim == 1) {
89 stokes_vec[0] = stokes_vec[0] * trans_mat(0, 0) +
90 (abs_vec_av.Kjj()[0] * rtp_planck_value + sca_vec_av[0]) /
91 ext_mat_av.Kjj()[0] * (1 - trans_mat(0, 0));
92 }
93
94 //--- Vector case: ---------------------------------------------------------
95
96 // We have here two cases, diagonal or non-diagonal ext_mat_gas
97 // For diagonal ext_mat_gas, we expect abs_vec_gas to only have a
98 // non-zero value in position 1.
99
100 //- Unpolarised
101// else if (extmat_case == 1 && unpol_abs_vec && unpol_sca_vec) {
102// const Numeric invK = 1.0 / ext_mat_av.Kjj()[0];
103// // Stokes dim 1
104// stokes_vec[0] = stokes_vec[0] * trans_mat(0, 0) +
105// (abs_vec_av.Kjj()[0] * rtp_planck_value + sca_vec_av[0]) *
106// invK * (1 - trans_mat(0, 0));
107//
108// // Stokes dims > 1
109// for (Index i = 1; i < stokes_dim; i++) {
110// stokes_vec[i] = stokes_vec[i] * trans_mat(i, i) +
111// sca_vec_av[i] * invK * (1 - trans_mat(i, i));
112// }
113// }
114
115 //- General case
116 else {
117 Matrix invK(stokes_dim, stokes_dim);
118 ext_mat_av.MatrixInverseAtPosition(invK);
119
120 Vector source{abs_vec_av.VectorAtPosition()};
121 source *= rtp_planck_value;
122
123 for (Index i = 0; i < stokes_dim; i++)
124 source[i] += sca_vec_av[i]; // b = abs_vec * B + sca_vec
125
126 // solve K^(-1)*b = x
127 Vector x(stokes_dim);
128 mult(x, invK, source);
129
130 Vector term1(stokes_dim);
131 Vector term2(stokes_dim);
132
133 Matrix ImT(stokes_dim, stokes_dim);
134 id_mat(ImT);
135 ImT -= trans_mat;
136 mult(term2, ImT, x); // term2: second term of the solution of the RTE with
137 //fixed scattered field
138
139 // term1: first term of solution of the RTE with fixed scattered field
140 mult(term1, trans_mat, stokes_vec);
141
142 for (Index i = 0; i < stokes_dim; i++)
143 stokes_vec[i] = term1[i] + term2[i]; // Compute the new Stokes Vector
144 }
145}
146
148 // Output and Input:
149 Tensor5View ext_mat_field,
150 Tensor4View abs_vec_field,
151 // Input:
152 const Agenda& spt_calc_agenda,
153 const Index& za_index,
154 const Index& aa_index,
155 const ArrayOfIndex& cloudbox_limits,
156 ConstTensor3View t_field,
157 ConstTensor4View pnd_field,
158 const Verbosity& verbosity) {
160
161 // Input variables are checked in the WSMs i_fieldUpdateSeqXXX, from
162 // where this function is called.
163
164 out3 << "Calculate scattering properties in cloudbox \n";
165
166 const Index atmosphere_dim = cloudbox_limits.nelem() / 2;
167 const Index N_se = pnd_field.nbooks();
168 const Index stokes_dim = ext_mat_field.ncols();
169
170 ARTS_ASSERT(atmosphere_dim == 1 || atmosphere_dim == 3);
171 ARTS_ASSERT(ext_mat_field.ncols() == ext_mat_field.nrows() &&
172 ext_mat_field.ncols() == abs_vec_field.ncols());
173
174 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
175
176 // If atmosohere_dim == 1
177 Index Nlat_cloud = 1;
178 Index Nlon_cloud = 1;
179
180 if (atmosphere_dim == 3) {
181 Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
182 Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
183 }
184
185 // Initialize ext_mat(_spt), abs_vec(_spt)
186 // Resize and initialize variables for storing optical properties
187 // of all scattering elements.
188 ArrayOfStokesVector abs_vec_spt_local(N_se);
189 for (auto& av : abs_vec_spt_local) {
190 av = StokesVector(1, stokes_dim);
191 av.SetZero();
192 }
193 ArrayOfPropagationMatrix ext_mat_spt_local(N_se);
194 for (auto& pm : ext_mat_spt_local) {
195 pm = PropagationMatrix(1, stokes_dim);
196 pm.SetZero();
197 }
198
199 StokesVector abs_vec_local;
200 PropagationMatrix ext_mat_local;
201 Numeric rtp_temperature_local;
202
203 // Calculate ext_mat, abs_vec for all points inside the cloudbox.
204 // sca_vec can be obtained from the workspace variable doit_scat_field.
205 // As we need the average for each layer, it makes sense to calculte
206 // the coefficients once and store them in an array instead of
207 // calculating at each point the coefficient of the point above and
208 // the point below.
209 // To use special interpolation functions for atmospheric fields we
210 // use ext_mat_field and abs_vec_field:
211
212 // Loop over all positions inside the cloudbox defined by the
213 // cloudbox_limits.
214 for (Index scat_p_index_local = 0; scat_p_index_local < Np_cloud;
215 scat_p_index_local++) {
216 for (Index scat_lat_index_local = 0; scat_lat_index_local < Nlat_cloud;
217 scat_lat_index_local++) {
218 for (Index scat_lon_index_local = 0; scat_lon_index_local < Nlon_cloud;
219 scat_lon_index_local++) {
220 if (atmosphere_dim == 1)
221 rtp_temperature_local =
222 t_field(scat_p_index_local + cloudbox_limits[0], 0, 0);
223 else
224 rtp_temperature_local =
225 t_field(scat_p_index_local + cloudbox_limits[0],
226 scat_lat_index_local + cloudbox_limits[2],
227 scat_lon_index_local + cloudbox_limits[4]);
228
229 //Calculate optical properties for individual scattering elements:
230 //( Execute agendas silently. )
232 ext_mat_spt_local,
233 abs_vec_spt_local,
234 scat_p_index_local,
235 scat_lat_index_local,
236 scat_lon_index_local,
237 rtp_temperature_local,
238 za_index,
239 aa_index,
240 spt_calc_agenda);
241 /*
242// so far missing here (accessed through workspace within agenda):
243// - scat_data
244// - za_grid, aa_grid
245// - f_index
246 opt_prop_sptFromScat_data(ext_mat_spt_local, abs_vec_spt_local,
247 scat_data, 1,
248 za_grid, aa_grid,
249 za_index, aa_index,
250 f_index,
251 rtp_temperature_local,
252 pnd_field,
253 scat_p_index_local,
254 scat_lat_index_local,
255 scat_lon_index_local,
256 verbosity);
257*/
258
259 opt_prop_bulkCalc(ext_mat_local,
260 abs_vec_local,
261 ext_mat_spt_local,
262 abs_vec_spt_local,
263 Tensor4{pnd_field},
264 scat_p_index_local,
265 scat_lat_index_local,
266 scat_lon_index_local,
267 verbosity);
268
269 // Store coefficients in arrays for the whole cloudbox.
270 abs_vec_field(scat_p_index_local,
271 scat_lat_index_local,
272 scat_lon_index_local,
273 joker) = abs_vec_local.VectorAtPosition();
274
275 ext_mat_local.MatrixAtPosition(ext_mat_field(scat_p_index_local,
276 scat_lat_index_local,
277 scat_lon_index_local,
278 joker,
279 joker));
280 }
281 }
282 }
283}
284
286 // Input and output
287 Tensor6View cloudbox_field_mono,
288 // ppath_step_agenda:
289 const Index& p_index,
290 const Index& za_index,
291 ConstVectorView za_grid,
292 const ArrayOfIndex& cloudbox_limits,
293 ConstTensor6View doit_scat_field,
294 // Calculate scalar gas absorption:
295 const Agenda& propmat_clearsky_agenda,
296 ConstTensor4View vmr_field,
297 // Propagation path calculation:
298 const Agenda& ppath_step_agenda,
299 const Numeric& ppath_lmax,
300 const Numeric& ppath_lraytrace,
301 ConstVectorView p_grid,
302 ConstTensor3View z_field,
303 ConstVectorView refellipsoid,
304 // Calculate thermal emission:
305 ConstTensor3View t_field,
306 ConstVectorView f_grid,
307 const Index& f_index,
308 //particle optical properties
309 ConstTensor5View ext_mat_field,
310 ConstTensor4View abs_vec_field,
311 const Agenda& surface_rtprop_agenda,
312 //const Agenda& surface_rtprop_agenda,
313 const Index& scat_za_interp,
314 const Verbosity& verbosity) {
315 Ppath ppath_step;
316 // Input variables are checked in the WSMs i_fieldUpdateSeqXXX, from
317 // where this function is called.
318
319 //Initialize ppath for 1D.
320 ppath_init_structure(ppath_step, 1, 1);
321 // See documentation of ppath_init_structure for understanding
322 // the parameters.
323
324 // Assign value to ppath.pos:
325 ppath_step.pos(0, 0) = z_field(p_index, 0, 0);
326 ppath_step.r[0] = refellipsoid[0] + z_field(p_index, 0, 0);
327
328 // Define the direction:
329 ppath_step.los(0, 0) = za_grid[za_index];
330
331 // Define the grid positions:
332 ppath_step.gp_p[0].idx = p_index;
333 ppath_step.gp_p[0].fd[0] = 0;
334 ppath_step.gp_p[0].fd[1] = 1;
335
336 // Call ppath_step_agenda:
338 ppath_step,
339 ppath_lmax,
340 ppath_lraytrace,
341 Vector(1, f_grid[f_index]),
342 ppath_step_agenda);
343
344 // Check whether the next point is inside or outside the
345 // cloudbox. Only if the next point lies inside the
346 // cloudbox a radiative transfer step caclulation has to
347 // be performed.
348
349 if ((cloudbox_limits[0] <= ppath_step.gp_p[1].idx &&
350 cloudbox_limits[1] > ppath_step.gp_p[1].idx) ||
351 (cloudbox_limits[1] == ppath_step.gp_p[1].idx &&
352 abs(ppath_step.gp_p[1].fd[0]) < 1e-6)) {
353 // Stokes dimension
354 const Index stokes_dim = cloudbox_field_mono.ncols();
355 // Number of species
356 const Index N_species = vmr_field.nbooks();
357
358 // Ppath_step normally has 2 points, the starting
359 // point and the intersection point.
360 // But there can be points in between, when a maximum
361 // lstep is given. We have to interpolate on all the
362 // points in the ppath_step.
363
364 // Initialize variables for interpolated values
365 Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.np, 0.);
366 Matrix abs_vec_int(stokes_dim, ppath_step.np, 0.);
367 Matrix sca_vec_int(stokes_dim, ppath_step.np, 0.);
368 Matrix cloudbox_field_mono_int(stokes_dim, ppath_step.np, 0.);
369 Vector t_int(ppath_step.np, 0.);
370 Matrix vmr_list_int(N_species, ppath_step.np, 0.);
371 Vector p_int(ppath_step.np, 0.);
372
373 interp_cloud_coeff1D(ext_mat_int,
374 abs_vec_int,
375 sca_vec_int,
376 cloudbox_field_mono_int,
377 t_int,
378 vmr_list_int,
379 p_int,
380 ext_mat_field,
381 abs_vec_field,
382 doit_scat_field,
383 cloudbox_field_mono,
384 t_field,
385 vmr_field,
386 p_grid,
387 ppath_step,
388 cloudbox_limits,
389 za_grid,
390 scat_za_interp,
391 verbosity);
392
393 // ppath_what_background(ppath_step) tells the
394 // radiative background. More information in the
395 // function get_iy_of_background.
396 // if there is no background we proceed the RT
397 Index bkgr = ppath_what_background(ppath_step);
398
399 // Radiative transfer from one layer to the next, starting
400 // at the intersection with the next layer and propagating
401 // to the considered point.
403 cloudbox_field_mono,
404 propmat_clearsky_agenda,
405 ppath_step,
406 t_int,
407 vmr_list_int,
408 ext_mat_int,
409 abs_vec_int,
410 sca_vec_int,
411 cloudbox_field_mono_int,
412 p_int,
413 cloudbox_limits,
414 f_grid,
415 f_index,
416 p_index,
417 0,
418 0,
419 za_index,
420 0,
421 verbosity);
422
423 // bkgr=2 indicates that the background is the surface
424 if (bkgr == 2) {
425 // cout << "hit surface "<< ppath_step.gp_p << endl;
427 cloudbox_field_mono,
428 surface_rtprop_agenda,
429 f_grid,
430 f_index,
431 stokes_dim,
432 ppath_step,
433 cloudbox_limits,
434 za_grid,
435 za_index);
436 }
437
438 } //end if inside cloudbox
439}
440
442 // Output
443 Tensor6View cloudbox_field_mono,
444 // ppath_step_agenda:
445 const Index& p_index,
446 const Index& za_index,
447 ConstVectorView za_grid,
448 const ArrayOfIndex& cloudbox_limits,
449 ConstTensor6View cloudbox_field_mono_old,
450 ConstTensor6View doit_scat_field,
451 // Calculate scalar gas absorption:
452 const Agenda& propmat_clearsky_agenda,
453 ConstTensor4View vmr_field,
454 // Propagation path calculation:
455 const Agenda& ppath_step_agenda,
456 const Numeric& ppath_lmax,
457 const Numeric& ppath_lraytrace,
458 ConstVectorView p_grid,
459 ConstTensor3View z_field,
460 ConstVectorView refellipsoid,
461 // Calculate thermal emission:
462 ConstTensor3View t_field,
463 ConstVectorView f_grid,
464 // used for surface ?
465 const Index& f_index,
466 //particle optical properties
467 ConstTensor5View ext_mat_field,
468 ConstTensor4View abs_vec_field,
469 const Agenda& surface_rtprop_agenda,
470 const Index& scat_za_interp,
471 const Verbosity& verbosity) {
472 Ppath ppath_step;
473 // Input variables are checked in the WSMs i_fieldUpdateSeqXXX, from
474 // where this function is called.
475
476 //Initialize ppath for 1D.
477 ppath_init_structure(ppath_step, 1, 1);
478 // See documentation of ppath_init_structure for understanding
479 // the parameters.
480
481 // Assign value to ppath.pos:
482 ppath_step.pos(0, 0) = z_field(p_index, 0, 0);
483 ppath_step.r[0] = refellipsoid[0] + z_field(p_index, 0, 0);
484
485 // Define the direction:
486 ppath_step.los(0, 0) = za_grid[za_index];
487
488 // Define the grid positions:
489 ppath_step.gp_p[0].idx = p_index;
490 ppath_step.gp_p[0].fd[0] = 0;
491 ppath_step.gp_p[0].fd[1] = 1;
492
493 // Call ppath_step_agenda:
495 ppath_step,
496 ppath_lmax,
497 ppath_lraytrace,
498 Vector(1, f_grid[f_index]),
499 ppath_step_agenda);
500
501 // Check whether the next point is inside or outside the
502 // cloudbox. Only if the next point lies inside the
503 // cloudbox a radiative transfer step caclulation has to
504 // be performed.
505
506 if ((cloudbox_limits[0] <= ppath_step.gp_p[1].idx &&
507 cloudbox_limits[1] > ppath_step.gp_p[1].idx) ||
508 (cloudbox_limits[1] == ppath_step.gp_p[1].idx &&
509 abs(ppath_step.gp_p[1].fd[0]) < 1e-6)) {
510 // Stokes dimension
511 const Index stokes_dim = cloudbox_field_mono.ncols();
512 // Number of species
513 const Index N_species = vmr_field.nbooks();
514
515 // Ppath_step normally has 2 points, the starting
516 // point and the intersection point.
517 // But there can be points in between, when a maximum
518 // lstep is given. We have to interpolate on all the
519 // points in the ppath_step.
520
521 // Initialize variables for interpolated values
522 Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.np, 0.);
523 Matrix abs_vec_int(stokes_dim, ppath_step.np, 0.);
524 Matrix sca_vec_int(stokes_dim, ppath_step.np, 0.);
525 Matrix cloudbox_field_mono_int(stokes_dim, ppath_step.np, 0.);
526 Vector t_int(ppath_step.np, 0.);
527 Matrix vmr_list_int(N_species, ppath_step.np, 0.);
528 Vector p_int(ppath_step.np, 0.);
529
530 interp_cloud_coeff1D(ext_mat_int,
531 abs_vec_int,
532 sca_vec_int,
533 cloudbox_field_mono_int,
534 t_int,
535 vmr_list_int,
536 p_int,
537 ext_mat_field,
538 abs_vec_field,
539 doit_scat_field,
540 cloudbox_field_mono_old,
541 t_field,
542 vmr_field,
543 p_grid,
544 ppath_step,
545 cloudbox_limits,
546 za_grid,
547 scat_za_interp,
548 verbosity);
549
550 // ppath_what_background(ppath_step) tells the
551 // radiative background. More information in the
552 // function get_iy_of_background.
553 // if there is no background we proceed the RT
554 Index bkgr = ppath_what_background(ppath_step);
555
556 // if 0, there is no background
557 // do this in any case. cause we need downwelling cloudbox_field_mono
558 // at the surface for calculating surface scattering
559
560 // Radiative transfer from one layer to the next, starting
561 // at the intersection with the next layer and propagating
562 // to the considered point.
564 cloudbox_field_mono,
565 propmat_clearsky_agenda,
566 ppath_step,
567 t_int,
568 vmr_list_int,
569 ext_mat_int,
570 abs_vec_int,
571 sca_vec_int,
572 cloudbox_field_mono_int,
573 p_int,
574 cloudbox_limits,
575 f_grid,
576 f_index,
577 p_index,
578 0,
579 0,
580 za_index,
581 0,
582 verbosity);
583
584 if (bkgr == 2) {
586 cloudbox_field_mono,
587 surface_rtprop_agenda,
588 f_grid,
589 f_index,
590 stokes_dim,
591 ppath_step,
592 cloudbox_limits,
593 za_grid,
594 za_index);
595 }
596 } //end if inside cloudbox
597}
598
600 Tensor6View cloudbox_field_mono,
601 const Index& p_index,
602 const Index& za_index,
603 ConstVectorView za_grid,
604 const ArrayOfIndex& cloudbox_limits,
605 ConstTensor6View doit_scat_field,
606 // Calculate scalar gas absorption:
607 const Agenda& propmat_clearsky_agenda,
608 ConstTensor4View vmr_field,
609 // Propagation path calculation:
610 ConstVectorView p_grid,
611 ConstTensor3View z_field,
612 // Calculate thermal emission:
613 ConstTensor3View t_field,
614 ConstVectorView f_grid,
615 const Index& f_index,
616 //particle opticla properties
617 ConstTensor5View ext_mat_field,
618 ConstTensor4View abs_vec_field,
619 const Verbosity& verbosity) {
621
622 const Index N_species = vmr_field.nbooks();
623 const Index stokes_dim = cloudbox_field_mono.ncols();
624 const Index atmosphere_dim = 1;
625 PropagationMatrix propmat_clearsky;
626 PropagationMatrix ext_mat;
627 StokesVector abs_vec;
628 Matrix matrix_tmp(stokes_dim, stokes_dim);
629 Vector vector_tmp(stokes_dim);
630 Vector rtp_vmr(N_species, 0.);
631 EnergyLevelMap rtp_nlte_dummy;
632 Vector sca_vec_av(stokes_dim, 0);
633
634 // Radiative transfer from one layer to the next, starting
635 // at the intersection with the next layer and propagating
636 // to the considered point.
637 Vector stokes_vec(stokes_dim);
638 Index bkgr;
639 if ((p_index == 0) && (za_grid[za_index] > 90)) {
640 bkgr = 2;
641 } else {
642 bkgr = 0;
643 }
644
645 // if 0, there is no background
646 if (bkgr == 0) {
647 if (za_grid[za_index] <= 90.0) {
648 stokes_vec = cloudbox_field_mono(
649 p_index - cloudbox_limits[0] + 1, 0, 0, za_index, 0, joker);
650 Numeric z_field_above = z_field(p_index + 1, 0, 0);
651 Numeric z_field_0 = z_field(p_index, 0, 0);
652
653 Numeric cos_theta, lstep;
654 if (za_grid[za_index] == 90.0) {
655 //cos_theta = cos((za_grid[za_index] -1)*PI/180.);
656 // cos_theta = cos(89.999999999999*PI/180.);
657 cos_theta = 1e-20;
658
659 } else {
660 cos_theta = cos(za_grid[za_index] * PI / 180.0);
661 }
662 Numeric dz = (z_field_above - z_field_0);
663
664 lstep = abs(dz / cos_theta);
665
666 // Average temperature
667 Numeric rtp_temperature =
668 0.5 * (t_field(p_index, 0, 0) + t_field(p_index + 1, 0, 0));
669 //
670 // Average pressure
671 Numeric rtp_pressure = 0.5 * (p_grid[p_index] + p_grid[p_index + 1]);
672
673 // Average vmrs
674 for (Index j = 0; j < N_species; j++)
675 rtp_vmr[j] = 0.5 * (vmr_field(j, p_index, 0, 0) +
676 vmr_field(j, p_index + 1, 0, 0));
677 //
678 // Calculate scalar gas absorption and add it to abs_vec
679 // and ext_mat.
680 //
681
682 const Vector rtp_mag_dummy(3, 0);
683 const Vector ppath_los_dummy;
684
685 StokesVector nlte_dummy; //FIXME: do this right?
686 ArrayOfPropagationMatrix partial_dummy; // This is right since there should be only clearsky partials
687 ArrayOfStokesVector partial_nlte_dummy; // This is right since there should be only clearsky partials
689 propmat_clearsky,
690 nlte_dummy,
691 partial_dummy,
692 partial_nlte_dummy,
694 {},
695 Vector{f_grid[Range(f_index, 1)]},
696 rtp_mag_dummy,
697 ppath_los_dummy,
698 rtp_pressure,
699 rtp_temperature,
700 rtp_nlte_dummy,
701 rtp_vmr,
702 propmat_clearsky_agenda);
703
704 opt_prop_sum_propmat_clearsky(ext_mat, abs_vec, propmat_clearsky);
705
706 for (Index k = 0; k < stokes_dim; k++) {
707 //
708 // Averaging of sca_vec:
709 //
710 sca_vec_av[k] +=
711 0.5 *
712 (doit_scat_field(
713 p_index - cloudbox_limits[0], 0, 0, za_index, 0, k) +
714 doit_scat_field(
715 p_index - cloudbox_limits[0] + 1, 0, 0, za_index, 0, k));
716 }
717
718 //
719 // Add average particle absorption to abs_vec.
720 //
721 abs_vec.AddAverageAtPosition(
722 abs_vec_field(p_index - cloudbox_limits[0], 0, 0, joker),
723 abs_vec_field(p_index - cloudbox_limits[0] + 1, 0, 0, joker));
724
725 //
726 // Add average particle extinction to ext_mat.
727 ext_mat.AddAverageAtPosition(
728 ext_mat_field(p_index - cloudbox_limits[0], 0, 0, joker, joker),
729 ext_mat_field(p_index - cloudbox_limits[0] + 1, 0, 0, joker, joker));
730
731 // Frequency
732 Numeric f = f_grid[f_index];
733 //
734 // Calculate Planck function
735 //
736 Numeric rte_planck_value = planck(f, rtp_temperature);
737
738 // Some messages:
739 if (out3.sufficient_priority()) {
740 vector_tmp = abs_vec.VectorAtPosition();
741 ext_mat.MatrixAtPosition(matrix_tmp);
742 out3 << "-----------------------------------------\n";
743 out3 << "Input for radiative transfer step \n"
744 << "calculation inside"
745 << " the cloudbox:"
746 << "\n";
747 out3 << "Stokes vector at intersection point: \n" << stokes_vec << "\n";
748 out3 << "lstep: ..." << lstep << "\n";
749 out3 << "------------------------------------------\n";
750 out3 << "Averaged coefficients: \n";
751 out3 << "Planck function: " << rte_planck_value << "\n";
752 out3 << "Scattering vector: " << sca_vec_av << "\n";
753 out3 << "Absorption vector: " << vector_tmp << "\n";
754 out3 << "Extinction matrix: " << matrix_tmp << "\n";
755
756 ARTS_ASSERT(!is_singular(matrix_tmp));
757 }
758
759 // Radiative transfer step calculation. The Stokes vector
760 // is updated until the considered point is reached.
761 rte_step_doit_replacement(stokes_vec,
762 Matrix(stokes_dim, stokes_dim),
763 ext_mat,
764 abs_vec,
765 sca_vec_av,
766 lstep,
767 rte_planck_value);
768
769 // Assign calculated Stokes Vector to cloudbox_field_mono.
770 cloudbox_field_mono(
771 p_index - cloudbox_limits[0], 0, 0, za_index, 0, joker) =
772 stokes_vec;
773 }
774 if (za_grid[za_index] > 90) {
775 stokes_vec = cloudbox_field_mono(
776 p_index - cloudbox_limits[0] - 1, 0, 0, za_index, 0, joker);
777 Numeric z_field_below = z_field(p_index - 1, 0, 0);
778 Numeric z_field_0 = z_field(p_index, 0, 0);
779
780 Numeric cos_theta, lstep;
781 if (za_grid[za_index] == 90.0) {
782 cos_theta = cos((za_grid[za_index] - 1) * PI / 180.);
783 //cos_theta = cos(90.00000001*PI/180.);
784 //cout<<"cos_theta"<<cos_theta<<endl;
785 } else {
786 cos_theta = cos(za_grid[za_index] * PI / 180.0);
787 }
788 Numeric dz = (z_field_0 - z_field_below);
789 lstep = abs(dz / cos_theta);
790
791 // Average temperature
792 Numeric rtp_temperature =
793 0.5 * (t_field(p_index, 0, 0) + t_field(p_index - 1, 0, 0));
794 //
795 // Average pressure
796 Numeric rtp_pressure = 0.5 * (p_grid[p_index] + p_grid[p_index - 1]);
797
798 //
799 // Average vmrs
800 for (Index k = 0; k < N_species; k++)
801 rtp_vmr[k] = 0.5 * (vmr_field(k, p_index, 0, 0) +
802 vmr_field(k, p_index - 1, 0, 0));
803 //
804 // Calculate scalar gas absorption and add it to abs_vec
805 // and ext_mat.
806 //
807
808 const Vector rtp_mag_dummy(3, 0);
809 const Vector ppath_los_dummy;
810 StokesVector nlte_dummy; //FIXME: do this right?
811 ArrayOfPropagationMatrix partial_dummy; // This is right since there should be only clearsky partials
812 ArrayOfStokesVector partial_nlte_dummy; // This is right since there should be only clearsky partials
814 propmat_clearsky,
815 nlte_dummy,
816 partial_dummy,
817 partial_nlte_dummy,
819 {},
820 Vector{f_grid[Range(f_index, 1)]},
821 rtp_mag_dummy,
822 ppath_los_dummy,
823 rtp_pressure,
824 rtp_temperature,
825 rtp_nlte_dummy,
826 rtp_vmr,
827 propmat_clearsky_agenda);
828
829 opt_prop_sum_propmat_clearsky(ext_mat, abs_vec, propmat_clearsky);
830
831 //
832 // Add average particle extinction to ext_mat.
833 //
834
835 // cout<<"cloudbox_limits[1]"<<cloudbox_limits[1]<<endl;
836 // cout<<"p_index - cloudbox_limits[0]"<<p_index - cloudbox_limits[0]<<endl;
837 for (Index k = 0; k < stokes_dim; k++) {
838 //
839 // Averaging of sca_vec:
840 //
841 sca_vec_av[k] +=
842 0.5 *
843 (doit_scat_field(
844 p_index - cloudbox_limits[0], 0, 0, za_index, 0, k) +
845 doit_scat_field(
846 p_index - cloudbox_limits[0] - 1, 0, 0, za_index, 0, k));
847 }
848 //
849 // Add average particle absorption to abs_vec.
850 //
851 abs_vec.AddAverageAtPosition(
852 abs_vec_field(p_index - cloudbox_limits[0], 0, 0, joker),
853 abs_vec_field(p_index - cloudbox_limits[0] - 1, 0, 0, joker));
854
855 //
856 // Add average particle extinction to ext_mat.
857 ext_mat.AddAverageAtPosition(
858 ext_mat_field(p_index - cloudbox_limits[0], 0, 0, joker, joker),
859 ext_mat_field(p_index - cloudbox_limits[0] + 1, 0, 0, joker, joker));
860
861 // Frequency
862 Numeric f = f_grid[f_index];
863 //
864 // Calculate Planck function
865 //
866 Numeric rte_planck_value = planck(f, rtp_temperature);
867
868 // Some messages:
869 if (out3.sufficient_priority()) {
870 vector_tmp = abs_vec.VectorAtPosition();
871 ext_mat.MatrixAtPosition(matrix_tmp);
872 out3 << "-----------------------------------------\n";
873 out3 << "Input for radiative transfer step \n"
874 << "calculation inside"
875 << " the cloudbox:"
876 << "\n";
877 out3 << "Stokes vector at intersection point: \n" << stokes_vec << "\n";
878 out3 << "lstep: ..." << lstep << "\n";
879 out3 << "------------------------------------------\n";
880 out3 << "Averaged coefficients: \n";
881 out3 << "Planck function: " << rte_planck_value << "\n";
882 out3 << "Scattering vector: " << sca_vec_av << "\n";
883 out3 << "Absorption vector: " << vector_tmp << "\n";
884 out3 << "Extinction matrix: " << matrix_tmp << "\n";
885
886 ARTS_ASSERT(!is_singular(matrix_tmp));
887 }
888
889 // Radiative transfer step calculation. The Stokes vector
890 // is updated until the considered point is reached.
891 rte_step_doit_replacement(stokes_vec,
892 Matrix(stokes_dim, stokes_dim),
893 ext_mat,
894 abs_vec,
895 sca_vec_av,
896 lstep,
897 rte_planck_value);
898
899 // Assign calculated Stokes Vector to cloudbox_field_mono.
900 cloudbox_field_mono(
901 p_index - cloudbox_limits[0], 0, 0, za_index, 0, joker) =
902 stokes_vec;
903 }
904
905 } // if loop end - for non_ground background
906
907 // bkgr=2 indicates that the background is the surface
908 else if (bkgr == 2) {
909 //Set rte_pos, rte_gp_p and rte_los to match the last point
910 //in ppath.
911 //pos
912 Vector rte_pos(atmosphere_dim);
913 Numeric z_field_0 = z_field(0, 0, 0);
914 rte_pos = z_field_0; //ppath_step.pos(np-1,Range(0,atmosphere_dim));
915 //los
916 Vector rte_los(1);
917 rte_los = za_grid[za_index]; //ppath_step.los(np-1,joker);
918 //gp_p
919 //GridPos rte_gp_p;
920 //rte_gp_p.idx = p_index;
921 //rte_gp_p.fd[0] = 0;
922 //rte_gp_p.fd[1] = 1;
923 //gridpos_copy( rte_gp_p, ppath_step.gp_p[np-1] );
924 // Executes the surface agenda
925 // FIXME: Convert to new agenda scheme before using
926 // surface_rtprop_agenda.execute();
927
929 "Surface reflections inside cloud box not yet handled.");
930 /*
931 See comment in function above
932 // Check returned variables
933 if( surface_emission.nrows() != f_grid.nelem() ||
934 surface_emission.ncols() != stokes_dim )
935 throw runtime_error(
936 "The size of the created *surface_emission* is not correct.");
937
938 Index nlos = surface_los.nrows();
939
940 // Define a local vector cloudbox_field_mono_sum which adds the
941 // products of groudnd_refl_coeffs with the downwelling
942 // radiation for each elements of surface_los
943 Vector cloudbox_field_mono_sum(stokes_dim,0);
944 // Loop over the surface_los elements
945 for( Index ilos=0; ilos < nlos; ilos++ )
946 {
947 if( stokes_dim == 1 )
948 {
949 cloudbox_field_mono_sum[0] += surface_refl_coeffs(ilos,f_index,0,0) *
950 cloudbox_field_mono(cloudbox_limits[0],
951 0, 0,
952 (za_grid.nelem() -1 - za_index), 0,
953 0);
954 }
955 else
956 {
957 Vector stokes_vec2(stokes_dim);
958 mult( stokes_vec2,
959 surface_refl_coeffs(ilos,0,joker,joker),
960 cloudbox_field_mono(cloudbox_limits[0],
961 0, 0,
962 (za_grid.nelem() -1 - za_index), 0,
963 joker));
964 for( Index is=0; is < stokes_dim; is++ )
965 {
966 cloudbox_field_mono_sum[is] += stokes_vec2[is];
967 }
968
969 }
970 }
971 // Copy from *cloudbox_field_mono_sum* to *cloudbox_field_mono*, and add the surface emission
972 for( Index is=0; is < stokes_dim; is++ )
973 {
974 cloudbox_field_mono (cloudbox_limits[0],
975 0, 0,
976 za_index, 0,
977 is) = cloudbox_field_mono_sum[is] + surface_emission(f_index,is);
978 }
979
980 //cout<<"za_grid"<<za_grid[za_index]<<endl;
981 //cout<<"p_index"<<p_index<<endl;
982 //cout<<"cloudboxlimit[0]"<<cloudbox_limits[0]<<endl;
983 // now the RT is done to the next point in the path.
984 //
985 Vector stokes_vec_local;
986 stokes_vec_local = cloudbox_field_mono (0,
987 0, 0,
988 za_index, 0,
989 joker);
990
991 Numeric z_field_above = z_field(p_index +1, 0, 0);
992 //Numeric z_field_0 = z_field(p_index, 0, 0);
993 Numeric cos_theta;
994 if (za_grid[za_index] == 90.0)
995 {
996 //cos_theta = cos((za_grid[za_index] -1)*PI/180.);
997 cos_theta = cos(90.00000001*PI/180.);
998 }
999 else
1000 {
1001 cos_theta = cos(za_grid[za_index]* PI/180.0);
1002 }
1003 Numeric dz = (z_field_above - z_field_0);
1004
1005 Numeric lstep = abs(dz/cos_theta) ;
1006
1007 // Average temperature
1008 Numeric rtp_temperature = 0.5 * (t_field(p_index,0,0)
1009 + t_field(p_index + 1,0,0));
1010
1011 //
1012 // Average pressure
1013 Numeric rtp_pressure = 0.5 * (p_grid[p_index]
1014 + p_grid[p_index + 1]);
1015
1016 //
1017 const Index N_species = vmr_field.nbooks();
1018 // Average vmrs
1019 for (Index k = 0; k < N_species; k++)
1020 {
1021 rtp_vmr[k] = 0.5 * (vmr_field(k,p_index,0,0) +
1022 vmr_field(k,p_index + 1,0,0));
1023 }
1024 //
1025 // Calculate scalar gas absorption and add it to abs_vec
1026 // and ext_mat.
1027 //
1028
1029 // FIXME: Convert to new agenda scheme before using
1030 //abs_scalar_gas_agenda.execute(p_index);
1031
1032 opt_prop_gas_agendaExecute(ext_mat, abs_vec, abs_scalar_gas,
1033 opt_prop_gas_agenda);
1034
1035 //
1036 // Add average particle extinction to ext_mat.
1037 //
1038 //cout<<"Reached Here ????????????????????????????????????????????????";
1039 for (Index k = 0; k < stokes_dim; k++)
1040 {
1041 for (Index j = 0; j < stokes_dim; j++)
1042 {
1043 ext_mat(0,k,j) += 0.5 *
1044 (ext_mat_field(p_index - cloudbox_limits[0],0,0,k,j) +
1045 ext_mat_field(p_index - cloudbox_limits[0]+ 1,0,0,k,j));
1046 }
1047
1048
1049 //
1050 //
1051 // Add average particle absorption to abs_vec.
1052 //
1053 abs_vec(0,k) += 0.5 *
1054 (abs_vec_field(p_index - cloudbox_limits[0],0,0,k) +
1055 abs_vec_field(p_index - cloudbox_limits[0]+1,0,0,k));
1056 //
1057 // Averaging of sca_vec:
1058 //
1059 sca_vec_av[k] = 0.5*( doit_scat_field(p_index- cloudbox_limits[0],
1060 0, 0, za_index, 0, k)
1061 + doit_scat_field(p_index- cloudbox_limits[0]+1,
1062 0, 0, za_index, 0, k)) ;
1063
1064 }
1065 // Frequency
1066 Numeric f = f_grid[f_index];
1067 //
1068 // Calculate Planck function
1069 //
1070 Numeric rte_planck_value = planck(f, rtp_temperature);
1071
1072 ARTS_ASSERT (!is_singular( ext_mat(0,joker,joker)));
1073
1074 // Radiative transfer step calculation. The Stokes vector
1075 // is updated until the considered point is reached.
1076 rte_step_doit(stokes_vec_local, ext_mat(0,joker,joker),
1077 abs_vec(0,joker),
1078 sca_vec_av, lstep, rte_planck_value);
1079 // Assign calculated Stokes Vector to cloudbox_field_mono.
1080 cloudbox_field_mono(p_index - cloudbox_limits[0],
1081 0, 0,
1082 za_index, 0,
1083 joker) = stokes_vec_local;
1084 */
1085 } //end else loop over surface
1086}
1087
1089 Tensor6View cloudbox_field_mono,
1090 // ppath_step_agenda:
1091 const Index& p_index,
1092 const Index& lat_index,
1093 const Index& lon_index,
1094 const Index& za_index,
1095 const Index& aa_index,
1096 ConstVectorView za_grid,
1097 ConstVectorView aa_grid,
1098 const ArrayOfIndex& cloudbox_limits,
1099 ConstTensor6View doit_scat_field,
1100 // Calculate scalar gas absorption:
1101 const Agenda& propmat_clearsky_agenda,
1102 ConstTensor4View vmr_field,
1103 // Propagation path calculation:
1104 const Agenda& ppath_step_agenda,
1105 const Numeric& ppath_lmax,
1106 const Numeric& ppath_lraytrace,
1107 ConstVectorView p_grid,
1108 ConstVectorView lat_grid,
1109 ConstVectorView lon_grid,
1110 ConstTensor3View z_field,
1111 ConstVectorView refellipsoid,
1112 // Calculate thermal emission:
1113 ConstTensor3View t_field,
1114 ConstVectorView f_grid,
1115 const Index& f_index,
1116 //particle optical properties
1117 ConstTensor5View ext_mat_field,
1118 ConstTensor4View abs_vec_field,
1119 const Index&, //scat_za_interp
1120 const Verbosity& verbosity) {
1122
1123 Ppath ppath_step;
1124 const Index stokes_dim = cloudbox_field_mono.ncols();
1125
1126 Vector sca_vec_av(stokes_dim, 0);
1127 Vector aa_g(aa_grid.nelem());
1128
1129 for (Index i = 0; i < aa_grid.nelem(); i++)
1130 aa_g[i] = aa_grid[i] - 180.;
1131
1132 //Initialize ppath for 3D.
1133 ppath_init_structure(ppath_step, 3, 1);
1134 // See documentation of ppath_init_structure for
1135 // understanding the parameters.
1136
1137 // The first dimension of pos are the points in
1138 // the propagation path.
1139 // Here we initialize the first point.
1140 // The second is: radius, latitude, longitude
1141
1142 ppath_step.pos(0, 2) = lon_grid[lon_index];
1143 ppath_step.pos(0, 1) = lat_grid[lat_index];
1144 ppath_step.pos(0, 0) = z_field(p_index, lat_index, lon_index);
1145 // As always on top of the lat. grid positions, OK to call refell2r:
1146 ppath_step.r[0] =
1147 refell2r(refellipsoid, ppath_step.pos(0, 1)) + ppath_step.pos(0, 0);
1148
1149 // Define the direction:
1150 ppath_step.los(0, 0) = za_grid[za_index];
1151 ppath_step.los(0, 1) = aa_g[aa_index];
1152
1153 // Define the grid positions:
1154 ppath_step.gp_p[0].idx = p_index;
1155 ppath_step.gp_p[0].fd[0] = 0.;
1156 ppath_step.gp_p[0].fd[1] = 1.;
1157
1158 ppath_step.gp_lat[0].idx = lat_index;
1159 ppath_step.gp_lat[0].fd[0] = 0.;
1160 ppath_step.gp_lat[0].fd[1] = 1.;
1161
1162 ppath_step.gp_lon[0].idx = lon_index;
1163 ppath_step.gp_lon[0].fd[0] = 0.;
1164 ppath_step.gp_lon[0].fd[1] = 1.;
1165
1166 // Call ppath_step_agenda:
1168 ppath_step,
1169 ppath_lmax,
1170 ppath_lraytrace,
1171 Vector(1, f_grid[f_index]),
1172 ppath_step_agenda);
1173
1174 // Check whether the next point is inside or outside the
1175 // cloudbox. Only if the next point lies inside the
1176 // cloudbox a radiative transfer step caclulation has to
1177 // be performed.
1178 if (is_inside_cloudbox(ppath_step, cloudbox_limits, true)) {
1179 // Gridpositions inside the cloudbox.
1180 // The optical properties are stored only inside the
1181 // cloudbox. For interpolation we use grids
1182 // inside the cloudbox.
1183
1184 ArrayOfGridPos cloud_gp_p = ppath_step.gp_p;
1185 ArrayOfGridPos cloud_gp_lat = ppath_step.gp_lat;
1186 ArrayOfGridPos cloud_gp_lon = ppath_step.gp_lon;
1187
1188 for (Index i = 0; i < ppath_step.np; i++) {
1189 cloud_gp_p[i].idx -= cloudbox_limits[0];
1190 cloud_gp_lat[i].idx -= cloudbox_limits[2];
1191 cloud_gp_lon[i].idx -= cloudbox_limits[4];
1192 }
1193 const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
1194 const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
1195 const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
1196 gridpos_upperend_check(cloud_gp_p[0], n1);
1197 gridpos_upperend_check(cloud_gp_p[ppath_step.np - 1], n1);
1198 gridpos_upperend_check(cloud_gp_lat[0], n2);
1199 gridpos_upperend_check(cloud_gp_lat[ppath_step.np - 1], n2);
1200 gridpos_upperend_check(cloud_gp_lon[0], n3);
1201 gridpos_upperend_check(cloud_gp_lon[ppath_step.np - 1], n3);
1202
1203 Matrix itw(ppath_step.np, 8);
1204 interpweights(itw, cloud_gp_p, cloud_gp_lat, cloud_gp_lon);
1205
1206 Matrix itw_p(ppath_step.np, 2);
1207 interpweights(itw_p, cloud_gp_p);
1208
1209 // The zenith angles and azimuth of the propagation path are
1210 // needed as we have to
1211 // interpolate the intensity field and the scattered field on the
1212 // right angles.
1213 VectorView los_grid_za = ppath_step.los(joker, 0);
1214 VectorView los_grid_aa = ppath_step.los(joker, 1);
1215
1216 for (Index i = 0; i < los_grid_aa.nelem(); i++)
1217 los_grid_aa[i] = los_grid_aa[i] + 180.;
1218
1219 ArrayOfGridPos gp_za(los_grid_za.nelem());
1220 gridpos(gp_za, za_grid, los_grid_za);
1221
1222 ArrayOfGridPos gp_aa(los_grid_aa.nelem());
1223 gridpos(gp_aa, aa_grid, los_grid_aa);
1224
1225 Matrix itw_p_za(ppath_step.np, 32);
1227 itw_p_za, cloud_gp_p, cloud_gp_lat, cloud_gp_lon, gp_za, gp_aa);
1228
1229 // Ppath_step normally has 2 points, the starting
1230 // point and the intersection point.
1231 // But there can be points in between, when a maximum
1232 // lstep is given. We have to interpolate on all the
1233 // points in the ppath_step.
1234
1235 Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.np);
1236 Matrix abs_vec_int(stokes_dim, ppath_step.np);
1237 Matrix sca_vec_int(stokes_dim, ppath_step.np, 0.);
1238 Matrix cloudbox_field_mono_int(stokes_dim, ppath_step.np, 0.);
1239 Vector t_int(ppath_step.np);
1240 Vector vmr_int(ppath_step.np);
1241 Vector p_int(ppath_step.np);
1242 Vector stokes_vec(stokes_dim);
1243 //Tensor3 ext_mat_gas(stokes_dim, stokes_dim, ppath_step.np);
1244 //Matrix abs_vec_gas(stokes_dim, ppath_step.np);
1245
1246 // Calculate the average of the coefficients for the layers
1247 // to be considered in the
1248 // radiative transfer calculation.
1249
1250 for (Index i = 0; i < stokes_dim; i++) {
1251 // Extinction matrix requires a second loop
1252 // over stokes_dim
1253 out3 << "Interpolate ext_mat:\n";
1254 for (Index j = 0; j < stokes_dim; j++) {
1255 //
1256 // Interpolation of ext_mat
1257 //
1258 interp(ext_mat_int(i, j, joker),
1259 itw,
1260 ext_mat_field(joker, joker, joker, i, j),
1261 cloud_gp_p,
1262 cloud_gp_lat,
1263 cloud_gp_lon);
1264 }
1265 // Absorption vector:
1266 //
1267 // Interpolation of abs_vec
1268 //
1269 interp(abs_vec_int(i, joker),
1270 itw,
1271 abs_vec_field(joker, joker, joker, i),
1272 cloud_gp_p,
1273 cloud_gp_lat,
1274 cloud_gp_lon);
1275 //
1276 // Scattered field:
1277 //
1278 // Interpolation of sca_vec:
1279 //
1280 out3 << "Interpolate doit_scat_field:\n";
1281 interp(sca_vec_int(i, joker),
1282 itw_p_za,
1283 doit_scat_field(joker, joker, joker, joker, joker, i),
1284 cloud_gp_p,
1285 cloud_gp_lat,
1286 cloud_gp_lon,
1287 gp_za,
1288 gp_aa);
1289 out3 << "Interpolate cloudbox_field_mono:\n";
1290 interp(cloudbox_field_mono_int(i, joker),
1291 itw_p_za,
1292 cloudbox_field_mono(joker, joker, joker, joker, joker, i),
1293 cloud_gp_p,
1294 cloud_gp_lat,
1295 cloud_gp_lon,
1296 gp_za,
1297 gp_aa);
1298 }
1299 //
1300 // Planck function
1301 //
1302 // Interpolate temperature field
1303 //
1304 out3 << "Interpolate temperature field\n";
1305 interp(t_int,
1306 itw,
1307 t_field(joker, joker, joker),
1308 ppath_step.gp_p,
1309 ppath_step.gp_lat,
1310 ppath_step.gp_lon);
1311
1312 //
1313 // The vmr_field is needed for the gaseous absorption
1314 // calculation.
1315 //
1316 const Index N_species = vmr_field.nbooks();
1317 //
1318 // Interpolated vmr_list, holds a vmr_list for each point in
1319 // ppath_step.
1320 //
1321 Matrix vmr_list_int(N_species, ppath_step.np);
1322
1323 for (Index i = 0; i < N_species; i++) {
1324 out3 << "Interpolate vmr field\n";
1325 interp(vmr_int,
1326 itw,
1327 vmr_field(i, joker, joker, joker),
1328 ppath_step.gp_p,
1329 ppath_step.gp_lat,
1330 ppath_step.gp_lon);
1331
1332 vmr_list_int(i, joker) = vmr_int;
1333 }
1334
1335 // Presssure (needed for the calculation of gas absorption)
1336 itw2p(p_int, p_grid, ppath_step.gp_p, itw_p);
1337
1338 out3 << "Calculate radiative transfer inside cloudbox.\n";
1340 cloudbox_field_mono,
1341 propmat_clearsky_agenda,
1342 ppath_step,
1343 t_int,
1344 vmr_list_int,
1345 ext_mat_int,
1346 abs_vec_int,
1347 sca_vec_int,
1348 cloudbox_field_mono_int,
1349 p_int,
1350 cloudbox_limits,
1351 f_grid,
1352 f_index,
1353 p_index,
1354 lat_index,
1355 lon_index,
1356 za_index,
1357 aa_index,
1358 verbosity);
1359 } //end if inside cloudbox
1360}
1361
1363 //Output
1364 Tensor6View cloudbox_field_mono,
1365 // Input
1366 const Agenda& propmat_clearsky_agenda,
1367 const Ppath& ppath_step,
1368 ConstVectorView t_int,
1369 ConstMatrixView vmr_list_int,
1370 ConstTensor3View ext_mat_int,
1371 ConstMatrixView abs_vec_int,
1372 ConstMatrixView sca_vec_int,
1373 ConstMatrixView cloudbox_field_mono_int,
1374 ConstVectorView p_int,
1375 const ArrayOfIndex& cloudbox_limits,
1376 ConstVectorView f_grid,
1377 const Index& f_index,
1378 const Index& p_index,
1379 const Index& lat_index,
1380 const Index& lon_index,
1381 const Index& za_index,
1382 const Index& aa_index,
1383 const Verbosity& verbosity) {
1385
1386 const Index N_species = vmr_list_int.nrows();
1387 const Index stokes_dim = cloudbox_field_mono.ncols();
1388 const Index atmosphere_dim = cloudbox_limits.nelem() / 2;
1389
1390 Vector sca_vec_av(stokes_dim, 0);
1391 Vector stokes_vec(stokes_dim, 0.);
1392 EnergyLevelMap rtp_nlte_dummy;
1393 Vector rtp_vmr_local(N_species, 0.);
1394
1395 // Two propmat_clearsky to average between
1396 PropagationMatrix cur_propmat_clearsky;
1397 PropagationMatrix prev_propmat_clearsky;
1398
1399 PropagationMatrix ext_mat_local;
1400 StokesVector abs_vec_local;
1401 Matrix matrix_tmp(stokes_dim, stokes_dim);
1402 Vector vector_tmp(stokes_dim);
1403
1404 // Incoming stokes vector
1405 stokes_vec = cloudbox_field_mono_int(joker, ppath_step.np - 1);
1406
1407 for (Index k = ppath_step.np - 1; k >= 0; k--) {
1408 // Save propmat_clearsky from previous level by
1409 // swapping it with current level
1410 swap(cur_propmat_clearsky, prev_propmat_clearsky);
1411
1412 //
1413 // Calculate scalar gas absorption
1414 //
1415 const Vector rtp_mag_dummy(3, 0);
1416 const Vector ppath_los_dummy;
1417
1418 StokesVector nlte_dummy; //FIXME: do this right?
1419 ArrayOfPropagationMatrix partial_dummy; // This is right since there should be only clearsky partials
1420 ArrayOfStokesVector partial_nlte_dummy; // This is right since there should be only clearsky partials
1422 cur_propmat_clearsky,
1423 nlte_dummy,
1424 partial_dummy,
1425 partial_nlte_dummy,
1427 {},
1428 Vector{f_grid[Range(f_index, 1)]},
1429 rtp_mag_dummy,
1430 ppath_los_dummy,
1431 p_int[k],
1432 t_int[k],
1433 rtp_nlte_dummy,
1434 Vector{vmr_list_int(joker, k)},
1435 propmat_clearsky_agenda);
1436
1437 // Skip any further calculations for the first point.
1438 // We need values at two ppath points before we can average.
1439 if (k == ppath_step.np - 1) continue;
1440
1441 // Average prev_propmat_clearsky with cur_propmat_clearsky
1442 prev_propmat_clearsky += cur_propmat_clearsky;
1443 prev_propmat_clearsky *= 0.5;
1444
1446 ext_mat_local, abs_vec_local, prev_propmat_clearsky);
1447
1448 for (Index i = 0; i < stokes_dim; i++) {
1449 //
1450 // Averaging of sca_vec:
1451 //
1452 sca_vec_av[i] = 0.5 * (sca_vec_int(i, k) + sca_vec_int(i, k + 1));
1453 }
1454
1455 //
1456 // Add average particle absorption to abs_vec.
1457 //
1458 abs_vec_local.AddAverageAtPosition(abs_vec_int(joker, k),
1459 abs_vec_int(joker, k + 1));
1460
1461 //
1462 // Add average particle extinction to ext_mat.
1463 //
1464 ext_mat_local.AddAverageAtPosition(ext_mat_int(joker, joker, k),
1465 ext_mat_int(joker, joker, k + 1));
1466
1467 // Frequency
1468 Numeric f = f_grid[f_index];
1469 //
1470 // Calculate Planck function
1471 //
1472 Numeric rte_planck_value = planck(f, 0.5 * (t_int[k] + t_int[k + 1]));
1473
1474 // Length of the path between the two layers.
1475 Numeric lstep = ppath_step.lstep[k];
1476
1477 // Some messages:
1478 if (out3.sufficient_priority()) {
1479 vector_tmp = abs_vec_local.VectorAtPosition();
1480 ext_mat_local.MatrixAtPosition(matrix_tmp);
1481 out3 << "-----------------------------------------\n";
1482 out3 << "Input for radiative transfer step \n"
1483 << "calculation inside"
1484 << " the cloudbox:"
1485 << "\n";
1486 out3 << "Stokes vector at intersection point: \n" << stokes_vec << "\n";
1487 out3 << "lstep: ..." << lstep << "\n";
1488 out3 << "------------------------------------------\n";
1489 out3 << "Averaged coefficients: \n";
1490 out3 << "Planck function: " << rte_planck_value << "\n";
1491 out3 << "Scattering vector: " << sca_vec_av << "\n";
1492 out3 << "Absorption vector: " << vector_tmp << "\n";
1493 out3 << "Extinction matrix: " << matrix_tmp << "\n";
1494
1495 ARTS_ASSERT(!is_singular(matrix_tmp));
1496 }
1497
1498 // Radiative transfer step calculation. The Stokes vector
1499 // is updated until the considered point is reached.
1500 rte_step_doit_replacement(stokes_vec,
1501 Matrix(stokes_dim, stokes_dim),
1502 ext_mat_local,
1503 abs_vec_local,
1504 sca_vec_av,
1505 lstep,
1506 rte_planck_value);
1507
1508 } // End of loop over ppath_step.
1509 // Assign calculated Stokes Vector to cloudbox_field_mono.
1510 if (atmosphere_dim == 1)
1511 cloudbox_field_mono(
1512 p_index - cloudbox_limits[0], 0, 0, za_index, 0, joker) =
1513 stokes_vec;
1514 else if (atmosphere_dim == 3)
1515 cloudbox_field_mono(p_index - cloudbox_limits[0],
1516 lat_index - cloudbox_limits[2],
1517 lon_index - cloudbox_limits[4],
1518 za_index,
1519 aa_index,
1520 joker) = stokes_vec;
1521}
1522
1524 //Output
1525 Tensor6View cloudbox_field_mono,
1526 //Input
1527 const Agenda& surface_rtprop_agenda,
1528 ConstVectorView f_grid,
1529 const Index& f_index,
1530 const Index& stokes_dim,
1531 const Ppath& ppath_step,
1532 const ArrayOfIndex& cloudbox_limits,
1533 ConstVectorView za_grid,
1534 const Index& za_index) {
1535 chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
1536
1537 Matrix iy;
1538
1539 // Local output of surface_rtprop_agenda.
1540 Numeric surface_skin_t;
1541 Matrix surface_emission;
1542 Matrix surface_los;
1543 Tensor4 surface_rmatrix;
1544
1545 //Set rte_pos and rte_los to match the last point in ppath.
1546
1547 Index np = ppath_step.np;
1548
1549 Vector rte_pos; // ppath_step.pos contains two columns for 1D
1550 rte_pos.resize(ppath_step.dim);
1551 rte_pos = ppath_step.pos(np - 1, Range(0, ppath_step.dim));
1552
1553 Vector rte_los;
1554 rte_los.resize(ppath_step.los.ncols());
1555 rte_los = ppath_step.los(np - 1, joker);
1556
1557 //Execute the surface_rtprop_agenda which gives the surface
1558 //parameters.
1560 surface_skin_t,
1561 surface_emission,
1562 surface_los,
1563 surface_rmatrix,
1564 Vector(1, f_grid[f_index]),
1565 rte_pos,
1566 rte_los,
1567 surface_rtprop_agenda);
1568
1569 iy = surface_emission;
1570
1571 Index nlos = surface_los.nrows();
1572
1573 if (nlos > 0) {
1574 Vector rtmp(stokes_dim); // Reflected Stokes vector for 1 frequency
1575
1576 for (Index ilos = 0; ilos < nlos; ilos++) {
1577 // Several things needs to be fixed here. As far as I understand it,
1578 // this works only for specular cases and if the lower cloudbox limit
1579 // is exactly at the surface (PE, 120828)
1580
1581 mult(rtmp,
1582 surface_rmatrix(ilos, 0, joker, joker),
1583 cloudbox_field_mono(cloudbox_limits[0],
1584 0,
1585 0,
1586 (za_grid.nelem() - 1 - za_index),
1587 0,
1588 joker));
1589 iy(0, joker) += rtmp;
1590 }
1591 }
1592 cloudbox_field_mono(cloudbox_limits[0], 0, 0, za_index, 0, joker) =
1593 iy(0, joker);
1594}
1595
1596void cloudbox_field_ngAcceleration(Tensor6& cloudbox_field_mono,
1597 const ArrayOfTensor6& acceleration_input,
1598 const Index& accelerated,
1599 const Verbosity&) {
1600 const Index N_p = cloudbox_field_mono.nvitrines();
1601 const Index N_za = cloudbox_field_mono.npages();
1602
1603 // Loop over 4 components of Stokes Vector
1604 for (Index i = 0; i < accelerated; ++i) {
1605 ConstMatrixView S1 = acceleration_input[0](joker, 0, 0, joker, 0, i);
1606 ConstMatrixView S2 = acceleration_input[1](joker, 0, 0, joker, 0, i);
1607 ConstMatrixView S3 = acceleration_input[2](joker, 0, 0, joker, 0, i);
1608 ConstMatrixView S4 = acceleration_input[3](joker, 0, 0, joker, 0, i);
1609
1610 ConstMatrixView J = S4;
1611 Matrix Q1;
1612 Matrix Q2;
1613 Matrix Q3;
1614 Numeric A1 = 0;
1615 Numeric A2B1 = 0;
1616 Numeric B2 = 0;
1617 Numeric C1 = 0;
1618 Numeric C2 = 0;
1619 Numeric NGA = 0;
1620 Numeric NGB = 0;
1621
1622 // Q1 = -2*S3 + S4 + S2
1623
1624 Q1 = S3;
1625 Q1 *= -2;
1626 Q1 += S4;
1627 Q1 += S2;
1628
1629 // Q2 = S4 - S3 - S2 + S1
1630 Q2 = S4;
1631 Q2 -= S3;
1632 Q2 -= S2;
1633 Q2 += S1;
1634
1635 //Q3 = S4 - S3
1636 Q3 = S4;
1637 Q3 -= S3;
1638
1639 for (Index p_index = 0; p_index < N_p; ++p_index) {
1640 for (Index za_index = 0; za_index < N_za; ++za_index) {
1641 A1 += Q1(p_index, za_index) * Q1(p_index, za_index) *
1642 J(p_index, za_index);
1643 A2B1 += Q2(p_index, za_index) * Q1(p_index, za_index) *
1644 J(p_index, za_index);
1645 B2 += Q2(p_index, za_index) * Q2(p_index, za_index) *
1646 J(p_index, za_index);
1647 C1 += Q1(p_index, za_index) * Q3(p_index, za_index) *
1648 J(p_index, za_index);
1649 C2 += Q2(p_index, za_index) * Q3(p_index, za_index) *
1650 J(p_index, za_index);
1651 }
1652 }
1653
1654 NGA = (C1 * B2 - C2 * A2B1) / (A1 * B2 - A2B1 * A2B1);
1655 NGB = (C2 * A1 - C1 * A2B1) / (A1 * B2 - A2B1 * A2B1);
1656
1657 if (!std::isnan(NGB) && !std::isnan(NGA)) {
1658 // Calculating the accelerated field
1659 for (Index p_index = 0; p_index < N_p; ++p_index) {
1660 for (Index za_index = 0; za_index < N_za; ++za_index) {
1661 Q1(p_index, za_index) = (1 - NGA - NGB) * S4(p_index, za_index) +
1662 NGA * S3(p_index, za_index) +
1663 NGB * S2(p_index, za_index);
1664 }
1665 }
1666 cloudbox_field_mono(joker, 0, 0, joker, 0, i) = Q1;
1667 }
1668 }
1669}
1670
1672 Tensor3View ext_mat_int,
1673 MatrixView abs_vec_int,
1674 MatrixView sca_vec_int,
1675 MatrixView cloudbox_field_mono_int,
1676 VectorView t_int,
1677 MatrixView vmr_list_int,
1678 VectorView p_int,
1679 //Input
1680 ConstTensor5View ext_mat_field,
1681 ConstTensor4View abs_vec_field,
1682 ConstTensor6View doit_scat_field,
1683 ConstTensor6View cloudbox_field_mono,
1684 ConstTensor3View t_field,
1685 ConstTensor4View vmr_field,
1686 ConstVectorView p_grid,
1687 const Ppath& ppath_step,
1688 const ArrayOfIndex& cloudbox_limits,
1689 ConstVectorView za_grid,
1690 const Index& scat_za_interp,
1691 const Verbosity& verbosity) {
1693
1694 // Stokes dimension
1695 const Index stokes_dim = cloudbox_field_mono.ncols();
1696
1697 // Gridpositions inside the cloudbox.
1698 // The optical properties are stored only inside the
1699 // cloudbox. For interpolation we use grids
1700 // inside the cloudbox.
1701 ArrayOfGridPos cloud_gp_p = ppath_step.gp_p;
1702
1703 for (Index i = 0; i < ppath_step.np; i++)
1704 cloud_gp_p[i].idx -= cloudbox_limits[0];
1705
1706 // Grid index for points at upper limit of cloudbox must be shifted
1707 const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
1708 gridpos_upperend_check(cloud_gp_p[0], n1);
1709 gridpos_upperend_check(cloud_gp_p[ppath_step.np - 1], n1);
1710
1711 Matrix itw(cloud_gp_p.nelem(), 2);
1712 interpweights(itw, cloud_gp_p);
1713
1714 // The zenith angles of the propagation path are needed as we have to
1715 // interpolate the intensity field and the scattered field on the
1716 // right angles.
1717 Vector los_grid{ppath_step.los(joker, 0)};
1718
1719 ArrayOfGridPos gp_za(los_grid.nelem());
1720 gridpos(gp_za, za_grid, los_grid);
1721
1722 Matrix itw_p_za(cloud_gp_p.nelem(), 4);
1723 interpweights(itw_p_za, cloud_gp_p, gp_za);
1724
1725 // Calculate the average of the coefficients for the layers
1726 // to be considered in the
1727 // radiative transfer calculation.
1728
1729 for (Index i = 0; i < stokes_dim; i++) {
1730 // Extinction matrix requires a second loop
1731 // over stokes_dim
1732 out3 << "Interpolate ext_mat:\n";
1733 for (Index j = 0; j < stokes_dim; j++) {
1734 //
1735 // Interpolation of ext_mat
1736 //
1737 interp(ext_mat_int(i, j, joker),
1738 itw,
1739 ext_mat_field(joker, 0, 0, i, j),
1740 cloud_gp_p);
1741 }
1742 // Particle absorption vector:
1743 //
1744 // Interpolation of abs_vec
1745 // //
1746 out3 << "Interpolate abs_vec:\n";
1747 interp(
1748 abs_vec_int(i, joker), itw, abs_vec_field(joker, 0, 0, i), cloud_gp_p);
1749 //
1750 // Scattered field:
1751 //
1752 //
1753
1754 out3 << "Interpolate doit_scat_field and cloudbox_field_mono:\n";
1755 if (scat_za_interp == 0) // linear interpolation
1756 {
1757 interp(sca_vec_int(i, joker),
1758 itw_p_za,
1759 doit_scat_field(joker, 0, 0, joker, 0, i),
1760 cloud_gp_p,
1761 gp_za);
1762 interp(cloudbox_field_mono_int(i, joker),
1763 itw_p_za,
1764 cloudbox_field_mono(joker, 0, 0, joker, 0, i),
1765 cloud_gp_p,
1766 gp_za);
1767 } else if (scat_za_interp == 1) //polynomial interpolation
1768 {
1769 // These intermediate variables are needed because polynomial
1770 // interpolation is not implemented as multidimensional
1771 // interpolation.
1772 Tensor3 sca_vec_int_za(
1773 stokes_dim, ppath_step.np, za_grid.nelem(), 0.);
1774 Tensor3 cloudbox_field_mono_int_za(
1775 stokes_dim, ppath_step.np, za_grid.nelem(), 0.);
1776 for (Index za = 0; za < za_grid.nelem(); za++) {
1777 interp(sca_vec_int_za(i, joker, za),
1778 itw,
1779 doit_scat_field(joker, 0, 0, za, 0, i),
1780 cloud_gp_p);
1781 out3 << "Interpolate cloudbox_field_mono:\n";
1782 interp(cloudbox_field_mono_int_za(i, joker, za),
1783 itw,
1784 cloudbox_field_mono(joker, 0, 0, za, 0, i),
1785 cloud_gp_p);
1786 }
1787 for (Index ip = 0; ip < ppath_step.np; ip++) {
1788 sca_vec_int(i, ip) = interp_poly(za_grid,
1789 sca_vec_int_za(i, ip, joker),
1790 los_grid[ip],
1791 gp_za[ip]);
1792 cloudbox_field_mono_int(i, ip) =
1793 interp_poly(za_grid,
1794 cloudbox_field_mono_int_za(i, ip, joker),
1795 los_grid[ip],
1796 gp_za[ip]);
1797 }
1798 }
1799 }
1800 //
1801 // Planck function
1802 //
1803 // Interpolate temperature field
1804 //
1805 out3 << "Interpolate temperature field\n";
1806 interp(t_int, itw, t_field(joker, 0, 0), ppath_step.gp_p);
1807 //
1808 // The vmr_field is needed for the gaseous absorption
1809 // calculation.
1810 //
1811 const Index N_species = vmr_field.nbooks();
1812 //
1813 // Interpolated vmr_list, holds a vmr_list for each point in
1814 // ppath_step.
1815 //
1816 Vector vmr_int(ppath_step.np);
1817
1818 for (Index i_sp = 0; i_sp < N_species; i_sp++) {
1819 out3 << "Interpolate vmr field\n";
1820 interp(vmr_int, itw, vmr_field(i_sp, joker, 0, 0), ppath_step.gp_p);
1821 vmr_list_int(i_sp, joker) = vmr_int;
1822 }
1823 //
1824 // Interpolate pressure
1825 //
1826 itw2p(p_int, p_grid, ppath_step.gp_p, itw);
1827}
1828
1829void za_gridOpt( //Output:
1830 Vector& za_grid_opt,
1831 Matrix& cloudbox_field_opt,
1832 // Input
1833 ConstVectorView za_grid_fine,
1834 ConstTensor6View cloudbox_field_mono,
1835 const Numeric& acc,
1836 const Index& scat_za_interp) {
1837 Index N_za = za_grid_fine.nelem();
1838
1839 ARTS_ASSERT(cloudbox_field_mono.npages() == N_za);
1840
1841 Index N_p = cloudbox_field_mono.nvitrines();
1842
1843 Vector i_approx_interp(N_za);
1844 Vector za_reduced(2);
1845
1846 ArrayOfIndex idx;
1847 idx.push_back(0);
1848 idx.push_back(N_za - 1);
1849 ArrayOfIndex idx_unsorted;
1850
1851 Numeric max_diff = 100;
1852
1853 ArrayOfGridPos gp_za(N_za);
1854 Matrix itw(za_grid_fine.nelem(), 2);
1855
1856 ArrayOfIndex i_sort;
1857 Vector diff_vec(N_za);
1858 Vector max_diff_za(N_p);
1859 ArrayOfIndex ind_za(N_p);
1860 Numeric max_diff_p;
1861 Index ind_p = 0;
1862
1863 while (max_diff > acc) {
1864 za_reduced.resize(idx.nelem());
1865 cloudbox_field_opt.resize(N_p, idx.nelem());
1866 max_diff_za = 0.;
1867 max_diff_p = 0.;
1868
1869 // Interpolate reduced intensity field on fine za_grid for
1870 // all pressure levels
1871 for (Index i_p = 0; i_p < N_p; i_p++) {
1872 for (Index i_za_red = 0; i_za_red < idx.nelem(); i_za_red++) {
1873 za_reduced[i_za_red] = za_grid_fine[idx[i_za_red]];
1874 cloudbox_field_opt(i_p, i_za_red) =
1875 cloudbox_field_mono(i_p, 0, 0, idx[i_za_red], 0, 0);
1876 }
1877 // Calculate grid positions
1878 gridpos(gp_za, za_reduced, za_grid_fine);
1879 //linear interpolation
1880 if (scat_za_interp == 0 || idx.nelem() < 3) {
1881 interpweights(itw, gp_za);
1882 interp(i_approx_interp, itw, cloudbox_field_opt(i_p, joker), gp_za);
1883 } else if (scat_za_interp == 1) {
1884 for (Index i_za = 0; i_za < N_za; i_za++) {
1885 i_approx_interp[i_za] = interp_poly(za_reduced,
1886 cloudbox_field_opt(i_p, joker),
1887 za_grid_fine[i_za],
1888 gp_za[i_za]);
1889 }
1890 } else
1891 // Interpolation method not defined
1892 ARTS_ASSERT(false);
1893
1894 // Calculate differences between approximated i-vector and
1895 // exact i_vector for the i_p pressure level
1896 for (Index i_za = 0; i_za < N_za; i_za++) {
1897 diff_vec[i_za] = abs(cloudbox_field_mono(i_p, 0, 0, i_za, 0, 0) -
1898 i_approx_interp[i_za]);
1899 if (diff_vec[i_za] > max_diff_za[i_p]) {
1900 max_diff_za[i_p] = diff_vec[i_za];
1901 ind_za[i_p] = i_za;
1902 }
1903 }
1904 // Take maximum value of max_diff_za
1905 if (max_diff_za[i_p] > max_diff_p) {
1906 max_diff_p = max_diff_za[i_p];
1907 ind_p = i_p;
1908 }
1909 }
1910
1911 //Transform in %
1912 max_diff =
1913 max_diff_p / cloudbox_field_mono(ind_p, 0, 0, ind_za[ind_p], 0, 0) * 100.;
1914
1915 idx.push_back(ind_za[ind_p]);
1916 idx_unsorted = idx;
1917
1918 i_sort.resize(idx_unsorted.nelem());
1919 get_sorted_indexes(i_sort, idx_unsorted);
1920
1921 for (Index i = 0; i < idx_unsorted.nelem(); i++)
1922 idx[i] = idx_unsorted[i_sort[i]];
1923
1924 za_reduced.resize(idx.nelem());
1925 }
1926
1927 za_grid_opt.resize(idx.nelem());
1928 cloudbox_field_opt.resize(N_p, idx.nelem());
1929 for (Index i = 0; i < idx.nelem(); i++) {
1930 za_grid_opt[i] = za_grid_fine[idx[i]];
1931 cloudbox_field_opt(joker, i) = cloudbox_field_mono(joker, 0, 0, idx[i], 0, 0);
1932 }
1933}
1934
1936 Tensor6& doit_scat_field,
1937 const Tensor6& cloudbox_field_mono,
1938 const ArrayOfIndex& cloudbox_limits,
1939 const Agenda& spt_calc_agenda,
1940 const Index& atmosphere_dim,
1941 const Vector& za_grid,
1942 const Vector& aa_grid,
1943 const Tensor4& pnd_field,
1944 const Tensor3& t_field,
1945 const Numeric& norm_error_threshold,
1946 const Index& norm_debug,
1947 const Verbosity& verbosity) {
1948 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
1949 "Only 1D is supported here for now");
1950
1953
1954 // Number of zenith angles.
1955 const Index Nza = za_grid.nelem();
1956
1957 ARTS_USER_ERROR_IF (za_grid[0] != 0. || za_grid[Nza - 1] != 180.,
1958 "The range of *za_grid* must [0 180].");
1959
1960 // Number of azimuth angles.
1961 const Index Naa = aa_grid.nelem();
1962
1963 ARTS_USER_ERROR_IF (Naa > 1 && (aa_grid[0] != 0. || aa_grid[Naa - 1] != 360.),
1964 "The range of *aa_grid* must [0 360].");
1965
1966 // Get stokes dimension from *doit_scat_field*:
1967 const Index stokes_dim = doit_scat_field.ncols();
1968 ARTS_ASSERT(stokes_dim > 0 || stokes_dim < 5);
1969
1970 // To use special interpolation functions for atmospheric fields we
1971 // use ext_mat_field and abs_vec_field:
1972 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1973 1,
1974 1,
1975 stokes_dim,
1976 stokes_dim,
1977 0.);
1978 Tensor4 abs_vec_field(
1979 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
1980
1981 const Index Np = doit_scat_field.nvitrines();
1982
1983 Tensor5 doit_scat_ext_field(doit_scat_field.nvitrines(),
1984 doit_scat_field.nshelves(),
1985 doit_scat_field.nbooks(),
1986 doit_scat_field.npages(),
1987 doit_scat_field.nrows(),
1988 0.);
1989
1990 Index aa_index_local = 0;
1991
1992 // Calculate scattering extinction field
1993 for (Index za_index_local = 0; za_index_local < Nza;
1994 za_index_local++) {
1995 // This function has to be called inside the angular loop, as
1996 // spt_calc_agenda takes *za_index* and *aa_index*
1997 // from the workspace.
1999 ext_mat_field,
2000 abs_vec_field,
2001 spt_calc_agenda,
2002 za_index_local,
2003 aa_index_local,
2004 cloudbox_limits,
2005 t_field,
2006 pnd_field,
2007 verbosity);
2008
2009 for (Index p_index = 0;
2010 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
2011 p_index++) {
2012 // For all in p_grid (in cloudbox):
2013 // I_ext = (ext_mat_field - abs_vec_field) * cloudbox_field_mono
2014 // equivalent to:
2015 // I_ext = I * (K11-a1) + Q * (K12 - a2) + U * (K13 - a3) + V * (K14 - a4)
2016 for (Index i = 0; i < stokes_dim; i++) {
2017 doit_scat_ext_field(p_index, 0, 0, za_index_local, 0) +=
2018 cloudbox_field_mono(p_index, 0, 0, za_index_local, 0, i) *
2019 (ext_mat_field(p_index, 0, 0, 0, i) -
2020 abs_vec_field(p_index, 0, 0, i));
2021 }
2022 }
2023 }
2024
2025 Numeric corr_max = .0;
2026 Index corr_max_p_index = -1;
2027
2028 for (Index p_index = 0; p_index < Np; p_index++) {
2029 // Calculate scattering integrals
2030 const Numeric scat_int = AngIntegrate_trapezoid(
2031 doit_scat_field(p_index, 0, 0, joker, 0, 0), za_grid);
2032
2033 const Numeric scat_ext_int = AngIntegrate_trapezoid(
2034 doit_scat_ext_field(p_index, 0, 0, joker, 0), za_grid);
2035
2036 // Calculate factor between scattered extinction field integral
2037 // and scattered field integral
2038 const Numeric corr_factor = scat_ext_int / scat_int;
2039
2040 // If no scattering is present, the correction factor can become
2041 // inf or nan. We just don't apply it for those cases.
2042 if (!std::isnan(corr_factor) && !std::isinf(corr_factor)) {
2043 if (abs(corr_factor) > abs(corr_max)) {
2044 corr_max = corr_factor;
2045 corr_max_p_index = p_index;
2046 }
2047 if (norm_debug) {
2048 out0 << " DOIT corr_factor: " << 1. - corr_factor
2049 << " ( scat_ext_int: " << scat_ext_int
2050 << ", scat_int: " << scat_int << ")"
2051 << " at p_index " << p_index << "\n";
2052 }
2053 ARTS_USER_ERROR_IF (abs(1. - corr_factor) > norm_error_threshold,
2054 "ERROR: DOIT correction factor exceeds threshold (=",
2055 norm_error_threshold, "): ", setprecision(4),
2056 1. - corr_factor, " at p_index ", p_index, "\n")
2057 if (abs(1. - corr_factor) > norm_error_threshold / 2.) {
2058 out0 << " WARNING: DOIT correction factor above threshold/2: "
2059 << 1. - corr_factor << " at p_index " << p_index << "\n";
2060 }
2061
2062 // Scale scattered field with correction factor
2063 doit_scat_field(p_index, 0, 0, joker, 0, joker) *= corr_factor;
2064 } else if (norm_debug) {
2065 out0 << " DOIT corr_factor ignored: " << 1. - corr_factor
2066 << " ( scat_ext_int: " << scat_ext_int << ", scat_int: " << scat_int
2067 << ")"
2068 << " at p_index " << p_index << "\n";
2069 }
2070 }
2071
2072 ostringstream os;
2073 if (corr_max_p_index != -1) {
2074 os << " Max. DOIT correction factor in this iteration: " << 1. - corr_max
2075 << " at p_index " << corr_max_p_index << "\n";
2076 } else {
2077 os << " No DOIT correction performed in this iteration.\n";
2078 }
2079
2080 if (norm_debug)
2081 out0 << os.str();
2082 else
2083 out2 << os.str();
2084}
Declarations for agendas.
This file contains the definition of Array.
Constants of physical expressions as constexpr.
Common ARTS conversions.
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &f_grid, const Agenda &input_agenda)
Definition auto_md.cc:27543
void spt_calc_agendaExecute(Workspace &ws, ArrayOfPropagationMatrix &ext_mat_spt, ArrayOfStokesVector &abs_vec_spt, const Index scat_p_index, const Index scat_lat_index, const Index scat_lon_index, const Numeric rtp_temperature, const Index za_index, const Index aa_index, const Agenda &input_agenda)
Definition auto_md.cc:27719
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 surface_rtprop_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition auto_md.cc:27763
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
The Agenda class.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Workspace class.
bool is_inside_cloudbox(const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, const bool include_boundaries)
Definition cloudbox.cc:523
Internal cloudbox functions.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#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
constexpr Numeric RAD2DEG
Definition doit.cc:42
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 cloud_RT_no_background(Workspace &ws, Tensor6View cloudbox_field_mono, const Agenda &propmat_clearsky_agenda, const Ppath &ppath_step, ConstVectorView t_int, ConstMatrixView vmr_list_int, ConstTensor3View ext_mat_int, ConstMatrixView abs_vec_int, ConstMatrixView sca_vec_int, ConstMatrixView cloudbox_field_mono_int, ConstVectorView p_int, const ArrayOfIndex &cloudbox_limits, ConstVectorView f_grid, const Index &f_index, const Index &p_index, const Index &lat_index, const Index &lon_index, const Index &za_index, const Index &aa_index, const Verbosity &verbosity)
Calculates RT in the cloudbox (1D)
Definition doit.cc:1362
void cloudbox_field_ngAcceleration(Tensor6 &cloudbox_field_mono, const ArrayOfTensor6 &acceleration_input, const Index &accelerated, const Verbosity &)
Convergence acceleration.
Definition doit.cc:1596
void interp_cloud_coeff1D(Tensor3View ext_mat_int, MatrixView abs_vec_int, MatrixView sca_vec_int, MatrixView cloudbox_field_mono_int, VectorView t_int, MatrixView vmr_list_int, VectorView p_int, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, ConstTensor6View doit_scat_field, ConstTensor6View cloudbox_field_mono, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView p_grid, const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, ConstVectorView za_grid, const Index &scat_za_interp, const Verbosity &verbosity)
Interpolate all inputs of the VRTE on a propagation path step.
Definition doit.cc:1671
void rte_step_doit_replacement(VectorView stokes_vec, MatrixView trans_mat, const PropagationMatrix &ext_mat_av, const StokesVector &abs_vec_av, ConstVectorView sca_vec_av, const Numeric &lstep, const Numeric &rtp_planck_value, const bool &trans_is_precalc)
Solves monochromatic VRTE for an atmospheric slab with constant conditions.
Definition doit.cc:46
void cloud_RT_surface(Workspace &ws, Tensor6View cloudbox_field_mono, const Agenda &surface_rtprop_agenda, ConstVectorView f_grid, const Index &f_index, const Index &stokes_dim, const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, ConstVectorView za_grid, const Index &za_index)
Calculates RT in the cloudbox.
Definition doit.cc:1523
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
constexpr Numeric PI
Definition doit.cc:41
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.
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition geodetic.cc:1248
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
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 opt_prop_bulkCalc(PropagationMatrix &ext_mat, StokesVector &abs_vec, const ArrayOfPropagationMatrix &ext_mat_spt, const ArrayOfStokesVector &abs_vec_spt, const Tensor4 &pnd_field, const Index &scat_p_index, const Index &scat_lat_index, const Index &scat_lon_index, const Verbosity &)
WORKSPACE METHOD: opt_prop_bulkCalc.
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Declarations having to do with the four output streams.
#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.
void opt_prop_sum_propmat_clearsky(PropagationMatrix &ext_mat, StokesVector &abs_vec, const PropagationMatrix &propmat_clearsky)
Get optical properties from propmat_clearsky.
Numeric planck(const Numeric &f, const Numeric &t)
planck
This file contains declerations of functions of physical character.
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition ppath.cc:1460
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
Initiates a Ppath structure to hold the given number of points.
Definition ppath.cc:1393
Propagation path structure and functions.
Declaration of functions in rte.cc.
Contains sorting routines.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition sorting.h:39
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Header file for special_interp.cc.
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Vector lstep
The length between ppath points.
Vector r
Radius of each ppath point.
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Index dim
Atmospheric dimensionality.
This file contains basic functions to handle XML data files.