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