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