ARTS  2.2.66
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 
31 /*===========================================================================
32  === External declarations
33  ===========================================================================*/
34 
35 #include <stdexcept>
36 #include <iostream>
37 #include <cstdlib>
38 #include <cmath>
39 #include "array.h"
40 #include "auto_md.h"
41 #include "matpackVII.h"
42 #include "ppath.h"
43 #include "agenda_class.h"
44 #include "physics_funcs.h"
45 #include "lin_alg.h"
46 #include "math_funcs.h"
47 #include "messages.h"
48 #include "xml_io.h"
49 #include "rte.h"
50 #include "special_interp.h"
51 #include "doit.h"
52 #include "logic.h"
53 #include "check_input.h"
54 #include "sorting.h"
55 #include "cloudbox.h"
56 #include "geodetic.h"
57 
58 extern const Numeric PI;
59 extern const Numeric RAD2DEG;
60 
61 
62 
64 
108 void rte_step_doit(//Output and Input:
109  VectorView stokes_vec,
110  MatrixView trans_mat,
111  //Input
112  ConstMatrixView ext_mat_av,
113  ConstVectorView abs_vec_av,
114  ConstVectorView sca_vec_av,
115  const Numeric& lstep,
116  const Numeric& rtp_planck_value,
117  const bool& trans_is_precalc )
118 {
119  //Stokes dimension:
120  Index stokes_dim = stokes_vec.nelem();
121 
122  //Check inputs:
123  assert(is_size(trans_mat, stokes_dim, stokes_dim));
124  assert(is_size(ext_mat_av, stokes_dim, stokes_dim));
125  assert(is_size(abs_vec_av, stokes_dim));
126  assert(is_size(sca_vec_av, stokes_dim));
127  assert( rtp_planck_value >= 0 );
128  assert( lstep >= 0 );
129  assert (!is_singular( ext_mat_av ));
130 
131  // Check, if only the first component of abs_vec is non-zero:
132  bool unpol_abs_vec = true;
133 
134  for (Index i = 1; i < stokes_dim; i++)
135  if (abs_vec_av[i] != 0)
136  unpol_abs_vec = false;
137 
138  bool unpol_sca_vec = true;
139 
140  for (Index i = 1; i < stokes_dim; i++)
141  if (sca_vec_av[i] != 0)
142  unpol_sca_vec = false;
143 
144  // Calculate transmission by general function, if not precalculated
145  Index extmat_case = 0;
146  if( !trans_is_precalc )
147  { ext2trans( trans_mat, extmat_case, ext_mat_av, lstep ); }
148 
149  //--- Scalar case: ---------------------------------------------------------
150  if( stokes_dim == 1 )
151  {
152  stokes_vec[0] = stokes_vec[0] * trans_mat(0,0) +
153  ( abs_vec_av[0] * rtp_planck_value + sca_vec_av[0] ) /
154  ext_mat_av(0,0) * (1 - trans_mat(0,0) );
155  }
156 
157 
158  //--- Vector case: ---------------------------------------------------------
159 
160  // We have here two cases, diagonal or non-diagonal ext_mat_gas
161  // For diagonal ext_mat_gas, we expect abs_vec_gas to only have a
162  // non-zero value in position 1.
163 
164  //- Unpolarised
165  else if( extmat_case==1 && unpol_abs_vec && unpol_sca_vec )
166  {
167  // Stokes dim 1
168  stokes_vec[0] = stokes_vec[0] * trans_mat(0,0) +
169  ( abs_vec_av[0] * rtp_planck_value + sca_vec_av[0] ) /
170  ext_mat_av(0,0) * (1 - trans_mat(0,0) );
171 
172  // Stokes dims > 1
173  for( Index i=1; i<stokes_dim; i++ )
174  {
175  stokes_vec[i] = stokes_vec[i] * trans_mat(i,i) + sca_vec_av[i] /
176  ext_mat_av(i,i) * (1 - trans_mat(i,i));
177  }
178  }
179 
180 
181  //- General case
182  else
183  {
184  //Initialize internal variables:
185 
186  // Matrix LU used for LU decompostion and as dummy variable:
187  Matrix LU(stokes_dim, stokes_dim);
188  ArrayOfIndex indx(stokes_dim); // index for pivoting information
189  Vector b(stokes_dim); // dummy variable
190  Vector x(stokes_dim); // solution vector for K^(-1)*b
191  Matrix I(stokes_dim, stokes_dim);
192 
193  Vector B_abs_vec(stokes_dim);
194  B_abs_vec = abs_vec_av;
195  B_abs_vec *= rtp_planck_value;
196 
197  for (Index i=0; i<stokes_dim; i++)
198  b[i] = B_abs_vec[i] + sca_vec_av[i]; // b = abs_vec * B + sca_vec
199 
200  // solve K^(-1)*b = x
201  ludcmp(LU, indx, ext_mat_av);
202  lubacksub(x, LU, b, indx);
203 
204  Matrix ext_mat_ds(stokes_dim, stokes_dim);
205  ext_mat_ds = ext_mat_av;
206  ext_mat_ds *= -lstep; // ext_mat_ds = -ext_mat*ds
207 
208  Vector term1(stokes_dim);
209  Vector term2(stokes_dim);
210 
211  id_mat(I);
212  for(Index i=0; i<stokes_dim; i++)
213  {
214  for(Index j=0; j<stokes_dim; j++)
215  LU(i,j) = I(i,j) - trans_mat(i,j); // take LU as dummy variable
216  }
217  mult(term2, LU, x); // term2: second term of the solution of the RTE with
218  //fixed scattered field
219 
220  // term1: first term of solution of the RTE with fixed scattered field
221  mult( term1, trans_mat, stokes_vec );
222 
223  for (Index i=0; i<stokes_dim; i++)
224  stokes_vec[i] = term1[i] + term2[i]; // Compute the new Stokes Vector
225  }
226 }
227 
228 
229 
230 
232 
258  // Output and Input:
259  Tensor5View ext_mat_field,
260  Tensor4View abs_vec_field,
261  // Input:
262  const Agenda& spt_calc_agenda,
263  const Agenda& opt_prop_part_agenda,
264  const Index& scat_za_index,
265  const Index& scat_aa_index,
266  const ArrayOfIndex& cloudbox_limits,
267  ConstTensor3View t_field,
268  ConstTensor4View pnd_field,
269  const Verbosity& verbosity)
270 {
271  CREATE_OUT3;
272 
273  // Input variables are checked in the WSMs i_fieldUpdateSeqXXX, from
274  // where this function is called.
275 
276  out3 << "Calculate scattering properties in cloudbox \n";
277 
278  const Index atmosphere_dim = cloudbox_limits.nelem()/2;
279  const Index N_pt = pnd_field.nbooks();
280  const Index stokes_dim = ext_mat_field.ncols();
281 
282  assert( atmosphere_dim == 1 || atmosphere_dim ==3 );
283  assert( ext_mat_field.ncols() == ext_mat_field.nrows() &&
284  ext_mat_field.ncols() == abs_vec_field.ncols());
285 
286  const Index Np_cloud = cloudbox_limits[1]-cloudbox_limits[0]+1;
287 
288  // If atmosohere_dim == 1
289  Index Nlat_cloud = 1;
290  Index Nlon_cloud = 1;
291 
292  if (atmosphere_dim == 3)
293  {
294  Nlat_cloud = cloudbox_limits[3]-cloudbox_limits[2]+1;
295  Nlon_cloud = cloudbox_limits[5]-cloudbox_limits[4]+1;
296  }
297 
298  // Initialize ext_mat(_spt), abs_vec(_spt)
299  // Resize and initialize variables for storing optical properties
300  // of cloud particles
301  Matrix abs_vec_spt_local(N_pt, stokes_dim, 0.);
302  Tensor3 ext_mat_spt_local(N_pt, stokes_dim, stokes_dim, 0.);
303  Matrix abs_vec_local;
304  Tensor3 ext_mat_local;
305  Numeric rtp_temperature_local;
306 
307  // Calculate ext_mat, abs_vec for all points inside the cloudbox.
308  // sca_vec can be obtained from the workspace variable doit_scat_field.
309  // As we need the average for each layer, it makes sense to calculte
310  // the coefficients once and store them in an array instead of
311  // calculating at each point the coefficient of the point above and
312  // the point below.
313  // To use special interpolation functions for atmospheric fields we
314  // use ext_mat_field and abs_vec_field:
315 
316  // Loop over all positions inside the cloudbox defined by the
317  // cloudbox_limits.
318  for(Index scat_p_index_local = 0; scat_p_index_local < Np_cloud;
319  scat_p_index_local ++)
320  {
321  for(Index scat_lat_index_local = 0; scat_lat_index_local < Nlat_cloud;
322  scat_lat_index_local ++)
323  {
324  for(Index scat_lon_index_local = 0;
325  scat_lon_index_local < Nlon_cloud;
326  scat_lon_index_local ++)
327  {
328  if (atmosphere_dim == 1)
329  rtp_temperature_local =
330  t_field(scat_p_index_local + cloudbox_limits[0], 0, 0);
331  else
332  rtp_temperature_local =
333  t_field(scat_p_index_local + cloudbox_limits[0],
334  scat_lat_index_local + cloudbox_limits[2],
335  scat_lon_index_local + cloudbox_limits[4]);
336 
337  //Calculate optical properties for single particle types:
338  //( Execute agendas silently. )
339  spt_calc_agendaExecute(ws, ext_mat_spt_local,
340  abs_vec_spt_local,
341  scat_p_index_local,
342  scat_lat_index_local,
343  scat_lon_index_local,
344  rtp_temperature_local,
345  scat_za_index,
346  scat_aa_index,
347  spt_calc_agenda);
348 
349  opt_prop_part_agendaExecute(ws, ext_mat_local, abs_vec_local,
350  ext_mat_spt_local,
351  abs_vec_spt_local,
352  scat_p_index_local,
353  scat_lat_index_local,
354  scat_lon_index_local,
355  opt_prop_part_agenda);
356 
357  // Store coefficients in arrays for the whole cloudbox.
358  abs_vec_field(scat_p_index_local, scat_lat_index_local,
359  scat_lon_index_local,
360  joker) = abs_vec_local(0, joker);
361 
362  ext_mat_field(scat_p_index_local, scat_lat_index_local,
363  scat_lon_index_local,
364  joker, joker) = ext_mat_local(0, joker, joker);
365  }
366  }
367  }
368 }
369 
370 
371 
372 
373 
375 
417  // Input and output
418  Tensor6View doit_i_field,
419  // ppath_step_agenda:
420  const Index& p_index,
421  const Index& scat_za_index,
422  ConstVectorView scat_za_grid,
423  const ArrayOfIndex& cloudbox_limits,
424  ConstTensor6View doit_scat_field,
425  // Calculate scalar gas absorption:
426  const Agenda& propmat_clearsky_agenda,
427  ConstTensor4View vmr_field,
428  // Propagation path calculation:
429  const Agenda& ppath_step_agenda,
430  const Numeric& ppath_lraytrace,
431  ConstVectorView p_grid,
432  ConstTensor3View z_field,
433  ConstVectorView refellipsoid,
434  // Calculate thermal emission:
435  ConstTensor3View t_field,
436  ConstVectorView f_grid,
437  const Index& f_index,
438  //particle optical properties
439  ConstTensor5View ext_mat_field,
440  ConstTensor4View abs_vec_field,
441  const Agenda& surface_rtprop_agenda,
442  //const Agenda& surface_rtprop_agenda,
443  const Index& scat_za_interp,
444  const Verbosity& verbosity)
445 {
446  Matrix iy;
447  Matrix surface_emission;
448  Matrix surface_los;
449  Tensor4 surface_rmatrix;
450  Ppath ppath_step;
451  // Input variables are checked in the WSMs i_fieldUpdateSeqXXX, from
452  // where this function is called.
453 
454  //Initialize ppath for 1D.
455  ppath_init_structure(ppath_step, 1, 1);
456  // See documentation of ppath_init_structure for understanding
457  // the parameters.
458 
459  // Assign value to ppath.pos:
460  ppath_step.pos(0,0) = z_field(p_index,0,0);
461  ppath_step.r[0] = refellipsoid[0] + z_field(p_index,0,0);
462 
463  // Define the direction:
464  ppath_step.los(0,0) = scat_za_grid[scat_za_index];
465 
466 
467  // Define the grid positions:
468  ppath_step.gp_p[0].idx = p_index;
469  ppath_step.gp_p[0].fd[0] = 0;
470  ppath_step.gp_p[0].fd[1] = 1;
471 
472  // Call ppath_step_agenda:
473  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field, z_field,
474  vmr_field,
475  Vector(1,f_grid[f_index]), ppath_step_agenda );
476 
477  // Check whether the next point is inside or outside the
478  // cloudbox. Only if the next point lies inside the
479  // cloudbox a radiative transfer step caclulation has to
480  // be performed.
481 
482  if((cloudbox_limits[0] <= ppath_step.gp_p[1].idx &&
483  cloudbox_limits[1] > ppath_step.gp_p[1].idx) ||
484  (cloudbox_limits[1] == ppath_step.gp_p[1].idx &&
485  abs(ppath_step.gp_p[1].fd[0]) < 1e-6))
486  {
487  // Stokes dimension
488  const Index stokes_dim = doit_i_field.ncols();
489  // Number of species
490  const Index N_species = vmr_field.nbooks();
491 
492  // Ppath_step normally has 2 points, the starting
493  // point and the intersection point.
494  // But there can be points in between, when a maximum
495  // lstep is given. We have to interpolate on all the
496  // points in the ppath_step.
497 
498  // Initialize variables for interpolated values
499  Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.np, 0.);
500  Matrix abs_vec_int(stokes_dim, ppath_step.np, 0.);
501  Matrix sca_vec_int(stokes_dim, ppath_step.np, 0.);
502  Matrix doit_i_field_int(stokes_dim, ppath_step.np, 0.);
503  Vector t_int(ppath_step.np, 0.);
504  Matrix vmr_list_int(N_species, ppath_step.np, 0.);
505  Vector p_int(ppath_step.np, 0.);
506 
507  interp_cloud_coeff1D(ext_mat_int, abs_vec_int, sca_vec_int,
508  doit_i_field_int, t_int, vmr_list_int, p_int,
509  ext_mat_field, abs_vec_field, doit_scat_field,
510  doit_i_field, t_field, vmr_field, p_grid,
511  ppath_step, cloudbox_limits, scat_za_grid,
512  scat_za_interp, verbosity);
513 
514 
515  // ppath_what_background(ppath_step) tells the
516  // radiative background. More information in the
517  // function get_iy_of_background.
518  // if there is no background we proceed the RT
519  Index bkgr = ppath_what_background(ppath_step);
520 
521  // Radiative transfer from one layer to the next, starting
522  // at the intersection with the next layer and propagating
523  // to the considered point.
524  cloud_RT_no_background(ws, doit_i_field,
525  propmat_clearsky_agenda,
526  ppath_step,
527  t_int, vmr_list_int,
528  ext_mat_int, abs_vec_int, sca_vec_int,
529  doit_i_field_int,
530  p_int, cloudbox_limits,
531  f_grid, f_index, p_index, 0, 0,
532  scat_za_index, 0, verbosity);
533 
534  // bkgr=2 indicates that the background is the surface
535  if (bkgr == 2)
536  {
537  // cout << "hit surface "<< ppath_step.gp_p << endl;
538  cloud_RT_surface(ws,
539  doit_i_field, surface_rtprop_agenda, f_grid,
540  f_index, stokes_dim, ppath_step, cloudbox_limits,
541  scat_za_grid, scat_za_index);
542 
543  }
544 
545  }//end if inside cloudbox
546 }
547 
549 /*
550  Basically the same as cloud_ppath_update1D, the only difference is that
551  i_field_old is always used as incoming Stokes vector.
552 
553  \author Claudia Emde
554  \date 2005-05-04
555 */
557  // Output
558  Tensor6View doit_i_field,
559  // ppath_step_agenda:
560  const Index& p_index,
561  const Index& scat_za_index,
562  ConstVectorView scat_za_grid,
563  const ArrayOfIndex& cloudbox_limits,
564  ConstTensor6View doit_i_field_old,
565  ConstTensor6View doit_scat_field,
566  // Calculate scalar gas absorption:
567  const Agenda& propmat_clearsky_agenda,
568  ConstTensor4View vmr_field,
569  // Propagation path calculation:
570  const Agenda& ppath_step_agenda,
571  const Numeric& ppath_lraytrace,
572  ConstVectorView p_grid,
573  ConstTensor3View z_field,
574  ConstVectorView refellipsoid,
575  // Calculate thermal emission:
576  ConstTensor3View t_field,
577  ConstVectorView f_grid,
578  // used for surface ?
579  const Index& f_index,
580  //particle optical properties
581  ConstTensor5View ext_mat_field,
582  ConstTensor4View abs_vec_field,
583  const Agenda& surface_rtprop_agenda,
584  const Index& scat_za_interp,
585  const Verbosity& verbosity)
586 {
587  Matrix iy;
588  Matrix surface_emission;
589  Matrix surface_los;
590  Tensor4 surface_rmatrix;
591  Ppath ppath_step;
592  // Input variables are checked in the WSMs i_fieldUpdateSeqXXX, from
593  // where this function is called.
594 
595  //Initialize ppath for 1D.
596  ppath_init_structure(ppath_step, 1, 1);
597  // See documentation of ppath_init_structure for understanding
598  // the parameters.
599 
600  // Assign value to ppath.pos:
601  ppath_step.pos(0,0) = z_field(p_index,0,0);
602  ppath_step.r[0] = refellipsoid[0] + z_field(p_index,0,0);
603 
604  // Define the direction:
605  ppath_step.los(0,0) = scat_za_grid[scat_za_index];
606 
607 
608  // Define the grid positions:
609  ppath_step.gp_p[0].idx = p_index;
610  ppath_step.gp_p[0].fd[0] = 0;
611  ppath_step.gp_p[0].fd[1] = 1;
612 
613  // Call ppath_step_agenda:
614  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field, z_field,
615  vmr_field,
616  Vector(1,f_grid[f_index]), ppath_step_agenda );
617 
618  // Check whether the next point is inside or outside the
619  // cloudbox. Only if the next point lies inside the
620  // cloudbox a radiative transfer step caclulation has to
621  // be performed.
622 
623  if((cloudbox_limits[0] <= ppath_step.gp_p[1].idx &&
624  cloudbox_limits[1] > ppath_step.gp_p[1].idx) ||
625  (cloudbox_limits[1] == ppath_step.gp_p[1].idx &&
626  abs(ppath_step.gp_p[1].fd[0]) < 1e-6))
627  {
628  // Stokes dimension
629  const Index stokes_dim = doit_i_field.ncols();
630  // Number of species
631  const Index N_species = vmr_field.nbooks();
632 
633  // Ppath_step normally has 2 points, the starting
634  // point and the intersection point.
635  // But there can be points in between, when a maximum
636  // lstep is given. We have to interpolate on all the
637  // points in the ppath_step.
638 
639  // Initialize variables for interpolated values
640  Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.np, 0.);
641  Matrix abs_vec_int(stokes_dim, ppath_step.np, 0.);
642  Matrix sca_vec_int(stokes_dim, ppath_step.np, 0.);
643  Matrix doit_i_field_int(stokes_dim, ppath_step.np, 0.);
644  Vector t_int(ppath_step.np, 0.);
645  Matrix vmr_list_int(N_species, ppath_step.np, 0.);
646  Vector p_int(ppath_step.np, 0.);
647 
648  interp_cloud_coeff1D(ext_mat_int, abs_vec_int, sca_vec_int,
649  doit_i_field_int, t_int, vmr_list_int, p_int,
650  ext_mat_field, abs_vec_field, doit_scat_field,
651  doit_i_field_old, t_field, vmr_field, p_grid,
652  ppath_step, cloudbox_limits, scat_za_grid,
653  scat_za_interp, verbosity);
654 
655  // ppath_what_background(ppath_step) tells the
656  // radiative background. More information in the
657  // function get_iy_of_background.
658  // if there is no background we proceed the RT
659  Index bkgr = ppath_what_background(ppath_step);
660 
661  // if 0, there is no background
662  // do this in any case. cause we need downwelling doit_i_field
663  // at the surface for calculating surface scattering
664 
665  // Radiative transfer from one layer to the next, starting
666  // at the intersection with the next layer and propagating
667  // to the considered point.
668  cloud_RT_no_background(ws, doit_i_field,
669  propmat_clearsky_agenda,
670  ppath_step,
671  t_int, vmr_list_int,
672  ext_mat_int, abs_vec_int, sca_vec_int,
673  doit_i_field_int,
674  p_int, cloudbox_limits,
675  f_grid, f_index, p_index, 0, 0,
676  scat_za_index, 0, verbosity);
677 
678  if (bkgr == 2)
679  {
680  cloud_RT_surface( ws,
681  doit_i_field, surface_rtprop_agenda, f_grid,
682  f_index, stokes_dim, ppath_step, cloudbox_limits,
683  scat_za_grid, scat_za_index);
684 
685  }
686  }//end if inside cloudbox
687 }
688 
689 
690 
691 
693 
739  Tensor6View doit_i_field,
740  // ppath_step_agenda:
741  const Index& p_index,
742  const Index& lat_index,
743  const Index& lon_index,
744  const Index& scat_za_index,
745  const Index& scat_aa_index,
746  ConstVectorView scat_za_grid,
747  ConstVectorView scat_aa_grid,
748  const ArrayOfIndex& cloudbox_limits,
749  ConstTensor6View doit_scat_field,
750  // Calculate scalar gas absorption:
751  const Agenda& propmat_clearsky_agenda,
752  ConstTensor4View vmr_field,
753  // Propagation path calculation:
754  const Agenda& ppath_step_agenda,
755  const Numeric& ppath_lraytrace,
756  ConstVectorView p_grid,
757  ConstVectorView lat_grid,
758  ConstVectorView lon_grid,
759  ConstTensor3View z_field,
760  ConstVectorView refellipsoid,
761  // Calculate thermal emission:
762  ConstTensor3View t_field,
763  ConstVectorView f_grid,
764  const Index& f_index,
765  //particle optical properties
766  ConstTensor5View ext_mat_field,
767  ConstTensor4View abs_vec_field,
768  const Index&, //scat_za_interp
769  const Verbosity& verbosity
770  )
771 {
772  CREATE_OUT3;
773 
774  Ppath ppath_step;
775  const Index stokes_dim = doit_i_field.ncols();
776 
777  Vector sca_vec_av(stokes_dim,0);
778  Vector aa_grid(scat_aa_grid.nelem());
779 
780  for(Index i = 0; i<scat_aa_grid.nelem(); i++)
781  aa_grid[i] = scat_aa_grid[i]-180.;
782 
783  //Initialize ppath for 3D.
784  ppath_init_structure(ppath_step, 3, 1);
785  // See documentation of ppath_init_structure for
786  // understanding the parameters.
787 
788  // The first dimension of pos are the points in
789  // the propagation path.
790  // Here we initialize the first point.
791  // The second is: radius, latitude, longitude
792 
793 
794  ppath_step.pos(0,2) = lon_grid[lon_index];
795  ppath_step.pos(0,1) = lat_grid[lat_index];
796  ppath_step.pos(0,0) = z_field( p_index, lat_index, lon_index );
797  // As always on top of the lat. grid positions, OK to call refell2r:
798  ppath_step.r[0] = refell2r( refellipsoid, ppath_step.pos(0,1) ) +
799  ppath_step.pos(0,0);
800 
801 
802  // Define the direction:
803  ppath_step.los(0,0) = scat_za_grid[scat_za_index];
804  ppath_step.los(0,1) = aa_grid[scat_aa_index];
805 
806  // Define the grid positions:
807  ppath_step.gp_p[0].idx = p_index;
808  ppath_step.gp_p[0].fd[0] = 0.;
809  ppath_step.gp_p[0].fd[1] = 1.;
810 
811  ppath_step.gp_lat[0].idx = lat_index;
812  ppath_step.gp_lat[0].fd[0] = 0.;
813  ppath_step.gp_lat[0].fd[1] = 1.;
814 
815  ppath_step.gp_lon[0].idx = lon_index;
816  ppath_step.gp_lon[0].fd[0] = 0.;
817  ppath_step.gp_lon[0].fd[1] = 1.;
818 
819  // Call ppath_step_agenda:
820  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field, z_field,
821  vmr_field,
822  Vector(1,f_grid[f_index]), ppath_step_agenda);
823 
824  // Check whether the next point is inside or outside the
825  // cloudbox. Only if the next point lies inside the
826  // cloudbox a radiative transfer step caclulation has to
827  // be performed.
828  if (is_inside_cloudbox(ppath_step, cloudbox_limits, true))
829  {
830  // Gridpositions inside the cloudbox.
831  // The optical properties are stored only inside the
832  // cloudbox. For interpolation we use grids
833  // inside the cloudbox.
834 
835  ArrayOfGridPos cloud_gp_p = ppath_step.gp_p;
836  ArrayOfGridPos cloud_gp_lat = ppath_step.gp_lat;
837  ArrayOfGridPos cloud_gp_lon = ppath_step.gp_lon;
838 
839  for(Index i = 0; i<ppath_step.np; i++ )
840  {
841  cloud_gp_p[i].idx -= cloudbox_limits[0];
842  cloud_gp_lat[i].idx -= cloudbox_limits[2];
843  cloud_gp_lon[i].idx -= cloudbox_limits[4];
844  }
845  const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
846  const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
847  const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
848  gridpos_upperend_check( cloud_gp_p[0], n1 );
849  gridpos_upperend_check( cloud_gp_p[ppath_step.np-1], n1 );
850  gridpos_upperend_check( cloud_gp_lat[0], n2 );
851  gridpos_upperend_check( cloud_gp_lat[ppath_step.np-1], n2);
852  gridpos_upperend_check( cloud_gp_lon[0], n3 );
853  gridpos_upperend_check( cloud_gp_lon[ppath_step.np-1], n3 );
854 
855  Matrix itw(ppath_step.np, 8);
856  interpweights(itw, cloud_gp_p, cloud_gp_lat, cloud_gp_lon);
857 
858  Matrix itw_p(ppath_step.np, 2);
859  interpweights(itw_p, cloud_gp_p);
860 
861  // The zenith angles and azimuth of the propagation path are
862  // needed as we have to
863  // interpolate the intensity field and the scattered field on the
864  // right angles.
865  VectorView los_grid_za = ppath_step.los(joker,0);
866  VectorView los_grid_aa = ppath_step.los(joker,1);
867 
868  for(Index i = 0; i<los_grid_aa.nelem(); i++)
869  los_grid_aa[i] = los_grid_aa[i] + 180.;
870 
871  ArrayOfGridPos gp_za(los_grid_za.nelem());
872  gridpos(gp_za, scat_za_grid, los_grid_za);
873 
874  ArrayOfGridPos gp_aa(los_grid_aa.nelem());
875  gridpos(gp_aa, scat_aa_grid, los_grid_aa);
876 
877  Matrix itw_p_za(ppath_step.np, 32);
878  interpweights(itw_p_za, cloud_gp_p, cloud_gp_lat, cloud_gp_lon,
879  gp_za, gp_aa);
880 
881  // Ppath_step normally has 2 points, the starting
882  // point and the intersection point.
883  // But there can be points in between, when a maximum
884  // lstep is given. We have to interpolate on all the
885  // points in the ppath_step.
886 
887  Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.np);
888  Matrix abs_vec_int(stokes_dim, ppath_step.np);
889  Matrix sca_vec_int(stokes_dim, ppath_step.np, 0.);
890  Matrix doit_i_field_int(stokes_dim, ppath_step.np, 0.);
891  Vector t_int(ppath_step.np);
892  Vector vmr_int(ppath_step.np);
893  Vector p_int(ppath_step.np);
894  Vector stokes_vec(stokes_dim);
895  //Tensor3 ext_mat_gas(stokes_dim, stokes_dim, ppath_step.np);
896  //Matrix abs_vec_gas(stokes_dim, ppath_step.np);
897 
898  // Calculate the average of the coefficients for the layers
899  // to be considered in the
900  // radiative transfer calculation.
901 
902  for (Index i = 0; i < stokes_dim; i++)
903  {
904  // Extinction matrix requires a second loop
905  // over stokes_dim
906  out3 << "Interpolate ext_mat:\n";
907  for (Index j = 0; j < stokes_dim; j++)
908  {
909  //
910  // Interpolation of ext_mat
911  //
912  interp( ext_mat_int(i, j, joker), itw,
913  ext_mat_field(joker, joker, joker, i, j), cloud_gp_p,
914  cloud_gp_lat, cloud_gp_lon);
915  }
916  // Absorption vector:
917  //
918  // Interpolation of abs_vec
919  //
920  interp( abs_vec_int(i,joker), itw,
921  abs_vec_field(joker, joker, joker, i),
922  cloud_gp_p, cloud_gp_lat, cloud_gp_lon);
923  //
924  // Scattered field:
925  //
926  // Interpolation of sca_vec:
927  //
928  out3 << "Interpolate doit_scat_field:\n";
929  interp( sca_vec_int(i, joker), itw_p_za,
930  doit_scat_field(joker, joker, joker, joker, joker, i),
931  cloud_gp_p,
932  cloud_gp_lat, cloud_gp_lon, gp_za, gp_aa);
933  out3 << "Interpolate doit_i_field:\n";
934  interp( doit_i_field_int(i, joker), itw_p_za,
935  doit_i_field(joker, joker, joker, joker, joker, i),
936  cloud_gp_p,
937  cloud_gp_lat, cloud_gp_lon, gp_za, gp_aa);
938  }
939  //
940  // Planck function
941  //
942  // Interpolate temperature field
943  //
944  out3 << "Interpolate temperature field\n";
945  interp( t_int, itw,
946  t_field(joker, joker, joker), ppath_step.gp_p,
947  ppath_step.gp_lat, ppath_step.gp_lon);
948 
949  //
950  // The vmr_field is needed for the gaseous absorption
951  // calculation.
952  //
953  const Index N_species = vmr_field.nbooks();
954  //
955  // Interpolated vmr_list, holds a vmr_list for each point in
956  // ppath_step.
957  //
958  Matrix vmr_list_int(N_species, ppath_step.np);
959 
960  for (Index i = 0; i < N_species; i++)
961  {
962  out3 << "Interpolate vmr field\n";
963  interp( vmr_int, itw,
964  vmr_field(i, joker, joker, joker), ppath_step.gp_p,
965  ppath_step.gp_lat, ppath_step.gp_lon );
966 
967  vmr_list_int(i, joker) = vmr_int;
968  }
969 
970  // Presssure (needed for the calculation of gas absorption)
971  itw2p( p_int, p_grid, ppath_step.gp_p, itw_p);
972 
973  out3 << "Calculate radiative transfer inside cloudbox.\n";
974  cloud_RT_no_background(ws, doit_i_field,
975  propmat_clearsky_agenda,
976  ppath_step,
977  t_int, vmr_list_int,
978  ext_mat_int, abs_vec_int, sca_vec_int,
979  doit_i_field_int,
980  p_int, cloudbox_limits,
981  f_grid, f_index, p_index, lat_index, lon_index,
982  scat_za_index, scat_aa_index, verbosity);
983  }//end if inside cloudbox
984 }
985 
987 /*
988  This function calculates RT in the cloudbox (1D) if the next intersected
989  level is an atmospheric level (in contrast to the surface).
990  It is used inside the functions cloud_ppath_update1DXXX.
991 
992  Output:
993  \param doit_i_field Radiation field in cloudbox. This variable is filled
994  with the updated values for a given zenith angle (scat_za_index) and
995  pressure (p_index).
996  Input:
997  \param propmat_clearsky_agenda Calculate gas absorption.
998  \param ppath_step Propagation path step from one pressure level to the next
999  (this can include several points)
1000  \param t_int Temperature values interpolated on propagation path points.
1001  \param vmr_list_int Interpolated volume mixing ratios.
1002  \param ext_mat_int Interpolated particle extinction matrix.
1003  \param abs_vec_int Interpolated particle absorption vector.
1004  \param sca_vec_int Interpolated particle scattering vector.
1005  \param doit_i_field_int Interpolated radiances.
1006  \param p_int Interpolated pressure values.
1007  \param cloudbox_limits Cloudbox_limits.
1008  \param f_grid Frequency grid.
1009  \param f_index Frequency index of (monochromatic) scattering calculation.
1010  \param p_index Pressure index in *doit_i_field*.
1011  \param scat_za_index Zenith angle index in *doit_i_field*.
1012 
1013  \author Claudia Emde
1014  \date 2005-05-13
1015 */
1017  //Output
1018  Tensor6View doit_i_field,
1019  // Input
1020  const Agenda& propmat_clearsky_agenda,
1021  const Ppath& ppath_step,
1022  ConstVectorView t_int,
1023  ConstMatrixView vmr_list_int,
1024  ConstTensor3View ext_mat_int,
1025  ConstMatrixView abs_vec_int,
1026  ConstMatrixView sca_vec_int,
1027  ConstMatrixView doit_i_field_int,
1028  ConstVectorView p_int,
1029  const ArrayOfIndex& cloudbox_limits,
1030  ConstVectorView f_grid,
1031  const Index& f_index,
1032  const Index& p_index,
1033  const Index& lat_index,
1034  const Index& lon_index,
1035  const Index& scat_za_index,
1036  const Index& scat_aa_index,
1037  const Verbosity& verbosity)
1038 {
1039  CREATE_OUT3;
1040 
1041  const Index N_species = vmr_list_int.nrows();
1042  const Index stokes_dim = doit_i_field.ncols();
1043  const Index atmosphere_dim = cloudbox_limits.nelem()/2;
1044 
1045  Vector sca_vec_av(stokes_dim,0);
1046  Vector stokes_vec(stokes_dim, 0.);
1047  Vector rtp_vmr_local(N_species,0.);
1048 
1049  // Two propmat_clearsky to average between
1050  Tensor4 cur_propmat_clearsky;
1051  Tensor4 prev_propmat_clearsky;
1052 
1053  Tensor3 ext_mat_local;
1054  Matrix abs_vec_local;
1055 
1056  // Incoming stokes vector
1057  stokes_vec = doit_i_field_int(joker, ppath_step.np-1);
1058 
1059  for( Index k = ppath_step.np-1; k >= 0; k--)
1060  {
1061  // Save propmat_clearsky from previous level by
1062  // swapping it with current level
1063  swap(cur_propmat_clearsky, prev_propmat_clearsky);
1064 
1065  //
1066  // Calculate scalar gas absorption
1067  //
1068  const Vector rtp_mag_dummy(3,0);
1069  const Vector ppath_los_dummy;
1070 
1071  propmat_clearsky_agendaExecute( ws, cur_propmat_clearsky,
1072  f_grid[Range(f_index, 1)],
1073  rtp_mag_dummy, ppath_los_dummy,
1074  p_int[k],
1075  t_int[k],
1076  vmr_list_int(joker,k),
1077  propmat_clearsky_agenda );
1078 
1079  // Skip any further calculations for the first point.
1080  // We need values at two ppath points before we can average.
1081  if (k == ppath_step.np-1)
1082  continue;
1083 
1084  // Average prev_propmat_clearsky with cur_propmat_clearsky
1085  prev_propmat_clearsky += cur_propmat_clearsky;
1086  prev_propmat_clearsky *= 0.5;
1087 
1088  opt_prop_sum_propmat_clearsky(ext_mat_local, abs_vec_local,
1089  prev_propmat_clearsky);
1090 
1091  //
1092  // Add average particle extinction to ext_mat.
1093  //
1094  for (Index i = 0; i < stokes_dim; i++)
1095  {
1096  for (Index j = 0; j < stokes_dim; j++)
1097  {
1098  ext_mat_local(0,i,j) += 0.5 *
1099  (ext_mat_int(i,j,k) + ext_mat_int(i,j,k+1));
1100  }
1101  //
1102  // Add average particle absorption to abs_vec.
1103  //
1104  abs_vec_local(0,i) += 0.5 *
1105  (abs_vec_int(i,k) + abs_vec_int(i,k+1));
1106 
1107  //
1108  // Averaging of sca_vec:
1109  //
1110  sca_vec_av[i] = 0.5 *
1111  (sca_vec_int(i, k) + sca_vec_int(i, k+1));
1112 
1113  }
1114  // Frequency
1115  Numeric f = f_grid[f_index];
1116  //
1117  // Calculate Planck function
1118  //
1119  Numeric rte_planck_value = planck(f, 0.5 * (t_int[k] + t_int[k+1]));
1120 
1121  // Length of the path between the two layers.
1122  Numeric lstep = ppath_step.lstep[k];
1123 
1124  // Some messages:
1125  out3 << "-----------------------------------------\n";
1126  out3 << "Input for radiative transfer step \n"
1127  << "calculation inside"
1128  << " the cloudbox:" << "\n";
1129  out3 << "Stokes vector at intersection point: \n"
1130  << stokes_vec
1131  << "\n";
1132  out3 << "lstep: ..." << lstep << "\n";
1133  out3 << "------------------------------------------\n";
1134  out3 << "Averaged coefficients: \n";
1135  out3 << "Planck function: " << rte_planck_value << "\n";
1136  out3 << "Scattering vector: " << sca_vec_av << "\n";
1137  out3 << "Absorption vector: " << abs_vec_local(0,joker) << "\n";
1138  out3 << "Extinction matrix: " << ext_mat_local(0,joker,joker) << "\n";
1139 
1140 
1141  assert (!is_singular( ext_mat_local(0,joker,joker)));
1142 
1143  // Radiative transfer step calculation. The Stokes vector
1144  // is updated until the considered point is reached.
1145  rte_step_doit(stokes_vec, Matrix(stokes_dim,stokes_dim),
1146  ext_mat_local(0,joker,joker), abs_vec_local(0,joker),
1147  sca_vec_av, lstep, rte_planck_value);
1148 
1149  }// End of loop over ppath_step.
1150  // Assign calculated Stokes Vector to doit_i_field.
1151  if (atmosphere_dim == 1)
1152  doit_i_field(p_index - cloudbox_limits[0], 0, 0, scat_za_index, 0, joker)
1153  = stokes_vec;
1154  else if (atmosphere_dim == 3)
1155  doit_i_field(p_index - cloudbox_limits[0],
1156  lat_index - cloudbox_limits[2],
1157  lon_index - cloudbox_limits[4],
1158  scat_za_index, scat_aa_index,
1159  joker) = stokes_vec;
1160 
1161 }
1162 
1163 
1164 
1166 /*
1167  This function calculates RT in the cloudbox if the next intersected
1168  level is the surface.
1169 
1170  CE (2006-05-29) Included surface_rtprop_agenda here.
1171 
1172  \author Claudia Emde
1173  \date 2005-05-13
1174 */
1176  //Output
1177  Tensor6View doit_i_field,
1178  //Input
1179  const Agenda& surface_rtprop_agenda,
1180  ConstVectorView f_grid,
1181  const Index& f_index,
1182  const Index& stokes_dim,
1183  const Ppath& ppath_step,
1184  const ArrayOfIndex& cloudbox_limits,
1185  ConstVectorView scat_za_grid,
1186  const Index& scat_za_index
1187  )
1188 {
1189  chk_not_empty( "surface_rtprop_agenda", surface_rtprop_agenda );
1190 
1191  Matrix iy;
1192 
1193  // Local output of surface_rtprop_agenda.
1194  Matrix surface_emission;
1195  Matrix surface_los;
1196  Tensor4 surface_rmatrix;
1197 
1198 
1199  //Set rte_pos and rte_los to match the last point in ppath.
1200 
1201  Index np = ppath_step.np;
1202 
1203  Vector rte_pos; // ppath_step.pos contains two columns for 1D
1204  rte_pos.resize( ppath_step.dim );
1205  rte_pos = ppath_step.pos(np-1,Range(0,ppath_step.dim));
1206 
1207  Vector rte_los;
1208  rte_los.resize( ppath_step.los.ncols() );
1209  rte_los = ppath_step.los(np-1,joker);
1210 
1211  //Execute the surface_rtprop_agenda which gives the surface
1212  //parameters.
1213 
1214  surface_rtprop_agendaExecute( ws, surface_emission, surface_los,
1215  surface_rmatrix, Vector(1,f_grid[f_index]),
1216  rte_pos, rte_los, surface_rtprop_agenda );
1217 
1218  iy = surface_emission;
1219 
1220  Index nlos = surface_los.nrows();
1221 
1222  if( nlos > 0 )
1223  {
1224  Vector rtmp(stokes_dim); // Reflected Stokes vector for 1 frequency
1225 
1226  for( Index ilos=0; ilos<nlos; ilos++ )
1227  {
1228  // Several things needs to be fixed here. As far as I understand it,
1229  // this works only for specular cases and if the lower cloudbox limit
1230  // is exactly at the surface (PE, 120828)
1231 
1232  mult( rtmp, surface_rmatrix(ilos,0,joker,joker),
1233  doit_i_field( cloudbox_limits[0], 0, 0,
1234  (scat_za_grid.nelem() -1 - scat_za_index), 0, joker) );
1235  iy(0,joker) += rtmp;
1236 
1237  doit_i_field( cloudbox_limits[0], 0, 0, scat_za_index, 0, joker ) =
1238  iy( 0, joker );
1239  }
1240  }
1241 }
1242 
1243 
1244 
1246 
1254  Tensor3View ext_mat_int,
1255  MatrixView abs_vec_int,
1256  MatrixView sca_vec_int,
1257  MatrixView doit_i_field_int,
1258  VectorView t_int,
1259  MatrixView vmr_list_int,
1260  VectorView p_int,
1261  //Input
1262  ConstTensor5View ext_mat_field,
1263  ConstTensor4View abs_vec_field,
1264  ConstTensor6View doit_scat_field,
1265  ConstTensor6View doit_i_field,
1266  ConstTensor3View t_field,
1267  ConstTensor4View vmr_field,
1268  ConstVectorView p_grid,
1269  const Ppath& ppath_step,
1270  const ArrayOfIndex& cloudbox_limits,
1271  ConstVectorView scat_za_grid,
1272  const Index& scat_za_interp,
1273  const Verbosity& verbosity)
1274 {
1275  CREATE_OUT3;
1276 
1277  // Stokes dimension
1278  const Index stokes_dim = doit_i_field.ncols();
1279 
1280  // Gridpositions inside the cloudbox.
1281  // The optical properties are stored only inside the
1282  // cloudbox. For interpolation we use grids
1283  // inside the cloudbox.
1284  ArrayOfGridPos cloud_gp_p = ppath_step.gp_p;
1285 
1286  for(Index i = 0; i < ppath_step.np; i++ )
1287  cloud_gp_p[i].idx -= cloudbox_limits[0];
1288 
1289  // Grid index for points at upper limit of cloudbox must be shifted
1290  const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
1291  gridpos_upperend_check( cloud_gp_p[0], n1 );
1292  gridpos_upperend_check( cloud_gp_p[ppath_step.np-1], n1 );
1293 
1294 
1295  Matrix itw(cloud_gp_p.nelem(),2);
1296  interpweights( itw, cloud_gp_p );
1297 
1298  // The zenith angles of the propagation path are needed as we have to
1299  // interpolate the intensity field and the scattered field on the
1300  // right angles.
1301  Vector los_grid = ppath_step.los(joker,0);
1302 
1303  ArrayOfGridPos gp_za(los_grid.nelem());
1304  gridpos(gp_za, scat_za_grid, los_grid);
1305 
1306  Matrix itw_p_za(cloud_gp_p.nelem(), 4);
1307  interpweights(itw_p_za, cloud_gp_p, gp_za );
1308 
1309  // Calculate the average of the coefficients for the layers
1310  // to be considered in the
1311  // radiative transfer calculation.
1312 
1313  for (Index i = 0; i < stokes_dim; i++)
1314  {
1315  // Extinction matrix requires a second loop
1316  // over stokes_dim
1317  out3 << "Interpolate ext_mat:\n";
1318  for (Index j = 0; j < stokes_dim; j++)
1319  {
1320  //
1321  // Interpolation of ext_mat
1322  //
1323  interp( ext_mat_int(i, j, joker), itw,
1324  ext_mat_field(joker, 0, 0, i, j), cloud_gp_p );
1325  }
1326  // Particle absorption vector:
1327  //
1328  // Interpolation of abs_vec
1329  // //
1330  out3 << "Interpolate abs_vec:\n";
1331  interp( abs_vec_int(i,joker), itw,
1332  abs_vec_field(joker, 0, 0, i), cloud_gp_p );
1333  //
1334  // Scattered field:
1335  //
1336  //
1337 
1338  out3 << "Interpolate doit_scat_field and doit_i_field:\n";
1339  if(scat_za_interp == 0) // linear interpolation
1340  {
1341  interp( sca_vec_int(i, joker), itw_p_za,
1342  doit_scat_field(joker, 0, 0, joker, 0, i),
1343  cloud_gp_p, gp_za);
1344  interp( doit_i_field_int(i, joker), itw_p_za,
1345  doit_i_field(joker, 0, 0, joker, 0, i), cloud_gp_p, gp_za);
1346  }
1347  else if (scat_za_interp == 1) //polynomial interpolation
1348  {
1349  // These intermediate variables are needed because polynomial
1350  // interpolation is not implemented as multidimensional
1351  // interpolation.
1352  Tensor3 sca_vec_int_za(stokes_dim, ppath_step.np,
1353  scat_za_grid.nelem(), 0.);
1354  Tensor3 doit_i_field_int_za(stokes_dim, ppath_step.np,
1355  scat_za_grid.nelem(), 0.);
1356  for (Index za = 0; za < scat_za_grid.nelem(); za++)
1357  {
1358  interp( sca_vec_int_za(i, joker, za), itw,
1359  doit_scat_field(joker, 0, 0, za, 0, i), cloud_gp_p);
1360  out3 << "Interpolate doit_i_field:\n";
1361  interp( doit_i_field_int_za(i, joker, za), itw,
1362  doit_i_field(joker, 0, 0, za, 0, i), cloud_gp_p);
1363  }
1364  for (Index ip = 0; ip < ppath_step.np; ip ++)
1365  {
1366  sca_vec_int(i, ip) =
1367  interp_poly(scat_za_grid, sca_vec_int_za(i, ip, joker),
1368  los_grid[ip], gp_za[ip]);
1369  doit_i_field_int(i, ip) =
1370  interp_poly(scat_za_grid, doit_i_field_int_za(i, ip, joker),
1371  los_grid[ip], gp_za[ip]);
1372  }
1373  }
1374  }
1375  //
1376  // Planck function
1377  //
1378  // Interpolate temperature field
1379  //
1380  out3 << "Interpolate temperature field\n";
1381  interp( t_int, itw, t_field(joker, 0, 0), ppath_step.gp_p );
1382  //
1383  // The vmr_field is needed for the gaseous absorption
1384  // calculation.
1385  //
1386  const Index N_species = vmr_field.nbooks();
1387  //
1388  // Interpolated vmr_list, holds a vmr_list for each point in
1389  // ppath_step.
1390  //
1391  Vector vmr_int(ppath_step.np);
1392 
1393  for (Index i_sp = 0; i_sp < N_species; i_sp ++)
1394  {
1395  out3 << "Interpolate vmr field\n";
1396  interp( vmr_int, itw, vmr_field(i_sp, joker, 0, 0), ppath_step.gp_p );
1397  vmr_list_int(i_sp, joker) = vmr_int;
1398  }
1399  //
1400  // Interpolate pressure
1401  //
1402  itw2p( p_int, p_grid, ppath_step.gp_p, itw);
1403 }
1404 
1405 
1406 
1408 
1447  Tensor6View doit_i_field,
1448  const Index& p_index,
1449  const Index& scat_za_index,
1450  ConstVectorView scat_za_grid,
1451  const ArrayOfIndex& cloudbox_limits,
1452  ConstTensor6View doit_scat_field,
1453  // Calculate scalar gas absorption:
1454  const Agenda& propmat_clearsky_agenda,
1455  ConstTensor4View vmr_field,
1456  // Propagation path calculation:
1457  ConstVectorView p_grid,
1458  ConstTensor3View z_field,
1459  // Calculate thermal emission:
1460  ConstTensor3View t_field,
1461  ConstVectorView f_grid,
1462  const Index& f_index,
1463  //particle opticla properties
1464  ConstTensor5View ext_mat_field,
1465  ConstTensor4View abs_vec_field,
1466  const Verbosity& verbosity)
1467 {
1468  CREATE_OUT3;
1469 
1470  const Index N_species = vmr_field.nbooks();
1471  const Index stokes_dim = doit_i_field.ncols();
1472  const Index atmosphere_dim = 1;
1473  Tensor4 propmat_clearsky;
1474  Tensor3 ext_mat;
1475  Matrix abs_vec;
1476  Vector rtp_vmr(N_species,0.);
1477  Vector sca_vec_av(stokes_dim,0);
1478 
1479  // Radiative transfer from one layer to the next, starting
1480  // at the intersection with the next layer and propagating
1481  // to the considered point.
1482  Vector stokes_vec(stokes_dim);
1483  Index bkgr;
1484  if((p_index == 0) && (scat_za_grid[scat_za_index] > 90))
1485  {
1486  bkgr = 2;
1487  }
1488  else
1489  {
1490  bkgr = 0;
1491  }
1492 
1493  // if 0, there is no background
1494  if (bkgr == 0)
1495  {
1496  if(scat_za_grid[scat_za_index] <= 90.0)
1497  {
1498  stokes_vec = doit_i_field(p_index-cloudbox_limits[0] +1, 0, 0, scat_za_index, 0, joker);
1499  Numeric z_field_above = z_field(p_index +1, 0, 0);
1500  Numeric z_field_0 = z_field(p_index, 0, 0);
1501 
1502  Numeric cos_theta, lstep;
1503  if (scat_za_grid[scat_za_index] == 90.0)
1504  {
1505  //cos_theta = cos((scat_za_grid[scat_za_index] -1)*PI/180.);
1506  // cos_theta = cos(89.999999999999*PI/180.);
1507  cos_theta = 1e-20;
1508 
1509  }
1510  else
1511  {
1512  cos_theta = cos(scat_za_grid[scat_za_index]* PI/180.0);
1513  }
1514  Numeric dz = (z_field_above - z_field_0);
1515 
1516  lstep = abs(dz/cos_theta) ;
1517 
1518  // Average temperature
1519  Numeric rtp_temperature = 0.5 * (t_field(p_index,0,0)
1520  + t_field(p_index + 1,0,0));
1521  //
1522  // Average pressure
1523  Numeric rtp_pressure = 0.5 * (p_grid[p_index]
1524  + p_grid[p_index + 1]);
1525 
1526  // Average vmrs
1527  for (Index j = 0; j < N_species; j++)
1528  rtp_vmr[j] = 0.5 * (vmr_field(j,p_index,0,0) +
1529  vmr_field(j,p_index + 1,0,0));
1530  //
1531  // Calculate scalar gas absorption and add it to abs_vec
1532  // and ext_mat.
1533  //
1534 
1535  const Vector rtp_mag_dummy(3,0);
1536  const Vector ppath_los_dummy;
1537 
1538  propmat_clearsky_agendaExecute(ws, propmat_clearsky,
1539  f_grid[Range(f_index, 1)],
1540  rtp_mag_dummy,ppath_los_dummy,
1541  rtp_pressure,
1542  rtp_temperature,
1543  rtp_vmr,
1544  propmat_clearsky_agenda);
1545 
1546  opt_prop_sum_propmat_clearsky(ext_mat, abs_vec, propmat_clearsky);
1547 
1548  //
1549  // Add average particle extinction to ext_mat.
1550 
1551  for (Index k = 0; k < stokes_dim; k++)
1552  {
1553  for (Index j = 0; j < stokes_dim; j++)
1554  {
1555 
1556  ext_mat(0,k,j) += 0.5 *
1557  (ext_mat_field(p_index - cloudbox_limits[0],0,0,k,j) +
1558  ext_mat_field(p_index - cloudbox_limits[0]+ 1,0,0,k,j));
1559  }
1560  //
1561  // Add average particle absorption to abs_vec.
1562  //
1563  abs_vec(0,k) += 0.5 *
1564  (abs_vec_field(p_index- cloudbox_limits[0],0,0,k) +
1565  abs_vec_field(p_index - cloudbox_limits[0]+ 1,0,0,k));
1566 
1567  //
1568  // Averaging of sca_vec:
1569  //
1570  sca_vec_av[k] += 0.5 *
1571  (doit_scat_field(p_index- cloudbox_limits[0],0,0,scat_za_index,0, k) +
1572  doit_scat_field(p_index - cloudbox_limits[0]+ 1,0,0,scat_za_index,0,k));
1573 
1574  }
1575  // Frequency
1576  Numeric f = f_grid[f_index];
1577  //
1578  // Calculate Planck function
1579  //
1580  Numeric rte_planck_value = planck(f, rtp_temperature);
1581 
1582  // Some messages:
1583  out3 << "-----------------------------------------\n";
1584  out3 << "Input for radiative transfer step \n"
1585  << "calculation inside"
1586  << " the cloudbox:" << "\n";
1587  out3 << "Stokes vector at intersection point: \n"
1588  << stokes_vec
1589  << "\n";
1590  out3 << "lstep: ..." << lstep << "\n";
1591  out3 << "------------------------------------------\n";
1592  out3 << "Averaged coefficients: \n";
1593  out3 << "Planck function: " << rte_planck_value << "\n";
1594  out3 << "Scattering vector: " << sca_vec_av << "\n";
1595  out3 << "Absorption vector: " << abs_vec(0,joker) << "\n";
1596  out3 << "Extinction matrix: " << ext_mat(0,joker,joker) << "\n";
1597 
1598 
1599  assert (!is_singular( ext_mat(0,joker,joker)));
1600 
1601  // Radiative transfer step calculation. The Stokes vector
1602  // is updated until the considered point is reached.
1603  rte_step_doit(stokes_vec, Matrix(stokes_dim,stokes_dim),
1604  ext_mat(0,joker,joker), abs_vec(0,joker),
1605  sca_vec_av, lstep, rte_planck_value);
1606 
1607  // Assign calculated Stokes Vector to doit_i_field.
1608  doit_i_field(p_index - cloudbox_limits[0],
1609  0, 0,
1610  scat_za_index, 0,
1611  joker) = stokes_vec;
1612  }
1613  if(scat_za_grid[scat_za_index] > 90)
1614  {
1615  stokes_vec = doit_i_field(p_index-cloudbox_limits[0] -1, 0, 0, scat_za_index, 0, joker);
1616  Numeric z_field_below = z_field(p_index -1, 0, 0);
1617  Numeric z_field_0 = z_field(p_index, 0, 0);
1618 
1619  Numeric cos_theta, lstep;
1620  if (scat_za_grid[scat_za_index] == 90.0)
1621  {
1622  cos_theta = cos((scat_za_grid[scat_za_index] -1)*PI/180.);
1623  //cos_theta = cos(90.00000001*PI/180.);
1624  //cout<<"cos_theta"<<cos_theta<<endl;
1625  }
1626  else
1627  {
1628  cos_theta = cos(scat_za_grid[scat_za_index]* PI/180.0);
1629  }
1630  Numeric dz = ( z_field_0 - z_field_below);
1631  lstep = abs(dz/cos_theta) ;
1632 
1633  // Average temperature
1634  Numeric rtp_temperature = 0.5 * (t_field(p_index,0,0)
1635  + t_field(p_index - 1,0,0));
1636  //
1637  // Average pressure
1638  Numeric rtp_pressure = 0.5 * (p_grid[p_index]
1639  + p_grid[p_index - 1]);
1640 
1641  //
1642  // Average vmrs
1643  for (Index k = 0; k < N_species; k++)
1644  rtp_vmr[k] = 0.5 * (vmr_field(k,p_index,0,0) +
1645  vmr_field(k,p_index - 1,0,0));
1646  //
1647  // Calculate scalar gas absorption and add it to abs_vec
1648  // and ext_mat.
1649  //
1650 
1651  const Vector rtp_mag_dummy(3,0);
1652  const Vector ppath_los_dummy;
1653 
1654  propmat_clearsky_agendaExecute( ws, propmat_clearsky,
1655  f_grid[Range(f_index, 1)],
1656  rtp_mag_dummy,ppath_los_dummy,
1657  rtp_pressure,
1658  rtp_temperature,
1659  rtp_vmr,
1660  propmat_clearsky_agenda );
1661 
1662  opt_prop_sum_propmat_clearsky(ext_mat, abs_vec, propmat_clearsky);
1663 
1664  //
1665  // Add average particle extinction to ext_mat.
1666  //
1667 
1668  // cout<<"cloudbox_limits[1]"<<cloudbox_limits[1]<<endl;
1669 // cout<<"p_index - cloudbox_limits[0]"<<p_index - cloudbox_limits[0]<<endl;
1670  for (Index k = 0; k < stokes_dim; k++)
1671  {
1672  for (Index j = 0; j < stokes_dim; j++)
1673  {
1674 
1675  ext_mat(0,k,j) += 0.5 *
1676  (ext_mat_field(p_index - cloudbox_limits[0],0,0,k,j) +
1677  ext_mat_field(p_index - cloudbox_limits[0]- 1,0,0,k,j));
1678  }
1679  //
1680  // Add average particle absorption to abs_vec.
1681  //
1682  abs_vec(0,k) += 0.5 *
1683  (abs_vec_field(p_index - cloudbox_limits[0],0,0,k) +
1684  abs_vec_field(p_index - cloudbox_limits[0]- 1,0,0,k));
1685 
1686  //
1687  // Averaging of sca_vec:
1688  //
1689  sca_vec_av[k] += 0.5 *
1690  (doit_scat_field(p_index- cloudbox_limits[0],0,0,scat_za_index,0, k) +
1691  doit_scat_field(p_index - cloudbox_limits[0]- 1,0,0,scat_za_index,0,k));
1692 
1693  }
1694  // Frequency
1695  Numeric f = f_grid[f_index];
1696  //
1697  // Calculate Planck function
1698  //
1699  Numeric rte_planck_value = planck(f, rtp_temperature);
1700 
1701  // Some messages:
1702  out3 << "-----------------------------------------\n";
1703  out3 << "Input for radiative transfer step \n"
1704  << "calculation inside"
1705  << " the cloudbox:" << "\n";
1706  out3 << "Stokes vector at intersection point: \n"
1707  << stokes_vec
1708  << "\n";
1709  out3 << "lstep: ..." << lstep << "\n";
1710  out3 << "------------------------------------------\n";
1711  out3 << "Averaged coefficients: \n";
1712  out3 << "Planck function: " << rte_planck_value << "\n";
1713  out3 << "Scattering vector: " << sca_vec_av << "\n";
1714  out3 << "Absorption vector: " << abs_vec(0,joker) << "\n";
1715  out3 << "Extinction matrix: " << ext_mat(0,joker,joker) << "\n";
1716 
1717 
1718  assert (!is_singular( ext_mat(0,joker,joker)));
1719 
1720  // Radiative transfer step calculation. The Stokes vector
1721  // is updated until the considered point is reached.
1722  rte_step_doit(stokes_vec, Matrix(stokes_dim,stokes_dim),
1723  ext_mat(0,joker,joker), abs_vec(0,joker),
1724  sca_vec_av, lstep, rte_planck_value);
1725 
1726  // Assign calculated Stokes Vector to doit_i_field.
1727  doit_i_field(p_index - cloudbox_limits[0],
1728  0, 0,
1729  scat_za_index, 0,
1730  joker) = stokes_vec;
1731  }
1732 
1733 
1734  }// if loop end - for non_ground background
1735 
1736  // bkgr=2 indicates that the background is the surface
1737  else if (bkgr == 2)
1738  {
1739  //Set rte_pos, rte_gp_p and rte_los to match the last point
1740  //in ppath.
1741  //pos
1742  Vector rte_pos( atmosphere_dim );
1743  Numeric z_field_0 = z_field(0, 0, 0);
1744  rte_pos = z_field_0; //ppath_step.pos(np-1,Range(0,atmosphere_dim));
1745  //los
1746  Vector rte_los(1);
1747  rte_los = scat_za_grid[scat_za_index];//ppath_step.los(np-1,joker);
1748  //gp_p
1749  //GridPos rte_gp_p;
1750  //rte_gp_p.idx = p_index;
1751  //rte_gp_p.fd[0] = 0;
1752  //rte_gp_p.fd[1] = 1;
1753  //gridpos_copy( rte_gp_p, ppath_step.gp_p[np-1] );
1754  // Executes the surface agenda
1755  // FIXME: Convert to new agenda scheme before using
1756  // surface_rtprop_agenda.execute();
1757 
1758  throw runtime_error(
1759  "Surface reflections inside cloud box not yet handled." );
1760  /*
1761  See comment in function above
1762  // Check returned variables
1763  if( surface_emission.nrows() != f_grid.nelem() ||
1764  surface_emission.ncols() != stokes_dim )
1765  throw runtime_error(
1766  "The size of the created *surface_emission* is not correct.");
1767 
1768  Index nlos = surface_los.nrows();
1769 
1770  // Define a local vector doit_i_field_sum which adds the
1771  // products of groudnd_refl_coeffs with the downwelling
1772  // radiation for each elements of surface_los
1773  Vector doit_i_field_sum(stokes_dim,0);
1774  // Loop over the surface_los elements
1775  for( Index ilos=0; ilos < nlos; ilos++ )
1776  {
1777  if( stokes_dim == 1 )
1778  {
1779  doit_i_field_sum[0] += surface_refl_coeffs(ilos,f_index,0,0) *
1780  doit_i_field(cloudbox_limits[0],
1781  0, 0,
1782  (scat_za_grid.nelem() -1 - scat_za_index), 0,
1783  0);
1784  }
1785  else
1786  {
1787  Vector stokes_vec2(stokes_dim);
1788  mult( stokes_vec2,
1789  surface_refl_coeffs(ilos,0,joker,joker),
1790  doit_i_field(cloudbox_limits[0],
1791  0, 0,
1792  (scat_za_grid.nelem() -1 - scat_za_index), 0,
1793  joker));
1794  for( Index is=0; is < stokes_dim; is++ )
1795  {
1796  doit_i_field_sum[is] += stokes_vec2[is];
1797  }
1798 
1799  }
1800  }
1801  // Copy from *doit_i_field_sum* to *doit_i_field*, and add the surface emission
1802  for( Index is=0; is < stokes_dim; is++ )
1803  {
1804  doit_i_field (cloudbox_limits[0],
1805  0, 0,
1806  scat_za_index, 0,
1807  is) = doit_i_field_sum[is] + surface_emission(f_index,is);
1808  }
1809 
1810  //cout<<"scat_za_grid"<<scat_za_grid[scat_za_index]<<endl;
1811  //cout<<"p_index"<<p_index<<endl;
1812  //cout<<"cloudboxlimit[0]"<<cloudbox_limits[0]<<endl;
1813  // now the RT is done to the next point in the path.
1814  //
1815  Vector stokes_vec_local;
1816  stokes_vec_local = doit_i_field (0,
1817  0, 0,
1818  scat_za_index, 0,
1819  joker);
1820 
1821  Numeric z_field_above = z_field(p_index +1, 0, 0);
1822  //Numeric z_field_0 = z_field(p_index, 0, 0);
1823  Numeric cos_theta;
1824  if (scat_za_grid[scat_za_index] == 90.0)
1825  {
1826  //cos_theta = cos((scat_za_grid[scat_za_index] -1)*PI/180.);
1827  cos_theta = cos(90.00000001*PI/180.);
1828  }
1829  else
1830  {
1831  cos_theta = cos(scat_za_grid[scat_za_index]* PI/180.0);
1832  }
1833  Numeric dz = (z_field_above - z_field_0);
1834 
1835  Numeric lstep = abs(dz/cos_theta) ;
1836 
1837  // Average temperature
1838  Numeric rtp_temperature = 0.5 * (t_field(p_index,0,0)
1839  + t_field(p_index + 1,0,0));
1840 
1841  //
1842  // Average pressure
1843  Numeric rtp_pressure = 0.5 * (p_grid[p_index]
1844  + p_grid[p_index + 1]);
1845 
1846  //
1847  const Index N_species = vmr_field.nbooks();
1848  // Average vmrs
1849  for (Index k = 0; k < N_species; k++)
1850  {
1851  rtp_vmr[k] = 0.5 * (vmr_field(k,p_index,0,0) +
1852  vmr_field(k,p_index + 1,0,0));
1853  }
1854  //
1855  // Calculate scalar gas absorption and add it to abs_vec
1856  // and ext_mat.
1857  //
1858 
1859  // FIXME: Convert to new agenda scheme before using
1860  //abs_scalar_gas_agenda.execute(p_index);
1861 
1862  opt_prop_gas_agendaExecute(ext_mat, abs_vec, abs_scalar_gas,
1863  opt_prop_gas_agenda);
1864 
1865  //
1866  // Add average particle extinction to ext_mat.
1867  //
1868  //cout<<"Reached Here ????????????????????????????????????????????????";
1869  for (Index k = 0; k < stokes_dim; k++)
1870  {
1871  for (Index j = 0; j < stokes_dim; j++)
1872  {
1873  ext_mat(0,k,j) += 0.5 *
1874  (ext_mat_field(p_index - cloudbox_limits[0],0,0,k,j) +
1875  ext_mat_field(p_index - cloudbox_limits[0]+ 1,0,0,k,j));
1876  }
1877 
1878 
1879  //
1880  //
1881  // Add average particle absorption to abs_vec.
1882  //
1883  abs_vec(0,k) += 0.5 *
1884  (abs_vec_field(p_index - cloudbox_limits[0],0,0,k) +
1885  abs_vec_field(p_index - cloudbox_limits[0]+1,0,0,k));
1886  //
1887  // Averaging of sca_vec:
1888  //
1889  sca_vec_av[k] = 0.5*( doit_scat_field(p_index- cloudbox_limits[0],
1890  0, 0, scat_za_index, 0, k)
1891  + doit_scat_field(p_index- cloudbox_limits[0]+1,
1892  0, 0, scat_za_index, 0, k)) ;
1893 
1894  }
1895  // Frequency
1896  Numeric f = f_grid[f_index];
1897  //
1898  // Calculate Planck function
1899  //
1900  Numeric rte_planck_value = planck(f, rtp_temperature);
1901 
1902  assert (!is_singular( ext_mat(0,joker,joker)));
1903 
1904  // Radiative transfer step calculation. The Stokes vector
1905  // is updated until the considered point is reached.
1906  rte_step_doit(stokes_vec_local, ext_mat(0,joker,joker),
1907  abs_vec(0,joker),
1908  sca_vec_av, lstep, rte_planck_value);
1909  // Assign calculated Stokes Vector to doit_i_field.
1910  doit_i_field(p_index - cloudbox_limits[0],
1911  0, 0,
1912  scat_za_index, 0,
1913  joker) = stokes_vec_local;
1914  */
1915  }//end else loop over surface
1916 }
1917 
1918 
1919 
1920 
1944 void za_gridOpt(//Output:
1945  Vector& za_grid_opt,
1946  Matrix& doit_i_field_opt,
1947  // Input
1948  ConstVectorView za_grid_fine,
1949  ConstTensor6View doit_i_field,
1950  const Numeric& acc,
1951  const Index& scat_za_interp)
1952 {
1953  Index N_za = za_grid_fine.nelem();
1954 
1955  assert(doit_i_field.npages() == N_za);
1956 
1957  Index N_p = doit_i_field.nvitrines();
1958 
1959  Vector i_approx_interp(N_za);
1960  Vector za_reduced(2);
1961 
1962  ArrayOfIndex idx;
1963  idx.push_back(0);
1964  idx.push_back(N_za-1);
1965  ArrayOfIndex idx_unsorted;
1966 
1967  Numeric max_diff = 100;
1968 
1969  ArrayOfGridPos gp_za(N_za);
1970  Matrix itw(za_grid_fine.nelem(), 2);
1971 
1972  ArrayOfIndex i_sort;
1973  Vector diff_vec(N_za);
1974  Vector max_diff_za(N_p);
1975  ArrayOfIndex ind_za(N_p);
1976  Numeric max_diff_p;
1977  Index ind_p=0;
1978 
1979  while( max_diff > acc )
1980  {
1981  za_reduced.resize(idx.nelem());
1982  doit_i_field_opt.resize(N_p, idx.nelem());
1983  max_diff_za = 0.;
1984  max_diff_p = 0.;
1985 
1986  // Interpolate reduced intensity field on fine za_grid for
1987  // all pressure levels
1988  for( Index i_p = 0; i_p < N_p; i_p++ )
1989  {
1990  for( Index i_za_red = 0; i_za_red < idx.nelem(); i_za_red ++)
1991  {
1992  za_reduced[i_za_red] = za_grid_fine[idx[i_za_red]];
1993  doit_i_field_opt(i_p, i_za_red) = doit_i_field(i_p, 0, 0, idx[i_za_red],
1994  0, 0);
1995  }
1996  // Calculate grid positions
1997  gridpos(gp_za, za_reduced, za_grid_fine);
1998  //linear interpolation
1999  if(scat_za_interp == 0 || idx.nelem() < 3)
2000  {
2001  interpweights(itw, gp_za);
2002  interp(i_approx_interp, itw, doit_i_field_opt(i_p, joker), gp_za);
2003  }
2004  else if(scat_za_interp == 1)
2005  {
2006  for(Index i_za = 0; i_za < N_za; i_za ++)
2007  {
2008  i_approx_interp[i_za] =
2009  interp_poly(za_reduced, doit_i_field_opt(i_p, joker),
2010  za_grid_fine[i_za],
2011  gp_za[i_za]);
2012  }
2013  }
2014  else
2015  // Interpolation method not defined
2016  assert(false);
2017 
2018  // Calculate differences between approximated i-vector and
2019  // exact i_vector for the i_p pressure level
2020  for (Index i_za = 0; i_za < N_za; i_za ++)
2021  {
2022  diff_vec[i_za] = abs( doit_i_field(i_p, 0, 0, i_za, 0 ,0)
2023  - i_approx_interp[i_za]);
2024  if( diff_vec[i_za] > max_diff_za[i_p] )
2025  {
2026  max_diff_za[i_p] = diff_vec[i_za];
2027  ind_za[i_p] = i_za;
2028  }
2029  }
2030  // Take maximum value of max_diff_za
2031  if( max_diff_za[i_p] > max_diff_p )
2032  {
2033  max_diff_p = max_diff_za[i_p];
2034  ind_p = i_p;
2035  }
2036  }
2037 
2038 
2039  //Transform in %
2040  max_diff = max_diff_p/doit_i_field(ind_p, 0, 0, ind_za[ind_p], 0, 0)*100.;
2041 
2042  idx.push_back(ind_za[ind_p]);
2043  idx_unsorted = idx;
2044 
2045  i_sort.resize(idx_unsorted.nelem());
2046  get_sorted_indexes(i_sort, idx_unsorted);
2047 
2048  for (Index i = 0; i<idx_unsorted.nelem(); i++)
2049  idx[i] = idx_unsorted[i_sort[i]];
2050 
2051  za_reduced.resize(idx.nelem());
2052  }
2053 
2054  za_grid_opt.resize(idx.nelem());
2055  doit_i_field_opt.resize(N_p, idx.nelem());
2056  for(Index i = 0; i<idx.nelem(); i++)
2057  {
2058  za_grid_opt[i] = za_grid_fine[idx[i]];
2059  doit_i_field_opt(joker, i) = doit_i_field(joker, 0, 0, idx[i], 0, 0);
2060  }
2061 }
2062 
2063 
2065 /*
2066  See WSM *iyInterpCloudboxField*.
2067 
2068  \param iy Out: As the WSV with same name.
2069  \param scat_i_p In: As the WSV with same name.
2070  \param scat_i_lat In: As the WSV with same name.
2071  \param scat_i_lon In: As the WSV with same name.
2072  \param doit_i_field1D_spectrum In: As the WSV with same name.
2073  \param rte_gp_p In: As the WSV with same name.
2074  \param rte_gp_lat In: As the WSV with same name.
2075  \param rte_gp_lon In: As the WSV with same name.
2076  \param rte_los In: As the WSV with same name.
2077  \param cloudbox_on In: As the WSV with same name.
2078  \param cloudbox_limits In: As the WSV with same name.
2079  \param atmosphere_dim In: As the WSV with same name.
2080  \param stokes_dim In: As the WSV with same name.
2081  \param scat_za_grid In: As the WSV with same name.
2082  \param scat_aa_grid In: As the WSV with same name.
2083  \param f_grid In: As the WSV with same name.
2084  \param p_grid In: As the WSV with same name.
2085  \param interpmeth Interpolation method. Can be "linear" or
2086  "polynomial".
2087  \param rigorous Fail if incoming field is not safely interpolable.
2088  \param maxratio Maximum allowed ratio of two radiances regarded as
2089  interpolable.
2090 
2091  \author Claudia Emde and Patrick Eriksson
2092  \date 2004-09-29
2093 */
2095  const Tensor7& scat_i_p,
2096  const Tensor7& scat_i_lat,
2097  const Tensor7& scat_i_lon,
2098  const Tensor4& doit_i_field1D_spectrum,
2099  const GridPos& rte_gp_p,
2100  const GridPos& rte_gp_lat,
2101  const GridPos& rte_gp_lon,
2102  const Vector& rte_los,
2103  const Index& cloudbox_on,
2104  const ArrayOfIndex& cloudbox_limits,
2105  const Index& atmosphere_dim,
2106  const Index& stokes_dim,
2107  const Vector& scat_za_grid,
2108  const Vector& scat_aa_grid,
2109  const Vector& f_grid,
2110  const String& interpmeth,
2111  const Index& rigorous,
2112  const Numeric& maxratio,
2113  const Verbosity& verbosity)
2114 {
2115  CREATE_OUT3;
2116 
2117  //--- Check input -----------------------------------------------------------
2118  if( !(atmosphere_dim == 1 || atmosphere_dim == 3) )
2119  throw runtime_error( "The atmospheric dimensionality must be 1 or 3.");
2120  if( !cloudbox_on )
2121  throw runtime_error( "The cloud box is not activated and no outgoing "
2122  "field can be returned." );
2123  if ( cloudbox_limits.nelem() != 2*atmosphere_dim )
2124  throw runtime_error(
2125  "*cloudbox_limits* is a vector which contains the upper and lower\n"
2126  "limit of the cloud for all atmospheric dimensions.\n"
2127  "So its length must be 2 x *atmosphere_dim*" );
2128  if( scat_za_grid.nelem() == 0 )
2129  throw runtime_error( "The variable *scat_za_grid* is empty. Are dummy "
2130  "values from *cloudboxOff used?" );
2131  if( !( interpmeth == "linear" || interpmeth == "polynomial" ) )
2132  throw runtime_error( "Unknown interpolation method. Possible choices are "
2133  "\"linear\" and \"polynomial\"." );
2134  if( interpmeth == "polynomial" && atmosphere_dim != 1 )
2135  throw runtime_error( "Polynomial interpolation method is only available "
2136  "for *atmosphere_dim* = 1." );
2137  if( scat_i_p.nlibraries() != f_grid.nelem() )
2138  throw runtime_error( "Inconsistency in size between f_grid and scat_i_p! "
2139  "(This method does not yet handle dispersion type calculations.)" );
2140  //---------------------------------------------------------------------------
2141 
2142 
2143  //--- Determine if at border or inside of cloudbox (or outside!)
2144  //
2145  // Let us introduce a number coding for cloud box borders.
2146  // Borders have the same number as position in *cloudbox_limits*.
2147  // Inside cloud box is coded as 99, and outside as > 100.
2148  Index border = 999;
2149  //
2150  //- Check if at any border
2151  if( is_gridpos_at_index_i( rte_gp_p, cloudbox_limits[0], false ) )
2152  { border = 0; }
2153  else if( is_gridpos_at_index_i( rte_gp_p, cloudbox_limits[1], false ) )
2154  { border = 1; }
2155  if( atmosphere_dim == 3 && border > 100 )
2156  {
2157  if( is_gridpos_at_index_i( rte_gp_lat, cloudbox_limits[2], false ) )
2158  { border = 2; }
2159  else if( is_gridpos_at_index_i( rte_gp_lat, cloudbox_limits[3], false ) )
2160  { border = 3; }
2161  else if( is_gridpos_at_index_i( rte_gp_lon, cloudbox_limits[4], false ) )
2162  { border = 4; }
2163  else if( is_gridpos_at_index_i( rte_gp_lon, cloudbox_limits[5], false ) )
2164  { border = 5; }
2165  }
2166 
2167  //
2168  //- Check if inside (till here border<100 means we are at some border)
2169  if( border > 100 )
2170  {
2171  // Assume inside as it is easiest to detect if outside (can be detected
2172  // check in one dimension at the time)
2173  bool inside = true;
2174  Numeric fgp;
2175 
2176  // Check in pressure dimension
2177  fgp = fractional_gp( rte_gp_p );
2178  if( fgp < Numeric(cloudbox_limits[0]) ||
2179  fgp > Numeric(cloudbox_limits[1]) )
2180  { inside = false; }
2181 
2182  // Check in lat and lon dimensions
2183  if( atmosphere_dim == 3 && inside )
2184  {
2185  fgp = fractional_gp( rte_gp_lat );
2186  if( fgp < Numeric(cloudbox_limits[2]) ||
2187  fgp > Numeric(cloudbox_limits[3]) )
2188  { inside = false; }
2189  fgp = fractional_gp( rte_gp_lon );
2190  if( fgp < Numeric(cloudbox_limits[4]) ||
2191  fgp > Numeric(cloudbox_limits[5]) )
2192  { inside = false; }
2193  }
2194 
2195  if( inside )
2196  { border = 99; }
2197  }
2198 
2199  // If outside, something is wrong
2200  if( border > 100 )
2201  {
2202  throw runtime_error(
2203  "Given position has been found to be outside the cloudbox." );
2204  }
2205 
2206  //- Sizes
2207  const Index nf = f_grid.nelem();
2208  DEBUG_ONLY (const Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1);
2209  const Index nza = scat_za_grid.nelem();
2210 
2211  //- Resize *iy*
2212  iy.resize( nf, stokes_dim );
2213 
2214 
2215  // Sensor points inside the cloudbox
2216  if( border == 99 )
2217  {
2218  if (atmosphere_dim == 3)
2219  {
2220  throw runtime_error(
2221  "3D DOIT calculations are not implemented\n"
2222  "for observations from inside the cloudbox.\n"
2223  );
2224  }
2225  else
2226  {
2227  assert(atmosphere_dim == 1);
2228 
2229  // *doit_i_field1D_spectra* is normally calculated internally:
2230  assert( is_size(doit_i_field1D_spectrum, nf, np, nza, stokes_dim) );
2231 
2232  out3 << " Interpolating outgoing field:\n";
2233  out3 << " zenith_angle: " << rte_los[0] << "\n";
2234  out3 << " Sensor inside cloudbox at position: " <<
2235  rte_gp_p << "\n";
2236 
2237  // Grid position in *scat_za_grid*
2238  GridPos gp_za;
2239  gridpos( gp_za, scat_za_grid, rte_los[0] );
2240 
2241  // Corresponding interpolation weights
2242  Vector itw_za(2);
2243  interpweights( itw_za, gp_za );
2244 
2245  // Grid position in *p_grid* (only cloudbox part because
2246  // doit_i_field1D_spectra is only defined inside the cloudbox
2247  GridPos gp_p;
2248  gp_p = rte_gp_p;
2249  gp_p.idx = rte_gp_p.idx - cloudbox_limits[0];
2250  gridpos_upperend_check( gp_p, cloudbox_limits[1] -
2251  cloudbox_limits[0] );
2252 
2253  // cout << gp_p << endl;
2254 
2255  Vector itw_p(2);
2256  interpweights( itw_p, gp_p );
2257 
2258  Vector iy_p(nza);
2259 
2260  if( interpmeth == "linear" )
2261  {
2262  for(Index is = 0; is < stokes_dim; is++ )
2263  {
2264  for(Index iv = 0; iv < nf; iv++ )
2265  {
2266  for (Index i_za = 0; i_za < nza; i_za++)
2267  {
2268  iy_p[i_za] = interp
2269  (itw_p, doit_i_field1D_spectrum(iv, joker, i_za, is),
2270  gp_p);
2271  }
2272  iy(iv,is) = interp( itw_za, iy_p, gp_za);
2273  }
2274  }
2275  }
2276  else
2277  {
2278  for(Index is = 0; is < stokes_dim; is++ )
2279  {
2280  for(Index iv = 0; iv < nf; iv++ )
2281  {
2282  for (Index i_za = 0; i_za < nza; i_za++)
2283  {
2284  iy_p[i_za] = interp
2285  (itw_p, doit_i_field1D_spectrum(iv, joker, i_za, is),
2286  gp_p);
2287  }
2288  iy(iv,is) = interp_poly( scat_za_grid, iy_p, rte_los[0],
2289  gp_za );
2290  }
2291  }
2292  }
2293 
2294  }
2295 
2296  }
2297 
2298  // Sensor outside the cloudbox
2299 
2300  // --- 1D ------------------------------------------------------------------
2301  else if( atmosphere_dim == 1 )
2302  {
2303  out3 << " Interpolating outgoing field:\n";
2304  out3 << " zenith_angle: " << rte_los[0] << "\n";
2305  if( border )
2306  out3 << " top side\n";
2307  else
2308  out3 << " bottom side\n";
2309 
2310  // Grid position in *scat_za_grid*
2311  GridPos gp;
2312  gridpos( gp, scat_za_grid, rte_los[0] );
2313 
2314  // Corresponding interpolation weights
2315  Vector itw(2);
2316  interpweights( itw, gp );
2317 
2318  if( interpmeth == "linear" )
2319  {
2320  for(Index is = 0; is < stokes_dim; is++ )
2321  {
2322  if( is==0 && rigorous )
2323  {
2324  for(Index iv = 0; iv < nf; iv++ )
2325  {
2326  // cout << scat_i_p(iv,border,0,0,gp.idx,0,is)/
2327  // scat_i_p(iv,border,0,0,gp.idx+1,0,is) << "\n";
2328  if( scat_i_p(iv,border,0,0,gp.idx,0,is)/
2329  scat_i_p(iv,border,0,0,gp.idx+1,0,is) > 1/maxratio &&
2330  scat_i_p(iv,border,0,0,gp.idx,0,is)/
2331  scat_i_p(iv,border,0,0,gp.idx+1,0,is) < maxratio )
2332  {
2333  iy(iv,is) = interp( itw,
2334  scat_i_p( iv, border, 0, 0, joker, 0, is ),
2335  gp );
2336  }
2337  else
2338  {
2339  ostringstream os;
2340  os << "ERROR: Radiance difference between interpolation\n"
2341  << "points is too large (factor " << maxratio << ") to\n"
2342  << "safely interpolate. This might be due to za_grid\n"
2343  << "being too coarse or the radiance field being a\n"
2344  << "step-like function.\n";
2345  os << "Happens at boundary " << border << " between zenith\n"
2346  << "angels " << scat_za_grid[gp.idx] << " and "
2347  << scat_za_grid[gp.idx+1] << "deg for frequency"
2348  << "#" << iv << ", where radiances are "
2349  << scat_i_p(iv,border,0,0,gp.idx,0,0)
2350  << " and " << scat_i_p(iv,border,0,0,gp.idx+1,0,0)
2351  << " W/(sr m2 Hz).";
2352  throw runtime_error(os.str());
2353  }
2354  }
2355  }
2356  else
2357  {
2358  for(Index iv = 0; iv < nf; iv++ )
2359  {
2360  iy(iv,is) = interp( itw,
2361  scat_i_p( iv, border, 0, 0, joker, 0, is ),
2362  gp );
2363  }
2364  }
2365  }
2366  }
2367  else
2368  {
2369  for(Index is = 0; is < stokes_dim; is++ )
2370  {
2371  for(Index iv = 0; iv < nf; iv++ )
2372  {
2373  iy(iv,is) = interp_poly( scat_za_grid,
2374  scat_i_p( iv, border, 0, 0, joker, 0, is ) , rte_los[0],
2375  gp );
2376  }
2377  }
2378  }
2379  }
2380 
2381  // --- 3D ------------------------------------------------------------------
2382  else
2383  {
2384  // Use asserts to check *scat_i_XXX* as these variables should to 99% be
2385  // calculated internally, and thus make it possible to avoid this check.
2386  assert ( is_size( scat_i_p, nf, 2, scat_i_p.nshelves(),
2387  scat_i_p.nbooks(), scat_za_grid.nelem(),
2388  scat_aa_grid.nelem(), stokes_dim ));
2389 
2390  assert ( is_size( scat_i_lat, nf, scat_i_lat.nvitrines(), 2,
2391  scat_i_p.nbooks(), scat_za_grid.nelem(),
2392  scat_aa_grid.nelem(), stokes_dim ));
2393  assert ( is_size( scat_i_lon, nf, scat_i_lat.nvitrines(),
2394  scat_i_p.nshelves(), 2, scat_za_grid.nelem(),
2395  scat_aa_grid.nelem(), stokes_dim ));
2396 
2397  out3 << " Interpolating outgoing field:\n";
2398  out3 << " zenith angle : " << rte_los[0] << "\n";
2399  out3 << " azimuth angle: " << rte_los[1]+180. << "\n";
2400 
2401 
2402  // Scattering angle grid positions
2403  GridPos gp_za, gp_aa;
2404  gridpos( gp_za, scat_za_grid, rte_los[0] );
2405  gridpos( gp_aa, scat_aa_grid, rte_los[1]+180. );
2406 
2407  // Interpolation weights (for 4D "red" interpolation)
2408  Vector itw(16);
2409 
2410  // Outgoing from pressure level
2411  if( border <= 1 )
2412  {
2413  // Lat and lon grid positions with respect to cloud box
2414  GridPos cb_gp_lat, cb_gp_lon;
2415  cb_gp_lat = rte_gp_lat;
2416  cb_gp_lon = rte_gp_lon;
2417  cb_gp_lat.idx -= cloudbox_limits[2];
2418  cb_gp_lon.idx -= cloudbox_limits[4];
2419  //
2420  gridpos_upperend_check( cb_gp_lat, cloudbox_limits[3] -
2421  cloudbox_limits[2] );
2422  gridpos_upperend_check( cb_gp_lon, cloudbox_limits[5] -
2423  cloudbox_limits[4] );
2424 
2425  interpweights( itw, cb_gp_lat, cb_gp_lon, gp_za, gp_aa );
2426 
2427  for(Index is = 0; is < stokes_dim; is++ )
2428  {
2429  for(Index iv = 0; iv < nf; iv++ )
2430  {
2431  iy(iv,is) = interp( itw,
2432  scat_i_p( iv, border, joker, joker, joker, joker, is ),
2433  cb_gp_lat, cb_gp_lon, gp_za, gp_aa );
2434  }
2435  }
2436  }
2437 
2438  // Outgoing from latitude level
2439  else if( border <= 3 )
2440  {
2441  // Pressure and lon grid positions with respect to cloud box
2442  GridPos cb_gp_p, cb_gp_lon;
2443  cb_gp_p = rte_gp_p;
2444  cb_gp_lon = rte_gp_lon;
2445  cb_gp_p.idx -= cloudbox_limits[0];
2446  cb_gp_lon.idx -= cloudbox_limits[4];
2447  //
2448  gridpos_upperend_check( cb_gp_p, cloudbox_limits[1] -
2449  cloudbox_limits[0] );
2450  gridpos_upperend_check( cb_gp_lon, cloudbox_limits[5] -
2451  cloudbox_limits[4] );
2452 
2453  interpweights( itw, cb_gp_p, cb_gp_lon, gp_za, gp_aa );
2454 
2455  for(Index is = 0; is < stokes_dim; is++ )
2456  {
2457  for(Index iv = 0; iv < nf; iv++ )
2458  {
2459  iy(iv,is) = interp( itw,
2460  scat_i_lat( iv, joker, border-2, joker, joker, joker, is ),
2461  cb_gp_p, cb_gp_lon, gp_za, gp_aa );
2462  }
2463  }
2464  }
2465 
2466  // Outgoing from longitude level
2467  else
2468  {
2469  // Pressure and lat grid positions with respect to cloud box
2470  GridPos cb_gp_p, cb_gp_lat;
2471  cb_gp_p = rte_gp_p;
2472  cb_gp_lat = rte_gp_lat;
2473  cb_gp_p.idx -= cloudbox_limits[0];
2474  cb_gp_lat.idx -= cloudbox_limits[2];
2475  //
2476  gridpos_upperend_check( cb_gp_p, cloudbox_limits[1] -
2477  cloudbox_limits[0] );
2478  gridpos_upperend_check( cb_gp_lat, cloudbox_limits[3] -
2479  cloudbox_limits[2] );
2480 
2481  interpweights( itw, cb_gp_p, cb_gp_lat, gp_za, gp_aa );
2482 
2483  for(Index is = 0; is < stokes_dim; is++ )
2484  {
2485  for(Index iv = 0; iv < nf; iv++ )
2486  {
2487  iy(iv,is) = interp( itw,
2488  scat_i_lon( iv, joker, joker, border-4, joker, joker, is ),
2489  cb_gp_p, cb_gp_lat, gp_za, gp_aa );
2490  }
2491  }
2492  }
2493  }
2494 }
2495 
2496 
2498 
2520 void
2522  Tensor6& doit_scat_field,
2523  const Tensor6& doit_i_field,
2524  const ArrayOfIndex& cloudbox_limits,
2525  const Agenda& spt_calc_agenda,
2526  const Index& atmosphere_dim,
2527  const Vector& scat_za_grid,
2528  const Vector& scat_aa_grid,
2529  const Tensor4& pnd_field,
2530  const Agenda& opt_prop_part_agenda,
2531  const Tensor3& t_field,
2532  const Numeric& norm_error_threshold,
2533  const Index& norm_debug,
2534  const Verbosity& verbosity)
2535 {
2536  if (atmosphere_dim != 1)
2537  throw runtime_error("Only 1D is supported here for now");
2538 
2539  CREATE_OUT0;
2540  CREATE_OUT2;
2541 
2542  // Number of zenith angles.
2543  const Index Nza = scat_za_grid.nelem();
2544 
2545  if (scat_za_grid[0] != 0. || scat_za_grid[Nza-1] != 180.)
2546  throw runtime_error("The range of *scat_za_grid* must [0 180].");
2547 
2548  // Number of azimuth angles.
2549  const Index Naa = scat_aa_grid.nelem();
2550 
2551  if (scat_aa_grid[0] != 0. || scat_aa_grid[Naa-1] != 360.)
2552  throw runtime_error("The range of *scat_aa_grid* must [0 360].");
2553 
2554  // Get stokes dimension from *doit_scat_field*:
2555  const Index stokes_dim = doit_scat_field.ncols();
2556  assert(stokes_dim > 0 || stokes_dim < 5);
2557 
2558  // To use special interpolation functions for atmospheric fields we
2559  // use ext_mat_field and abs_vec_field:
2560  Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
2561  stokes_dim, stokes_dim, 0.);
2562  Tensor4 abs_vec_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
2563  stokes_dim, 0.);
2564 
2565  const Index Np = doit_scat_field.nvitrines();
2566 
2567  Tensor5 doit_scat_ext_field(doit_scat_field.nvitrines(),
2568  doit_scat_field.nshelves(),
2569  doit_scat_field.nbooks(),
2570  doit_scat_field.npages(),
2571  doit_scat_field.nrows(),
2572  0.);
2573 
2574  Index scat_aa_index_local = 0;
2575 
2576  // Calculate scattering extinction field
2577  for(Index scat_za_index_local = 0; scat_za_index_local < Nza;
2578  scat_za_index_local ++)
2579  {
2580  // This function has to be called inside the angular loop, as
2581  // spt_calc_agenda takes *scat_za_index* and *scat_aa_index*
2582  // from the workspace.
2583  // *scat_p_index* is needed for communication with agenda
2584  // *opt_prop_part_agenda*.
2585  cloud_fieldsCalc(ws, ext_mat_field, abs_vec_field,
2586  spt_calc_agenda, opt_prop_part_agenda,
2587  scat_za_index_local, scat_aa_index_local,
2588  cloudbox_limits, t_field, pnd_field, verbosity);
2589 
2590  for(Index p_index = 0;
2591  p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
2592  p_index ++)
2593  {
2594  // For all in p_grid (in cloudbox):
2595  // I_ext = (ext_mat_field - abs_vec_field) * doit_i_field
2596  // equivalent to:
2597  // I_ext = I * (K11-a1) + Q * (K12 - a2) + U * (K13 - a3) + V * (K14 - a4)
2598  for (Index i = 0; i < stokes_dim; i++)
2599  {
2600  doit_scat_ext_field(p_index, 0, 0, scat_za_index_local, 0)
2601  += doit_i_field(p_index, 0, 0, scat_za_index_local, 0, i)
2602  * (ext_mat_field(p_index, 0, 0, 0, i) - abs_vec_field(p_index, 0, 0, i));
2603  }
2604  }
2605  }
2606 
2607  Numeric corr_max = .0;
2608  Index corr_max_p_index = -1;
2609 
2610  for (Index p_index = 0; p_index < Np; p_index++)
2611  {
2612  // Calculate scattering integrals
2613  const Numeric scat_int
2614  = AngIntegrate_trapezoid(doit_scat_field(p_index, 0, 0, joker, 0, 0),
2615  scat_za_grid);
2616 
2617  const Numeric scat_ext_int
2618  = AngIntegrate_trapezoid(doit_scat_ext_field(p_index, 0, 0, joker, 0),
2619  scat_za_grid);
2620 
2621  // Calculate factor between scattered extinction field integral
2622  // and scattered field integral
2623  const Numeric corr_factor = scat_ext_int / scat_int;
2624 
2625  // If no scattering is present, the correction factor can become
2626  // inf or nan. We just don't apply it for those cases.
2627  if (!isnan(corr_factor) && !isinf(corr_factor))
2628  {
2629  if (abs(corr_factor) > abs(corr_max))
2630  {
2631  corr_max = corr_factor;
2632  corr_max_p_index = p_index;
2633  }
2634  if (abs(1.-corr_factor) > norm_error_threshold)
2635  {
2636  ostringstream os;
2637  os << "ERROR: DOIT correction factor exceeds threshold: "
2638  << setprecision(4) << 1.-corr_factor << " at p_index " << p_index << "\n";
2639  throw runtime_error(os.str());
2640  }
2641  else if (abs(1.-corr_factor) > norm_error_threshold/2.)
2642  {
2643  out0 << " WARNING: DOIT correction factor above threshold/2: "
2644  << 1.-corr_factor << " at p_index " << p_index << "\n";
2645  }
2646 
2647  // Scale scattered field with correction factor
2648  doit_scat_field(p_index, 0, 0, joker, 0, joker) *= corr_factor;
2649  }
2650  }
2651 
2652 
2653  ArtsOut& norm_out = out2;
2654  if (norm_debug) norm_out = out0;
2655 
2656  if (corr_max_p_index != -1)
2657  {
2658  norm_out << " Max. DOIT correction factor in this iteration: " << 1.-corr_max
2659  << " at p_index " << corr_max_p_index << "\n";
2660  }
2661  else
2662  {
2663  norm_out << " No DOIT correction performed in this iteration.\n";
2664  }
2665 }
2666 
2667 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void opt_prop_part_agendaExecute(Workspace &ws, Tensor3 &ext_mat, Matrix &abs_vec, const Tensor3 &ext_mat_spt, const Matrix &abs_vec_spt, const Index scat_p_index, const Index scat_lat_index, const Index scat_lon_index, const Agenda &input_agenda)
Definition: auto_md.cc:14459
ArrayOfGridPos gp_lat
Definition: ppath.h:77
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
void swap(Vector &v1, Vector &v2)
Swaps two objects.
Definition: matpackI.cc:813
The VectorView class.
Definition: matpackI.h:372
Index nrows() const
Returns the number of rows.
Definition: matpackV.cc:53
The Tensor4View class.
Definition: matpackIV.h:243
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:142
Declarations having to do with the four output streams.
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
void doit_scat_fieldNormalize(Workspace &ws, Tensor6 &doit_scat_field, const Tensor6 &doit_i_field, const ArrayOfIndex &cloudbox_limits, const Agenda &spt_calc_agenda, const Index &atmosphere_dim, const Vector &scat_za_grid, const Vector &scat_aa_grid, const Tensor4 &pnd_field, const Agenda &opt_prop_part_agenda, const Tensor3 &t_field, const Numeric &norm_error_threshold, const Index &norm_debug, const Verbosity &verbosity)
Normalization of scattered field.
Definition: doit.cc:2521
Matrix los
Definition: ppath.h:68
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVI.cc:32
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Index dim
Definition: ppath.h:60
A constant view of a Tensor6.
Definition: matpackVI.h:159
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lraytrace, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:14720
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
Vector lstep
Definition: ppath.h:70
Linear algebra functions.
void cloud_RT_no_background(Workspace &ws, Tensor6View doit_i_field, 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 doit_i_field_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 &scat_za_index, const Index &scat_aa_index, const Verbosity &verbosity)
cloud_RT_no_background
Definition: doit.cc:1016
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Matrix pos
Definition: ppath.h:67
The Tensor6View class.
Definition: matpackVI.h:449
This file contains basic functions to handle XML data files.
The Tensor7 class.
Definition: matpackVII.h:1931
void iy_interp_cloudbox_field(Matrix &iy, const Tensor7 &scat_i_p, const Tensor7 &scat_i_lat, const Tensor7 &scat_i_lon, const Tensor4 &doit_i_field1D_spectrum, const GridPos &rte_gp_p, const GridPos &rte_gp_lat, const GridPos &rte_gp_lon, const Vector &rte_los, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &scat_za_grid, const Vector &scat_aa_grid, const Vector &f_grid, const String &interpmeth, const Index &rigorous, const Numeric &maxratio, const Verbosity &verbosity)
Interpolation of cloud box intensity field.
Definition: doit.cc:2094
Vector r
Definition: ppath.h:69
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Contains sorting routines.
Structure to store a grid position.
Definition: interpolation.h:74
A constant view of a Tensor4.
Definition: matpackIV.h:141
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
This file contains the definition of Array.
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:768
void cloud_RT_surface(Workspace &ws, Tensor6View doit_i_field, 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 scat_za_grid, const Index &scat_za_index)
cloud_RT_surface
Definition: doit.cc:1175
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
#define DEBUG_ONLY(...)
Definition: debug.h:37
void cloud_ppath_update3D(Workspace &ws, Tensor6View doit_i_field, const Index &p_index, const Index &lat_index, const Index &lon_index, const Index &scat_za_index, const Index &scat_aa_index, ConstVectorView scat_za_grid, ConstVectorView scat_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_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:738
Declarations for agendas.
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
ppath_init_structure
Definition: ppath.cc:1726
The Tensor3View class.
Definition: matpackIII.h:232
Index nrows() const
Returns the number of rows.
Definition: matpackVI.cc:56
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void cloud_ppath_update1D_noseq(Workspace &ws, Tensor6View doit_i_field, const Index &p_index, const Index &scat_za_index, ConstVectorView scat_za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_i_field_old, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, 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)
cloud_ppath_update1D_noseq
Definition: doit.cc:556
bool is_inside_cloudbox(const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, const bool include_boundaries)
Definition: cloudbox.cc:892
void propmat_clearsky_agendaExecute(Workspace &ws, Tensor4 &propmat_clearsky, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:13039
void ludcmp(MatrixView LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:54
A constant view of a Tensor5.
Definition: matpackV.h:152
#define abs(x)
Definition: continua.cc:20458
const Joker joker
Index idx
Definition: interpolation.h:75
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:32
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1199
Header file for special_interp.cc.
Propagation path structure and functions.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:79
The Tensor5View class.
Definition: matpackV.h:276
Header file for logic.cc.
Radiative transfer in cloudbox.
This can be used to make arrays out of anything.
Definition: array.h:40
Index nshelves() const
Returns the number of shelves.
Definition: matpackVI.cc:38
bool is_singular(ConstMatrixView A)
Checks if a square matrix is singular.
Definition: logic.cc:377
Index ncols() const
Returns the number of columns.
Definition: matpackVI.cc:62
void cloud_fieldsCalc(Workspace &ws, Tensor5View ext_mat_field, Tensor4View abs_vec_field, const Agenda &spt_calc_agenda, const Agenda &opt_prop_part_agenda, const Index &scat_za_index, const Index &scat_aa_index, const ArrayOfIndex &cloudbox_limits, ConstTensor3View t_field, ConstTensor4View pnd_field, const Verbosity &verbosity)
cloud_fieldsCalc
Definition: doit.cc:257
void interp_cloud_coeff1D(Tensor3View ext_mat_int, MatrixView abs_vec_int, MatrixView sca_vec_int, MatrixView doit_i_field_int, VectorView t_int, MatrixView vmr_list_int, VectorView p_int, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, ConstTensor6View doit_scat_field, ConstTensor6View doit_i_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView p_grid, const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, ConstVectorView scat_za_grid, const Index &scat_za_interp, const Verbosity &verbosity)
interp_cloud_coeff1D
Definition: doit.cc:1253
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
void cloud_ppath_update1D(Workspace &ws, Tensor6View doit_i_field, const Index &p_index, const Index &scat_za_index, ConstVectorView scat_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_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)
cloud_ppath_update1D
Definition: doit.cc:416
A constant view of a Tensor3.
Definition: matpackIII.h:139
The Tensor6 class.
Definition: matpackVI.h:950
void cloud_ppath_update1D_planeparallel(Workspace &ws, Tensor6View doit_i_field, const Index &p_index, const Index &scat_za_index, ConstVectorView scat_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:1446
A constant view of a Vector.
Definition: matpackI.h:292
Index np
Definition: ppath.h:61
ArrayOfGridPos gp_lon
Definition: ppath.h:78
Index npages() const
Returns the number of pages.
Definition: matpackVI.cc:50
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:44
const Numeric PI
A constant view of a Matrix.
Definition: matpackI.h:596
Numeric planck(const Numeric &f, const Numeric &t)
planck
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
const Numeric RAD2DEG
void spt_calc_agendaExecute(Workspace &ws, Tensor3 &ext_mat_spt, Matrix &abs_vec_spt, const Index scat_p_index, const Index scat_lat_index, const Index scat_lon_index, const Numeric rtp_temperature, const Index scat_za_index, const Index scat_aa_index, const Agenda &input_agenda)
Definition: auto_md.cc:14960
Index ncols() const
Returns the number of columns.
Definition: matpackV.cc:59
#define CREATE_OUT0
Definition: messages.h:211
#define CREATE_OUT3
Definition: messages.h:214
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:91
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVII.cc:38
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
Internal cloudbox functions.
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void za_gridOpt(Vector &za_grid_opt, Matrix &doit_i_field_opt, ConstVectorView za_grid_fine, ConstTensor6View doit_i_field, const Numeric &acc, const Index &scat_za_interp)
Definition: doit.cc:1944
#define CREATE_OUT2
Definition: messages.h:213
void surface_rtprop_agendaExecute(Workspace &ws, 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:15046
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:50
ArrayOfGridPos gp_p
Definition: ppath.h:76
void rte_step_doit(VectorView stokes_vec, MatrixView trans_mat, ConstMatrixView ext_mat_av, ConstVectorView abs_vec_av, ConstVectorView sca_vec_av, const Numeric &lstep, const Numeric &rtp_planck_value, const bool &trans_is_precalc)
rte_step_doit
Definition: doit.cc:108
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Definition: rte.cc:903
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:327
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
itw2p
The Tensor5 class.
Definition: matpackV.h:451
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
void opt_prop_sum_propmat_clearsky(Tensor3 &ext_mat, Matrix &abs_vec, const Tensor4 propmat_clearsky)
Get optical properties from propmat_clearsky.
Index nbooks() const
Returns the number of books.
Definition: matpackVI.cc:44