ARTS  2.0.49
m_rte.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2010
2  Patrick Eriksson <patrick.eriksson@chalmers.se>
3  Stefan Buehler <sbuehler@ltu.se>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
20 
21 
22 /*===========================================================================
23  === File description
24  ===========================================================================*/
25 
39 /*===========================================================================
40  === External declarations
41  ===========================================================================*/
42 
43 #include <cmath>
44 #include <stdexcept>
45 #include "arts.h"
46 #include "arts_omp.h"
47 #include "auto_md.h"
48 #include "check_input.h"
49 #include "jacobian.h"
50 #include "logic.h"
51 #include "math_funcs.h"
52 #include "messages.h"
53 #include "montecarlo.h"
54 #include "physics_funcs.h"
55 #include "ppath.h"
56 #include "rte.h"
57 #include "special_interp.h"
58 
59 extern const String ABSSPECIES_MAINTAG;
60 extern const String TEMPERATURE_MAINTAG;
61 
62 
63 // A macro to loop analytical jacobian quantities
64 #define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do) \
65  for( Index iq=0; iq<jacobian_quantities.nelem(); iq++ ) \
66  { \
67  if( jacobian_quantities[iq].Analytical() ) \
68  { what_to_do } \
69  }
70 
71 
72 
73 
74 /*===========================================================================
75  === Special sub-functions for analytical jacobians (in alphabetical order)
76  === (Main sub-functions get_iy_of_background and iyb_calc placed separately)
77  ===========================================================================*/
78 
79 
81 
95 // A small help function, to make the code below cleaner
96 void from_dpath_to_dx(
97  MatrixView diy_dx,
98  ConstMatrixView diy_dq,
99  const Numeric& w )
100 {
101  for( Index irow=0; irow<diy_dx.nrows(); irow++ )
102  {
103  for( Index icol=0; icol<diy_dx.ncols(); icol++ )
104  { diy_dx(irow,icol) += w * diy_dq(irow,icol); }
105  }
106 }
107 // Kept as a local function as long as just used here.
108 // We trust here gridpos, and "extrapolation points" are identified simply
109 // by a fd outside [0,1].
110 void add_extrap( ArrayOfGridPos& gp )
111 {
112  for( Index i=0; i<gp.nelem(); i++ )
113  {
114  if( gp[i].fd[0] < 0 )
115  {
116  gp[i].fd[0] = 0;
117  gp[i].fd[1] = 1;
118  }
119  else if( gp[i].fd[0] > 1 )
120  {
121  gp[i].fd[0] = 1;
122  gp[i].fd[1] = 0;
123  }
124  }
125 }
126 //
127 void diy_from_path_to_rgrids(
128  Tensor3View diy_dx,
129  const RetrievalQuantity& jacobian_quantity,
130  ConstTensor3View diy_dpath,
131  const Index& atmosphere_dim,
132  const Ppath& ppath,
133  ConstVectorView ppath_p )
134 {
135  // We want here an extrapolation to infinity ->
136  // extremly high extrapolation factor
137  const Numeric extpolfac = 1.0e99;
138 
139  // Retrieval grid of interest
140  Vector r_grid;
141 
142  if( ppath.np > 1 ) // Otherwise nothing to do here
143  {
144  // Pressure
145  r_grid = jacobian_quantity.Grids()[0];
146  Index nr1 = r_grid.nelem();
147  ArrayOfGridPos gp_p(ppath.np);
148  p2gridpos( gp_p, r_grid, ppath_p, extpolfac );
149  add_extrap( gp_p );
150 
151  // Latitude
152  Index nr2 = 1;
153  ArrayOfGridPos gp_lat;
154  if( atmosphere_dim > 1 )
155  {
156  gp_lat.resize(ppath.np);
157  r_grid = jacobian_quantity.Grids()[1];
158  nr2 = r_grid.nelem();
159  gridpos( gp_lat, r_grid, ppath.pos(joker,1), extpolfac );
160  add_extrap( gp_lat );
161  }
162 
163  // Longitude
164  ArrayOfGridPos gp_lon;
165  if( atmosphere_dim > 2 )
166  {
167  gp_lon.resize(ppath.np);
168  r_grid = jacobian_quantity.Grids()[2];
169  gridpos( gp_lon, r_grid, ppath.pos(joker,2), extpolfac );
170  add_extrap( gp_lon );
171  }
172 
173  //- 1D
174  if( atmosphere_dim == 1 )
175  {
176  for( Index ip=0; ip<ppath.np; ip++ )
177  {
178  if( gp_p[ip].fd[0] < 1 )
179  {
180  from_dpath_to_dx( diy_dx(gp_p[ip].idx,joker,joker),
181  diy_dpath(ip,joker,joker), gp_p[ip].fd[1] );
182  }
183  if( gp_p[ip].fd[0] > 0 )
184  {
185  from_dpath_to_dx( diy_dx(gp_p[ip].idx+1,joker,joker),
186  diy_dpath(ip,joker,joker), gp_p[ip].fd[0] );
187  }
188  }
189  }
190 
191  //- 2D
192  else if( atmosphere_dim == 2 )
193  {
194  for( Index ip=0; ip<ppath.np; ip++ )
195  {
196  Index ix = nr1*gp_lat[ip].idx + gp_p[ip].idx;
197  // Low lat, low p
198  from_dpath_to_dx( diy_dx(ix,joker,joker),
199  diy_dpath(ip,joker,joker),
200  gp_lat[ip].fd[1]*gp_p[ip].fd[1] );
201  // Low lat, high p
202  from_dpath_to_dx( diy_dx(ix+1,joker,joker),
203  diy_dpath(ip,joker,joker),
204  gp_lat[ip].fd[1]*gp_p[ip].fd[0] );
205  // High lat, low p
206  from_dpath_to_dx( diy_dx(ix+nr1,joker,joker),
207  diy_dpath(ip,joker,joker),
208  gp_lat[ip].fd[0]*gp_p[ip].fd[1] );
209  // High lat, high p
210  from_dpath_to_dx( diy_dx(ix+nr1+1,joker,joker),
211  diy_dpath(ip,joker,joker),
212  gp_lat[ip].fd[0]*gp_p[ip].fd[0] );
213  }
214  }
215 
216  //- 3D
217  else if( atmosphere_dim == 3 )
218  {
219  for( Index ip=0; ip<ppath.np; ip++ )
220  {
221  Index ix = nr2*nr1*gp_lon[ip].idx +
222  nr1*gp_lat[ip].idx + gp_p[ip].idx;
223  // Low lon, low lat, low p
224  from_dpath_to_dx( diy_dx(ix,joker,joker),
225  diy_dpath(ip,joker,joker),
226  gp_lon[ip].fd[1]*gp_lat[ip].fd[1]*gp_p[ip].fd[1]);
227  // Low lon, low lat, high p
228  from_dpath_to_dx( diy_dx(ix+1,joker,joker),
229  diy_dpath(ip,joker,joker),
230  gp_lon[ip].fd[1]*gp_lat[ip].fd[1]*gp_p[ip].fd[0]);
231  // Low lon, high lat, low p
232  from_dpath_to_dx( diy_dx(ix+nr1,joker,joker),
233  diy_dpath(ip,joker,joker),
234  gp_lon[ip].fd[1]*gp_lat[ip].fd[0]*gp_p[ip].fd[1]);
235  // Low lon, high lat, high p
236  from_dpath_to_dx( diy_dx(ix+nr1+1,joker,joker),
237  diy_dpath(ip,joker,joker),
238  gp_lon[ip].fd[1]*gp_lat[ip].fd[0]*gp_p[ip].fd[0]);
239 
240  // Increase *ix* (to be valid for high lon level)
241  ix += nr2*nr1;
242 
243  // High lon, low lat, low p
244  from_dpath_to_dx( diy_dx(ix,joker,joker),
245  diy_dpath(ip,joker,joker),
246  gp_lon[ip].fd[0]*gp_lat[ip].fd[1]*gp_p[ip].fd[1]);
247  // High lon, low lat, high p
248  from_dpath_to_dx( diy_dx(ix+1,joker,joker),
249  diy_dpath(ip,joker,joker),
250  gp_lon[ip].fd[0]*gp_lat[ip].fd[1]*gp_p[ip].fd[0]);
251  // High lon, high lat, low p
252  from_dpath_to_dx( diy_dx(ix+nr1,joker,joker),
253  diy_dpath(ip,joker,joker),
254  gp_lon[ip].fd[0]*gp_lat[ip].fd[0]*gp_p[ip].fd[1]);
255  // High lon, high lat, high p
256  from_dpath_to_dx( diy_dx(ix+nr1+1,joker,joker),
257  diy_dpath(ip,joker,joker),
258  gp_lon[ip].fd[0]*gp_lat[ip].fd[0]*gp_p[ip].fd[0]);
259  }
260  }
261  }
262 }
263 
264 
265 
267 
286 void get_pointers_for_analytical_jacobians(
287  ArrayOfIndex& abs_species_i,
288  ArrayOfIndex& is_t,
289  const ArrayOfRetrievalQuantity& jacobian_quantities,
290  const ArrayOfArrayOfSpeciesTag& abs_species )
291 {
292 
293  FOR_ANALYTICAL_JACOBIANS_DO(
294  //
295  if( jacobian_quantities[iq].MainTag() == TEMPERATURE_MAINTAG )
296  { is_t[iq] = 1; }
297  else
298  { is_t[iq] = 0; }
299  //
300  if( jacobian_quantities[iq].MainTag() == ABSSPECIES_MAINTAG )
301  {
302  ArrayOfSpeciesTag atag;
303  array_species_tag_from_string( atag, jacobian_quantities[iq].Subtag() );
304  abs_species_i[iq] = chk_contains( "abs_species", abs_species, atag );
305  }
306  else
307  { abs_species_i[iq] = -1; }
308  )
309 }
310 
311 
312 
313 
314 
315 /*===========================================================================
316  === Main help functions (get_iy_of_background and iyb_calc)
317  ===========================================================================*/
318 
320 
360 void get_iy_of_background(
361  Workspace& ws,
362  Matrix& iy,
363  Matrix& iy_error,
364  Index& iy_error_type,
365  Matrix& iy_aux,
366  ArrayOfTensor3& diy_dx,
367  ConstTensor3View iy_transmission,
368  const Index& jacobian_do,
369  const Ppath& ppath,
370  const Index& atmosphere_dim,
371  ConstVectorView p_grid,
372  ConstVectorView lat_grid,
373  ConstVectorView lon_grid,
374  ConstTensor3View t_field,
375  ConstTensor3View z_field,
376  ConstTensor4View vmr_field,
377  const Matrix& r_geoid,
378  const Matrix& z_surface,
379  const Index& cloudbox_on,
380  const Index& stokes_dim,
381  ConstVectorView f_grid,
382  const Agenda& iy_clearsky_agenda,
383  const Agenda& iy_space_agenda,
384  const Agenda& surface_prop_agenda,
385  const Agenda& iy_cloudbox_agenda,
386  const Verbosity& verbosity)
387 {
388  CREATE_OUT3
389 
390  // Some sizes
391  const Index nf = f_grid.nelem();
392  const Index np = ppath.np;
393 
394  // Set rte_pos, rte_gp_XXX and rte_los to match the last point in ppath.
395  // The agendas below use different combinations of these variables.
396  //
397  // Note that the Ppath positions (ppath.pos) for 1D have one column more
398  // than expected by most functions. Only the first atmosphere_dim values
399  // shall be copied.
400  //
401  Vector rte_pos;
402  Vector rte_los;
403  GridPos rte_gp_p;
404  GridPos rte_gp_lat;
405  GridPos rte_gp_lon;
406  rte_pos.resize( atmosphere_dim );
407  rte_pos = ppath.pos(np-1,Range(0,atmosphere_dim));
408  rte_los.resize( ppath.los.ncols() );
409  rte_los = ppath.los(np-1,joker);
410  gridpos_copy( rte_gp_p, ppath.gp_p[np-1] );
411  if( atmosphere_dim > 1 )
412  { gridpos_copy( rte_gp_lat, ppath.gp_lat[np-1] ); }
413  if( atmosphere_dim > 2 )
414  { gridpos_copy( rte_gp_lon, ppath.gp_lon[np-1] ); }
415 
416  out3 << "Radiative background: " << ppath.background << "\n";
417 
418 
419  // Handle the different background cases
420  //
421  switch( ppath_what_background( ppath ) )
422  {
423 
424  case 1: //--- Space ----------------------------------------------------
425  {
426  iy_space_agendaExecute( ws, iy, rte_pos, rte_los, iy_space_agenda );
427 
428  if( iy.ncols() != stokes_dim || iy.nrows() != nf )
429  {
430  ostringstream os;
431  os << "expected size = [" << stokes_dim << "," << nf << "]\n"
432  << "size of iy = [" << iy.nrows() << "," << iy.ncols()<< "]\n"
433  << "The size of *iy* returned from *iy_space_agenda* is "
434  << "not correct.";
435  throw runtime_error( os.str() );
436  }
437  }
438  break;
439 
440  case 2: //--- The surface -----------------------------------------------
441  {
442  // Call *surface_prop_agenda*
443  //
444  Matrix surface_los;
445  Tensor4 surface_rmatrix;
446  Matrix surface_emission;
447  //
448  surface_prop_agendaExecute( ws, surface_emission, surface_los,
449  surface_rmatrix, rte_pos, rte_los,
450  rte_gp_p, rte_gp_lat, rte_gp_lon,
451  surface_prop_agenda );
452 
453  // Check output of *surface_prop_agenda*
454  //
455  const Index nlos = surface_los.nrows();
456  //
457  if( nlos ) // if 0, blackbody ground and some checks are not needed
458  {
459  if( surface_los.ncols() != rte_los.nelem() )
460  throw runtime_error(
461  "Number of columns in *surface_los* is not correct." );
462  if( nlos != surface_rmatrix.nbooks() )
463  throw runtime_error(
464  "Mismatch in size of *surface_los* and *surface_rmatrix*." );
465  if( surface_rmatrix.npages() != nf )
466  throw runtime_error(
467  "Mismatch in size of *surface_rmatrix* and *f_grid*." );
468  if( surface_rmatrix.nrows() != stokes_dim ||
469  surface_rmatrix.ncols() != stokes_dim )
470  throw runtime_error(
471  "Mismatch between size of *surface_rmatrix* and *stokes_dim*." );
472  }
473  if( surface_emission.ncols() != stokes_dim )
474  throw runtime_error(
475  "Mismatch between size of *surface_emission* and *stokes_dim*." );
476  if( surface_emission.nrows() != nf )
477  throw runtime_error(
478  "Mismatch in size of *surface_emission* and f_grid*." );
479  //---------------------------------------------------------------------
480 
481  // Variable to hold down-welling radiation
482  Tensor3 I( nlos, nf, stokes_dim );
483 
484  // Loop *surface_los*-es. If no such LOS, we are ready.
485  if( nlos > 0 )
486  {
487  for( Index ilos=0; ilos<nlos; ilos++ )
488  {
489  Vector los = surface_los(ilos,joker);
490 
491  // Include surface reflection matrix in *iy_transmission*
492  // If iy_transmission is empty, this is interpreted as the
493  // variable is not needed.
494  //
495  Tensor3 iy_trans_new;
496  //
497  if( iy_transmission.npages() )
498  {
499  iy_transmission_mult( iy_trans_new, iy_transmission,
500  surface_rmatrix(ilos,joker,joker,joker) );
501  }
502 
503  // Calculate angular tilt of the surface
504  // (sign(los[0]) to handle negative za for 2D)
505  const Numeric atilt = surfacetilt( atmosphere_dim, lat_grid,
506  lon_grid, r_geoid, z_surface,
507  rte_gp_lat, rte_gp_lat, los );
508  const Numeric zamax = 90 - sign(los[0])*atilt;
509 
510  // I considered to add a check that surface_los is above the
511  // horizon, but that would force e.g. surface_specular_los to
512  // actually calculate the surface tilt, which causes
513  // unnecessary calculation overhead. The angles are now moved
514  // to be just above the horizon, which should be acceptable.
515 
516  // Make sure that we have some margin to the "horizon"
517  // (otherwise numerical problems can create an infinite loop)
518  if( atmosphere_dim == 2 && los[0]<0 ) //2D with za<0
519  { los[0] = max( -zamax+0.1, los[0] ); }
520  else
521  { los[0] = min( zamax-0.1, los[0] ); }
522 
523  // Calculate downwelling radiation for LOS ilos
524  //
525  // The variable iy_clearsky_agenda can in fact be
526  // iy_clearsky_BASIC_agenda
527  //
528  if( iy_clearsky_agenda.name() == "iy_clearsky_basic_agenda" )
529  {
530  iy_clearsky_basic_agendaExecute( ws, iy, rte_pos, los,
531  cloudbox_on, iy_clearsky_agenda);
532  }
533  else
534  {
535  iy_clearsky_agendaExecute( ws, iy, iy_error, iy_error_type,
536  iy_aux, diy_dx, 0, iy_trans_new,
537  rte_pos, los, cloudbox_on, jacobian_do,
538  p_grid, lat_grid, lon_grid, t_field,
539  z_field, vmr_field, iy_clearsky_agenda );
540  }
541 
542  I(ilos,joker,joker) = iy;
543  }
544  }
545 
546  // Add up
547  surface_calc( iy, I, surface_los, surface_rmatrix, surface_emission );
548  }
549  break;
550 
551 
552  case 3: //--- Cloudbox boundary or interior ------------------------------
553  case 4:
554  {
555  iy_cloudbox_agendaExecute( ws, iy, iy_error, iy_error_type,
556  iy_aux, diy_dx, iy_transmission,
557  rte_pos, rte_los, rte_gp_p, rte_gp_lat,
558  rte_gp_lon, iy_cloudbox_agenda );
559 
560  if( iy.nrows() != nf || iy.ncols() != stokes_dim )
561  {
562  CREATE_OUT1
563  out1 << "expected size = [" << nf << "," << stokes_dim << "]\n";
564  out1 << "iy size = [" << iy.nrows() << "," << iy.ncols()<< "]\n";
565  throw runtime_error( "The size of *iy* returned from "
566  "*iy_cloudbox_agenda* is not correct." );
567  }
568  }
569  break;
570 
571  default: //--- ????? ----------------------------------------------------
572  // Are we here, the coding is wrong somewhere
573  assert( false );
574  }
575 }
576 
577 
578 
580 
588 void iyb_calc(
589  Workspace& ws,
590  Vector& iyb,
591  Vector& iyb_error,
592  Index& iy_error_type,
593  Vector& iyb_aux,
594  ArrayOfMatrix& diyb_dx,
595  const Index& imblock,
596  const Index& atmosphere_dim,
597  ConstVectorView p_grid,
598  ConstVectorView lat_grid,
599  ConstVectorView lon_grid,
600  ConstTensor3View t_field,
601  ConstTensor3View z_field,
602  ConstTensor4View vmr_field,
603  const Index& cloudbox_on,
604  const Index& stokes_dim,
605  ConstVectorView f_grid,
606  ConstMatrixView sensor_pos,
607  ConstMatrixView sensor_los,
608  ConstVectorView mblock_za_grid,
609  ConstVectorView mblock_aa_grid,
610  const Index& antenna_dim,
611  const Agenda& iy_clearsky_agenda,
612  const String& y_unit,
613  const Index& j_analytical_do,
614  const ArrayOfRetrievalQuantity& jacobian_quantities,
615  const ArrayOfArrayOfIndex& jacobian_indices,
616  const Verbosity& verbosity )
617 {
618  // Sizes
619  const Index nf = f_grid.nelem();
620  const Index nza = mblock_za_grid.nelem();
621  Index naa = mblock_aa_grid.nelem();
622  if( antenna_dim == 1 )
623  { naa = 1; }
624  const Index niyb = nf * nza * naa * stokes_dim;
625 
626  // Set up size of containers for data of 1 measurement block.
627  // (can not be made below due to parallalisation)
628  iyb.resize( niyb );
629  iyb_error.resize( niyb );
630  iyb_aux.resize( niyb );
631  Index aux_set = 0;
632  //
633  iyb_error = 0;
634  //
635  if( j_analytical_do )
636  {
637  diyb_dx.resize( jacobian_indices.nelem() );
638  FOR_ANALYTICAL_JACOBIANS_DO(
639  diyb_dx[iq].resize( niyb, jacobian_indices[iq][1] -
640  jacobian_indices[iq][0] + 1 );
641  )
642  }
643  else
644  { diyb_dx.resize( 0 ); }
645 
646  // Polarisation index variable
647  // If some kind of modified Stokes vector will be introduced, this must be
648  // made to a WSV
649  ArrayOfIndex i_pol(stokes_dim);
650  for( Index is=0; is<stokes_dim; is++ )
651  { i_pol[is] = is + 1; }
652 
653  // We have to make a local copy of the Workspace and the agendas because
654  // only non-reference types can be declared firstprivate in OpenMP
655  Workspace l_ws (ws);
656  Agenda l_iy_clearsky_agenda (iy_clearsky_agenda);
657 
658  // Start of actual calculations
659 /*#pragma omp parallel for \
660  if(!arts_omp_in_parallel() && nza>1) \
661  default(none) \
662  firstprivate(l_ws, l_iy_clearsky_agenda) \
663  shared(sensor_los, mblock_za_grid, mblock_aa_grid, vmr_field, \
664  t_field, lon_grid, lat_grid, p_grid, f_grid, sensor_pos, \
665  joker, naa) */
666 #pragma omp parallel for \
667  if(!arts_omp_in_parallel() && nza>1) \
668  firstprivate(l_ws, l_iy_clearsky_agenda)
669  for( Index iza=0; iza<nza; iza++ )
670  {
671  // The try block here is necessary to correctly handle
672  // exceptions inside the parallel region.
673  try
674  {
675  for( Index iaa=0; iaa<naa; iaa++ )
676  {
677  //--- LOS of interest
678  //
679  Vector los( sensor_los.ncols() );
680  //
681  los = sensor_los( imblock, joker );
682  los[0] += mblock_za_grid[iza];
683 
684  // Handle za/aa_grid "out-of-bounds" and mapping effects
685  //
686  if( antenna_dim == 2 )
687  { map_daa( los[0], los[1], los[0], los[1],
688  mblock_aa_grid[iaa] ); }
689  else if( atmosphere_dim == 1 && abs(los[0]-90) > 90 )
690  { if( los[0] < 0 ) { los[0] = -los[0]; }
691  else if( los[0] > 180 ) { los[0] = 360 - los[0]; } }
692  else if( atmosphere_dim == 2 && abs(los[0]) > 180 )
693  { los[0] = -sign(los[0])*360 + los[0]; }
694  else if( atmosphere_dim == 3 && abs(los[0]-90) > 90 )
695  { map_daa( los[0], los[1], los[0], los[1], 0 ); }
696 
697 
698  // Calculate iy and associated variables
699  //
700  Matrix iy(0,0), iy_error(0,0), iy_aux(0,0);
701  Tensor3 iy_transmission(0,0,0);
702  ArrayOfTensor3 diy_dx(0);
703  //
704  iy_error_type = 0;
705  //
706  iy_clearsky_agendaExecute( l_ws, iy, iy_error, iy_error_type,
707  iy_aux, diy_dx, 1, iy_transmission,
708  sensor_pos(imblock,joker), los, cloudbox_on,
709  j_analytical_do, p_grid, lat_grid, lon_grid,
710  t_field, z_field, vmr_field, l_iy_clearsky_agenda );
711 
712  // Start row in iyb etc. for present LOS
713  //
714  const Index row0 = ( iza*naa + iaa ) * nf * stokes_dim;
715 
716  // Jacobian part (must be converted to Tb before iy for PlanckBT)
717  //
718  if( j_analytical_do )
719  {
720  FOR_ANALYTICAL_JACOBIANS_DO(
721  //
722  apply_y_unit2( diy_dx[iq], iy, y_unit, f_grid, i_pol );
723  //
724  for( Index ip=0; ip<jacobian_indices[iq][1] -
725  jacobian_indices[iq][0]+1; ip++ )
726  {
727  for( Index is=0; is<stokes_dim; is++ )
728  {
729  diyb_dx[iq](Range(row0+is,nf,stokes_dim),ip)=
730  diy_dx[iq](ip,joker,is);
731  }
732  }
733  )
734  }
735 
736  // iy : unit conversion and copy to iyb
737  // iy_error : unit conversion and copy to iyb_error
738  // iy_aux : copy to iyb_aux (if iy_aux filled)
739  //
740  if( iy_error_type > 0 )
741  { apply_y_unit2( iy_error, iy, y_unit, f_grid, i_pol ); }
742  //
743  apply_y_unit( iy, y_unit, f_grid, i_pol );
744  //
745  for( Index is=0; is<stokes_dim; is++ )
746  {
747  iyb[Range(row0+is,nf,stokes_dim)] = iy(joker,is);
748  //
749  if( iy_error_type > 0 )
750  {
751  iyb_error[Range(row0+is,nf,stokes_dim)] =
752  iy_error(joker,is);
753  }
754  //
755  if( iy_aux.nrows() )
756  {
757  iyb_aux[Range(row0+is,nf,stokes_dim)] = iy_aux(joker,is);
758  aux_set = 1;
759  }
760  }
761  } // End aa loop
762  } // End try
763 
764  catch (runtime_error e)
765  {
766  CREATE_OUT0
767  exit_or_rethrow(e.what(), out0);
768  }
769  } // End za loop
770 
771  // If no aux, set to size 0 to flag this
772  if( !aux_set )
773  { iyb_aux.resize(0); }
774 }
775 
776 
777 
778 
779 
780 
781 
782 
783 
784 /*===========================================================================
785  === The workspace methods (in alphabetical order)
786  ===========================================================================*/
787 
788 /* Workspace method: Doxygen documentation will be auto-generated */
789 void iyBeerLambertStandardClearsky(
790  Workspace& ws,
791  Matrix& iy,
792  Matrix& iy_error,
793  Index& iy_error_type,
794  Matrix& iy_aux,
795  ArrayOfTensor3& diy_dx,
796  const Index& iy_agenda_call1,
797  const Tensor3& iy_transmission,
798  const Vector& rte_pos,
799  const Vector& rte_los,
800  const Index& jacobian_do,
801  const Index& atmosphere_dim,
802  const Vector& p_grid,
803  const Vector& lat_grid,
804  const Vector& lon_grid,
805  const Tensor3& z_field,
806  const Tensor3& t_field,
807  const Tensor4& vmr_field,
808  const Tensor3& wind_u_field,
809  const Tensor3& wind_v_field,
810  const Tensor3& wind_w_field,
811  const Matrix& r_geoid,
812  const Matrix& z_surface,
813  const Index& cloudbox_on,
814  const ArrayOfIndex& cloudbox_limits,
815  const Index& stokes_dim,
816  const Vector& f_grid,
817  const ArrayOfArrayOfSpeciesTag& abs_species,
818  const Agenda& ppath_step_agenda,
819  const Agenda& abs_scalar_gas_agenda,
820  const Agenda& iy_clearsky_agenda,
821  const Agenda& iy_space_agenda,
822  const Agenda& surface_prop_agenda,
823  const Agenda& iy_cloudbox_agenda,
824  const ArrayOfRetrievalQuantity& jacobian_quantities,
825  const ArrayOfArrayOfIndex& jacobian_indices,
826  const Verbosity& verbosity )
827 {
828  // See initial comments of iyEmissionStandardClearsky
829 
830  // Sizes
831  //
832  const Index nf = f_grid.nelem();
833  const Index nq = jacobian_quantities.nelem();
834 
835  // Determine if there are any analytical jacobians to handle, and if primary
836  // call, resize diy_dx.
837  //
838  Index j_analytical_do = 0;
839  //
840  if( jacobian_do ) { FOR_ANALYTICAL_JACOBIANS_DO( j_analytical_do = 1; ) }
841  //
842  if( iy_agenda_call1 && j_analytical_do )
843  {
844  diy_dx.resize( jacobian_indices.nelem() );
845  //
846  FOR_ANALYTICAL_JACOBIANS_DO(
847  diy_dx[iq].resize( jacobian_indices[iq][1] - jacobian_indices[iq][0] +
848  1, nf, stokes_dim );
849  diy_dx[iq] = 0.0;
850  )
851  }
852 
853  //- Determine propagation path
854  //
855  Ppath ppath;
856  //
857  ppath_calc( ws, ppath, ppath_step_agenda, atmosphere_dim, p_grid,
858  lat_grid, lon_grid, z_field, r_geoid, z_surface,
859  cloudbox_on, cloudbox_limits, rte_pos, rte_los, 1,
860  verbosity );
861 
862  // Get atmospheric and RT quantities for each ppath point/step
863  //
864  // If np = 1, we only need to determine the radiative background
865  //
866  // "atmvars"
867  Vector ppath_p, ppath_t, ppath_wind_u, ppath_wind_v, ppath_wind_w;
868  Matrix ppath_vmr;
869  // "rtvars"
870  Vector total_tau;
871  Matrix emission_dummy, ppath_tau;
872  Tensor3 ppath_abs_scalar, iy_trans_new;
873  Agenda agenda_dummy;
874  //
875  const Index np = ppath.np;
876  //
877  if( np > 1 )
878  {
879  // Get pressure, temperature and VMRs
880  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
881  ppath_wind_u, ppath_wind_v, ppath_wind_w,
882  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
883  wind_u_field, wind_v_field, wind_w_field );
884 
885  // Absorption and optical thickness for each step
886  get_ppath_rtvars( ws, ppath_abs_scalar, ppath_tau, total_tau,
887  emission_dummy, abs_scalar_gas_agenda, agenda_dummy,
888  ppath, ppath_p, ppath_t, ppath_vmr, ppath_wind_u,
889  ppath_wind_v, ppath_wind_w, f_grid, atmosphere_dim, 0 );
890  }
891  else // To handle cases inside the cloudbox, or totally outside the atmosphere
892  {
893  total_tau.resize( nf );
894  total_tau = 0;
895  }
896 
897  // iy_transmission
898  //
899  if( iy_agenda_call1 )
900  { iy_transmission_for_scalar_tau( iy_trans_new, stokes_dim, total_tau ); }
901  else
902  { iy_transmission_mult_scalar_tau( iy_trans_new, iy_transmission,
903  total_tau ); }
904 
905  // Radiative background
906  //
907  get_iy_of_background( ws, iy, iy_error, iy_error_type, iy_aux, diy_dx,
908  iy_trans_new, jacobian_do, ppath, atmosphere_dim,
909  p_grid, lat_grid, lon_grid, t_field, z_field,
910  vmr_field, r_geoid, z_surface, cloudbox_on, stokes_dim,
911  f_grid, iy_clearsky_agenda, iy_space_agenda,
912  surface_prop_agenda, iy_cloudbox_agenda, verbosity );
913 
914  // Do RT calculations
915  //
916  if( np > 1 )
917  {
918  // Create container for the derivatives with respect to changes
919  // at each ppath point. And find matching abs_species-index and
920  // "temperature flag" (for analytical jacobians).
921  //
922  ArrayOfTensor3 diy_dpath;
923  ArrayOfIndex abs_species_i, is_t;
924  //
925  const Numeric dt = 0.1;
926  Tensor3 ppath_as2;
927  //
928  if( j_analytical_do )
929  {
930  diy_dpath.resize( nq );
931  abs_species_i.resize( nq );
932  is_t.resize( nq );
933  //
934  FOR_ANALYTICAL_JACOBIANS_DO(
935  diy_dpath[iq].resize( np, nf, stokes_dim );
936  diy_dpath[iq] = 0.0;
937  )
938  get_pointers_for_analytical_jacobians( abs_species_i, is_t,
939  jacobian_quantities, abs_species );
940  //
941  // Determine if temperature is among the analytical jac. quantities.
942  // If yes, get emission and absorption for disturbed temperature
943  Index do_t=0;
944  for( Index iq=0; iq<is_t.nelem() && do_t==0; iq++ )
945  { if( is_t[iq] ) { do_t = 1; } }
946  if( do_t )
947  {
948  Matrix tau_dummy;
949  Vector total_tau_dummy, t2 = ppath_t;
950  t2 += dt;
951  get_ppath_rtvars( ws, ppath_as2, tau_dummy, total_tau_dummy,
952  emission_dummy, abs_scalar_gas_agenda,
953  agenda_dummy, ppath, ppath_p, t2, ppath_vmr,
954  ppath_wind_u, ppath_wind_v, ppath_wind_w,
955  f_grid, atmosphere_dim, 0 );
956  }
957  }
958 
959  // Loop ppath steps
960  for( Index ip=np-2; ip>=0; ip-- )
961  {
962 
963  // Jacobian quantities
964  if( j_analytical_do )
965  {
966  // Common terms introduced for efficiency and clarity
967  Vector X(nf); // See AUG
968  //
969  for( Index iv=0; iv<nf; iv++ )
970  {
971  X[iv] = 0.5 * ppath.l_step[ip] * exp( -total_tau[iv] );
972  }
973 
974  // Loop quantities
975  // Loop quantities
976  for( Index iq=0; iq<nq; iq++ )
977  {
978  if( jacobian_quantities[iq].Analytical() )
979  {
980  // Absorbing species
981  const Index isp = abs_species_i[iq];
982  if( isp >= 0 )
983  {
984  // Scaling factors to handle retrieval unit
985  // (gives the cross-section to apply)
986  Numeric unitscf1, unitscf2;
987  vmrunitscf( unitscf1,
988  jacobian_quantities[iq].Mode(),
989  ppath_vmr(isp,ip), ppath_p[ip],
990  ppath_t[ip] );
991  vmrunitscf( unitscf2,
992  jacobian_quantities[iq].Mode(),
993  ppath_vmr(isp,ip+1), ppath_p[ip+1],
994  ppath_t[ip+1] );
995  //
996  for( Index iv=0; iv<nf; iv++ )
997  {
998  // All Stokes components equal
999  for( Index is=0; is<stokes_dim; is++ )
1000  {
1001  const Numeric Z = -X[iv] * iy(iv,is);
1002  diy_dpath[iq](ip ,iv,is) += Z * unitscf1 *
1003  ppath_abs_scalar(iv,isp,ip);
1004  diy_dpath[iq](ip+1,iv,is) += Z * unitscf2 *
1005  ppath_abs_scalar(iv,isp,ip+1);
1006  }
1007  }
1008  }
1009 
1010  // Temperature
1011  else if( is_t[iq] )
1012  {
1013  for( Index iv=0; iv<nf; iv++ )
1014  {
1015  // The terms associated with Dtau/Dk:
1016  const Numeric k1 =
1017  ppath_abs_scalar(iv,joker,ip ).sum();
1018  const Numeric k2 =
1019  ppath_abs_scalar(iv,joker,ip+1).sum();
1020  const Numeric dkdt1 =
1021  ( ppath_as2(iv,joker,ip ).sum() - k1 ) / dt;
1022  const Numeric dkdt2 =
1023  ( ppath_as2(iv,joker,ip+1).sum() - k2 ) / dt;
1024  for( Index is=0; is<stokes_dim; is++ )
1025  {
1026  const Numeric Z = -X[iv] * iy(iv,is);
1027  diy_dpath[iq](ip ,iv,is) += Z * dkdt1;
1028  diy_dpath[iq](ip+1,iv,is) += Z * dkdt2;
1029  }
1030  //
1031  // The terms associated with Delta-s:
1032  if( jacobian_quantities[iq].Subtag() == "HSE on" )
1033  {
1034  const Numeric kbar = 0.5 * ( k1 + k2 );
1035  for( Index is=0; is<stokes_dim; is++ )
1036  {
1037  const Numeric Z = -X[iv] * iy(iv,is);
1038  diy_dpath[iq](ip ,iv,is) += Z * kbar /
1039  ppath_t[ip];
1040  diy_dpath[iq](ip+1,iv,is) += Z * kbar /
1041  ppath_t[ip+1];
1042  }
1043  } //hse
1044  } // frequency
1045  } // if is_t
1046  } // if analytical
1047  } // for iq
1048 
1049  // Update total_tau (not used for spectrum calculation)
1050  for( Index iv=0; iv<nf; iv++ )
1051  { total_tau[iv] -= ppath_tau(iv,ip); }
1052  }
1053 
1054  // Spectrum
1055  //
1056  for( Index iv=0; iv<nf; iv++ )
1057  {
1058  const Numeric step_tr = exp( -ppath_tau(iv,ip) );
1059  //
1060  for( Index is=0; is<stokes_dim; is++ )
1061  { iy(iv,is) *= step_tr; }
1062  }
1063  }
1064 
1065  // Map jacobians from ppath to retrieval grids
1066  // (this operation corresponds to the term Dx_i/Dx)
1067  if( j_analytical_do )
1068  {
1069  // Weight with iy_transmission
1070  if( !iy_agenda_call1 )
1071  {
1072  Matrix X, Y(stokes_dim,diy_dpath[0].npages());
1073  //
1074  FOR_ANALYTICAL_JACOBIANS_DO(
1075  for( Index iv=0; iv<nf; iv++ )
1076  {
1077  X = transpose( diy_dpath[iq](joker,iv,joker) );
1078  mult( Y, iy_transmission(iv,joker,joker), X );
1079  diy_dpath[iq](joker,iv,joker) = transpose( Y );
1080  }
1081  )
1082  }
1083 
1084  // Map to retrieval grids
1085  FOR_ANALYTICAL_JACOBIANS_DO(
1086  diy_from_path_to_rgrids( diy_dx[iq], jacobian_quantities[iq],
1087  diy_dpath[iq], atmosphere_dim, ppath, ppath_p );
1088  )
1089  }
1090  }
1091 }
1092 
1093 
1094 
1095 
1096 /* Workspace method: Doxygen documentation will be auto-generated */
1097 void iyBeerLambertStandardCloudbox(
1098  Workspace& ws,
1099  Matrix& iy,
1100  Matrix& iy_error,
1101  Index& iy_error_type,
1102  Matrix& iy_aux,
1103  ArrayOfTensor3& diy_dx,
1104  const Tensor3& iy_transmission,
1105  const Vector& rte_pos,
1106  const Vector& rte_los,
1107  const Index& jacobian_do,
1108  const Index& atmosphere_dim,
1109  const Vector& p_grid,
1110  const Vector& lat_grid,
1111  const Vector& lon_grid,
1112  const Tensor3& z_field,
1113  const Tensor3& t_field,
1114  const Tensor4& vmr_field,
1115  const Matrix& r_geoid,
1116  const Matrix& z_surface,
1117  const Index& cloudbox_on,
1118  const ArrayOfIndex& cloudbox_limits,
1119  const Index& stokes_dim,
1120  const Vector& f_grid,
1121  const Agenda& ppath_step_agenda,
1122  const Agenda& abs_scalar_gas_agenda,
1123  const Agenda& iy_clearsky_agenda,
1124  const Tensor4& pnd_field,
1125  const Index& use_mean_scat_data,
1126  const ArrayOfSingleScatteringData& scat_data_raw,
1127  const Agenda& opt_prop_gas_agenda,
1128  const Verbosity& verbosity)
1129 {
1130  // Input checks
1131  if( !cloudbox_on )
1132  throw runtime_error( "The cloudbox must be defined to use this method." );
1133  if( jacobian_do )
1134  throw runtime_error(
1135  "This method does not provide any jacobians (jacobian_do must be 0)." );
1136 
1137  // Sizes
1138  //
1139  const Index nf = f_grid.nelem();
1140 
1141  // Determine ppath through the cloudbox
1142  //
1143  Ppath ppath;
1144  //
1145  ppath_calc( ws, ppath, ppath_step_agenda, atmosphere_dim, p_grid,
1146  lat_grid, lon_grid, z_field, r_geoid, z_surface,
1147  cloudbox_on, cloudbox_limits, rte_pos, rte_los, 0,
1148  verbosity );
1149 
1150  // Check radiative background
1151  const Index bkgr = ppath_what_background( ppath );
1152  if( bkgr == 2 )
1153  throw runtime_error( "Observations where (unscattered) propagation path "
1154  "hits the surface inside the cloudbox are not yet "
1155  "handled by this method." );
1156  assert( bkgr == 3 );
1157 
1158  // Get atmospheric and RT quantities for each ppath point/step (inside box)
1159  //
1160  // If np = 1, we only need to determine the radiative background
1161  //
1162  // "atmvars"
1163  Vector ppath_p, ppath_t, ppath_wind_u, ppath_wind_v, ppath_wind_w;
1164  Matrix ppath_vmr, ppath_pnd;
1165  // "rtvars"
1166  Matrix emission_dummy, ppath_tau;
1167  Tensor3 wind_field_dummy(0,0,0), iy_trans_new;
1168  Tensor3 ppath_asp_abs_vec, ppath_pnd_abs_vec, total_transmission;
1169  Tensor4 ppath_asp_ext_mat, ppath_pnd_ext_mat, ppath_transmission;
1170  Agenda agenda_dummy;
1171  Array<ArrayOfSingleScatteringData> scat_data;
1172  //
1173  const Index np = ppath.np;
1174  //
1175  if( np > 1 )
1176  {
1177  // Get pressure, temperature and VMRs
1178  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
1179  ppath_wind_u, ppath_wind_v, ppath_wind_w,
1180  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
1181  wind_field_dummy, wind_field_dummy, wind_field_dummy );
1182 
1183  // Particle number densities
1184  get_ppath_pnd( ppath_pnd, ppath, atmosphere_dim, cloudbox_limits,
1185  pnd_field );
1186 
1187  // Absorption and optical thickness for each step
1188  get_ppath_cloudrtvars( ws, ppath_asp_abs_vec, ppath_asp_ext_mat,
1189  ppath_pnd_abs_vec, ppath_pnd_ext_mat,
1190  ppath_transmission, total_transmission,
1191  emission_dummy, scat_data, abs_scalar_gas_agenda,
1192  agenda_dummy, opt_prop_gas_agenda,
1193  ppath, ppath_p, ppath_t, ppath_vmr, ppath_wind_u,
1194  ppath_wind_v, ppath_wind_w, ppath_pnd,
1195  use_mean_scat_data, scat_data_raw, stokes_dim,
1196  f_grid, atmosphere_dim, 0, verbosity );
1197 
1198  // *scat_data* not used. Free the memory.
1199  scat_data.resize( 0 );
1200  }
1201  else // Just in case, should not happen
1202  { assert( 0 ); }
1203 
1204  // iy_transmission
1205  //
1206  iy_transmission_mult( iy_trans_new, iy_transmission, total_transmission );
1207 
1208  // Get iy for unscattered direction
1209  //
1210  // Note that the Ppath positions (ppath.pos) for 1D have one column more
1211  // than expected by most functions. Only the first atmosphere_dim values
1212  // shall be copied.
1213  //
1214  {
1215  Vector rte_pos2, rte_los2;
1216  rte_pos2 = ppath.pos(ppath.np-1,Range(0,atmosphere_dim));
1217  rte_los2 = ppath.los(ppath.np-1,joker);
1218  //
1219  iy_clearsky_agendaExecute( ws, iy, iy_error, iy_error_type,
1220  iy_aux, diy_dx, 0, iy_trans_new,
1221  rte_pos2, rte_los2, 0, jacobian_do,
1222  p_grid, lat_grid, lon_grid, t_field,
1223  z_field, vmr_field, iy_clearsky_agenda );
1224  }
1225 
1226  // Without jacobians the complete RT calculations are just multiplication
1227  // with *total_transmission*
1228  for( Index iv=0; iv<nf; iv++ )
1229  {
1230  const Vector tmp = iy(iv,joker);
1231  mult( iy(iv,joker), total_transmission(iv,joker,joker), tmp );
1232  }
1233 }
1234 
1235 
1236 
1237 
1238 
1239 /* Workspace method: Doxygen documentation will be auto-generated */
1240 void iyEmissionStandardClearsky(
1241  Workspace& ws,
1242  Matrix& iy,
1243  Matrix& iy_error,
1244  Index& iy_error_type,
1245  Matrix& iy_aux,
1246  ArrayOfTensor3& diy_dx,
1247  const Index& iy_agenda_call1,
1248  const Tensor3& iy_transmission,
1249  const Vector& rte_pos,
1250  const Vector& rte_los,
1251  const Index& jacobian_do,
1252  const Index& atmosphere_dim,
1253  const Vector& p_grid,
1254  const Vector& lat_grid,
1255  const Vector& lon_grid,
1256  const Tensor3& z_field,
1257  const Tensor3& t_field,
1258  const Tensor4& vmr_field,
1259  const Tensor3& wind_u_field,
1260  const Tensor3& wind_v_field,
1261  const Tensor3& wind_w_field,
1262  const Matrix& r_geoid,
1263  const Matrix& z_surface,
1264  const Index& cloudbox_on,
1265  const ArrayOfIndex& cloudbox_limits,
1266  const Index& stokes_dim,
1267  const Vector& f_grid,
1268  const ArrayOfArrayOfSpeciesTag& abs_species,
1269  const Agenda& ppath_step_agenda,
1270  const Agenda& emission_agenda,
1271  const Agenda& abs_scalar_gas_agenda,
1272  const Agenda& iy_clearsky_agenda,
1273  const Agenda& iy_space_agenda,
1274  const Agenda& surface_prop_agenda,
1275  const Agenda& iy_cloudbox_agenda,
1276  const ArrayOfRetrievalQuantity& jacobian_quantities,
1277  const ArrayOfArrayOfIndex& jacobian_indices,
1278  const Verbosity& verbosity )
1279 {
1280  // The method can in principle be used "stand-alone", but for efficiency
1281  // reasons we skip all checks needed to handle such usage. Those checks are
1282  // done by yCalc.
1283 
1284  // iy_error, iy_error_type, iy_aux and diy_dx are initiated by iyb_calc
1285  // to be empty/zero. Accordingly, if any radiative background wants
1286  // to flag an error, iy_error must be resized. No subsequent scaling of the
1287  // error is made. Thus the iy_transmission should be considered when adding
1288  // a term to the error. A setting of iy_aux must also include a resizing.
1289 
1290  // The WSV *iy_transmission* holds the transmission between the sensor
1291  // and the end of the propagation path.
1292 
1293  // Sizes
1294  //
1295  const Index nf = f_grid.nelem();
1296  const Index nq = jacobian_quantities.nelem();
1297 
1298  // Determine if there are any analytical jacobians to handle, and if primary
1299  // call, resize diy_dx.
1300  //
1301  Index j_analytical_do = 0;
1302  //
1303  if( jacobian_do ) { FOR_ANALYTICAL_JACOBIANS_DO( j_analytical_do = 1; ) }
1304  //
1305  if( iy_agenda_call1 && j_analytical_do )
1306  {
1307  diy_dx.resize( nq );
1308  //
1309  FOR_ANALYTICAL_JACOBIANS_DO(
1310  diy_dx[iq].resize( jacobian_indices[iq][1] - jacobian_indices[iq][0] +
1311  1, nf, stokes_dim );
1312  diy_dx[iq] = 0.0;
1313  )
1314  }
1315 
1316  //- Determine propagation path
1317  //
1318  Ppath ppath;
1319  //
1320  ppath_calc( ws, ppath, ppath_step_agenda, atmosphere_dim, p_grid,
1321  lat_grid, lon_grid, z_field, r_geoid, z_surface,
1322  cloudbox_on, cloudbox_limits, rte_pos, rte_los, 1,
1323  verbosity );
1324 
1325  // Get atmospheric and RT quantities for each ppath point/step
1326  //
1327  // If np = 1, we only need to determine the radiative background
1328  //
1329  // "atmvars"
1330  Vector ppath_p, ppath_t, ppath_wind_u, ppath_wind_v, ppath_wind_w;
1331  Matrix ppath_vmr;
1332  // "rtvars"
1333  Vector total_tau;
1334  Matrix ppath_emission, ppath_tau;
1335  Tensor3 ppath_abs_scalar, iy_trans_new;
1336  //
1337  const Index np = ppath.np;
1338  //
1339  if( np > 1 )
1340  {
1341  // Get pressure, temperature and VMRs
1342  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
1343  ppath_wind_u, ppath_wind_v, ppath_wind_w,
1344  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
1345  wind_u_field, wind_v_field, wind_w_field );
1346 
1347  // Get emission, absorption and optical thickness for each step
1348  get_ppath_rtvars( ws, ppath_abs_scalar, ppath_tau, total_tau,
1349  ppath_emission, abs_scalar_gas_agenda, emission_agenda,
1350  ppath, ppath_p, ppath_t, ppath_vmr, ppath_wind_u,
1351  ppath_wind_v, ppath_wind_w, f_grid, atmosphere_dim, 1 );
1352  }
1353  else // To handle cases inside the cloudbox, or totally outside the atmosphere
1354  {
1355  total_tau.resize( nf );
1356  total_tau = 0;
1357  }
1358 
1359  // iy_transmission
1360  //
1361  if( iy_agenda_call1 )
1362  { iy_transmission_for_scalar_tau( iy_trans_new, stokes_dim, total_tau ); }
1363  else
1364  { iy_transmission_mult_scalar_tau( iy_trans_new, iy_transmission,
1365  total_tau ); }
1366 
1367  // Radiative background
1368  //
1369  get_iy_of_background( ws, iy, iy_error, iy_error_type, iy_aux, diy_dx,
1370  iy_trans_new, jacobian_do, ppath, atmosphere_dim,
1371  p_grid, lat_grid, lon_grid, t_field, z_field,
1372  vmr_field, r_geoid, z_surface, cloudbox_on, stokes_dim,
1373  f_grid, iy_clearsky_agenda, iy_space_agenda,
1374  surface_prop_agenda, iy_cloudbox_agenda, verbosity);
1375 
1376  // Do RT calculations
1377  //
1378  if( np > 1 )
1379  {
1380  // Create container for the derivatives with respect to changes
1381  // at each ppath point. And find matching abs_species-index and
1382  // "temperature flag" (for analytical jacobians).
1383  //
1384  ArrayOfTensor3 diy_dpath;
1385  ArrayOfIndex abs_species_i, is_t;
1386  //
1387  const Numeric dt = 0.1;
1388  Matrix ppath_e2;
1389  Tensor3 ppath_as2;
1390  //
1391  if( j_analytical_do )
1392  {
1393  diy_dpath.resize( nq );
1394  abs_species_i.resize( nq );
1395  is_t.resize( nq );
1396  //
1397  FOR_ANALYTICAL_JACOBIANS_DO(
1398  diy_dpath[iq].resize( np, nf, stokes_dim );
1399  diy_dpath[iq] = 0.0;
1400  )
1401  get_pointers_for_analytical_jacobians( abs_species_i, is_t,
1402  jacobian_quantities, abs_species );
1403  //
1404  // Determine if temperature is among the analytical jac. quantities.
1405  // If yes, get emission and absorption for disturbed temperature
1406  Index do_t = 0;
1407  for( Index iq=0; iq<is_t.nelem() && do_t==0; iq++ )
1408  { if( is_t[iq] ) { do_t = 1; } }
1409  if( do_t )
1410  {
1411  Matrix tau_dummy;
1412  Vector total_tau_dummy, t2 = ppath_t;
1413  t2 += dt;
1414  get_ppath_rtvars( ws, ppath_as2, tau_dummy, total_tau_dummy,
1415  ppath_e2, abs_scalar_gas_agenda,
1416  emission_agenda, ppath, ppath_p, t2, ppath_vmr,
1417  ppath_wind_u, ppath_wind_v, ppath_wind_w,
1418  f_grid, atmosphere_dim, 1 );
1419  }
1420  }
1421 
1422  // Loop ppath steps
1423  for( Index ip=np-2; ip>=0; ip-- )
1424  {
1425  Vector esource(nf); // B-bar
1426  Matrix iy_new(nf,stokes_dim); // Intensity at end of step
1427 
1428  // Spectrum at end of ppath step
1429  // (iy updated to iy_new at end of ip-loop)
1430  for( Index iv=0; iv<nf; iv++ )
1431  {
1432  const Numeric step_tr = exp( -ppath_tau(iv,ip) );
1433  //
1434  // First Stokes element:
1435  esource[iv] = 0.5 * ( ppath_emission(iv,ip) +
1436  ppath_emission(iv,ip+1) );
1437  iy_new(iv,0) = iy(iv,0) * step_tr + esource[iv] * (1-step_tr);
1438  //
1439  // Higher Stokes:
1440  for( Index is=1; is<stokes_dim; is++ )
1441  { iy_new(iv,is) = step_tr * iy(iv,is); }
1442  }
1443 
1444  // Jacobian quantities
1445  if( j_analytical_do )
1446  {
1447  // Common terms introduced for efficiency and clarity
1448  Vector X(nf), Y(nf); // See AUG
1449  //
1450  for( Index iv=0; iv<nf; iv++ )
1451  {
1452  X[iv] = 0.5 * ppath.l_step[ip] * exp( -total_tau[iv] );
1453  Y[iv] = X[iv] * ( esource[iv] - iy(iv,0) );
1454  }
1455 
1456  // Loop quantities
1457  for( Index iq=0; iq<nq; iq++ )
1458  {
1459  if( jacobian_quantities[iq].Analytical() )
1460  {
1461  // Absorbing species
1462  const Index isp = abs_species_i[iq];
1463  if( isp >= 0 )
1464  {
1465  // Scaling factors to handle retrieval unit
1466  // (gives the cross-section to apply)
1467  Numeric unitscf1, unitscf2;
1468  vmrunitscf( unitscf1,
1469  jacobian_quantities[iq].Mode(),
1470  ppath_vmr(isp,ip), ppath_p[ip],
1471  ppath_t[ip] );
1472  vmrunitscf( unitscf2,
1473  jacobian_quantities[iq].Mode(),
1474  ppath_vmr(isp,ip+1), ppath_p[ip+1],
1475  ppath_t[ip+1] );
1476  //
1477  for( Index iv=0; iv<nf; iv++ )
1478  {
1479  // Stokes component 1
1480  diy_dpath[iq](ip ,iv,0) += Y[iv] * unitscf1 *
1481  ppath_abs_scalar(iv,isp,ip);
1482  diy_dpath[iq](ip+1,iv,0) += Y[iv] * unitscf2 *
1483  ppath_abs_scalar(iv,isp,ip+1);
1484  // Higher stokes components
1485  for( Index is=1; is<stokes_dim; is++ )
1486  {
1487  const Numeric Z = -X[iv] * iy(iv,is);
1488  diy_dpath[iq](ip ,iv,is) += Z * unitscf1 *
1489  ppath_abs_scalar(iv,isp,ip);
1490  diy_dpath[iq](ip+1,iv,is) += Z * unitscf2 *
1491  ppath_abs_scalar(iv,isp,ip+1);
1492  }
1493  }
1494  }
1495 
1496  // Temperature
1497  else if( is_t[iq] )
1498  {
1499  for( Index iv=0; iv<nf; iv++ )
1500  {
1501  // The terms associated with Dtau/Dk:
1502  const Numeric k1 =
1503  ppath_abs_scalar(iv,joker,ip ).sum();
1504  const Numeric k2 =
1505  ppath_abs_scalar(iv,joker,ip+1).sum();
1506  const Numeric dkdt1 =
1507  ( ppath_as2(iv,joker,ip ).sum() - k1 ) / dt;
1508  const Numeric dkdt2 =
1509  ( ppath_as2(iv,joker,ip+1).sum() - k2 ) / dt;
1510  // Stokes 1:
1511  diy_dpath[iq](ip ,iv,0) += Y[iv] * dkdt1;
1512  diy_dpath[iq](ip+1,iv,0) += Y[iv] * dkdt2;
1513  // Higher Stokes
1514  for( Index is=1; is<stokes_dim; is++ )
1515  {
1516  const Numeric Z = -X[iv] * iy(iv,is);
1517  diy_dpath[iq](ip ,iv,is) += Z * dkdt1;
1518  diy_dpath[iq](ip+1,iv,is) += Z * dkdt2;
1519  }
1520  //
1521  // The terms associated with B-bar:
1522  const Numeric V = 0.5 *
1523  exp(-(total_tau[iv]-ppath_tau(iv,ip))) *
1524  ( 1.0 - exp(-ppath_tau(iv,ip)) );
1525  diy_dpath[iq](ip ,iv,0) += V *
1526  ( ppath_e2(iv,ip) -
1527  ppath_emission(iv,ip) ) / dt;
1528  diy_dpath[iq](ip+1,iv,0) += V *
1529  ( ppath_e2(iv,ip+1) -
1530  ppath_emission(iv,ip+1) ) / dt;
1531  // Zero for higher Stokes
1532  //
1533  // The terms associated with Delta-s:
1534  if( jacobian_quantities[iq].Subtag() == "HSE on" )
1535  {
1536  // Stokes 1:
1537  const Numeric kbar = 0.5 * ( k1 + k2 );
1538  diy_dpath[iq](ip ,iv,0) += Y[iv] * kbar /
1539  ppath_t[ip];
1540  diy_dpath[iq](ip+1,iv,0) += Y[iv] * kbar /
1541  ppath_t[ip+1];
1542  // Higher Stokes
1543  for( Index is=1; is<stokes_dim; is++ )
1544  {
1545  const Numeric Z = -X[iv] * iy(iv,is);
1546  diy_dpath[iq](ip ,iv,is) += Z * kbar /
1547  ppath_t[ip];
1548  diy_dpath[iq](ip+1,iv,is) += Z * kbar /
1549  ppath_t[ip+1];
1550  }
1551  } //hse
1552  } // frequency
1553  } // if is_t
1554  } // if analytical
1555  } // for iq
1556 
1557  // Update total_tau (not used for spectrum calculation)
1558  for( Index iv=0; iv<nf; iv++ )
1559  { total_tau[iv] -= ppath_tau(iv,ip); }
1560  }
1561 
1562  // Update iy
1563  iy = iy_new;
1564  }
1565 
1566  // Map jacobians from ppath to retrieval grids
1567  // (this operation corresponds to the term Dx_i/Dx)
1568  if( j_analytical_do )
1569  {
1570  // Weight with iy_transmission
1571  if( !iy_agenda_call1 )
1572  {
1573  Matrix X, Y(stokes_dim,diy_dpath[0].npages());
1574  //
1575  FOR_ANALYTICAL_JACOBIANS_DO(
1576  for( Index iv=0; iv<nf; iv++ )
1577  {
1578  X = transpose( diy_dpath[iq](joker,iv,joker) );
1579  mult( Y, iy_transmission(iv,joker,joker), X );
1580  diy_dpath[iq](joker,iv,joker) = transpose( Y );
1581  }
1582  )
1583  }
1584 
1585  // Map to retrieval grids
1586  FOR_ANALYTICAL_JACOBIANS_DO(
1587  diy_from_path_to_rgrids( diy_dx[iq], jacobian_quantities[iq],
1588  diy_dpath[iq], atmosphere_dim, ppath, ppath_p );
1589  )
1590  }
1591  } // if np>1
1592 
1593  // Set iy_aux, if background is TOA or surface
1594  //
1595  if( ppath_what_background( ppath ) <= 2 )
1596  {
1597  iy_aux.resize( nf, stokes_dim );
1598  //
1599  for( Index iv=0; iv<nf; iv++ )
1600  {
1601  for( Index is=0; is<stokes_dim; is++ )
1602  {
1603  iy_aux( iv, is ) = iy_trans_new( iv, is, is );
1604  }
1605  }
1606  }
1607 }
1608 
1609 
1610 
1611 /* Workspace method: Doxygen documentation will be auto-generated */
1612 void iyEmissionStandardClearskyBasic(
1613  Workspace& ws,
1614  Matrix& iy,
1615  const Vector& rte_pos,
1616  const Vector& rte_los,
1617  const Index& jacobian_do,
1618  const Index& atmosphere_dim,
1619  const Vector& p_grid,
1620  const Vector& lat_grid,
1621  const Vector& lon_grid,
1622  const Tensor3& z_field,
1623  const Tensor3& t_field,
1624  const Tensor4& vmr_field,
1625  const Tensor3& wind_u_field,
1626  const Tensor3& wind_v_field,
1627  const Tensor3& wind_w_field,
1628  const Matrix& r_geoid,
1629  const Matrix& z_surface,
1630  const Index& cloudbox_on,
1631  const ArrayOfIndex& cloudbox_limits,
1632  const Index& stokes_dim,
1633  const Vector& f_grid,
1634  const Agenda& ppath_step_agenda,
1635  const Agenda& emission_agenda,
1636  const Agenda& abs_scalar_gas_agenda,
1637  const Agenda& iy_clearsky_basic_agenda,
1638  const Agenda& iy_space_agenda,
1639  const Agenda& surface_prop_agenda,
1640  const Agenda& iy_cloudbox_agenda,
1641  const Verbosity& verbosity )
1642 {
1643  // See initial comments of iyEmissionStandardClearsky
1644 
1645  if( jacobian_do )
1646  throw runtime_error(
1647  "This method does not provide any jacobians (jacobian_do must be 0)" );
1648 
1649  // Sizes
1650  //
1651  const Index nf = f_grid.nelem();
1652 
1653  //- Determine propagation path
1654  //
1655  Ppath ppath;
1656  //
1657  ppath_calc( ws, ppath, ppath_step_agenda, atmosphere_dim, p_grid,
1658  lat_grid, lon_grid, z_field, r_geoid, z_surface,
1659  cloudbox_on, cloudbox_limits, rte_pos, rte_los, 1,
1660  verbosity );
1661 
1662  // Get quantities required for each ppath point
1663  //
1664  // If np = 1, we only need to determine the radiative background
1665  // "atmvars"
1666  Vector ppath_p, ppath_t, ppath_wind_u, ppath_wind_v, ppath_wind_w;
1667  Matrix ppath_vmr;
1668  // "rtvars"
1669  Vector total_tau;
1670  Matrix ppath_emission, ppath_tau;
1671  Tensor3 ppath_abs_scalar, iy_trans_new;
1672  //
1673  const Index np = ppath.np;
1674  //
1675  if( np > 1 )
1676  {
1677  // Get pressure, temperature and VMRs
1678  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
1679  ppath_wind_u, ppath_wind_v, ppath_wind_w,
1680  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
1681  wind_u_field, wind_v_field, wind_w_field );
1682 
1683  // Get emission, absorption and optical thickness for each step
1684  get_ppath_rtvars( ws, ppath_abs_scalar, ppath_tau, total_tau,
1685  ppath_emission, abs_scalar_gas_agenda, emission_agenda,
1686  ppath, ppath_p, ppath_t, ppath_vmr, ppath_wind_u,
1687  ppath_wind_v, ppath_wind_w, f_grid, atmosphere_dim, 1 );
1688  }
1689 
1690  // Radiative background
1691  //
1692  Matrix dummy1, dummy3;
1693  Index dummy2;
1694  ArrayOfTensor3 dummy4;
1695  Tensor3 dummy5;
1696  //
1697  get_iy_of_background( ws, iy, dummy1, dummy2, dummy3, dummy4, dummy5,
1698  0, ppath, atmosphere_dim, p_grid, lat_grid, lon_grid,
1699  t_field, z_field, vmr_field, r_geoid, z_surface,
1700  cloudbox_on, stokes_dim, f_grid,
1701  iy_clearsky_basic_agenda, iy_space_agenda,
1702  surface_prop_agenda, iy_cloudbox_agenda, verbosity );
1703 
1704  // Do RT calculations
1705  //
1706  if( np > 1 )
1707  {
1708  // Loop ppath steps
1709  for( Index ip=np-2; ip>=0; ip-- )
1710  {
1711  // Loop frequencies
1712  for( Index iv=0; iv<nf; iv++ )
1713  {
1714  const Numeric step_tr = exp( -ppath_tau(iv,ip) );
1715  //
1716  iy(iv,0) = iy(iv,0) * step_tr + 0.5 * (1-step_tr) *
1717  ( ppath_emission(iv,ip) + ppath_emission(iv,ip+1) );
1718  //
1719  for( Index is=1; is<stokes_dim; is++ )
1720  { iy(iv,is) *= step_tr; }
1721  }
1722  }
1723  }
1724 }
1725 
1726 
1727 
1728 
1729 
1730 /* Workspace method: Doxygen documentation will be auto-generated */
1731 void iyMC(
1732  Workspace& ws,
1733  Matrix& iy,
1734  Matrix& iy_error,
1735  Index& iy_error_type,
1736  Matrix& iy_aux,
1737  ArrayOfTensor3& diy_dx,
1738  const Index& iy_agenda_call1,
1739  const Tensor3& iy_transmission,
1740  const Vector& rte_pos,
1741  const Vector& rte_los,
1742  const Index& jacobian_do,
1743  const Index& atmosphere_dim,
1744  const Vector& p_grid,
1745  const Vector& lat_grid,
1746  const Vector& lon_grid,
1747  const Tensor3& z_field,
1748  const Tensor3& t_field,
1749  const Tensor4& vmr_field,
1750  const Matrix& r_geoid,
1751  const Matrix& z_surface,
1752  const Index& cloudbox_on,
1753  const ArrayOfIndex& cloudbox_limits,
1754  const Index& cloudbox_checked,
1755  const Index& stokes_dim,
1756  const Vector& f_grid,
1757  const ArrayOfSingleScatteringData& scat_data_raw,
1758  const Agenda& iy_space_agenda,
1759  const Agenda& surface_prop_agenda,
1760  const Agenda& abs_scalar_gas_agenda,
1761  const Agenda& opt_prop_gas_agenda,
1762  const Tensor4& pnd_field,
1763  const String& y_unit,
1764  const Numeric& mc_std_err,
1765  const Index& mc_max_time,
1766  const Index& mc_max_iter,
1767  const Verbosity& verbosity)
1768 {
1769  // Throw error if unsupported features are requested
1770  if( atmosphere_dim != 3 )
1771  throw runtime_error(
1772  "Only 3D atmospheres are allowed (atmosphere_dim must be 3)" );
1773  if( !cloudbox_on )
1774  throw runtime_error(
1775  "The cloudbox must be activated (cloudbox_on must be 1)" );
1776  if( jacobian_do )
1777  throw runtime_error(
1778  "This method does not provide any jacobians (jacobian_do must be 0)" );
1779  if( !iy_agenda_call1 )
1780  throw runtime_error(
1781  "Recursive usage not possible (iy_agenda_call1 must be 1)" );
1782  if( iy_transmission.ncols() )
1783  throw runtime_error( "*iy_transmission* must be empty" );
1784 
1785  // See initial comments of iyEmissionStandardClearsky
1786 
1787  // Size output variables
1788  //
1789  const Index nf = f_grid.nelem();
1790  //
1791  iy.resize( nf, stokes_dim );
1792  iy_error.resize( nf, stokes_dim );
1793  iy_error_type = 1;
1794 
1795  // To avoid warning about unused variables
1796  iy_aux.resize(0,0);
1797  diy_dx.resize(0);
1798 
1799  // Some MC variables are only local here
1800  Tensor3 mc_points;
1801  Index mc_iteration_count, mc_seed;
1802  //
1803  MCAntenna mc_antenna;
1804  mc_antenna.set_pencil_beam();
1805 
1806  // Pos and los must be matrices
1807  Matrix pos(1,3), los(1,2);
1808  //
1809  pos(0,joker) = rte_pos;
1810  los(0,joker) = rte_los;
1811 
1812  for( Index f_index=0; f_index<nf; f_index++ )
1813  {
1814  ArrayOfSingleScatteringData scat_data_mono;
1815 
1816  scat_data_monoCalc( scat_data_mono, scat_data_raw,
1817  f_grid, f_index, verbosity );
1818 
1819  // Seed reset for each loop. If not done, the errors
1820  // appear to be highly correlated.
1821  MCSetSeedFromTime( mc_seed, verbosity );
1822 
1823  Vector y, mc_error;
1824 
1825  MCGeneral( ws, y, mc_iteration_count, mc_error, mc_points, mc_antenna,
1826  f_grid, f_index, pos, los, stokes_dim, atmosphere_dim,
1827  iy_space_agenda, surface_prop_agenda, opt_prop_gas_agenda,
1828  abs_scalar_gas_agenda, p_grid, lat_grid, lon_grid, z_field,
1829  r_geoid, z_surface, t_field, vmr_field,
1830  cloudbox_on, cloudbox_limits,
1831  pnd_field, scat_data_mono, 1, cloudbox_checked,
1832  mc_seed, y_unit, mc_std_err, mc_max_time, mc_max_iter,
1833  verbosity);
1834  // GH 2011-06-17, mc_z_field_is_1D);
1835 
1836  assert( y.nelem() == stokes_dim );
1837 
1838  // Data returned can not be in Tb
1839  if ( y_unit == "RJBT" )
1840  {
1841  const Numeric scfac = invrayjean( 1, f_grid[f_index] );
1842  y /= scfac;
1843  mc_error /= scfac;
1844  }
1845 
1846  iy(f_index,joker) = y;
1847  iy_error(f_index,joker) = mc_error;
1848  }
1849 }
1850 
1851 
1852 
1853 
1854 
1855 /* Workspace method: Doxygen documentation will be auto-generated */
1856 void yCalc(
1857  Workspace& ws,
1858  Vector& y,
1859  Vector& y_f,
1860  ArrayOfIndex& y_pol,
1861  Matrix& y_pos,
1862  Matrix& y_los,
1863  Vector& y_error,
1864  Vector& y_aux,
1865  Matrix& jacobian,
1866  const Index& basics_checked,
1867  const Index& atmosphere_dim,
1868  const Vector& p_grid,
1869  const Vector& lat_grid,
1870  const Vector& lon_grid,
1871  const Tensor3& t_field,
1872  const Tensor3& z_field,
1873  const Tensor4& vmr_field,
1874  const Index& cloudbox_on,
1875  const Index& cloudbox_checked,
1876  const Index& stokes_dim,
1877  const Vector& f_grid,
1878  const Matrix& sensor_pos,
1879  const Matrix& sensor_los,
1880  const Vector& mblock_za_grid,
1881  const Vector& mblock_aa_grid,
1882  const Index& antenna_dim,
1883  const Sparse& sensor_response,
1884  const Vector& sensor_response_f,
1885  const ArrayOfIndex& sensor_response_pol,
1886  const Vector& sensor_response_za,
1887  const Vector& sensor_response_aa,
1888  const Agenda& iy_clearsky_agenda,
1889  const String& y_unit,
1890  const Agenda& jacobian_agenda,
1891  const Index& jacobian_do,
1892  const ArrayOfRetrievalQuantity& jacobian_quantities,
1893  const ArrayOfArrayOfIndex& jacobian_indices,
1894  const Verbosity& verbosity )
1895 {
1896  // Some sizes
1897  const Index nf = f_grid.nelem();
1898  const Index nza = mblock_za_grid.nelem();
1899  Index naa = mblock_aa_grid.nelem();
1900  if( antenna_dim == 1 )
1901  { naa = 1; }
1902  const Index n1y = sensor_response.nrows();
1903  const Index nmblock = sensor_pos.nrows();
1904  const Index niyb = nf * nza * naa * stokes_dim;
1905 
1906 
1907  //---------------------------------------------------------------------------
1908  // Input checks
1909  //---------------------------------------------------------------------------
1910 
1911  // Basics and cloudbox OK?
1912  //
1913  if( !basics_checked )
1914  throw runtime_error( "The atmosphere and basic control varaibles must be "
1915  "flagged to have passed a consistency check (basics_checked=1)." );
1916  if( !cloudbox_checked )
1917  throw runtime_error( "The cloudbox must be flagged to have passed a "
1918  "consistency check (cloudbox_checked=1)." );
1919 
1920  // Sensor position and LOS.
1921  //
1922  if( sensor_pos.ncols() != atmosphere_dim )
1923  throw runtime_error( "The number of columns of sensor_pos must be "
1924  "equal to the atmospheric dimensionality." );
1925  if( atmosphere_dim <= 2 && sensor_los.ncols() != 1 )
1926  throw runtime_error( "For 1D and 2D, sensor_los shall have one column." );
1927  if( atmosphere_dim == 3 && sensor_los.ncols() != 2 )
1928  throw runtime_error( "For 3D, sensor_los shall have two columns." );
1929  if( sensor_los.nrows() != nmblock )
1930  {
1931  ostringstream os;
1932  os << "The number of rows of sensor_pos and sensor_los must be "
1933  << "identical, but sensor_pos has " << nmblock << " rows,\n"
1934  << "while sensor_los has " << sensor_los.nrows() << " rows.";
1935  throw runtime_error( os.str() );
1936  }
1937  if( max( sensor_los(joker,0) ) > 180 )
1938  throw runtime_error(
1939  "First column of *sensor_los* is not allowed to have values above 180." );
1940  if( atmosphere_dim == 2 )
1941  {
1942  if( min( sensor_los(joker,0) ) < -180 )
1943  throw runtime_error( "For atmosphere_dim = 2, first column of "
1944  "*sensor_los* is not allowed to have values below -180." );
1945  }
1946  else
1947  {
1948  if( min( sensor_los(joker,0) ) < 0 )
1949  throw runtime_error( "For atmosphere_dim != 2, first column of "
1950  "*sensor_los* is not allowed to have values below 0." );
1951  }
1952  if( atmosphere_dim == 3 && max( sensor_los(joker,1) ) > 180 )
1953  throw runtime_error(
1954  "Second column of *sensor_los* is not allowed to have values above 180." );
1955 
1956  // Antenna
1957  //
1958  chk_if_in_range( "antenna_dim", antenna_dim, 1, 2 );
1959  //
1960  if( nza == 0 )
1961  throw runtime_error( "The measurement block zenith angle grid is empty." );
1962  chk_if_increasing( "mblock_za_grid", mblock_za_grid );
1963  //
1964  if( antenna_dim == 1 )
1965  {
1966  if( mblock_aa_grid.nelem() != 0 )
1967  throw runtime_error(
1968  "For antenna_dim = 1, the azimuthal angle grid must be empty." );
1969  }
1970  else
1971  {
1972  if( atmosphere_dim < 3 )
1973  throw runtime_error( "2D antennas (antenna_dim=2) can only be "
1974  "used with 3D atmospheres." );
1975  if( mblock_aa_grid.nelem() == 0 )
1976  throw runtime_error(
1977  "The measurement block azimuthal angle grid is empty." );
1978  chk_if_increasing( "mblock_aa_grid", mblock_aa_grid );
1979  }
1980 
1981  // Sensor
1982  //
1983  if( sensor_response.ncols() != niyb )
1984  {
1985  ostringstream os;
1986  os << "The *sensor_response* matrix does not have the right size,\n"
1987  << "either the method *sensor_responseInit* has not been run or some\n"
1988  << "of the other sensor response methods has not been correctly\n"
1989  << "configured.";
1990  throw runtime_error( os.str() );
1991  }
1992 
1993  // Sensor aux variables
1994  //
1995  if( n1y != sensor_response_f.nelem() || n1y != sensor_response_pol.nelem() ||
1996  n1y != sensor_response_za.nelem() || n1y != sensor_response_za.nelem() )
1997  {
1998  ostringstream os;
1999  os << "Sensor auxiliary variables do not have the correct size.\n"
2000  << "The following variables should all have same size:\n"
2001  << "length(y) for one block : " << n1y << "\n"
2002  << "sensor_response_f.nelem() : " << sensor_response_f.nelem()
2003  << "\nsensor_response_pol.nelem(): " << sensor_response_pol.nelem()
2004  << "\nsensor_response_za.nelem() : " << sensor_response_za.nelem()
2005  << "\nsensor_response_aa.nelem() : " << sensor_response_za.nelem()
2006  << "\n";
2007  throw runtime_error( os.str() );
2008  }
2009 
2010 
2011  //---------------------------------------------------------------------------
2012  // Allocations and resizing
2013  //---------------------------------------------------------------------------
2014 
2015  // Resize *y* and *y_XXX*
2016  //
2017  y.resize( nmblock*n1y );
2018  y_f.resize( nmblock*n1y );
2019  y_pol.resize( nmblock*n1y );
2020  y_pos.resize( nmblock*n1y, sensor_pos.ncols() );
2021  y_los.resize( nmblock*n1y, sensor_los.ncols() );
2022  y_error.resize( nmblock*n1y );
2023  y_aux.resize( nmblock*n1y );
2024  //
2025  y_error = 0;
2026  y_aux = 999;
2027 
2028  // Jacobian variables
2029  //
2030  Index j_analytical_do = 0;
2031  //
2032  if( jacobian_do )
2033  {
2034  jacobian.resize( nmblock*n1y,
2035  jacobian_indices[jacobian_indices.nelem()-1][1]+1 );
2036  jacobian = 0;
2037  //
2038  FOR_ANALYTICAL_JACOBIANS_DO(
2039  j_analytical_do = 1;
2040  )
2041  }
2042  else
2043  { jacobian.resize(0, 0); }
2044 
2045 
2046  //---------------------------------------------------------------------------
2047  // The calculations
2048  //---------------------------------------------------------------------------
2049 
2050  // We have to make a local copy of the Workspace and the agendas because
2051  // only non-reference types can be declared firstprivate in OpenMP
2052  Workspace l_ws (ws);
2053  Agenda l_jacobian_agenda (jacobian_agenda);
2054  Agenda l_iy_clearsky_agenda (iy_clearsky_agenda);
2055 
2056 /*#pragma omp parallel for \
2057  if(!arts_omp_in_parallel() && nmblock>1 && nmblock>=nza) \
2058  default(none) \
2059  firstprivate(l_ws, l_jacobian_agenda, l_iy_clearsky_agenda) \
2060  shared(j_analytical_do, sensor_los, mblock_za_grid, mblock_aa_grid, \
2061  vmr_field, t_field, lon_grid, lat_grid, p_grid, f_grid, \
2062  sensor_pos, joker, naa)*/
2063 #pragma omp parallel for \
2064  if(!arts_omp_in_parallel() && nmblock>1 && nmblock>=nza) \
2065  firstprivate(l_ws, l_jacobian_agenda, l_iy_clearsky_agenda)
2066  for( Index imblock=0; imblock<nmblock; imblock++ )
2067  {
2068  // Calculate monochromatic pencil beam data for 1 measurement block
2069  //
2070  Vector iyb, iyb_aux, iyb_error, yb(n1y);
2071  Index iy_error_type;
2072  ArrayOfMatrix diyb_dx;
2073  //
2074  iyb_calc( l_ws, iyb, iyb_error, iy_error_type, iyb_aux, diyb_dx,
2075  imblock, atmosphere_dim, p_grid, lat_grid, lon_grid, t_field,
2076  z_field, vmr_field, cloudbox_on, stokes_dim, f_grid, sensor_pos,
2077  sensor_los, mblock_za_grid, mblock_aa_grid, antenna_dim,
2078  l_iy_clearsky_agenda, y_unit, j_analytical_do,
2079  jacobian_quantities, jacobian_indices, verbosity );
2080 
2081 
2082  // Apply sensor response matrix on iyb, and put into y
2083  //
2084  const Range rowind = get_rowindex_for_mblock( sensor_response, imblock );
2085  const Index row0 = rowind.get_start();
2086  //
2087  mult( yb, sensor_response, iyb );
2088  //
2089  y[rowind] = yb; // *yb* also used below, as input to jacobian_agenda
2090 
2091  // Apply sensor response matrix on iyb_error, and put into y_error
2092  // (nothing to do if type is 0).
2093  if( iy_error_type == 1 )
2094  {
2095  for( Index i=0; i<n1y; i++ )
2096  {
2097  const Index ithis = row0+i;
2098  y_error[ithis] = 0;
2099  for( Index icol=0; icol<niyb; icol++ )
2100  {
2101  y_error[ithis] += pow(
2102  sensor_response(i,icol)*iyb_error[icol], (Numeric)2.0 );
2103  }
2104  y_error[ithis] = sqrt( y_error[ithis] );
2105  }
2106  }
2107  else if( iy_error_type == 2 )
2108  { mult( y_error[rowind], sensor_response, iyb_error ); }
2109 
2110  // Apply sensor response matrix on iyb_aux, and put into y_aux
2111  // (if filled)
2112  if( iyb_aux.nelem() )
2113  { mult( y_aux[rowind], sensor_response, iyb_aux ); }
2114 
2115  // Fill information variables
2116  //
2117  for( Index i=0; i<n1y; i++ )
2118  {
2119  y_f[row0+i] = sensor_response_f[i];
2120  y_pol[row0+i] = sensor_response_pol[i];
2121  y_pos(row0+i,joker) = sensor_pos(imblock,joker);
2122  y_los(row0+i,0) = sensor_los(imblock,0) + sensor_response_za[i];
2123  if( sensor_response_aa.nelem() )
2124  {
2125  y_los(row0+i,1) = sensor_los(imblock,0) + sensor_response_aa[i];
2126  }
2127  }
2128 
2129  // Apply sensor response matrix on diyb_dx, and put into jacobian
2130  // (that is, analytical jacobian part)
2131  //
2132  if( j_analytical_do )
2133  {
2134  FOR_ANALYTICAL_JACOBIANS_DO(
2135  mult( jacobian(rowind, Range(jacobian_indices[iq][0],
2136  jacobian_indices[iq][1]-jacobian_indices[iq][0]+1)),
2137  sensor_response, diyb_dx[iq] );
2138  )
2139  }
2140 
2141  // Rest of *jacobian*
2142  //
2143  if( jacobian_do )
2144  {
2145  jacobian_agendaExecute( l_ws, jacobian, imblock, iyb, yb,
2146  l_jacobian_agenda );
2147  }
2148  } // End mblock loop
2149 }
2150 
2151 
2152 
2153 
2154 
2155 /* Workspace method: Doxygen documentation will be auto-generated */
2156 void y_unitApply(
2157  Vector& y,
2158  Vector& y_error,
2159  Matrix& jacobian,
2160  const Vector& y_f,
2161  const ArrayOfIndex& y_pol,
2162  const String& y_unit,
2163  const Verbosity&)
2164 {
2165  if( y_unit == "1" )
2166  { throw runtime_error(
2167  "No need to use this method with *y_unit* = \"1\"." ); }
2168 
2169  if( max(y) > 1e-3 )
2170  {
2171  ostringstream os;
2172  os << "The spectrum vector *y* is required to have original radiance\n"
2173  << "unit, but this seems not to be the case. This as a value above\n"
2174  << "1e-3 is found in *y*.";
2175  throw runtime_error( os.str() );
2176  }
2177 
2178  // Are jacobian and y_error set?
2179  //
2180  const Index ny = y.nelem();
2181  //
2182  const bool do_j = jacobian.nrows() == ny;
2183  const bool do_e = y_error.nelem() == ny;
2184 
2185  // Some jacobian quantities can not be handled
2186  if( do_j && max(jacobian) > 1e-3 )
2187  {
2188  ostringstream os;
2189  os << "The method can not be used with jacobian quantities that are not\n"
2190  << "obtained through radiative transfer calculations. One example on\n"
2191  << "quantity that can not be handled is *jacobianAddPolyfit*.\n"
2192  << "The maximum value of *jacobian* indicates that one or several\n"
2193  << "such jacobian quantities are included.";
2194  throw runtime_error( os.str() );
2195  }
2196 
2197  // All conversions could be handled as Planck-BT, but before everything is
2198  // tested properly and it is clear which version that is most efficient, both
2199  // versions are kept
2200 
2201  // Planck-Tb
2202  //--------------------------------------------------------------------------
2203  if( y_unit == "PlanckBT" )
2204  {
2205  // Hard to use telescoping here as the data are sorted differently in y
2206  // and jacobian, than what is expected apply_y_unit. Copy to temporary
2207  // variables instead.
2208 
2209  // Handle the elements in "frequency chunks"
2210 
2211  Index i0 = 0; // Index of first element for present chunk
2212  //
2213  while( i0 < ny )
2214  {
2215  // Find number values for this chunk
2216  Index n = 1;
2217  //
2218  while( i0+n < ny && y_f[i0] == y_f[i0+n] )
2219  { n++; }
2220 
2221  Matrix yv(1,n);
2222  ArrayOfIndex i_pol(n);
2223  bool any_quv = false;
2224  //
2225  for( Index i=0; i<n; i++ )
2226  {
2227  const Index ix=i0+i;
2228  yv(0,i) = y[ix];
2229  i_pol[i] = y_pol[ix];
2230  if( i_pol[i] > 1 && i_pol[i] < 5 )
2231  { any_quv = true; }
2232  }
2233 
2234  // Index of elements to convert
2235  Range ii( i0, n );
2236 
2237  if( do_e || do_j )
2238  {
2239  if( any_quv && i_pol[0] != 1 )
2240  {
2241  ostringstream os;
2242  os << "The conversion to PlanckBT, of the Jacobian and "
2243  << "errors for Q, U and V, requires that I (first Stokes "
2244  << "element) is at hand and that the data are sorted in "
2245  << "such way that I comes first for each frequency.";
2246  throw runtime_error( os.str() );
2247  }
2248 
2249  // y_error
2250  if( do_e )
2251  {
2252  Tensor3 e(1,1,n);
2253  e(0,0,joker) = y_error[ii];
2254  apply_y_unit2( e, yv, y_unit, y_f[i0], i_pol );
2255  y_error[ii] = e(0,0,joker);
2256  }
2257 
2258  // Jacobian
2259  if( do_j )
2260  {
2261  Tensor3 J(jacobian.ncols(),1,n);
2262  J(joker,0,joker) = transpose( jacobian(ii,joker) );
2263  apply_y_unit2( J, yv, y_unit, y_f[i0], i_pol );
2264  jacobian(ii,joker) = transpose( J(joker,0,joker) );
2265  }
2266  }
2267 
2268  // y (must be done last)
2269  apply_y_unit( yv, y_unit, y_f[i0], i_pol );
2270  y[ii] = yv(0,joker);
2271 
2272  i0 += n;
2273  }
2274  }
2275 
2276 
2277  // Other conversions
2278  //--------------------------------------------------------------------------
2279  else
2280  {
2281  // Here we take each element of y separately.
2282 
2283  Matrix yv(1,1);
2284  ArrayOfIndex i_pol(1);
2285 
2286  for( Index i=0; i<ny; i++ )
2287  {
2288  yv(0,0) = y[i]; // To avoid repeated telescoping
2289  i_pol[0] = y_pol[i];
2290 
2291  // y_error
2292  if( do_e )
2293  {
2294  apply_y_unit2( MatrixView(VectorView( y_error[i] )), yv,
2295  y_unit, y_f[i], i_pol );
2296  }
2297 
2298  // Jacobian
2299  if( do_j )
2300  {
2301  apply_y_unit2( MatrixView( jacobian(i,joker) ), yv,
2302  y_unit, y_f[i], i_pol );
2303  }
2304 
2305  // y (must be done last)
2306  apply_y_unit( yv, y_unit, y_f[i], i_pol );
2307  y[i] = yv(0,0);
2308  }
2309  }
2310 }
2311 
2312 
2313 
2314 
2315 
2316 
2317 /* Workspace method: Doxygen documentation will be auto-generated */
2318 void yCalc2(
2319  Workspace& ws,
2320  Vector& y,
2321  Vector& y_f,
2322  ArrayOfIndex& y_pol,
2323  Matrix& y_pos,
2324  Matrix& y_los,
2325  Vector& y_error,
2326  Vector& y_aux,
2327  Matrix& jacobian,
2328  const Index& basics_checked,
2329  const Index& atmosphere_dim,
2330  const Vector& p_grid,
2331  const Vector& lat_grid,
2332  const Vector& lon_grid,
2333  const Tensor3& t_field,
2334  const Tensor3& z_field,
2335  const Tensor4& vmr_field,
2336  const Index& cloudbox_on,
2337  const Index& cloudbox_checked,
2338  const Index& stokes_dim,
2339  const Vector& f_grid,
2340  const Matrix& sensor_pos,
2341  const Matrix& sensor_los,
2342  const Vector& mblock_za_grid,
2343  const Vector& mblock_aa_grid,
2344  const Index& antenna_dim,
2345  const Agenda& sensor_response_agenda,
2346  const Agenda& iy_clearsky_agenda,
2347  const String& y_unit,
2348  const Agenda& jacobian_agenda,
2349  const Index& jacobian_do,
2350  const ArrayOfRetrievalQuantity& jacobian_quantities,
2351  const ArrayOfArrayOfIndex& jacobian_indices,
2352  const Verbosity& verbosity )
2353 {
2354  // Some sizes
2355  const Index nf = f_grid.nelem();
2356  const Index nza = mblock_za_grid.nelem();
2357  Index naa = mblock_aa_grid.nelem();
2358  if( antenna_dim == 1 )
2359  { naa = 1; }
2360  const Index nmblock = sensor_pos.nrows();
2361  const Index niyb = nf * nza * naa * stokes_dim;
2362 
2363 
2364  //---------------------------------------------------------------------------
2365  // Input checks
2366  //---------------------------------------------------------------------------
2367 
2368  // Basics and cloudbox OK?
2369  //
2370  if( !basics_checked )
2371  throw runtime_error( "The atmosphere and basic control varaibles must be "
2372  "flagged to have passed a consistency check (basics_checked=1)." );
2373  if( !cloudbox_checked )
2374  throw runtime_error( "The cloudbox must be flagged to have passed a "
2375  "consistency check (cloudbox_checked=1)." );
2376 
2377  // Sensor position and LOS.
2378  //
2379  if( sensor_pos.ncols() != atmosphere_dim )
2380  throw runtime_error( "The number of columns of sensor_pos must be "
2381  "equal to the atmospheric dimensionality." );
2382  if( atmosphere_dim <= 2 && sensor_los.ncols() != 1 )
2383  throw runtime_error( "For 1D and 2D, sensor_los shall have one column." );
2384  if( atmosphere_dim == 3 && sensor_los.ncols() != 2 )
2385  throw runtime_error( "For 3D, sensor_los shall have two columns." );
2386  if( sensor_los.nrows() != nmblock )
2387  {
2388  ostringstream os;
2389  os << "The number of rows of sensor_pos and sensor_los must be "
2390  << "identical, but sensor_pos has " << nmblock << " rows,\n"
2391  << "while sensor_los has " << sensor_los.nrows() << " rows.";
2392  throw runtime_error( os.str() );
2393  }
2394  if( max( sensor_los(joker,0) ) > 180 )
2395  throw runtime_error(
2396  "First column of *sensor_los* is not allowed to have values above 180." );
2397  if( atmosphere_dim == 2 )
2398  {
2399  if( min( sensor_los(joker,0) ) < -180 )
2400  throw runtime_error( "For atmosphere_dim = 2, first column of "
2401  "*sensor_los* is not allowed to have values below -180." );
2402  }
2403  else
2404  {
2405  if( min( sensor_los(joker,0) ) < 0 )
2406  throw runtime_error( "For atmosphere_dim != 2, first column of "
2407  "*sensor_los* is not allowed to have values below 0." );
2408  }
2409  if( atmosphere_dim == 3 && max( sensor_los(joker,1) ) > 180 )
2410  throw runtime_error(
2411  "Second column of *sensor_los* is not allowed to have values above 180." );
2412 
2413  // Antenna
2414  //
2415  chk_if_in_range( "antenna_dim", antenna_dim, 1, 2 );
2416  //
2417  if( nza == 0 )
2418  throw runtime_error( "The measurement block zenith angle grid is empty." );
2419  chk_if_increasing( "mblock_za_grid", mblock_za_grid );
2420  //
2421  if( antenna_dim == 1 )
2422  {
2423  if( mblock_aa_grid.nelem() != 0 )
2424  throw runtime_error(
2425  "For antenna_dim = 1, the azimuthal angle grid must be empty." );
2426  }
2427  else
2428  {
2429  if( atmosphere_dim < 3 )
2430  throw runtime_error( "2D antennas (antenna_dim=2) can only be "
2431  "used with 3D atmospheres." );
2432  if( mblock_aa_grid.nelem() == 0 )
2433  throw runtime_error(
2434  "The measurement block azimuthal angle grid is empty." );
2435  chk_if_increasing( "mblock_aa_grid", mblock_aa_grid );
2436  }
2437 
2438 
2439  //---------------------------------------------------------------------------
2440  // Allocations
2441  //---------------------------------------------------------------------------
2442 
2443  // Temporary storage for *y* etc.
2444  //
2445  ArrayOfVector ya(nmblock), yf(nmblock), yerror(nmblock), yaux(nmblock);
2446  ArrayOfMatrix ypos(nmblock), ylos(nmblock), ja(nmblock);
2447  ArrayOfArrayOfIndex ypol(nmblock);
2448 
2449  // Jacobian variables
2450  //
2451  Index j_analytical_do = 0;
2452  //
2453  if( jacobian_do )
2454  {
2455  FOR_ANALYTICAL_JACOBIANS_DO(
2456  j_analytical_do = 1;
2457  )
2458  }
2459 
2460 
2461  //---------------------------------------------------------------------------
2462  // The calculations
2463  //---------------------------------------------------------------------------
2464 
2465  // We have to make a local copy of the Workspace and the agendas because
2466  // only non-reference types can be declared firstprivate in OpenMP
2467  Workspace l_ws (ws);
2468  Agenda l_jacobian_agenda (jacobian_agenda);
2469  Agenda l_iy_clearsky_agenda (iy_clearsky_agenda);
2470  Agenda l_sensor_response_agenda (sensor_response_agenda);
2471 
2472 /*#pragma omp parallel for \
2473  if(!arts_omp_in_parallel() && nmblock>1 && nmblock>=nza) \
2474  default(none) \
2475  firstprivate(l_ws, l_jacobian_agenda, l_iy_clearsky_agenda, \
2476  l_sensor_response_agenda) \
2477  shared(j_analytical_do, sensor_los, mblock_za_grid, mblock_aa_grid, \
2478  vmr_field, t_field, lon_grid, lat_grid, p_grid, f_grid, \
2479  sensor_pos, joker, naa)*/
2480 #pragma omp parallel for \
2481  if(!arts_omp_in_parallel() && nmblock>1 && nmblock>=nza) \
2482  firstprivate(l_ws, l_jacobian_agenda, l_iy_clearsky_agenda)
2483  for( Index imblock=0; imblock<nmblock; imblock++ )
2484  {
2485  // Calculate monochromatic pencil beam data for 1 measurement block
2486  //
2487  Vector iyb, iyb_aux, iyb_error;
2488  Index iy_error_type;
2489  ArrayOfMatrix diyb_dx;
2490  //
2491  iyb_calc( l_ws, iyb, iyb_error, iy_error_type, iyb_aux, diyb_dx,
2492  imblock, atmosphere_dim, p_grid, lat_grid, lon_grid, t_field,
2493  z_field, vmr_field, cloudbox_on, stokes_dim, f_grid, sensor_pos,
2494  sensor_los, mblock_za_grid, mblock_aa_grid, antenna_dim,
2495  l_iy_clearsky_agenda, y_unit, j_analytical_do,
2496  jacobian_quantities, jacobian_indices, verbosity );
2497 
2498  // Obtain sensor response data
2499  //
2500  Sparse sensor_response;
2501  Vector sensor_response_f;
2502  ArrayOfIndex sensor_response_pol;
2503  Vector sensor_response_za;
2504  Vector sensor_response_aa;
2505  //
2506  sensor_response_agendaExecute( l_ws, sensor_response, sensor_response_f,
2507  sensor_response_pol, sensor_response_za,
2508  sensor_response_aa,
2509  imblock, l_sensor_response_agenda );
2510 
2511  if( sensor_response.ncols() != niyb )
2512  {
2513  ostringstream os;
2514  os << "The number of columns in *sensor_response* and the length\n"
2515  << "of *iyb* do not match.";
2516  throw runtime_error( os.str() );
2517  }
2518  //
2519  const Index yl = sensor_response.nrows(); // Length of this part of y
2520  //
2521  if( yl!=sensor_response_f.nelem() || yl!=sensor_response_pol.nelem() ||
2522  yl!= sensor_response_za.nelem() || yl!=sensor_response_za.nelem() )
2523  {
2524  ostringstream os;
2525  os << "Sensor auxiliary variables do not have the correct size.\n"
2526  << "The following variables should all have same size:\n"
2527  << "length(y) for one block : " << yl << "\n"
2528  << "sensor_response_f.nelem() : " << sensor_response_f.nelem()
2529  << "\nsensor_response_pol.nelem(): " << sensor_response_pol.nelem()
2530  << "\nsensor_response_za.nelem() : " << sensor_response_za.nelem()
2531  << "\nsensor_response_aa.nelem() : " << sensor_response_za.nelem()
2532  << "\n";
2533  throw runtime_error( os.str() );
2534  }
2535 
2536  // Set sizes of ya and associated variables
2537  //
2538  ya[imblock].resize( yl );
2539  yf[imblock].resize( yl );
2540  ypol[imblock].resize( yl );
2541  ypos[imblock].resize( yl, sensor_pos.ncols() );
2542  ylos[imblock].resize( yl, sensor_los.ncols() );
2543  yerror[imblock].resize( yl );
2544  yaux[imblock].resize( yl );
2545  if( jacobian_do )
2546  { ja[imblock].resize( yl,
2547  jacobian_indices[jacobian_indices.nelem()-1][1]+1 ); }
2548 
2549  // Apply sensor response matrix on iyb, and put into ya
2550  //
2551  mult( ya[imblock], sensor_response, iyb );
2552 
2553  // Apply sensor response matrix on iyb_error, and put into yerror
2554  //
2555  if( iy_error_type == 0 )
2556  { yerror[imblock] = 0; }
2557  else if( iy_error_type == 1 )
2558  {
2559  for( Index i=0; i<yl; i++ )
2560  {
2561  yerror[imblock][i] = 0;
2562  for( Index icol=0; icol<niyb; icol++ )
2563  {
2564  yerror[imblock][i] += pow(
2565  sensor_response(i,icol)*iyb_error[icol], (Numeric)2.0 );
2566  }
2567  yerror[imblock][i] = sqrt( yerror[imblock][i] );
2568  }
2569  }
2570  else if( iy_error_type == 2 )
2571  { mult( yerror[imblock], sensor_response, iyb_error ); }
2572 
2573  // Apply sensor response matrix on iyb_aux, and put into y_aux
2574  // (if filled)
2575  if( iyb_aux.nelem() )
2576  { mult( yaux[imblock], sensor_response, iyb_aux ); }
2577  else
2578  { yaux[imblock] = 999; }
2579 
2580  // Fill information variables
2581  //
2582  for( Index i=0; i<yl; i++ )
2583  {
2584  yf[imblock][i] = sensor_response_f[i];
2585  ypol[imblock][i] = sensor_response_pol[i];
2586  ypos[imblock](i,joker) = sensor_pos(imblock,joker);
2587  ylos[imblock](i,0) = sensor_los(imblock,0)+sensor_response_za[i];
2588  if( sensor_response_aa.nelem() )
2589  { ylos[imblock](i,1) = sensor_los(imblock,0)+sensor_response_aa[i];}
2590  }
2591 
2592  // Apply sensor response matrix on diyb_dx, and put into jacobian
2593  // (that is, analytical jacobian part)
2594  //
2595  if( j_analytical_do )
2596  {
2597  FOR_ANALYTICAL_JACOBIANS_DO(
2598  mult( ja[imblock](joker, Range(jacobian_indices[iq][0],
2599  jacobian_indices[iq][1]-jacobian_indices[iq][0]+1)),
2600  sensor_response, diyb_dx[iq] );
2601  )
2602  }
2603 
2604  // Rest of *jacobian*
2605  //
2606  if( jacobian_do )
2607  {
2608  jacobian_agendaExecute( l_ws, ja[imblock], imblock, iyb, ya[imblock],
2609  l_jacobian_agenda );
2610  }
2611  } // End mblock loop
2612 
2613 
2614  // Combine data from the mblocks to form output variables
2615  //
2616  Index i0 =0 , n = 0;
2617  //
2618  for( Index m=0; m<nmblock; m++ )
2619  { n += ya[m].nelem(); }
2620  //
2621  y.resize( n );
2622  y_f.resize( n );
2623  y_pol.resize( n );
2624  y_pos.resize( n, sensor_pos.ncols() );
2625  y_los.resize( n, sensor_los.ncols() );
2626  y_error.resize( n );
2627  y_aux.resize( n );
2628  if( jacobian_do )
2629  { jacobian.resize( n, jacobian_indices[jacobian_indices.nelem()-1][1]+1 );}
2630  //
2631  for( Index m=0; m<nmblock; m++ )
2632  {
2633  for( Index i=0; i<ya[m].nelem(); i++ )
2634  {
2635  const Index r = i0 + i;
2636  y[r] = ya[m][i];
2637  y_f[r] = yf[m][i];
2638  y_pol[r] = ypol[m][i];
2639  y_pos(r,joker) = ypos[m](i,joker);
2640  y_los(r,joker) = ylos[m](i,joker);
2641  y_error[r] = yerror[m][i];
2642  y_aux[r] = yaux[m][i];
2643  if( jacobian_do )
2644  { jacobian(r,joker) = ja[m](i,joker); }
2645  }
2646  i0 += ya[m].nelem();
2647  }
2648 }
Matrix
The Matrix class.
Definition: matpackI.h:767
ppath_calc
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const bool &outside_cloudbox, const Verbosity &verbosity)
ppath_calc
Definition: ppath.cc:6171
ABSSPECIES_MAINTAG
const String ABSSPECIES_MAINTAG
apply_y_unit2
void apply_y_unit2(Tensor3View J, ConstMatrixView iy, const String &y_unit, ConstVectorView f_grid, const ArrayOfIndex &i_pol)
apply_y_unit2
Definition: rte.cc:184
iy_clearsky_agendaExecute
void iy_clearsky_agendaExecute(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Index cloudbox_on, const Index jacobian_do, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Agenda &input_agenda)
Definition: auto_md.cc:9483
iyEmissionStandardClearskyBasic
void iyEmissionStandardClearskyBasic(Workspace &ws, Matrix &iy, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const Agenda &ppath_step_agenda, const Agenda &emission_agenda, const Agenda &abs_scalar_gas_agenda, const Agenda &iy_clearsky_basic_agenda, const Agenda &iy_space_agenda, const Agenda &surface_prop_agenda, const Agenda &iy_cloudbox_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandardClearskyBasic.
Definition: m_rte.cc:1612
MatrixView
The MatrixView class.
Definition: matpackI.h:668
iy_transmission_mult
void iy_transmission_mult(Tensor3 &iy_trans_new, ConstTensor3View iy_transmission, ConstTensor3View trans)
iy_transmission_mult
Definition: rte.cc:929
auto_md.h
gridpos_copy
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Definition: interpolation.cc:482
Sparse::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:59
array_species_tag_from_string
void array_species_tag_from_string(ArrayOfSpeciesTag &tags, const String &names)
Converts a String to ArrayOfSpeciesTag.
Definition: abs_species_tags.cc:509
RetrievalQuantity::Grids
const ArrayOfVector & Grids() const
Grids.
Definition: jacobian.h:102
chk_if_increasing
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:126
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:167
Tensor3
The Tensor3 class.
Definition: matpackIII.h:340
yCalc2
void yCalc2(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Vector &y_error, Vector &y_aux, Matrix &jacobian, const Index &basics_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Agenda &sensor_response_agenda, const Agenda &iy_clearsky_agenda, const String &y_unit, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: yCalc2.
Definition: m_rte.cc:2318
sign
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:473
p2gridpos
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
p2gridpos
Definition: special_interp.cc:810
joker
const Joker joker
get_ppath_rtvars
void get_ppath_rtvars(Workspace &ws, Tensor3 &ppath_abs_scalar, Matrix &ppath_tau, Vector &total_tau, Matrix &ppath_emission, const Agenda &abs_scalar_gas_agenda, const Agenda &emission_agenda, const Ppath &ppath, ConstVectorView ppath_p, ConstVectorView ppath_t, ConstMatrixView ppath_vmr, ConstVectorView ppath_wind_u, ConstVectorView ppath_wind_v, ConstVectorView ppath_wind_w, ConstVectorView f_grid, const Index &atmosphere_dim, const Index &emission_do)
get_ppath_rtvars
Definition: rte.cc:531
iyb_calc
void iyb_calc(Workspace &ws, Vector &iyb, Vector &iyb_error, Index &iy_error_type, Vector &iyb_aux, ArrayOfMatrix &diyb_dx, const Index &imblock, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, const Index &cloudbox_on, const Index &stokes_dim, ConstVectorView f_grid, ConstMatrixView sensor_pos, ConstMatrixView sensor_los, ConstVectorView mblock_za_grid, ConstVectorView mblock_aa_grid, const Index &antenna_dim, const Agenda &iy_clearsky_agenda, const String &y_unit, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
iyb_calc
Definition: m_rte.cc:588
get_ppath_cloudrtvars
void get_ppath_cloudrtvars(Workspace &ws, Tensor3 &ppath_asp_abs_vec, Tensor4 &ppath_asp_ext_mat, Tensor3 &ppath_pnd_abs_vec, Tensor4 &ppath_pnd_ext_mat, Tensor4 &ppath_transmission, Tensor3 &total_transmission, Matrix &ppath_emission, Array< ArrayOfSingleScatteringData > &scat_data, const Agenda &abs_scalar_gas_agenda, const Agenda &emission_agenda, const Agenda &opt_prop_gas_agenda, const Ppath &ppath, ConstVectorView ppath_p, ConstVectorView ppath_t, ConstMatrixView ppath_vmr, ConstVectorView ppath_wind_u, ConstVectorView ppath_wind_v, ConstVectorView ppath_wind_w, ConstMatrixView ppath_pnd, const Index &use_mean_scat_data, const ArrayOfSingleScatteringData &scat_data_raw, const Index &stokes_dim, ConstVectorView f_grid, const Index &atmosphere_dim, const Index &emission_do, const Verbosity &verbosity)
get_ppath_cloudrtvars
Definition: rte.cc:672
get_ppath_atmvars
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, Matrix &ppath_vmr, Vector &ppath_wind_u, Vector &ppath_wind_v, Vector &ppath_wind_w, const Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstTensor3View wind_u_field, ConstTensor3View wind_v_field, ConstTensor3View wind_w_field)
get_ppath_atmvars
Definition: rte.cc:383
Ppath::background
String background
Definition: ppath.h:70
Sparse
The Sparse class.
Definition: matpackII.h:55
vmrunitscf
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
vmrunitscf
Definition: jacobian.cc:561
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: m_rte.cc:64
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:771
yCalc
void yCalc(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Vector &y_error, Vector &y_aux, Matrix &jacobian, const Index &basics_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Vector &sensor_response_za, const Vector &sensor_response_aa, const Agenda &iy_clearsky_agenda, const String &y_unit, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: yCalc.
Definition: m_rte.cc:1856
exit_or_rethrow
void exit_or_rethrow(const String &m, ArtsOut &out)
Exit ARTS or re-throw error.
Definition: arts.cc:98
iy_transmission_for_scalar_tau
void iy_transmission_for_scalar_tau(Tensor3 &iy_transmission, const Index &stokes_dim, ConstVectorView tau)
iy_transmission_for_scalar_tau
Definition: rte.cc:885
iy_cloudbox_agendaExecute
void iy_cloudbox_agendaExecute(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const GridPos &rte_gp_p, const GridPos &rte_gp_lat, const GridPos &rte_gp_lon, const Agenda &input_agenda)
Definition: auto_md.cc:9656
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:796
add_extrap
void add_extrap(ArrayOfGridPos &gp)
Definition: m_rte.cc:110
MCAntenna
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:56
surface_prop_agendaExecute
void surface_prop_agendaExecute(Workspace &ws, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &rte_pos, const Vector &rte_los, const GridPos &rte_gp_p, const GridPos &rte_gp_lat, const GridPos &rte_gp_lon, const Agenda &input_agenda)
Definition: auto_md.cc:10590
iyEmissionStandardClearsky
void iyEmissionStandardClearsky(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Agenda &ppath_step_agenda, const Agenda &emission_agenda, const Agenda &abs_scalar_gas_agenda, const Agenda &iy_clearsky_agenda, const Agenda &iy_space_agenda, const Agenda &surface_prop_agenda, const Agenda &iy_cloudbox_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandardClearsky.
Definition: m_rte.cc:1240
montecarlo.h
Tensor4
The Tensor4 class.
Definition: matpackIV.h:375
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
Ppath::gp_p
ArrayOfGridPos gp_p
Definition: ppath.h:66
CREATE_OUT1
#define CREATE_OUT1
Definition: messages.h:206
Agenda
The Agenda class.
Definition: agenda_class.h:44
apply_y_unit
void apply_y_unit(MatrixView iy, const String &y_unit, ConstVectorView f_grid, const ArrayOfIndex &i_pol)
apply_y_unit
Definition: rte.cc:82
get_rowindex_for_mblock
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &imblock)
get_rowindex_for_mblock
Definition: rte.cc:861
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:151
iy_clearsky_basic_agendaExecute
void iy_clearsky_basic_agendaExecute(Workspace &ws, Matrix &iy, const Vector &rte_pos, const Vector &rte_los, const Index cloudbox_on, const Agenda &input_agenda)
Definition: auto_md.cc:9589
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:149
invrayjean
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
Definition: physics_funcs.cc:170
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:234
Array< GridPos >
TEMPERATURE_MAINTAG
const String TEMPERATURE_MAINTAG
mult
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1607
y_unitApply
void y_unitApply(Vector &y, Vector &y_error, Matrix &jacobian, const Vector &y_f, const ArrayOfIndex &y_pol, const String &y_unit, const Verbosity &)
WORKSPACE METHOD: y_unitApply.
Definition: m_rte.cc:2156
diy_from_path_to_rgrids
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Definition: m_rte.cc:127
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:208
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:205
Ppath::gp_lat
ArrayOfGridPos gp_lat
Definition: ppath.h:67
messages.h
Declarations having to do with the four output streams.
MCSetSeedFromTime
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
Definition: m_montecarlo.cc:817
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:802
my_basic_string
The implementation for String, the ARTS string class.
Definition: mystring.h:62
abs
#define abs(x)
Definition: continua.cc:13094
Ppath::los
Matrix los
Definition: ppath.h:69
Ppath::gp_lon
ArrayOfGridPos gp_lon
Definition: ppath.h:68
Ppath::np
Index np
Definition: ppath.h:61
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:175
VectorView
The VectorView class.
Definition: matpackI.h:373
Agenda::name
String name() const
Agenda name.
Definition: agenda_class.cc:614
physics_funcs.h
ConstTensor4View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:78
MCGeneral
void MCGeneral(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &iy_space_agenda, const Agenda &surface_prop_agenda, const Agenda &opt_prop_gas_agenda, const Agenda &abs_scalar_gas_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Matrix &r_geoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_mono, const Index &basics_checked, const Index &cloudbox_checked, const Index &mc_seed, const String &y_unit, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
Definition: m_montecarlo.cc:185
jacobian.h
Declarations required for the calculation of jacobians.
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:50
Ppath::l_step
Vector l_step
Definition: ppath.h:65
ConstTensor4View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:66
map_daa
void map_daa(double &za, double &aa, const double &za0, const double &aa0, const double &aa_grid)
Definition: ppath.cc:967
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:60
math_funcs.h
transpose
ConstMatrixView transpose(ConstMatrixView m)
Const version of transpose.
Definition: matpackI.cc:1689
get_iy_of_background
void get_iy_of_background(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, ConstTensor3View iy_transmission, const Index &jacobian_do, const Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const Index &stokes_dim, ConstVectorView f_grid, const Agenda &iy_clearsky_agenda, const Agenda &iy_space_agenda, const Agenda &surface_prop_agenda, const Agenda &iy_cloudbox_agenda, const Verbosity &verbosity)
get_iy_of_background
Definition: m_rte.cc:360
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1549
iy_space_agendaExecute
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &rte_pos, const Vector &rte_los, const Agenda &input_agenda)
Definition: auto_md.cc:9744
scat_data_monoCalc
void scat_data_monoCalc(ArrayOfSingleScatteringData &scat_data_mono, const ArrayOfSingleScatteringData &scat_data_raw, const Vector &f_grid, const Index &f_index, const Verbosity &)
WORKSPACE METHOD: scat_data_monoCalc.
Definition: m_optproperties.cc:1035
max
#define max(a, b)
Definition: continua.cc:13097
iyBeerLambertStandardCloudbox
void iyBeerLambertStandardCloudbox(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const Agenda &ppath_step_agenda, const Agenda &abs_scalar_gas_agenda, const Agenda &iy_clearsky_agenda, const Tensor4 &pnd_field, const Index &use_mean_scat_data, const ArrayOfSingleScatteringData &scat_data_raw, const Agenda &opt_prop_gas_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: iyBeerLambertStandardCloudbox.
Definition: m_rte.cc:1097
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:591
jacobian_agendaExecute
void jacobian_agendaExecute(Workspace &ws, Matrix &jacobian, const Index imblock, const Vector &iyb, const Vector &yb, const Agenda &input_agenda)
Definition: auto_md.cc:9808
ConstTensor4View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:72
Range
The range class.
Definition: matpackI.h:148
GridPos
Structure to store a grid position.
Definition: interpolation.h:74
get_ppath_pnd
void get_ppath_pnd(Matrix &ppath_pnd, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, ConstTensor4View pnd_field)
get_ppath_pnd
Definition: rte.cc:479
ppath.h
Propagation path structure and functions.
surface_calc
void surface_calc(Matrix &iy, const Tensor3 &I, const Matrix &surface_los, const Tensor4 &surface_rmatrix, const Matrix &surface_emission)
surface_calc
Definition: rte.cc:1245
logic.h
Header file for logic.cc.
Sparse::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:53
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:157
min
#define min(a, b)
Definition: continua.cc:13096
MCAntenna::set_pencil_beam
void set_pencil_beam(void)
makes the antenna pattern a pencil beam
Definition: mc_antenna.cc:88
from_dpath_to_dx
void from_dpath_to_dx(MatrixView diy_dx, ConstMatrixView diy_dq, const Numeric &w)
diy_from_path_to_rgrids
Definition: m_rte.cc:96
rte.h
Declaration of functions in rte.cc.
Workspace
Workspace class.
Definition: workspace_ng.h:47
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:147
special_interp.h
Header file for special_interp.cc.
iyBeerLambertStandardClearsky
void iyBeerLambertStandardClearsky(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Agenda &ppath_step_agenda, const Agenda &abs_scalar_gas_agenda, const Agenda &iy_clearsky_agenda, const Agenda &iy_space_agenda, const Agenda &surface_prop_agenda, const Agenda &iy_cloudbox_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: iyBeerLambertStandardClearsky.
Definition: m_rte.cc:789
chk_contains
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition: check_input.h:193
Range::get_start
Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:196
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:95
RetrievalQuantity
Contains the data for one retrieval quantity.
Definition: jacobian.h:45
sensor_response_agendaExecute
void sensor_response_agendaExecute(Workspace &ws, Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Vector &sensor_response_za, Vector &sensor_response_aa, const Index imblock, const Agenda &input_agenda)
Definition: auto_md.cc:10438
check_input.h
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:2561
Vector
The Vector class.
Definition: matpackI.h:555
arts_omp.h
Header file for helper functions for OpenMP.
iy_transmission_mult_scalar_tau
void iy_transmission_mult_scalar_tau(Tensor3 &iy_trans_new, ConstTensor3View iy_transmission, ConstVectorView tau)
iy_transmission_mult_scalar_tau
Definition: rte.cc:981
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:300
Ppath::pos
Matrix pos
Definition: ppath.h:63
iyMC
void iyMC(Workspace &ws, Matrix &iy, Matrix &iy_error, Index &iy_error_type, Matrix &iy_aux, ArrayOfTensor3 &diy_dx, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Matrix &r_geoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &cloudbox_checked, const Index &stokes_dim, const Vector &f_grid, const ArrayOfSingleScatteringData &scat_data_raw, const Agenda &iy_space_agenda, const Agenda &surface_prop_agenda, const Agenda &abs_scalar_gas_agenda, const Agenda &opt_prop_gas_agenda, const Tensor4 &pnd_field, const String &y_unit, const Numeric &mc_std_err, const Index &mc_max_time, const Index &mc_max_iter, const Verbosity &verbosity)
WORKSPACE METHOD: iyMC.
Definition: m_rte.cc:1731
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:172
get_pointers_for_analytical_jacobians
void get_pointers_for_analytical_jacobians(ArrayOfIndex &abs_species_i, ArrayOfIndex &is_t, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: m_rte.cc:286
surfacetilt
double surfacetilt(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstMatrixView r_geoid, ConstMatrixView z_surface, const GridPos &gp_lat, const GridPos &gp_lon, ConstVectorView los)
surfacetilt
Definition: ppath.cc:1348
arts.h
The global header file for ARTS.