ARTS  2.2.66
m_rte.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012
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 "geodetic.h"
50 #include "jacobian.h"
51 #include "logic.h"
52 #include "math_funcs.h"
53 #include "messages.h"
54 #include "montecarlo.h"
55 #include "physics_funcs.h"
56 #include "ppath.h"
57 #include "rte.h"
58 #include "special_interp.h"
59 
60 extern const Numeric PI;
61 extern const Numeric SPEED_OF_LIGHT;
62 extern const String ABSSPECIES_MAINTAG;
63 extern const String TEMPERATURE_MAINTAG;
64 extern const String WIND_MAINTAG;
65 
66 
67 
68 /*===========================================================================
69  === Workspace methods
70  ===========================================================================*/
71 
72 /* Workspace method: Doxygen documentation will be auto-generated */
73 void iyApplyUnit(
74  Matrix& iy,
75  ArrayOfTensor4& iy_aux,
76  const Index& stokes_dim,
77  const Vector& f_grid,
78  const ArrayOfString& iy_aux_vars,
79  const String& iy_unit,
80  const Verbosity&)
81 {
82  if( iy_unit == "1" )
83  throw runtime_error( "No need to use this method with *iy_unit* = \"1\"." );
84 
85  if( max(iy(joker,0)) > 1e-3 )
86  {
87  ostringstream os;
88  os << "The spectrum matrix *iy* is required to have original radiance\n"
89  << "unit, but this seems not to be the case. This as a value above\n"
90  << "1e-3 is found in *iy*.";
91  throw runtime_error( os.str() );
92  }
93 
94  // Polarisation index variable
95  ArrayOfIndex i_pol(stokes_dim);
96  for( Index is=0; is<stokes_dim; is++ )
97  { i_pol[is] = is + 1; }
98 
99  apply_iy_unit( iy, iy_unit, f_grid, 1, i_pol );
100 
101  for( Index i=0; i<iy_aux_vars.nelem(); i++ )
102  {
103  if( iy_aux_vars[i] == "iy" || iy_aux_vars[i] == "Error" ||
104  iy_aux_vars[i] == "Error (uncorrelated)" )
105  {
106  if( iy_aux[i].nrows() > 1 )
107  throw runtime_error( "Data marked as \"iy\" or \"Error\" "
108  "have incorrect size." );
109  for( Index j=0; j<iy_aux[i].ncols(); j++ )
110  { apply_iy_unit( iy_aux[i](joker,joker,0,j), iy_unit, f_grid, 1,
111  i_pol ); }
112  }
113  }
114 }
115 
116 
117 
118 
119 
120 /* Workspace method: Doxygen documentation will be auto-generated */
121 void iyCalc(
122  Workspace& ws,
123  Matrix& iy,
124  ArrayOfTensor4& iy_aux,
125  Ppath& ppath,
126  const Index& atmfields_checked,
127  const Index& atmgeom_checked,
128  const ArrayOfString& iy_aux_vars,
129  const Vector& f_grid,
130  const Tensor3& t_field,
131  const Tensor3& z_field,
132  const Tensor4& vmr_field,
133  const Index& cloudbox_on,
134  const Index& cloudbox_checked,
135  const Vector& rte_pos,
136  const Vector& rte_los,
137  const Vector& rte_pos2,
138  const Agenda& iy_main_agenda,
139  const Verbosity& )
140 {
141  // Basics
142  //
143  if( atmfields_checked != 1 )
144  throw runtime_error( "The atmospheric fields must be flagged to have "
145  "passed a consistency check (atmfields_checked=1)." );
146  if( atmgeom_checked != 1 )
147  throw runtime_error( "The atmospheric geometry must be flagged to have "
148  "passed a consistency check (atmgeom_checked=1)." );
149  if( cloudbox_checked != 1 )
150  throw runtime_error( "The cloudbox must be flagged to have "
151  "passed a consistency check (cloudbox_checked=1)." );
152 
153 
154  // iy_transmission is just input and can be left empty for first call
155  Tensor3 iy_transmission(0,0,0);
156 
157  ArrayOfTensor3 diy_dx;
158 
159  iy_main_agendaExecute( ws, iy, iy_aux, ppath, diy_dx, 1, iy_transmission,
160  iy_aux_vars, cloudbox_on, 0, t_field,
161  z_field, vmr_field, f_grid, rte_pos, rte_los, rte_pos2,
162  iy_main_agenda );
163 }
164 
165 
166 
167 
168 
169 /* Workspace method: Doxygen documentation will be auto-generated */
170 void iyEmissionStandard(
171  Workspace& ws,
172  Matrix& iy,
173  ArrayOfTensor4& iy_aux,
174  Ppath& ppath,
175  ArrayOfTensor3& diy_dx,
176  const Index& stokes_dim,
177  const Vector& f_grid,
178  const Index& atmosphere_dim,
179  const Vector& p_grid,
180  const Tensor3& z_field,
181  const Tensor3& t_field,
182  const Tensor4& vmr_field,
183  const ArrayOfArrayOfSpeciesTag& abs_species,
184  const Tensor3& wind_u_field,
185  const Tensor3& wind_v_field,
186  const Tensor3& wind_w_field,
187  const Tensor3& mag_u_field,
188  const Tensor3& mag_v_field,
189  const Tensor3& mag_w_field,
190  const Index& cloudbox_on,
191  const String& iy_unit,
192  const ArrayOfString& iy_aux_vars,
193  const Index& jacobian_do,
194  const ArrayOfRetrievalQuantity& jacobian_quantities,
195  const ArrayOfArrayOfIndex& jacobian_indices,
196  const Agenda& ppath_agenda,
197  const Agenda& blackbody_radiation_agenda,
198  const Agenda& propmat_clearsky_agenda,
199  const Agenda& iy_main_agenda,
200  const Agenda& iy_space_agenda,
201  const Agenda& iy_surface_agenda,
202  const Agenda& iy_cloudbox_agenda,
203  const Index& iy_agenda_call1,
204  const Tensor3& iy_transmission,
205  const Vector& rte_pos,
206  const Vector& rte_los,
207  const Vector& rte_pos2,
208  const Numeric& rte_alonglos_v,
209  const Numeric& ppath_lraytrace,
210  const Verbosity& verbosity )
211 {
212  // Determine propagation path
213  //
214  ppath_agendaExecute( ws, ppath, ppath_lraytrace, rte_pos, rte_los, rte_pos2,
215  cloudbox_on, 0, t_field, z_field, vmr_field, f_grid,
216  ppath_agenda );
217 
218  // Some basic sizes
219  //
220  const Index nf = f_grid.nelem();
221  const Index ns = stokes_dim;
222  const Index np = ppath.np;
223  const Index nq = jacobian_quantities.nelem();
224 
225 
226  //### jacobian part #########################################################
227  // Initialise analytical jacobians (diy_dx and help variables)
228  //
229  Index j_analytical_do = 0;
230  ArrayOfTensor3 diy_dpath;
231  ArrayOfIndex jac_species_i(0), jac_is_t(0), jac_wind_i(0);
232  //
233  if( jacobian_do ) { FOR_ANALYTICAL_JACOBIANS_DO( j_analytical_do = 1; ) }
234  //
235  if( !j_analytical_do )
236  { diy_dx.resize( 0 ); }
237  else
238  {
239  diy_dpath.resize( nq );
240  jac_species_i.resize( nq );
241  jac_is_t.resize( nq );
242  jac_wind_i.resize( nq );
243  //
244  FOR_ANALYTICAL_JACOBIANS_DO(
245  diy_dpath[iq].resize( np, nf, ns );
246  diy_dpath[iq] = 0.0;
247  )
248  get_pointers_for_analytical_jacobians( jac_species_i, jac_is_t,
249  jac_wind_i, jacobian_quantities, abs_species );
250  if( iy_agenda_call1 )
251  {
252  diy_dx.resize( nq );
253  //
254  FOR_ANALYTICAL_JACOBIANS_DO( diy_dx[iq].resize(
255  jacobian_indices[iq][1]-jacobian_indices[iq][0]+1, nf, ns );
256  diy_dx[iq] = 0.0;
257  )
258  }
259  }
260  //###########################################################################
261 
262 
263  // Set up variable with index of species where we want abs_per_species.
264  // This variable can below be extended in iy_aux part.
265  //
266  ArrayOfIndex iaps(0);
267  //
268  for( Index i=0; i<jac_species_i.nelem(); i++ )
269  {
270  if( jac_species_i[i] >= 0 )
271  { iaps.push_back( jac_species_i[i] ); }
272  }
273 
274 
275  //=== iy_aux part ===========================================================
276  Index auxPressure = -1,
277  auxTemperature = -1,
278  auxAbsSum = -1,
279  auxBackground = -1,
280  auxIy = -1,
281  auxTrans = -1,
282  auxOptDepth = -1;
283  ArrayOfIndex auxAbsSpecies(0), auxAbsIsp(0);
284  ArrayOfIndex auxVmrSpecies(0), auxVmrIsp(0);
285  //
286  if( !iy_agenda_call1 )
287  { iy_aux.resize( 0 ); }
288  else
289  {
290  const Index naux = iy_aux_vars.nelem();
291  iy_aux.resize( naux );
292  //
293  for( Index i=0; i<naux; i++ )
294  {
295  if( iy_aux_vars[i] == "Pressure" )
296  { auxPressure = i; iy_aux[i].resize( 1, 1, 1, np ); }
297  else if( iy_aux_vars[i] == "Temperature" )
298  { auxTemperature = i; iy_aux[i].resize( 1, 1, 1, np ); }
299  else if( iy_aux_vars[i].substr(0,13) == "VMR, species " )
300  {
301  Index ispecies;
302  istringstream is(iy_aux_vars[i].substr(13,2));
303  is >> ispecies;
304  if( ispecies < 0 || ispecies>=abs_species.nelem() )
305  {
306  ostringstream os;
307  os << "You have selected VMR of species with index "
308  << ispecies << ".\nThis species does not exist!";
309  throw runtime_error( os.str() );
310  }
311  auxVmrSpecies.push_back(i);
312  auxVmrIsp.push_back(ispecies);
313  iy_aux[i].resize( 1, 1, 1, np );
314  }
315  else if( iy_aux_vars[i] == "Absorption, summed" )
316  { auxAbsSum = i; iy_aux[i].resize( nf, ns, ns, np ); }
317  else if( iy_aux_vars[i].substr(0,20) == "Absorption, species " )
318  {
319  Index ispecies;
320  istringstream is(iy_aux_vars[i].substr(20,2));
321  is >> ispecies;
322  if( ispecies < 0 || ispecies>=abs_species.nelem() )
323  {
324  ostringstream os;
325  os << "You have selected absorption species with index "
326  << ispecies << ".\nThis species does not exist!";
327  throw runtime_error( os.str() );
328  }
329  auxAbsSpecies.push_back(i);
330  const Index ihit = find_first( iaps, ispecies );
331  if( ihit >= 0 )
332  { auxAbsIsp.push_back( ihit ); }
333  else
334  {
335  iaps.push_back(ispecies);
336  auxAbsIsp.push_back( iaps.nelem()-1 );
337  }
338  iy_aux[i].resize( nf, ns, ns, np );
339  }
340  else if( iy_aux_vars[i] == "Radiative background" )
341  { auxBackground = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
342  else if( iy_aux_vars[i] == "iy" && auxIy < 0 )
343  { auxIy = i; iy_aux[i].resize( nf, ns, 1, np ); }
344  else if( iy_aux_vars[i] == "Transmission" && auxTrans < 0 )
345  { auxTrans = i; iy_aux[i].resize( nf, ns, ns, np ); }
346  else if( iy_aux_vars[i] == "Optical depth" )
347  { auxOptDepth = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
348  else if( iy_aux_vars[i].substr(0,14) == "Mass content, " )
349  { iy_aux[i].resize( 0, 0, 0, 0 ); }
350  else if( iy_aux_vars[i].substr(0,10) == "PND, type " )
351  { iy_aux[i].resize( 0, 0, 0, 0 ); }
352  else
353  {
354  ostringstream os;
355  os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
356  << "\"\nThis choice is not recognised.";
357  throw runtime_error( os.str() );
358  }
359  }
360  }
361  //===========================================================================
362 
363 
364  // Get atmospheric and attenuation quantities for each ppath point/step
365  //
366  // "atmvars"
367  Vector ppath_p, ppath_t;
368  Matrix ppath_vmr, ppath_wind, ppath_mag, ppath_f;
369  // Attenuation vars
370  Tensor4 ppath_abs;
371  Tensor5 abs_per_species;
372  Tensor4 trans_partial, trans_cumulat;
373  Matrix ppath_blackrad;
374  Vector scalar_tau;
375  ArrayOfArrayOfIndex extmat_case;
376  //
377  if( np > 1 )
378  {
379  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
380  ppath_wind, ppath_mag,
381  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
382  wind_u_field, wind_v_field, wind_w_field,
383  mag_u_field, mag_v_field, mag_w_field );
384  get_ppath_f( ppath_f, ppath, f_grid, atmosphere_dim,
385  rte_alonglos_v, ppath_wind );
386  get_ppath_abs( ws, ppath_abs, abs_per_species,
387  propmat_clearsky_agenda, ppath,
388  ppath_p, ppath_t, ppath_vmr, ppath_f,
389  ppath_mag, f_grid, stokes_dim, iaps );
390  get_ppath_trans( trans_partial, extmat_case, trans_cumulat,
391  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim );
392  get_ppath_blackrad( ws, ppath_blackrad, blackbody_radiation_agenda,
393  ppath, ppath_t, ppath_f );
394  }
395  else // For cases inside the cloudbox, or totally outside the atmosphere,
396  { // set zero optical thickness and unit transmission
397  scalar_tau.resize( nf );
398  scalar_tau = 0;
399  trans_cumulat.resize( nf, ns, ns, np );
400  for( Index iv=0; iv<nf; iv++ )
401  { id_mat( trans_cumulat(iv,joker,joker,np-1) ); }
402  }
403 
404  // iy_transmission
405  //
406  Tensor3 iy_trans_new;
407  //
408  if( iy_agenda_call1 )
409  { iy_trans_new = trans_cumulat(joker,joker,joker,np-1); }
410  else
411  { iy_transmission_mult( iy_trans_new, iy_transmission,
412  trans_cumulat(joker,joker,joker,np-1) ); }
413 
414  // Radiative background
415  //
416  get_iy_of_background( ws, iy, diy_dx,
417  iy_trans_new, jacobian_do, ppath, rte_pos2,
418  atmosphere_dim, t_field, z_field, vmr_field,
419  cloudbox_on, stokes_dim, f_grid, iy_main_agenda,
420  iy_space_agenda, iy_surface_agenda, iy_cloudbox_agenda,
421  verbosity );
422 
423 
424  //=== iy_aux part ===========================================================
425  // Fill parts of iy_aux that are defined even for np=1.
426  // Radiative background
427  if( auxBackground >= 0 )
428  { iy_aux[auxBackground](joker,0,0,0) = (Numeric)min( (Index)2,
429  ppath_what_background(ppath)-1); }
430  // Radiance
431  if( auxIy >= 0 )
432  { iy_aux[auxIy](joker,joker,0,np-1) = iy; }
433  // Transmission variables
434  if( auxTrans >= 0 )
435  {
436  if( np == 1 )
437  { for( Index iv=0; iv<nf; iv++ ) {
438  id_mat( iy_aux[auxTrans](iv,joker,joker,0) ); } }
439  else
440  { iy_aux[auxTrans] = trans_cumulat; }
441  }
442  if( auxOptDepth >= 0 )
443  { iy_aux[auxOptDepth](joker,0,0,0) = scalar_tau; }
444  //===========================================================================
445 
446 
447  // Do RT calculations
448  //
449  if( np > 1 )
450  {
451  //### jacobian part #####################################################
452  // Create container for the derivatives with respect to changes
453  // at each ppath point. And find matching abs_species-index and
454  // "temperature flag" (for analytical jacobians).
455  //
456  //
457  const Numeric dt = 0.1;
458  const Numeric dw = 5;
459  Tensor4 ppath_at2, ppath_awu, ppath_awv, ppath_aww;
460  Matrix ppath_bt2;
461  //
462  if( j_analytical_do )
463  {
464  // Determine if temperature is among the analytical jac. quantities.
465  // If yes, get emission and absorption for disturbed temperature
466  // Same for wind, but disturb only absorption
467  for( Index iq=0; iq<jac_is_t.nelem(); iq++ )
468  {
469  if( jac_is_t[iq] )
470  {
471  Tensor5 dummy_abs_per_species;
472  Vector t2 = ppath_t; t2 += dt;
473  get_ppath_abs( ws, ppath_at2, dummy_abs_per_species,
474  propmat_clearsky_agenda, ppath, ppath_p,
475  t2, ppath_vmr, ppath_f, ppath_mag, f_grid,
476  stokes_dim, ArrayOfIndex(0) );
477  get_ppath_blackrad( ws, ppath_bt2, blackbody_radiation_agenda,
478  ppath, t2, ppath_f );
479  }
480  else if( jac_wind_i[iq] )
481  {
482  if( jac_wind_i[iq] == 1 )
483  {
484  Tensor5 dummy_abs_per_species;
485  Matrix f2, w2 = ppath_wind; w2(0,joker) += dw;
486  get_ppath_f( f2, ppath, f_grid, atmosphere_dim,
487  rte_alonglos_v, w2 );
488  get_ppath_abs( ws, ppath_awu, dummy_abs_per_species,
489  propmat_clearsky_agenda, ppath, ppath_p,
490  ppath_t, ppath_vmr, f2, ppath_mag, f_grid,
491  stokes_dim, ArrayOfIndex(0) );
492  }
493  else if( jac_wind_i[iq] == 2 )
494  {
495  Tensor5 dummy_abs_per_species;
496  Matrix f2, w2 = ppath_wind; w2(1,joker) += dw;
497  get_ppath_f( f2, ppath, f_grid, atmosphere_dim,
498  rte_alonglos_v, w2 );
499  get_ppath_abs( ws, ppath_awv, dummy_abs_per_species,
500  propmat_clearsky_agenda, ppath, ppath_p,
501  ppath_t, ppath_vmr, f2, ppath_mag, f_grid,
502  stokes_dim, ArrayOfIndex(0) );
503  }
504  else if( jac_wind_i[iq] == 3 )
505  {
506  Tensor5 dummy_abs_per_species;
507  Matrix f2, w2 = ppath_wind; w2(2,joker) += dw;
508  get_ppath_f( f2, ppath, f_grid, atmosphere_dim,
509  rte_alonglos_v, w2 );
510  get_ppath_abs( ws, ppath_aww, dummy_abs_per_species,
511  propmat_clearsky_agenda, ppath, ppath_p,
512  ppath_t, ppath_vmr, f2, ppath_mag, f_grid,
513  stokes_dim, ArrayOfIndex(0) );
514  }
515  }
516  }
517  }
518  //#######################################################################
519 
520  //=== iy_aux part =======================================================
521  // iy_aux for point np-1:
522  // Pressure
523  if( auxPressure >= 0 )
524  { iy_aux[auxPressure](0,0,0,np-1) = ppath_p[np-1]; }
525  // Temperature
526  if( auxTemperature >= 0 )
527  { iy_aux[auxTemperature](0,0,0,np-1) = ppath_t[np-1]; }
528  // VMR
529  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
530  { iy_aux[auxVmrSpecies[j]](0,0,0,np-1) = ppath_vmr(auxVmrIsp[j],np-1); }
531  // Absorption
532  if( auxAbsSum >= 0 )
533  { for( Index iv=0; iv<nf; iv++ ) {
534  for( Index is1=0; is1<ns; is1++ ){
535  for( Index is2=0; is2<ns; is2++ ){
536  iy_aux[auxAbsSum](iv,is1,is2,np-1) =
537  ppath_abs(iv,is1,is2,np-1); } } } }
538  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
539  { for( Index iv=0; iv<nf; iv++ ) {
540  for( Index is1=0; is1<ns; is1++ ){
541  for( Index is2=0; is2<ns; is2++ ){
542  iy_aux[auxAbsSpecies[j]](iv,is1,is2,np-1) =
543  abs_per_species(auxAbsIsp[j],iv,is1,is2,np-1); } } } }
544  // Radiance for this point is handled above
545  //=======================================================================
546 
547 
548  // Loop ppath steps
549  for( Index ip=np-2; ip>=0; ip-- )
550  {
551  // Path step average of B: Bbar
552  //
553  Vector bbar(nf);
554  //
555  for( Index iv=0; iv<nf; iv++ )
556  { bbar[iv] = 0.5 * ( ppath_blackrad(iv,ip) +
557  ppath_blackrad(iv,ip+1) ); }
558 
559  //### jacobian part #################################################
560  if( j_analytical_do )
561  {
562  // Calculate (si-bi)
563  Matrix sibi(nf,ns);
564  for( Index iv=0; iv<nf; iv++ )
565  {
566  sibi(iv,0) = iy(iv,0) - bbar[iv];
567  for( Index is=1; is<ns; is++ )
568  { sibi(iv,is) = iy(iv,is); }
569  }
570 
571  // Loop quantities
572  //
573  Index iiaps = -1; // Index with respect to abs_per_species.
574  // Jacobian species stored first and in order
575  for( Index iq=0; iq<nq; iq++ )
576  {
577  if( jacobian_quantities[iq].Analytical() )
578  {
579  //- Absorbing species -----------------------------------
580  if( jac_species_i[iq] >= 0 )
581  {
582  // Index with respect to abs_per_species and ppath_vmr
583  iiaps += 1;
584  const Index isp = jac_species_i[iq];
585 
586  // Scaling factors to handle retrieval unit
587  Numeric unitscf1, unitscf2;
588  vmrunitscf( unitscf1,
589  jacobian_quantities[iq].Mode(),
590  ppath_vmr(isp,ip), ppath_p[ip],
591  ppath_t[ip] );
592  vmrunitscf( unitscf2,
593  jacobian_quantities[iq].Mode(),
594  ppath_vmr(isp,ip+1), ppath_p[ip+1],
595  ppath_t[ip+1] );
596 
597  for( Index iv=0; iv<nf; iv++ )
598  {
599  // Diagonal transmission matrix
600  if( extmat_case[ip][iv] == 1 )
601  {
602  const Numeric x = -0.5 * ppath.lstep[ip] *
603  trans_cumulat(iv,0,0,ip+1);
604  const Numeric y = x * sibi(iv,0);
605  // Stokes component 1
606  diy_dpath[iq](ip ,iv,0) += y * unitscf1 *
607  abs_per_species(iiaps,iv,0,0,ip );
608  diy_dpath[iq](ip+1,iv,0) += y * unitscf2 *
609  abs_per_species(iiaps,iv,0,0,ip+1);
610  // Higher stokes components
611  for( Index is=1; is<ns; is++ )
612  {
613  const Numeric z = x * iy(iv,is);
614  diy_dpath[iq](ip ,iv,is) += z * unitscf1*
615  abs_per_species(iiaps,iv,0,0,ip );
616  diy_dpath[iq](ip+1,iv,is) += z * unitscf2*
617  abs_per_species(iiaps,iv,0,0,ip+1);
618  }
619  }
620 
621  // General case
622  else
623  {
624  // Size of disturbance, a relative number
625  const Numeric dd = 1e-3;
626  // Disturb for ip
627  Matrix ext_mat(ns,ns), dtdx(ns,ns);
628  //
629  for( Index is1=0; is1<ns; is1++ ) {
630  for( Index is2=0; is2<ns; is2++ ) {
631  ext_mat(is1,is2) = 0.5 * ( dd *
632  abs_per_species(iiaps,iv,is1,is2,ip ) +
633  ppath_abs(iv,is1,is2,ip) +
634  ppath_abs(iv,is1,is2,ip+1) );
635  } }
636  ext2trans( dtdx, extmat_case[ip][iv],
637  ext_mat, ppath.lstep[ip] );
638  for( Index is1=0; is1<ns; is1++ ) {
639  for( Index is2=0; is2<ns; is2++ ) {
640  dtdx(is1,is2) = (unitscf1/dd) *
641  ( dtdx(is1,is2) -
642  trans_partial(iv,is1,is2,ip) );
643  } }
644  Vector x(ns), y(ns);
645  mult( x, dtdx, sibi(iv,joker) );
646  mult( y, trans_cumulat(iv,joker,joker,ip),
647  x );
648  diy_dpath[iq](ip,iv,joker) += y;
649  //
650  // Disturb for ip+1
651  for( Index is1=0; is1<ns; is1++ ) {
652  for( Index is2=0; is2<ns; is2++ ) {
653  ext_mat(is1,is2) = 0.5 * ( dd *
654  abs_per_species(iiaps,iv,is1,is2,ip+1) +
655  ppath_abs(iv,is1,is2,ip) +
656  ppath_abs(iv,is1,is2,ip+1) );
657  } }
658  ext2trans( dtdx, extmat_case[ip][iv],
659  ext_mat, ppath.lstep[ip] );
660  for( Index is1=0; is1<ns; is1++ ) {
661  for( Index is2=0; is2<ns; is2++ ) {
662  dtdx(is1,is2) = (unitscf2/dd) *
663  ( dtdx(is1,is2) -
664  trans_partial(iv,is1,is2,ip) );
665  } }
666  mult( x, dtdx, sibi(iv,joker) );
667  mult( y, trans_cumulat(iv,joker,joker,ip),
668  x );
669  diy_dpath[iq](ip+1,iv,joker) += y;
670  }
671  }
672  }
673 
674  //- Winds -----------------------------------------------
675  else if( jac_wind_i[iq] )
676  {
677  for( Index iv=0; iv<nf; iv++ )
678  {
679  // Create pointer to relevant disturbed
680  // absorption to use. V-component is first guess.
681  Tensor4* ppath_awx = &ppath_awv;
682  if( jac_wind_i[iq] == 1 )
683  { ppath_awx = &ppath_awu; }
684  else if( jac_wind_i[iq] == 3 )
685  { ppath_awx = &ppath_aww; }
686 
687  // Diagonal transmission matrix
688  if( extmat_case[ip][iv] == 1 )
689  {
690  const Numeric dkdx1 = (1/dw) * (
691  (*ppath_awx)(iv,0,0,ip ) -
692  ppath_abs(iv,0,0,ip ) );
693  const Numeric dkdx2 = (1/dw) * (
694  (*ppath_awx)(iv,0,0,ip+1) -
695  ppath_abs(iv,0,0,ip+1) );
696  const Numeric x = -0.5 * ppath.lstep[ip] *
697  trans_cumulat(iv,0,0,ip+1);
698  const Numeric y = x * sibi(iv,0);
699  // Stokes component 1
700  diy_dpath[iq](ip ,iv,0) += y * dkdx1;
701  diy_dpath[iq](ip+1,iv,0) += y * dkdx2;
702  // Higher stokes components
703  for( Index is=1; is<ns; is++ )
704  {
705  const Numeric z = x * iy(iv,is);
706  diy_dpath[iq](ip ,iv,is) += z * dkdx1;
707  diy_dpath[iq](ip+1,iv,is) += z * dkdx2;
708  }
709  }
710 
711  // General case
712  else
713  {
714  // Disturb for ip
715  Matrix ext_mat(ns,ns), dtdx(ns,ns);
716  //
717  for( Index is1=0; is1<ns; is1++ ) {
718  for( Index is2=0; is2<ns; is2++ ) {
719  ext_mat(is1,is2) = 0.5 * (
720  (*ppath_awx)(iv,is1,is2,ip ) +
721  ppath_abs(iv,is1,is2,ip+1) );
722  } }
723  ext2trans( dtdx, extmat_case[ip][iv],
724  ext_mat, ppath.lstep[ip] );
725  for( Index is1=0; is1<ns; is1++ ) {
726  for( Index is2=0; is2<ns; is2++ ) {
727  dtdx(is1,is2) = (1/dw) *
728  ( dtdx(is1,is2) -
729  trans_partial(iv,is1,is2,ip) );
730  } }
731  Vector x(ns), y(ns);
732  mult( x, dtdx, sibi(iv,joker) );
733  mult( y, trans_cumulat(iv,joker,joker,ip),
734  x );
735  diy_dpath[iq](ip,iv,joker) += y;
736  //
737  // Disturb for ip+1
738  for( Index is1=0; is1<ns; is1++ ) {
739  for( Index is2=0; is2<ns; is2++ ) {
740  ext_mat(is1,is2) = 0.5 * (
741  ppath_abs(iv,is1,is2,ip ) +
742  (*ppath_awx)(iv,is1,is2,ip+1) );
743  } }
744  ext2trans( dtdx, extmat_case[ip][iv],
745  ext_mat, ppath.lstep[ip] );
746  for( Index is1=0; is1<ns; is1++ ) {
747  for( Index is2=0; is2<ns; is2++ ) {
748  dtdx(is1,is2) = (1/dw) *
749  ( dtdx(is1,is2) -
750  trans_partial(iv,is1,is2,ip) );
751  } }
752  mult( x, dtdx, sibi(iv,joker) );
753  mult( y, trans_cumulat(iv,joker,joker,ip),
754  x );
755  diy_dpath[iq](ip+1,iv,joker) += y;
756  }
757  }
758  }
759 
760  //- Temperature -----------------------------------------
761  else if( jac_is_t[iq] )
762  {
763  for( Index iv=0; iv<nf; iv++ )
764  {
765  // Diagonal transmission matrix
766  if( extmat_case[ip][iv] == 1 )
767  {
768  const Numeric dkdt1 = 1/dt * (
769  ppath_at2(iv,0,0,ip ) -
770  ppath_abs(iv,0,0,ip ) );
771  const Numeric dkdt2 = 1/dt * (
772  ppath_at2(iv,0,0,ip+1) -
773  ppath_abs(iv,0,0,ip+1) );
774  const Numeric x = -0.5 * ppath.lstep[ip] *
775  trans_cumulat(iv,0,0,ip+1);
776  const Numeric y = x * sibi(iv,0);
777  // Stokes 1:
778  diy_dpath[iq](ip ,iv,0) += y * dkdt1;
779  diy_dpath[iq](ip+1,iv,0) += y * dkdt2;
780  // Higher Stokes
781  for( Index is=1; is<ns; is++ )
782  {
783  const Numeric z = x * iy(iv,is);
784  diy_dpath[iq](ip ,iv,is) += z * dkdt1;
785  diy_dpath[iq](ip+1,iv,is) += z * dkdt2;
786  }
787  //
788  // The terms associated with B-bar:
789  const Numeric v = 0.5 *
790  trans_cumulat(iv,0,0,ip) *
791  ( 1.0 - trans_partial(iv,0,0,ip));
792  diy_dpath[iq](ip ,iv,0) += v/dt * (
793  ppath_bt2(iv,ip) -
794  ppath_blackrad(iv,ip) );
795  diy_dpath[iq](ip+1,iv,0) += v/dt * (
796  ppath_bt2(iv,ip+1) -
797  ppath_blackrad(iv,ip+1) );
798  // Zero for higher Stokes
799  //
800  // The terms associated with Delta-s:
801  if( jacobian_quantities[iq].Subtag() ==
802  "HSE on" )
803  {
804  // Stokes 1:
805  const Numeric kbar = 0.5 * (
806  ppath_abs(iv,0,0,ip ) +
807  ppath_abs(iv,0,0,ip+1) );
808  diy_dpath[iq](ip ,iv,0) += y * kbar /
809  ppath_t[ip];
810  diy_dpath[iq](ip+1,iv,0) += y * kbar /
811  ppath_t[ip+1];
812  // Higher Stokes
813  for( Index is=1; is<ns; is++ )
814  {
815  const Numeric z = x * iy(iv,is);
816  diy_dpath[iq](ip ,iv,is) +=
817  z * kbar / ppath_t[ip];
818  diy_dpath[iq](ip+1,iv,is) +=
819  z * kbar / ppath_t[ip+1];
820  }
821  } //hse
822  }
823  // General case
824  else
825  {
826  // Disturb for ip
827  Matrix ext_mat(ns,ns), dtdx(ns,ns);
828  //
829  for( Index is1=0; is1<ns; is1++ ) {
830  for( Index is2=0; is2<ns; is2++ ) {
831  ext_mat(is1,is2) = 0.5 * (
832  ppath_at2(iv,is1,is2,ip ) +
833  ppath_abs(iv,is1,is2,ip+1) );
834  } }
835  ext2trans( dtdx, extmat_case[ip][iv],
836  ext_mat, ppath.lstep[ip] );
837  for( Index is1=0; is1<ns; is1++ ) {
838  for( Index is2=0; is2<ns; is2++ ) {
839  dtdx(is1,is2) = (1/dt) *
840  ( dtdx(is1,is2) -
841  trans_partial(iv,is1,is2,ip) );
842  } }
843  Vector x(ns), y(ns);
844  mult( x, dtdx, sibi(iv,joker) );
845  mult( y, trans_cumulat(iv,joker,joker,ip),
846  x );
847  diy_dpath[iq](ip,iv,joker) += y;
848  //
849  // Disturb for ip+1
850  for( Index is1=0; is1<ns; is1++ ) {
851  for( Index is2=0; is2<ns; is2++ ) {
852  ext_mat(is1,is2) = 0.5 * (
853  ppath_abs(iv,is1,is2,ip ) +
854  ppath_at2(iv,is1,is2,ip+1) );
855  } }
856  ext2trans( dtdx, extmat_case[ip][iv],
857  ext_mat, ppath.lstep[ip] );
858  for( Index is1=0; is1<ns; is1++ ) {
859  for( Index is2=0; is2<ns; is2++ ) {
860  dtdx(is1,is2) = (1/dt) *
861  ( dtdx(is1,is2) -
862  trans_partial(iv,is1,is2,ip) );
863  } }
864  mult( x, dtdx, sibi(iv,joker) );
865  mult( y, trans_cumulat(iv,joker,joker,ip),
866  x );
867  diy_dpath[iq](ip+1,iv,joker) += y;
868  //
869  // The terms associated with B-bar:
870  const Numeric v =
871  ( 1.0 - trans_partial(iv,0,0,ip));
872  Numeric dbdt = 0.5/dt * ( ppath_bt2(iv,ip) -
873  ppath_blackrad(iv,ip) );
874  x[0] = v * dbdt;
875  for( Index is=1; is<ns; is++ )
876  { x[is] = -trans_partial(iv,is,0,ip)*dbdt; }
877  mult( y, trans_cumulat(iv,joker,joker,ip),
878  x );
879  diy_dpath[iq](ip,iv,joker) += y;
880  // Some for ip+1
881  dbdt = 0.5/dt * ( ppath_bt2(iv,ip+1) -
882  ppath_blackrad(iv,ip+1) );
883  x[0] = v * dbdt;
884  for( Index is=1; is<ns; is++ )
885  { x[is] = -trans_partial(iv,is,0,ip)*dbdt; }
886  mult( y, trans_cumulat(iv,joker,joker,ip),
887  x );
888  diy_dpath[iq](ip+1,iv,joker) += y;
889  //
890  // The terms associated with Delta-s:
891  if( jacobian_quantities[iq].Subtag() ==
892  "HSE on" )
893  {
894  for( Index is1=0; is1<ns; is1++ ) {
895  for( Index is2=0; is2<ns; is2++ ) {
896  ext_mat(is1,is2) = 0.5 * (
897  ppath_abs(iv,is1,is2,ip ) +
898  ppath_abs(iv,is1,is2,ip+1) );
899  } }
900  // dl for disturbed tbar
901  const Numeric tbar = 0.5 * (
902  ppath_t[ip] + ppath_t[ip+1] );
903  const Numeric dl = ppath.lstep[ip] *
904  ( 1 + dt/tbar );
905  ext2trans( dtdx, extmat_case[ip][iv],
906  ext_mat, dl );
907  for( Index is1=0; is1<ns; is1++ ) {
908  for( Index is2=0; is2<ns; is2++ ) {
909  dtdx(is1,is2) = (1/dt) *
910  ( dtdx(is1,is2) -
911  trans_partial(iv,is1,is2,ip) );
912  } }
913  mult( x, dtdx, sibi(iv,joker) );
914  mult( y, trans_cumulat(iv,joker,joker,ip),
915  x );
916  // Contribution shared between the two
917  // points and is proportional to 1/t
918  // See also AUG.
919  for( Index is=0; is<ns; is++ )
920  {
921  diy_dpath[iq](ip ,iv,is) += y[is] *
922  0.5 * tbar / ppath_t[ip];
923  diy_dpath[iq](ip+1,iv,is) += y[is] *
924  0.5 * tbar / ppath_t[ip+1];
925  }
926  } // HSE
927  } // General case
928  } // Frequency loop
929  } // Temperature
930  } // if this analytical
931  } // for iq
932  } // if any analytical
933  //###################################################################
934 
935 
936  // Spectrum at end of ppath step
937  emission_rtstep( iy, stokes_dim, bbar, extmat_case[ip],
938  trans_partial(joker,joker,joker,ip) );
939 
940 
941  //=== iy_aux part ===================================================
942  // Pressure
943  if( auxPressure >= 0 )
944  { iy_aux[auxPressure](0,0,0,ip) = ppath_p[ip]; }
945  // Temperature
946  if( auxTemperature >= 0 )
947  { iy_aux[auxTemperature](0,0,0,ip) = ppath_t[ip]; }
948  // VMR
949  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
950  { iy_aux[auxVmrSpecies[j]](0,0,0,ip) = ppath_vmr(auxVmrIsp[j],ip);}
951  // Absorption
952  if( auxAbsSum >= 0 )
953  { for( Index iv=0; iv<nf; iv++ ) {
954  for( Index is1=0; is1<ns; is1++ ){
955  for( Index is2=0; is2<ns; is2++ ){
956  iy_aux[auxAbsSum](iv,is1,is2,ip) =
957  ppath_abs(iv,is1,is2,ip); } } } }
958  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
959  { for( Index iv=0; iv<nf; iv++ ) {
960  for( Index is1=0; is1<ns; is1++ ){
961  for( Index is2=0; is2<ns; is2++ ){
962  iy_aux[auxAbsSpecies[j]](iv,is1,is2,ip) =
963  abs_per_species(auxAbsIsp[j],iv,is1,is2,ip); } } } }
964  // Radiance
965  if( auxIy >= 0 )
966  { iy_aux[auxIy](joker,joker,0,ip) = iy; }
967  //===================================================================
968  }
969 
970  //### jacobian part #####################################################
971  // Map jacobians from ppath to retrieval grids
972  // (this operation corresponds to the term Dx_i/Dx)
973  if( j_analytical_do )
974  {
975  // Weight with iy_transmission
976  if( !iy_agenda_call1 )
977  {
978  Matrix X, Y(ns,diy_dpath[0].npages());
979  //
980  FOR_ANALYTICAL_JACOBIANS_DO(
981  for( Index iv=0; iv<nf; iv++ )
982  {
983  X = transpose( diy_dpath[iq](joker,iv,joker) );
984  mult( Y, iy_transmission(iv,joker,joker), X );
985  diy_dpath[iq](joker,iv,joker) = transpose( Y );
986  }
987  )
988  }
989 
990  // Map to retrieval grids
991  FOR_ANALYTICAL_JACOBIANS_DO(
992  diy_from_path_to_rgrids( diy_dx[iq], jacobian_quantities[iq],
993  diy_dpath[iq], atmosphere_dim, ppath, ppath_p );
994  )
995  }
996  //#######################################################################
997  } // if np>1
998 
999 
1000  // Unit conversions
1001  if( iy_agenda_call1 )
1002  {
1003  // If any conversion, check that standard form of Planck used
1004  if( !chk_if_std_blackbody_agenda( ws, blackbody_radiation_agenda ) )
1005  {
1006  ostringstream os;
1007  os << "When any unit conversion is applied, "
1008  << "*blackbody_radiation_agenda\nmust use "
1009  << "*blackbody_radiationPlanck* or a corresponding WSM.\nA test "
1010  << "call of the agenda indicates that this is not the case.";
1011  throw runtime_error( os.str() );
1012  }
1013 
1014  // Determine refractive index to use for the n2 radiance law
1015  Numeric n = 1.0; // First guess is that sensor is in space
1016  //
1017  if( ppath.end_lstep == 0 ) // If true, sensor inside the atmosphere
1018  { n = ppath.nreal[np-1]; }
1019 
1020  // Polarisation index variable
1021  ArrayOfIndex i_pol(ns);
1022  for( Index is=0; is<ns; is++ )
1023  { i_pol[is] = is + 1; }
1024 
1025  // Jacobian part (must be converted to Tb before iy for PlanckBT)
1026  //
1027  if( j_analytical_do )
1028  {
1029  FOR_ANALYTICAL_JACOBIANS_DO( apply_iy_unit2( diy_dx[iq], iy, iy_unit,
1030  f_grid, n, i_pol ); )
1031  }
1032 
1033  // iy
1034  apply_iy_unit( iy, iy_unit, f_grid, n, i_pol );
1035 
1036  // iy_aux
1037  for( Index q=0; q<iy_aux.nelem(); q++ )
1038  {
1039  if( iy_aux_vars[q] == "iy")
1040  {
1041  for( Index ip=0; ip<ppath.np; ip++ )
1042  { apply_iy_unit( iy_aux[q](joker,joker,0,ip), iy_unit, f_grid,
1043  ppath.nreal[ip], i_pol ); }
1044  }
1045  }
1046  }
1047 }
1048 
1049 
1050 
1051 
1052 
1053 /* Workspace method: Doxygen documentation will be auto-generated */
1054 void iyLoopFrequencies(
1055  Workspace& ws,
1056  Matrix& iy,
1057  ArrayOfTensor4& iy_aux,
1058  Ppath& ppath,
1059  ArrayOfTensor3& diy_dx,
1060  const ArrayOfString& iy_aux_vars,
1061  const Index& stokes_dim,
1062  const Vector& f_grid,
1063  const Tensor3& t_field,
1064  const Tensor3& z_field,
1065  const Tensor4& vmr_field,
1066  const Index& cloudbox_on,
1067  const Index& iy_agenda_call1,
1068  const Tensor3& iy_transmission,
1069  const Vector& rte_pos,
1070  const Vector& rte_los,
1071  const Vector& rte_pos2,
1072  const Index& jacobian_do,
1073  const Agenda& iy_sub_agenda,
1074  const Verbosity& )
1075 {
1076  // Throw error if unsupported features are requested
1077  if( !iy_agenda_call1 )
1078  throw runtime_error(
1079  "Recursive usage not possible (iy_agenda_call1 must be 1)" );
1080  if( iy_transmission.ncols() )
1081  throw runtime_error( "*iy_transmission* must be empty" );
1082 
1083  const Index nf = f_grid.nelem();
1084 
1085  for( Index i=0; i<nf; i++ )
1086  {
1087  // Variables for 1 frequency
1088  Matrix iy1;
1089  ArrayOfTensor4 iy_aux1;
1090  ArrayOfTensor3 diy_dx1;
1091 
1092  iy_sub_agendaExecute( ws, iy1, iy_aux1, ppath, diy_dx1,
1093  1, iy_transmission, iy_aux_vars, cloudbox_on,
1094  jacobian_do, t_field, z_field, vmr_field,
1095  Vector(1,f_grid[i]),
1096  rte_pos, rte_los, rte_pos2, iy_sub_agenda );
1097 
1098  // After first frequency, give output its size
1099  if( i == 0 )
1100  {
1101  iy.resize( nf, stokes_dim );
1102  //
1103  iy_aux.resize( iy_aux1.nelem() );
1104  for( Index q=0; q<iy_aux1.nelem(); q++ )
1105  {
1106  if( iy_aux1[q].ncols() > 1 )
1107  throw runtime_error( "When using this method, *iy_aux_vars* "
1108  "is not allowed to include along-the-path variables." );
1109  iy_aux[q].resize(nf,iy_aux1[q].npages(),iy_aux1[q].nrows(),1);
1110  }
1111  //
1112  diy_dx.resize( diy_dx1.nelem() );
1113  for( Index q=0; q<diy_dx1.nelem(); q++ )
1114  { diy_dx[q].resize( diy_dx1[q].npages(), nf, stokes_dim ); }
1115  }
1116 
1117  // Copy to output variables
1118  iy(i,joker) = iy1(0,joker);
1119  for( Index q=0; q<iy_aux1.nelem(); q++ )
1120  { iy_aux[q](i,joker,joker,0) = iy_aux1[q](0,joker,joker,0); }
1121  for( Index q=0; q<diy_dx1.nelem(); q++ )
1122  { diy_dx[q](joker,i,joker) = diy_dx1[q](joker,0,joker); }
1123  }
1124 }
1125 
1126 
1127 
1128 
1129 
1130 /* Workspace method: Doxygen documentation will be auto-generated */
1131 void iyMC(
1132  Workspace& ws,
1133  Matrix& iy,
1134  ArrayOfTensor4& iy_aux,
1135  ArrayOfTensor3& diy_dx,
1136  const Index& iy_agenda_call1,
1137  const Tensor3& iy_transmission,
1138  const Vector& rte_pos,
1139  const Vector& rte_los,
1140  const ArrayOfString& iy_aux_vars,
1141  const Index& jacobian_do,
1142  const Index& atmosphere_dim,
1143  const Vector& p_grid,
1144  const Vector& lat_grid,
1145  const Vector& lon_grid,
1146  const Tensor3& z_field,
1147  const Tensor3& t_field,
1148  const Tensor4& vmr_field,
1149  const Vector& refellipsoid,
1150  const Matrix& z_surface,
1151  const Index& cloudbox_on,
1152  const ArrayOfIndex& cloudbox_limits,
1153  const Index& stokes_dim,
1154  const Vector& f_grid,
1155  const ArrayOfSingleScatteringData& scat_data_array,
1156  const Agenda& iy_space_agenda,
1157  const Agenda& surface_rtprop_agenda,
1158  const Agenda& propmat_clearsky_agenda,
1159  const Agenda& ppath_step_agenda,
1160  const Numeric& ppath_lraytrace,
1161  const Tensor4& pnd_field,
1162  const String& iy_unit,
1163  const Numeric& mc_std_err,
1164  const Index& mc_max_time,
1165  const Index& mc_max_iter,
1166  const Index& mc_min_iter,
1167  const Verbosity& verbosity)
1168 {
1169  // Throw error if unsupported features are requested
1170  if( atmosphere_dim != 3 )
1171  throw runtime_error(
1172  "Only 3D atmospheres are allowed (atmosphere_dim must be 3)" );
1173  if( !cloudbox_on )
1174  throw runtime_error(
1175  "The cloudbox must be activated (cloudbox_on must be 1)" );
1176  if( jacobian_do )
1177  throw runtime_error(
1178  "This method does not provide any jacobians (jacobian_do must be 0)" );
1179  if( !iy_agenda_call1 )
1180  throw runtime_error(
1181  "Recursive usage not possible (iy_agenda_call1 must be 1)" );
1182  if( iy_transmission.ncols() )
1183  throw runtime_error( "*iy_transmission* must be empty" );
1184 
1185 
1186  // Size output variables
1187  //
1188  const Index nf = f_grid.nelem();
1189  //
1190  iy.resize( nf, stokes_dim );
1191  diy_dx.resize(0);
1192  //
1193  //=== iy_aux part ===========================================================
1194  Index auxError = -1;
1195  {
1196  const Index naux = iy_aux_vars.nelem();
1197  iy_aux.resize( naux );
1198  //
1199  for( Index i=0; i<naux; i++ )
1200  {
1201  if( iy_aux_vars[i] == "Error (uncorrelated)" )
1202  { auxError = i; iy_aux[i].resize( nf, stokes_dim, 1, 1 ); }
1203  else
1204  {
1205  ostringstream os;
1206  os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
1207  << "\"\nThis choice is not recognised.";
1208  throw runtime_error( os.str() );
1209  }
1210  }
1211  }
1212  //===========================================================================
1213 
1214  // Some MC variables are only local here
1215  //
1216  MCAntenna mc_antenna;
1217  mc_antenna.set_pencil_beam();
1218 
1219  // Pos and los must be matrices
1220  Matrix pos(1,3), los(1,2);
1221  //
1222  pos(0,joker) = rte_pos;
1223  los(0,joker) = rte_los;
1224 
1225  Workspace l_ws (ws);
1226  Agenda l_ppath_step_agenda (ppath_step_agenda);
1227  Agenda l_iy_space_agenda (iy_space_agenda);
1228  Agenda l_propmat_clearsky_agenda (propmat_clearsky_agenda);
1229  Agenda l_surface_rtprop_agenda (surface_rtprop_agenda);
1230 
1231  String fail_msg;
1232  bool failed = false;
1233 
1234  if (nf)
1235 #pragma omp parallel for \
1236  if (!arts_omp_in_parallel() && nf > 1) \
1237  firstprivate(l_ws, l_ppath_step_agenda, l_iy_space_agenda, \
1238  l_propmat_clearsky_agenda, l_surface_rtprop_agenda)
1239  for( Index f_index=0; f_index<nf; f_index++ )
1240  {
1241  if (failed) continue;
1242 
1243  try {
1244  ArrayOfSingleScatteringData scat_data_array_mono;
1245 
1246  scat_data_array_monoCalc( scat_data_array_mono, scat_data_array,
1247  f_grid, f_index, verbosity );
1248 
1249  // Seed reset for each loop. If not done, the errors
1250  // appear to be highly correlated.
1251  Index mc_seed;
1252  MCSetSeedFromTime( mc_seed, verbosity );
1253 
1254  Vector y, mc_error;
1255  Index mc_iteration_count;
1256  Tensor3 mc_points;
1257 
1258  MCGeneral( l_ws, y, mc_iteration_count, mc_error, mc_points, mc_antenna,
1259  f_grid, f_index, pos, los, stokes_dim, atmosphere_dim,
1260  l_ppath_step_agenda, ppath_lraytrace, l_iy_space_agenda,
1261  l_surface_rtprop_agenda, l_propmat_clearsky_agenda,
1262  p_grid, lat_grid, lon_grid, z_field,
1263  refellipsoid, z_surface, t_field, vmr_field,
1264  cloudbox_on, cloudbox_limits,
1265  pnd_field, scat_data_array_mono, 1, 1, 1,
1266  mc_seed, iy_unit, mc_std_err, mc_max_time, mc_max_iter,
1267  mc_min_iter, verbosity);
1268  //cout << "Error: " << mc_error << endl;
1269  //cout << "N photons: " << mc_iteration_count << endl;
1270  assert( y.nelem() == stokes_dim );
1271 
1272  iy(f_index,joker) = y;
1273 
1274  if( auxError >= 0 )
1275  { iy_aux[auxError](f_index,joker,0,0) = mc_error; }
1276  } catch (runtime_error e) {
1277  ostringstream os;
1278  os << "Error for f_index = " << f_index << " (" << f_grid[f_index]
1279  << ")" << endl << e.what();
1280 #pragma omp critical (iyMC_fail)
1281  { failed = true; fail_msg = os.str(); }
1282  continue;
1283  }
1284  }
1285 
1286  if (failed)
1287  throw runtime_error(fail_msg);
1288 }
1289 
1290 
1291 
1292 
1293 
1294 /* Workspace method: Doxygen documentation will be auto-generated */
1295 void iyReplaceFromAux(
1296  Matrix& iy,
1297  const ArrayOfTensor4& iy_aux,
1298  const ArrayOfString& iy_aux_vars,
1299  const Index& jacobian_do,
1300  const String& aux_var,
1301  const Verbosity& )
1302 {
1303  if( iy_aux.nelem() != iy_aux_vars.nelem() )
1304  throw runtime_error( "*iy_aux* and *iy_aux_vars* must have the same "
1305  "number of elements." );
1306 
1307  if( jacobian_do )
1308  throw runtime_error( "This method can not provide any jacobians and "
1309  "*jacobian_do* must be 0." );
1310 
1311  bool ready = false;
1312 
1313  for( Index i=0; i<iy_aux.nelem() && !ready; i++ )
1314  {
1315  if( iy_aux_vars[i] == aux_var )
1316  {
1317  if( iy_aux[i].nrows() > 1 || iy_aux[i].ncols() > 1 )
1318  {
1319  throw runtime_error( "If an auxiliary variable shall be inserted "
1320  "in *iy*, its row and page dimensions must have size 1." );
1321  }
1322  if( iy_aux[i].nbooks() != iy.nrows() )
1323  {
1324  throw runtime_error( "If an auxiliary variable shall be inserted "
1325  "in *iy*, its frequency dimension must match"
1326  "the length of existing *iy*." );
1327  }
1328 
1329  iy = 0;
1330 
1331  for( Index iv=0; iv<iy.nrows(); iv++ )
1332  {
1333  for( Index is=0; is<iy_aux[i].npages(); is++ )
1334  {
1335  iy(iv,is) = iy_aux[i](iv,is,0,0);
1336  }
1337  }
1338 
1339  ready = true;
1340  }
1341  }
1342 
1343  if( !ready )
1344  throw runtime_error( "The selected auxiliary variable to insert in *iy* "
1345  "is either not defined at all or is not set." );
1346 
1347 }
1348 
1349 
1350 
1351 
1352 
1353 /* Workspace method: Doxygen documentation will be auto-generated */
1354 void iy_auxFillParticleVariables(
1355  ArrayOfTensor4& iy_aux,
1356  const Index& atmfields_checked,
1357  const Index& cloudbox_checked,
1358  const Index& atmosphere_dim,
1359  const Index& cloudbox_on,
1360  const ArrayOfIndex& cloudbox_limits,
1361  const Tensor4& pnd_field,
1362  const Matrix& particle_masses,
1363  const Ppath& ppath,
1364  const ArrayOfString& iy_aux_vars,
1365  const Verbosity& )
1366 {
1367  // Some sizes
1368  const Index np = ppath.np;
1369  const Index naux = iy_aux_vars.nelem();
1370 
1371  // Input checks
1372  if( atmfields_checked != 1 )
1373  throw runtime_error( "The atmospheric fields must be flagged to have "
1374  "passed a consistency check (atmfields_checked=1)." );
1375  if( cloudbox_checked != 1 )
1376  throw runtime_error( "The cloudbox must be flagged to have "
1377  "passed a consistency check (cloudbox_checked=1)." );
1378  if( !cloudbox_on )
1379  throw runtime_error(
1380  "The cloudbox must be activated (cloudbox_on must be 1)" );
1381  if( iy_aux.nelem() != naux )
1382  throw runtime_error( "*iy_aux_vars* and *iy_aux* must have the same array "
1383  "length. (You can not call this WSM before the main "
1384  "iy-WSM.)" );
1385 
1386  // Analayse iy_aux_vars
1387  ArrayOfIndex auxPartCont(0), auxPartContI(0);
1388  ArrayOfIndex auxPartField(0), auxPartFieldI(0);
1389  //
1390  for( Index i=0; i<naux; i++ )
1391  {
1392  if( iy_aux_vars[i].substr(0,14) == "Mass content, " )
1393  {
1394  Index icont;
1395  istringstream is(iy_aux_vars[i].substr(14,2));
1396  is >> icont;
1397  if( icont < 0 || icont>=particle_masses.ncols() )
1398  {
1399  ostringstream os;
1400  os << "You have selected particle mass content category with "
1401  << "index " << icont << ".\nThis category is not defined!";
1402  throw runtime_error( os.str() );
1403  }
1404  auxPartCont.push_back(i);
1405  auxPartContI.push_back(icont);
1406  iy_aux[i].resize( 1, 1, 1, np );
1407  }
1408  else if( iy_aux_vars[i].substr(0,10) == "PND, type " )
1409  {
1410  Index ip;
1411  istringstream is(iy_aux_vars[i].substr(10,2));
1412  is >> ip;
1413  if( ip < 0 || ip>=pnd_field.nbooks() )
1414  {
1415  ostringstream os;
1416  os << "You have selected particle number density field with "
1417  << "index " << ip << ".\nThis field is not defined!";
1418  throw runtime_error( os.str() );
1419  }
1420  auxPartField.push_back(i);
1421  auxPartFieldI.push_back(ip);
1422  iy_aux[i].resize( 1, 1, 1, np );
1423  }
1424  }
1425 
1426  if( auxPartCont.nelem() + auxPartField.nelem() > 0 )
1427  {
1428  // PND along the ppath
1429  Matrix ppath_pnd( pnd_field.nbooks(), np, 0 );
1430  //
1431  for( Index ip=0; ip<np; ip++ )
1432  {
1433  Matrix itw( 1, Index(pow(2.0,Numeric(atmosphere_dim))) );
1434 
1435  ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
1436  GridPos gp_lat, gp_lon;
1437  if( atmosphere_dim >= 2 ) { gridpos_copy( gp_lat, ppath.gp_lat[ip] );}
1438  if( atmosphere_dim == 3 ) { gridpos_copy( gp_lon, ppath.gp_lon[ip] );}
1439  if( is_gp_inside_cloudbox( ppath.gp_p[ip], gp_lat, gp_lon,
1440  cloudbox_limits, true, atmosphere_dim ) )
1441  {
1442  interp_cloudfield_gp2itw( itw(0,joker),
1443  gpc_p[0], gpc_lat[0], gpc_lon[0],
1444  ppath.gp_p[ip], gp_lat, gp_lon,
1445  atmosphere_dim, cloudbox_limits );
1446  for( Index i=0; i<pnd_field.nbooks(); i++ )
1447  {
1448  interp_atmfield_by_itw( ppath_pnd(i,ip), atmosphere_dim,
1449  pnd_field(i,joker,joker,joker),
1450  gpc_p, gpc_lat, gpc_lon, itw );
1451  }
1452  }
1453  }
1454 
1455  // Loop ppath steps
1456  for( Index ip=0; ip<np; ip++ )
1457  {
1458  // Particle mass content
1459  for( Index j=0; j<auxPartCont.nelem(); j++ )
1460  { iy_aux[auxPartCont[j]](0,0,0,ip) = ppath_pnd(joker,ip) *
1461  particle_masses(joker,auxPartContI[j]); }
1462  // Particle field
1463  for( Index j=0; j<auxPartField.nelem(); j++ )
1464  { iy_aux[auxPartField[j]](0,0,0,ip) =
1465  ppath_pnd(auxPartFieldI[j],ip); }
1466  }
1467  }
1468 }
1469 
1470 
1471 void yCalc_mblock_loop_body(
1472  bool& failed,
1473  String& fail_msg,
1474  ArrayOfArrayOfVector& iyb_aux_array,
1475  Workspace& ws,
1476  Vector& y,
1477  Vector& y_f,
1478  ArrayOfIndex& y_pol,
1479  Matrix& y_pos,
1480  Matrix& y_los,
1481  Matrix& jacobian,
1482  const Index& atmosphere_dim,
1483  const Tensor3& t_field,
1484  const Tensor3& z_field,
1485  const Tensor4& vmr_field,
1486  const Index& cloudbox_on,
1487  const Index& stokes_dim,
1488  const Vector& f_grid,
1489  const Matrix& sensor_pos,
1490  const Matrix& sensor_los,
1491  const Matrix& transmitter_pos,
1492  const Vector& mblock_za_grid,
1493  const Vector& mblock_aa_grid,
1494  const Index& antenna_dim,
1495  const Sparse& sensor_response,
1496  const Vector& sensor_response_f,
1497  const ArrayOfIndex& sensor_response_pol,
1498  const Vector& sensor_response_za,
1499  const Vector& sensor_response_aa,
1500  const Agenda& iy_main_agenda,
1501  const Agenda& jacobian_agenda,
1502  const Index& jacobian_do,
1503  const ArrayOfRetrievalQuantity& jacobian_quantities,
1504  const ArrayOfArrayOfIndex& jacobian_indices,
1505  const ArrayOfString& iy_aux_vars,
1506  const Verbosity& verbosity,
1507  const Index& mblock_index,
1508  const Index& n1y,
1509  const Index& j_analytical_do)
1510 {
1511  try
1512  {
1513  // Calculate monochromatic pencil beam data for 1 measurement block
1514  //
1515  Vector iyb, iyb_error, yb(n1y);
1516  ArrayOfMatrix diyb_dx;
1517  //
1518  iyb_calc(ws, iyb, iyb_aux_array[mblock_index], diyb_dx,
1519  mblock_index, atmosphere_dim, t_field, z_field, vmr_field,
1520  cloudbox_on, stokes_dim, f_grid, sensor_pos, sensor_los,
1521  transmitter_pos, mblock_za_grid, mblock_aa_grid, antenna_dim,
1522  iy_main_agenda, j_analytical_do, jacobian_quantities,
1523  jacobian_indices, iy_aux_vars, verbosity);
1524 
1525 
1526  // Apply sensor response matrix on iyb, and put into y
1527  //
1528  const Range rowind = get_rowindex_for_mblock(sensor_response,
1529  mblock_index);
1530  const Index row0 = rowind.get_start();
1531  //
1532  mult( yb, sensor_response, iyb );
1533  //
1534  y[rowind] = yb; // *yb* also used below, as input to jacobian_agenda
1535 
1536  // Fill information variables
1537  //
1538  for( Index i=0; i<n1y; i++ )
1539  {
1540  y_f[row0+i] = sensor_response_f[i];
1541  y_pol[row0+i] = sensor_response_pol[i];
1542  y_pos(row0+i,joker) = sensor_pos(mblock_index,joker);
1543  y_los(row0+i,0) = sensor_los(mblock_index,0) +
1544  sensor_response_za[i];
1545  if( sensor_response_aa.nelem() )
1546  {
1547  y_los(row0+i,1) = sensor_los(mblock_index,0) +
1548  sensor_response_aa[i];
1549  }
1550  }
1551 
1552  // Apply sensor response matrix on diyb_dx, and put into jacobian
1553  // (that is, analytical jacobian part)
1554  //
1555  if( j_analytical_do )
1556  {
1557  FOR_ANALYTICAL_JACOBIANS_DO(
1558  mult(jacobian(rowind,
1559  Range(jacobian_indices[iq][0],
1560  jacobian_indices[iq][1]-jacobian_indices[iq][0]+1)),
1561  sensor_response, diyb_dx[iq] );
1562  )
1563  }
1564 
1565  // Rest of *jacobian*
1566  //
1567  if( jacobian_do )
1568  {
1569  jacobian_agendaExecute( ws, jacobian, mblock_index, iyb, yb,
1570  jacobian_agenda );
1571  }
1572  }
1573  catch (runtime_error e)
1574  {
1575 #pragma omp critical (yCalc_fail)
1576  { fail_msg = e.what(); failed = true; }
1577  }
1578 }
1579 
1580 
1581 
1582 /* Workspace method: Doxygen documentation will be auto-generated */
1583 void yCalc(
1584  Workspace& ws,
1585  Vector& y,
1586  Vector& y_f,
1587  ArrayOfIndex& y_pol,
1588  Matrix& y_pos,
1589  Matrix& y_los,
1590  ArrayOfVector& y_aux,
1591  Matrix& jacobian,
1592  const Index& atmfields_checked,
1593  const Index& atmgeom_checked,
1594  const Index& atmosphere_dim,
1595  const Tensor3& t_field,
1596  const Tensor3& z_field,
1597  const Tensor4& vmr_field,
1598  const Index& cloudbox_on,
1599  const Index& cloudbox_checked,
1600  const Index& sensor_checked,
1601  const Index& stokes_dim,
1602  const Vector& f_grid,
1603  const Matrix& sensor_pos,
1604  const Matrix& sensor_los,
1605  const Matrix& transmitter_pos,
1606  const Vector& mblock_za_grid,
1607  const Vector& mblock_aa_grid,
1608  const Index& antenna_dim,
1609  const Sparse& sensor_response,
1610  const Vector& sensor_response_f,
1611  const ArrayOfIndex& sensor_response_pol,
1612  const Vector& sensor_response_za,
1613  const Vector& sensor_response_aa,
1614  const Agenda& iy_main_agenda,
1615  const Agenda& jacobian_agenda,
1616  const Index& jacobian_do,
1617  const ArrayOfRetrievalQuantity& jacobian_quantities,
1618  const ArrayOfArrayOfIndex& jacobian_indices,
1619  const ArrayOfString& iy_aux_vars,
1620  const Verbosity& verbosity )
1621 {
1622  CREATE_OUT3;
1623 
1624  // Basics
1625  //
1626  chk_if_in_range( "stokes_dim", stokes_dim, 1, 4 );
1627  //
1628  if ( f_grid.nelem() == 0 )
1629  { throw runtime_error ( "The frequency grid is empty." ); }
1630  chk_if_increasing ( "f_grid", f_grid );
1631  if( f_grid[0] <= 0)
1632  { throw runtime_error( "All frequencies in *f_grid* must be > 0." ); }
1633  //
1634  if( atmfields_checked != 1 )
1635  throw runtime_error( "The atmospheric fields must be flagged to have "
1636  "passed a consistency check (atmfields_checked=1)." );
1637  if( atmgeom_checked != 1 )
1638  throw runtime_error( "The atmospheric geometry must be flagged to have "
1639  "passed a consistency check (atmgeom_checked=1)." );
1640  if( cloudbox_checked != 1 )
1641  throw runtime_error( "The cloudbox must be flagged to have "
1642  "passed a consistency check (cloudbox_checked=1)." );
1643  if( sensor_checked != 1 )
1644  throw runtime_error( "The sensor variables must be flagged to have "
1645  "passed a consistency check (sensor_checked=1)." );
1646 
1647 
1648  // Some sizes
1649  const Index nf = f_grid.nelem();
1650  const Index nza = mblock_za_grid.nelem();
1651  Index naa = mblock_aa_grid.nelem();
1652  if( antenna_dim == 1 )
1653  { naa = 1; }
1654  const Index n1y = sensor_response.nrows();
1655  const Index nmblock = sensor_pos.nrows();
1656  const Index niyb = nf * nza * naa * stokes_dim;
1657 
1658 
1659  //---------------------------------------------------------------------------
1660  // Allocations and resizing
1661  //---------------------------------------------------------------------------
1662 
1663  // Resize *y* and *y_XXX*
1664  //
1665  y.resize( nmblock*n1y );
1666  y_f.resize( nmblock*n1y );
1667  y_pol.resize( nmblock*n1y );
1668  y_pos.resize( nmblock*n1y, sensor_pos.ncols() );
1669  y_los.resize( nmblock*n1y, sensor_los.ncols() );
1670 
1671  // For y_aux we don't know the number of quantities, and we need to
1672  // store all output
1673  ArrayOfArrayOfVector iyb_aux_array( nmblock );
1674 
1675  // Jacobian variables
1676  //
1677  Index j_analytical_do = 0;
1678  //
1679  if( jacobian_do )
1680  {
1681  jacobian.resize( nmblock*n1y,
1682  jacobian_indices[jacobian_indices.nelem()-1][1]+1 );
1683  jacobian = 0;
1684  //
1685  FOR_ANALYTICAL_JACOBIANS_DO(
1686  j_analytical_do = 1;
1687  )
1688  }
1689  else
1690  { jacobian.resize( 0, 0 ); }
1691 
1692 
1693  //---------------------------------------------------------------------------
1694  // The calculations
1695  //---------------------------------------------------------------------------
1696 
1697  String fail_msg;
1698  bool failed = false;
1699 
1700  if (nmblock >= arts_omp_get_max_threads()
1701  || (nf <= nmblock && nmblock >= nza))
1702  {
1703  out3 << " Parallelizing mblock loop (" << nmblock << " iterations)\n";
1704 
1705  // We have to make a local copy of the Workspace and the agendas because
1706  // only non-reference types can be declared firstprivate in OpenMP
1707  Workspace l_ws (ws);
1708  Agenda l_jacobian_agenda (jacobian_agenda);
1709  Agenda l_iy_main_agenda (iy_main_agenda);
1710 
1711 #pragma omp parallel for \
1712 firstprivate(l_ws, l_jacobian_agenda, l_iy_main_agenda)
1713  for( Index mblock_index=0; mblock_index<nmblock; mblock_index++ )
1714  {
1715  // Skip remaining iterations if an error occurred
1716  if (failed) continue;
1717 
1718  yCalc_mblock_loop_body(failed,
1719  fail_msg,
1720  iyb_aux_array,
1721  l_ws,
1722  y,
1723  y_f,
1724  y_pol,
1725  y_pos,
1726  y_los,
1727  jacobian,
1728  atmosphere_dim,
1729  t_field,
1730  z_field,
1731  vmr_field,
1732  cloudbox_on,
1733  stokes_dim,
1734  f_grid,
1735  sensor_pos,
1736  sensor_los,
1737  transmitter_pos,
1738  mblock_za_grid,
1739  mblock_aa_grid,
1740  antenna_dim,
1741  sensor_response,
1742  sensor_response_f,
1743  sensor_response_pol,
1744  sensor_response_za,
1745  sensor_response_aa,
1746  l_iy_main_agenda,
1747  l_jacobian_agenda,
1748  jacobian_do,
1749  jacobian_quantities,
1750  jacobian_indices,
1751  iy_aux_vars,
1752  verbosity,
1753  mblock_index,
1754  n1y,
1755  j_analytical_do);
1756 
1757  } // End mblock loop
1758  }
1759  else
1760  {
1761  out3 << " Not parallelizing mblock loop (" << nmblock << " iterations)\n";
1762 
1763  for( Index mblock_index=0; mblock_index<nmblock; mblock_index++ )
1764  {
1765  // Skip remaining iterations if an error occurred
1766  if (failed) continue;
1767 
1768  yCalc_mblock_loop_body(failed,
1769  fail_msg,
1770  iyb_aux_array,
1771  ws,
1772  y,
1773  y_f,
1774  y_pol,
1775  y_pos,
1776  y_los,
1777  jacobian,
1778  atmosphere_dim,
1779  t_field,
1780  z_field,
1781  vmr_field,
1782  cloudbox_on,
1783  stokes_dim,
1784  f_grid,
1785  sensor_pos,
1786  sensor_los,
1787  transmitter_pos,
1788  mblock_za_grid,
1789  mblock_aa_grid,
1790  antenna_dim,
1791  sensor_response,
1792  sensor_response_f,
1793  sensor_response_pol,
1794  sensor_response_za,
1795  sensor_response_aa,
1796  iy_main_agenda,
1797  jacobian_agenda,
1798  jacobian_do,
1799  jacobian_quantities,
1800  jacobian_indices,
1801  iy_aux_vars,
1802  verbosity,
1803  mblock_index,
1804  n1y,
1805  j_analytical_do);
1806 
1807  } // End mblock loop
1808  }
1809 
1810  // Rethrow exception if a runtime error occurred in the mblock loop
1811  if (failed) throw runtime_error(fail_msg);
1812 
1813  // Compile y_aux
1814  //
1815  const Index nq = iyb_aux_array[0].nelem();
1816  y_aux.resize( nq );
1817  //
1818  for( Index q=0; q<nq; q++ )
1819  {
1820  y_aux[q].resize( nmblock*n1y );
1821  //
1822  for( Index mblock_index=0; mblock_index<nmblock; mblock_index++ )
1823  {
1824  const Range rowind = get_rowindex_for_mblock( sensor_response,
1825  mblock_index );
1826  const Index row0 = rowind.get_start();
1827 
1828  // The sensor response must be applied in a special way for
1829  // uncorrelated errors. Schematically: sqrt( H.^2 * y.^2 )
1830  if( iy_aux_vars[q] == "Error (uncorrelated)" )
1831  {
1832  for( Index i=0; i<n1y; i++ )
1833  {
1834  const Index row = row0+i;
1835  y_aux[q][row] = 0;
1836  for( Index j=0; j<niyb; j++ )
1837  { y_aux[q][row] += pow( sensor_response(i,j) *
1838  iyb_aux_array[mblock_index][q][j], (Numeric)2.0 ); }
1839  y_aux[q][row] = sqrt( y_aux[q][row] );
1840  }
1841  }
1842  else
1843  { mult( y_aux[q][rowind], sensor_response,
1844  iyb_aux_array[mblock_index][q] ); }
1845  }
1846  }
1847 }
1848 
1849 
1850 
1851 /* Workspace method: Doxygen documentation will be auto-generated */
1852 void yCalcAppend(
1853  Workspace& ws,
1854  Vector& y,
1855  Vector& y_f,
1856  ArrayOfIndex& y_pol,
1857  Matrix& y_pos,
1858  Matrix& y_los,
1859  ArrayOfVector& y_aux,
1860  Matrix& jacobian,
1861  ArrayOfRetrievalQuantity& jacobian_quantities,
1862  ArrayOfArrayOfIndex& jacobian_indices,
1863  const Index& atmfields_checked,
1864  const Index& atmgeom_checked,
1865  const Index& atmosphere_dim,
1866  const Tensor3& t_field,
1867  const Tensor3& z_field,
1868  const Tensor4& vmr_field,
1869  const Index& cloudbox_on,
1870  const Index& cloudbox_checked,
1871  const Index& sensor_checked,
1872  const Index& stokes_dim,
1873  const Vector& f_grid,
1874  const Matrix& sensor_pos,
1875  const Matrix& sensor_los,
1876  const Matrix& transmitter_pos,
1877  const Vector& mblock_za_grid,
1878  const Vector& mblock_aa_grid,
1879  const Index& antenna_dim,
1880  const Sparse& sensor_response,
1881  const Vector& sensor_response_f,
1882  const ArrayOfIndex& sensor_response_pol,
1883  const Vector& sensor_response_za,
1884  const Vector& sensor_response_aa,
1885  const Agenda& iy_main_agenda,
1886  const Agenda& jacobian_agenda,
1887  const Index& jacobian_do,
1888  const ArrayOfString& iy_aux_vars,
1889  const ArrayOfRetrievalQuantity& jacobian_quantities1,
1890  const ArrayOfArrayOfIndex& jacobian_indices1,
1891  const Index& append_instrument_wfs,
1892  const Verbosity& verbosity )
1893 {
1894  // Some initial checks of old measurement
1895  const Index n1 = y.nelem();
1896  Index nrq1 = 0;
1897  if( n1 == 0 )
1898  throw runtime_error( "Input *y* is empty. Use *yCalc*" );
1899  if( y_f.nelem() != n1 )
1900  throw runtime_error( "Lengths of input *y* and *y_f* are inconsistent." );
1901  if( y_pol.nelem() != n1 )
1902  throw runtime_error( "Lengths of input *y* and *y_pol* are inconsistent." );
1903  if( y_pos.nrows() != n1 )
1904  throw runtime_error( "Sizes of input *y* and *y_pos* are inconsistent." );
1905  if( y_los.nrows() != n1 )
1906  throw runtime_error( "Sizes of input *y* and *y_los* are inconsistent." );
1907  if( jacobian_do )
1908  {
1909  nrq1 = jacobian_quantities1.nelem();
1910  if( jacobian.nrows() != n1 )
1911  throw runtime_error( "Sizes of *y* and *jacobian* are inconsistent." );
1912  if( jacobian_indices1.nelem() != nrq1 )
1913  throw runtime_error( "Lengths of *jacobian_quantities_copy* and "
1914  "*jacobian_indices_copy* are inconsistent." );
1915  if( jacobian.ncols() != jacobian_indices1[nrq1-1][1]+1 )
1916  throw runtime_error( "Size of input *jacobian* and max value in "
1917  "*jacobian_indices_copy* are inconsistent." );
1918  }
1919 
1920  // Calculate new measurement
1921  //
1922  Vector y2, y_f2;
1923  Matrix y_pos2, y_los2, jacobian2;
1924  ArrayOfIndex y_pol2;
1925  ArrayOfVector y_aux2;
1926  //
1927  yCalc( ws, y2, y_f2, y_pol2, y_pos2, y_los2, y_aux2, jacobian2,
1928  atmfields_checked, atmgeom_checked, atmosphere_dim, t_field,
1929  z_field, vmr_field, cloudbox_on, cloudbox_checked, sensor_checked,
1930  stokes_dim, f_grid, sensor_pos, sensor_los, transmitter_pos,
1931  mblock_za_grid, mblock_aa_grid, antenna_dim, sensor_response,
1932  sensor_response_f, sensor_response_pol, sensor_response_za,
1933  sensor_response_aa, iy_main_agenda, jacobian_agenda, jacobian_do,
1934  jacobian_quantities, jacobian_indices, iy_aux_vars, verbosity );
1935 
1936  // Consistency checks
1937  if( y_pos.ncols() != y_pos2.ncols() )
1938  throw runtime_error(
1939  "Different number of columns in *y_pos* between the measurements." );
1940  if( y_los.ncols() != y_los2.ncols() )
1941  throw runtime_error(
1942  "Different number of columns in *y_pos* between the measurements." );
1943 
1944 
1945  // y and y_XXX
1946  //
1947  const Index n2 = y2.nelem();
1948  //
1949  {
1950  // Make copy of old measurement
1951  const Vector y1=y, y_f1=y_f;
1952  const Matrix y_pos1=y_pos, y_los1=y_los;
1953  const ArrayOfIndex y_pol1=y_pol;
1954  const ArrayOfVector y_aux1=y_aux;
1955  //
1956  y.resize( n1+n2 );
1957  y[Range(0,n1)] = y1; y[Range(n1,n2)] = y2;
1958  //
1959  y_f.resize( n1+n2 );
1960  y_f[Range(0,n1)] = y_f1; y_f[Range(n1,n2)] = y_f2;
1961  //
1962  y_pos.resize( n1+n2, y_pos1.ncols() );
1963  y_pos(Range(0,n1),joker) = y_pos1; y_pos(Range(n1,n2),joker) = y_pos2;
1964  //
1965  y_los.resize( n1+n2, y_los1.ncols() );
1966  y_los(Range(0,n1),joker) = y_los1; y_los(Range(n1,n2),joker) = y_los2;
1967  //
1968  y_pol.resize( n1+n2 );
1969  for( Index i=0; i<n1; i++ )
1970  { y_pol[i] = y_pol1[i]; }
1971  for( Index i=0; i<n2; i++ )
1972  { y_pol[n1+i] = y_pol2[i]; }
1973 
1974  // y_aux
1975  const Index na1 = y_aux1.nelem();
1976  const Index na2 = y_aux2.nelem();
1977  const Index na = max(na1,na2);
1978  //
1979  y_aux.resize( na );
1980  //
1981  for( Index a=0; a<na; a++ )
1982  {
1983  y_aux[a].resize( n1+n2 );
1984  if( a < na1 )
1985  { y_aux[a][Range(0,n1)] = y_aux1[a]; }
1986  else
1987  { y_aux[a][Range(0,n1)] = 0; }
1988  if( a < na2 )
1989  { y_aux[a][Range(n1,n2)] = y_aux2[a]; }
1990  else
1991  { y_aux[a][Range(n1,n2)] = 0; }
1992  }
1993  }
1994 
1995  // Jacobian and friends
1996  if( jacobian_do )
1997  {
1998  // Put in old jacobian_quantities and jacobian_indices as first part in
1999  // new version of these variables
2000  ArrayOfRetrievalQuantity jacobian_quantities2 = jacobian_quantities;
2001  ArrayOfArrayOfIndex jacobian_indices2 = jacobian_indices;
2002  //
2003  jacobian_quantities = jacobian_quantities1;
2004  jacobian_indices = jacobian_indices1;
2005 
2006  // Loop new jacobian_quantities to determine how new jacobian data shall
2007  // be inserted
2008  //
2009  const Index nrq2 = jacobian_quantities2.nelem();
2010  ArrayOfIndex map_table(nrq2);
2011  //
2012  for( Index q2=0; q2<nrq2; q2++ )
2013  {
2014 
2015  Index pos = -1;
2016 
2017  // Compare to old quantities, to determine if append shall be
2018  // considered. Some special checks performed here, grids checked later
2019  if( jacobian_quantities2[q2].MainTag() == ABSSPECIES_MAINTAG ||
2020  jacobian_quantities2[q2].MainTag() == TEMPERATURE_MAINTAG ||
2021  jacobian_quantities2[q2].MainTag() == WIND_MAINTAG ||
2022  append_instrument_wfs )
2023  {
2024  for( Index q1=0; q1<nrq1; q1++ && pos < 0 )
2025  {
2026  if( jacobian_quantities2[q2].MainTag() ==
2027  jacobian_quantities1[q1].MainTag() )
2028  {
2029  // Absorption species
2030  if( jacobian_quantities2[q2].MainTag() ==
2031  ABSSPECIES_MAINTAG )
2032  {
2033  if( jacobian_quantities2[q2].Subtag() ==
2034  jacobian_quantities1[q1].Subtag() )
2035  {
2036  if( jacobian_quantities2[q2].Mode() ==
2037  jacobian_quantities1[q1].Mode() )
2038  { pos = q1; }
2039  else
2040  {
2041  ostringstream os;
2042  os << "Jacobians for "
2043  << jacobian_quantities2[q2].MainTag()<<"/"
2044  << jacobian_quantities2[q2].Subtag()
2045  << " shall be appended.\nThis requires "
2046  << "that the same retrieval unit is used "
2047  << "but it seems that this requirement is "
2048  << "not met.";
2049  throw runtime_error(os.str());
2050  }
2051  }
2052  }
2053  // Temperature
2054  else if( jacobian_quantities2[q2].MainTag() ==
2055  TEMPERATURE_MAINTAG )
2056  {
2057  if( jacobian_quantities2[q2].Subtag() ==
2058  jacobian_quantities1[q1].Subtag() )
2059  { pos = q1; }
2060  else
2061  {
2062  ostringstream os;
2063  os << "Jacobians for "
2064  << jacobian_quantities2[q2].MainTag()<<"/"
2065  << jacobian_quantities2[q2].Subtag()
2066  << " shall be appended.\nThis requires "
2067  << "that HSE is either ON or OFF for both "
2068  << "parts but it seems that this requirement "
2069  << "is not met.";
2070  throw runtime_error(os.str());
2071  }
2072  }
2073  // Other
2074  else if( jacobian_quantities2[q2].Subtag() ==
2075  jacobian_quantities1[q1].Subtag() )
2076  { pos = q1; }
2077  }
2078  }
2079  }
2080 
2081 
2082  // New quantity
2083  if( pos < 0 )
2084  {
2085  map_table[q2] = jacobian_quantities.nelem();
2086  jacobian_quantities.push_back( jacobian_quantities2[q2] );
2087  ArrayOfIndex indices(2);
2088  indices[0] = jacobian_indices[jacobian_indices.nelem()-1][1]+1;
2089  indices[1] = indices[0] + jacobian_indices2[q2][1] -
2090  jacobian_indices2[q2][0];
2091  jacobian_indices.push_back( indices );
2092  }
2093  // Existing quantity
2094  else
2095  {
2096  map_table[q2] = pos;
2097  // Check if grids are equal
2098  ArrayOfVector grids1 = jacobian_quantities1[pos].Grids();
2099  ArrayOfVector grids2 = jacobian_quantities2[q2].Grids();
2100  bool any_wrong = false;
2101  if( grids1.nelem() != grids2.nelem() )
2102  { any_wrong = true; }
2103  else
2104  {
2105  for( Index g=0; g<grids1.nelem(); g++ )
2106  {
2107  if( grids1[g].nelem() != grids2[g].nelem() )
2108  { any_wrong = true; }
2109  else
2110  {
2111  for( Index e=0; e<grids1[g].nelem(); e++ )
2112  {
2113  const Numeric v1 = grids1[g][e];
2114  const Numeric v2 = grids2[g][e];
2115  if( ( v1 == 0 && abs(v2) > 1e-9 ) ||
2116  abs(v1-v2)/v1 > 1e-6 )
2117  { any_wrong = true; }
2118  }
2119  }
2120  }
2121  }
2122  if( any_wrong )
2123  {
2124  ostringstream os;
2125  os << "Jacobians for "
2126  << jacobian_quantities2[q2].MainTag()
2127  << "/" << jacobian_quantities2[q2].Subtag()
2128  << " shall be appended.\nThis requires that the "
2129  << "same grids are used for both measurements,\nbut "
2130  << "it seems that this requirement is not met.";
2131  throw runtime_error(os.str());
2132  }
2133  }
2134  }
2135 
2136  // Create and fill *jacobian*
2137  //
2138  const Index nrq = jacobian_quantities.nelem();
2139  const Matrix jacobian1 = jacobian;
2140  //
2141  jacobian.resize( n1+n2, jacobian_indices[nrq-1][1]+1 );
2142  jacobian = 0;
2143  //
2144  // Put in old part in top-left corner
2145  jacobian(Range(0,n1),Range(0,jacobian_indices1[nrq1-1][1]+1)) = jacobian1;
2146  // New parts
2147  for( Index q2=0; q2<nrq2; q2++ )
2148  {
2149  jacobian(Range(n1,n2),Range(jacobian_indices[map_table[q2]][0],
2150  jacobian_indices[map_table[q2]][1]-
2151  jacobian_indices[map_table[q2]][0]+1) ) =
2152  jacobian2(joker,Range(jacobian_indices2[q2][0],
2153  jacobian_indices2[q2][1]-
2154  jacobian_indices2[q2][0]+1) );
2155  }
2156  }
2157 }
2158 
2159 
2160 
2161 /* Workspace method: Doxygen documentation will be auto-generated */
2162 void yApplyUnit(
2163  Vector& y,
2164  Matrix& jacobian,
2165  const Vector& y_f,
2166  const ArrayOfIndex& y_pol,
2167  const String& iy_unit,
2168  const Verbosity&)
2169 {
2170  if( iy_unit == "1" )
2171  { throw runtime_error(
2172  "No need to use this method with *iy_unit* = \"1\"." ); }
2173 
2174  if( max(y) > 1e-3 )
2175  {
2176  ostringstream os;
2177  os << "The spectrum vector *y* is required to have original radiance\n"
2178  << "unit, but this seems not to be the case. This as a value above\n"
2179  << "1e-3 is found in *y*.";
2180  throw runtime_error( os.str() );
2181  }
2182 
2183  // Is jacobian set?
2184  //
2185  const Index ny = y.nelem();
2186  //
2187  const bool do_j = jacobian.nrows() == ny;
2188 
2189  // Some jacobian quantities can not be handled
2190  if( do_j && max(jacobian) > 1e-3 )
2191  {
2192  ostringstream os;
2193  os << "The method can not be used with jacobian quantities that are not\n"
2194  << "obtained through radiative transfer calculations. One example on\n"
2195  << "quantity that can not be handled is *jacobianAddPolyfit*.\n"
2196  << "The maximum value of *jacobian* indicates that one or several\n"
2197  << "such jacobian quantities are included.";
2198  throw runtime_error( os.str() );
2199  }
2200 
2201  // Planck-Tb
2202  //--------------------------------------------------------------------------
2203  if( iy_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_iy_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 of 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_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  // Jacobian
2250  if( do_j )
2251  {
2252  Tensor3 J(jacobian.ncols(),1,n);
2253  J(joker,0,joker) = transpose( jacobian(ii,joker) );
2254  apply_iy_unit2( J, yv, iy_unit, y_f[i0], 1, i_pol );
2255  jacobian(ii,joker) = transpose( J(joker,0,joker) );
2256  }
2257  }
2258 
2259  // y (must be done last)
2260  apply_iy_unit( yv, iy_unit, y_f[i0], 1, i_pol );
2261  y[ii] = yv(0,joker);
2262 
2263  i0 += n;
2264  }
2265  }
2266 
2267 
2268  // Other conversions
2269  //--------------------------------------------------------------------------
2270  else
2271  {
2272  // Here we take each element of y separately.
2273 
2274  Matrix yv(1,1);
2275  ArrayOfIndex i_pol(1);
2276 
2277  for( Index i=0; i<ny; i++ )
2278  {
2279  yv(0,0) = y[i];
2280  i_pol[0] = y_pol[i];
2281 
2282  // Jacobian
2283  if( do_j )
2284  { apply_iy_unit2( MatrixView( jacobian(i,joker) ), yv,
2285  iy_unit, y_f[i], 1, i_pol ); }
2286 
2287  // y (must be done last)
2288  apply_iy_unit( yv, iy_unit, y_f[i], 1, i_pol );
2289  y[i] = yv(0,0);
2290  }
2291  }
2292 }
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
Range::get_start
Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:196
ppath_agendaExecute
void ppath_agendaExecute(Workspace &ws, Ppath &ppath, const Numeric ppath_lraytrace, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index cloudbox_on, const Index ppath_inside_cloudbox_do, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:14625
Ppath::gp_lat
ArrayOfGridPos gp_lat
Definition: ppath.h:77
find_first
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:281
MCAntenna::set_pencil_beam
void set_pencil_beam(void)
makes the antenna pattern a pencil beam
Definition: mc_antenna.cc:88
chk_if_std_blackbody_agenda
bool chk_if_std_blackbody_agenda(Workspace &ws, const Agenda &blackbody_radiation_agenda)
Checks if blackbody_radiation_agenda returns frequency based radiance.
Definition: check_input.cc:2149
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 &ppath_step_agenda, const Numeric &ppath_lraytrace, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, 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_array_mono, const Index &atmfields_checked, const Index &atmgeom_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 Index &min_iter, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
Definition: m_montecarlo.cc:103
get_ppath_f
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, ConstVectorView f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, ConstMatrixView ppath_wind)
get_ppath_f
Definition: rte.cc:1690
iyMC
void iyMC(Workspace &ws, Matrix &iy, ArrayOfTensor4 &iy_aux, ArrayOfTensor3 &diy_dx, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const ArrayOfString &iy_aux_vars, 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 Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const ArrayOfSingleScatteringData &scat_data_array, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &ppath_step_agenda, const Numeric &ppath_lraytrace, const Tensor4 &pnd_field, const String &iy_unit, const Numeric &mc_std_err, const Index &mc_max_time, const Index &mc_max_iter, const Index &mc_min_iter, const Verbosity &verbosity)
WORKSPACE METHOD: iyMC.
Definition: m_rte.cc:1131
Agenda
The Agenda class.
Definition: agenda_class.h:44
transpose
ConstMatrixView transpose(ConstMatrixView m)
Const version of transpose.
Definition: matpackI.cc:1799
iy_transmission_mult
void iy_transmission_mult(Tensor3 &iy_trans_total, ConstTensor3View iy_trans_old, ConstTensor3View iy_trans_new)
iy_transmission_mult
Definition: rte.cc:2309
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:176
iy_auxFillParticleVariables
void iy_auxFillParticleVariables(ArrayOfTensor4 &iy_aux, const Index &atmfields_checked, const Index &cloudbox_checked, const Index &atmosphere_dim, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Matrix &particle_masses, const Ppath &ppath, const ArrayOfString &iy_aux_vars, const Verbosity &)
WORKSPACE METHOD: iy_auxFillParticleVariables.
Definition: m_rte.cc:1354
iy_sub_agendaExecute
void iy_sub_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfTensor4 &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmission, const ArrayOfString &iy_aux_vars, const Index cloudbox_on, const Index jacobian_do, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition: auto_md.cc:13960
messages.h
Declarations having to do with the four output streams.
q
#define q
Definition: continua.cc:21469
iyb_calc
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, const Index &mblock_index, const Index &atmosphere_dim, 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, ConstMatrixView transmitter_pos, ConstVectorView mblock_za_grid, ConstVectorView mblock_aa_grid, const Index &antenna_dim, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
iyb_calc
Definition: rte.cc:2101
chk_if_increasing
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:132
jacobian_agendaExecute
void jacobian_agendaExecute(Workspace &ws, Matrix &jacobian, const Index mblock_index, const Vector &iyb, const Vector &yb, const Agenda &input_agenda)
Definition: auto_md.cc:14245
jacobian.h
Declarations required for the calculation of jacobians.
Vector
The Vector class.
Definition: matpackI.h:556
MatrixView
The MatrixView class.
Definition: matpackI.h:679
apply_iy_unit
void apply_iy_unit(MatrixView iy, const String &y_unit, ConstVectorView f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
apply_iy_unit
Definition: rte.cc:127
Sparse
The Sparse class.
Definition: matpackII.h:55
Tensor4
The Tensor4 class.
Definition: matpackIV.h:383
Range
The range class.
Definition: matpackI.h:148
Ppath::lstep
Vector lstep
Definition: ppath.h:70
arts_omp_get_max_threads
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:47
get_ppath_blackrad
void get_ppath_blackrad(Workspace &ws, Matrix &ppath_blackrad, const Agenda &blackbody_radiation_agenda, const Ppath &ppath, ConstVectorView ppath_t, ConstMatrixView ppath_f)
get_ppath_blackrad
Definition: rte.cc:1487
WIND_MAINTAG
const String WIND_MAINTAG
get_ppath_abs
void get_ppath_abs(Workspace &ws, Tensor4 &ppath_abs, Tensor5 &abs_per_species, const Agenda &propmat_clearsky_agenda, const Ppath &ppath, ConstVectorView ppath_p, ConstVectorView ppath_t, ConstMatrixView ppath_vmr, ConstMatrixView ppath_f, ConstMatrixView ppath_mag, ConstVectorView f_grid, const Index &stokes_dim, const ArrayOfIndex &ispecies)
get_ppath_abs
Definition: rte.cc:1351
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:114
TEMPERATURE_MAINTAG
const String TEMPERATURE_MAINTAG
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
GridPos
Structure to store a grid position.
Definition: interpolation.h:74
ArrayOfIndex
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:40
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
Ppath::end_lstep
Numeric end_lstep
Definition: ppath.h:73
iy_main_agendaExecute
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfTensor4 &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmission, const ArrayOfString &iy_aux_vars, const Index cloudbox_on, const Index jacobian_do, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition: auto_md.cc:13776
get_ppath_atmvars
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, 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, ConstTensor3View mag_u_field, ConstTensor3View mag_v_field, ConstTensor3View mag_w_field)
get_ppath_atmvars
Definition: rte.cc:1233
mult
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
my_basic_string
The implementation for String, the ARTS string class.
Definition: mystring.h:63
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
get_pointers_for_analytical_jacobians
void get_pointers_for_analytical_jacobians(ArrayOfIndex &abs_species_i, ArrayOfIndex &is_t, ArrayOfIndex &wind_i, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:292
Tensor3
The Tensor3 class.
Definition: matpackIII.h:348
arts.h
The global header file for ARTS.
vmrunitscf
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
vmrunitscf
Definition: jacobian.cc:853
PI
const Numeric PI
scat_data_array_monoCalc
void scat_data_array_monoCalc(ArrayOfSingleScatteringData &scat_data_array_mono, const ArrayOfSingleScatteringData &scat_data_array, const Vector &f_grid, const Index &f_index, const Verbosity &)
WORKSPACE METHOD: scat_data_array_monoCalc.
Definition: m_optproperties.cc:1095
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
max
#define max(a, b)
Definition: continua.cc:20461
math_funcs.h
interp_cloudfield_gp2itw
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
interp_cloudfield_gp2itw
Definition: special_interp.cc:301
Sparse::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:56
iyReplaceFromAux
void iyReplaceFromAux(Matrix &iy, const ArrayOfTensor4 &iy_aux, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const String &aux_var, const Verbosity &)
WORKSPACE METHOD: iyReplaceFromAux.
Definition: m_rte.cc:1295
iyCalc
void iyCalc(Workspace &ws, Matrix &iy, ArrayOfTensor4 &iy_aux, Ppath &ppath, const Index &atmfields_checked, const Index &atmgeom_checked, const ArrayOfString &iy_aux_vars, const Vector &f_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &iy_main_agenda, const Verbosity &)
WORKSPACE METHOD: iyCalc.
Definition: m_rte.cc:121
iyApplyUnit
void iyApplyUnit(Matrix &iy, ArrayOfTensor4 &iy_aux, const Index &stokes_dim, const Vector &f_grid, const ArrayOfString &iy_aux_vars, const String &iy_unit, const Verbosity &)
WORKSPACE METHOD: iyApplyUnit.
Definition: m_rte.cc:73
yCalc_mblock_loop_body
void yCalc_mblock_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfVector &iyb_aux_array, Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Matrix &jacobian, const Index &atmosphere_dim, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, 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_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity, const Index &mblock_index, const Index &n1y, const Index &j_analytical_do)
Definition: m_rte.cc:1471
abs
#define abs(x)
Definition: continua.cc:20458
joker
const Joker joker
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Matrix
The Matrix class.
Definition: matpackI.h:788
Ppath::nreal
Vector nreal
Definition: ppath.h:74
auto_md.h
MCAntenna
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:56
gridpos_copy
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Definition: interpolation.cc:517
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
special_interp.h
Header file for special_interp.cc.
ppath.h
Propagation path structure and functions.
q1
#define q1
Definition: continua.cc:21197
logic.h
Header file for logic.cc.
emission_rtstep
void emission_rtstep(Matrix &iy, const Index &stokes_dim, ConstVectorView bbar, ArrayOfIndex &extmat_case, ConstTensor3View t)
emission_rtstep
Definition: rte.cc:822
Array
This can be used to make arrays out of anything.
Definition: array.h:40
MCSetSeedFromTime
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
Definition: m_montecarlo.cc:448
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:101
yApplyUnit
void yApplyUnit(Vector &y, Matrix &jacobian, const Vector &y_f, const ArrayOfIndex &y_pol, const String &iy_unit, const Verbosity &)
WORKSPACE METHOD: yApplyUnit.
Definition: m_rte.cc:2162
yCalc
void yCalc(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &jacobian, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, 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_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
WORKSPACE METHOD: yCalc.
Definition: m_rte.cc:1583
Vector::resize
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
ns
#define ns
Definition: continua.cc:21931
min
#define min(a, b)
Definition: continua.cc:20460
apply_iy_unit2
void apply_iy_unit2(Tensor3View J, ConstMatrixView iy, const String &y_unit, ConstVectorView f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
apply_iy_unit2
Definition: rte.cc:232
Ppath::np
Index np
Definition: ppath.h:61
Ppath::gp_lon
ArrayOfGridPos gp_lon
Definition: ppath.h:78
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: jacobian.cc:110
get_rowindex_for_mblock
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
get_rowindex_for_mblock
Definition: rte.cc:1960
Workspace
Workspace class.
Definition: workspace_ng.h:47
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:214
get_iy_of_background
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, ConstTensor3View iy_transmission, const Index &jacobian_do, const Ppath &ppath, ConstVectorView rte_pos2, const Index &atmosphere_dim, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, const Index &cloudbox_on, const Index &stokes_dim, ConstVectorView f_grid, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Verbosity &verbosity)
get_iy_of_background
Definition: rte.cc:1106
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
arts_omp.h
Header file for helper functions for OpenMP.
iyLoopFrequencies
void iyLoopFrequencies(Workspace &ws, Matrix &iy, ArrayOfTensor4 &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const ArrayOfString &iy_aux_vars, const Index &stokes_dim, const Vector &f_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index &jacobian_do, const Agenda &iy_sub_agenda, const Verbosity &)
WORKSPACE METHOD: iyLoopFrequencies.
Definition: m_rte.cc:1054
is_gp_inside_cloudbox
bool is_gp_inside_cloudbox(const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, const bool &include_boundaries, const Index &atmosphere_dim)
Definition: cloudbox.cc:808
get_ppath_trans
void get_ppath_trans(Tensor4 &trans_partial, ArrayOfArrayOfIndex &extmat_case, Tensor4 &trans_cumulat, Vector &scalar_tau, const Ppath &ppath, ConstTensor4View &ppath_abs, ConstVectorView f_grid, const Index &stokes_dim)
get_ppath_trans
Definition: rte.cc:1770
iyEmissionStandard
void iyEmissionStandard(Workspace &ws, Matrix &iy, ArrayOfTensor4 &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Agenda &ppath_agenda, const Agenda &blackbody_radiation_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Numeric &rte_alonglos_v, const Numeric &ppath_lraytrace, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandard.
Definition: m_rte.cc:170
physics_funcs.h
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
ABSSPECIES_MAINTAG
const String ABSSPECIES_MAINTAG
Ppath::gp_p
ArrayOfGridPos gp_p
Definition: ppath.h:76
interp_atmfield_by_itw
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
interp_atmfield_by_itw
Definition: special_interp.cc:159
Tensor4::resize
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1403
ext2trans
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Definition: rte.cc:903
Tensor5
The Tensor5 class.
Definition: matpackV.h:451
rte.h
Declaration of functions in rte.cc.
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
yCalcAppend
void yCalcAppend(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &jacobian, ArrayOfRetrievalQuantity &jacobian_quantities, ArrayOfArrayOfIndex &jacobian_indices, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, 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_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfString &iy_aux_vars, const ArrayOfRetrievalQuantity &jacobian_quantities1, const ArrayOfArrayOfIndex &jacobian_indices1, const Index &append_instrument_wfs, const Verbosity &verbosity)
WORKSPACE METHOD: yCalcAppend.
Definition: m_rte.cc:1852
check_input.h