ARTS  2.2.66
m_fos.cc
Go to the documentation of this file.
1 /* Copyright (C) 2013
2  Patrick Eriksson <patrick.eriksson@chalmers.se>
3 
4  This program is free software; you can redistribute it and/or modify it
5  under the terms of the GNU General Public License as published by the
6  Free Software Foundation; either version 2, or (at your option) any
7  later version.
8 
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with this program; if not, write to the Free Software
16  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17  USA. */
18 
19 
20 
21 /*===========================================================================
22  === File description
23  ===========================================================================*/
24 
38 /*===========================================================================
39  === External declarations
40  ===========================================================================*/
41 
42 #include <cmath>
43 #include <stdexcept>
44 #include "arts.h"
45 #include "arts_omp.h"
46 #include "auto_md.h"
47 #include "doit.h"
48 #include "math_funcs.h"
49 #include "montecarlo.h"
50 #include "rte.h"
51 
52 extern const Numeric DEG2RAD;
53 extern const Numeric RAD2DEG;
54 extern const Numeric PI;
55 
56 
57 // FOS implemented as an internal function, to allow an recursive algorithm
58 void fos(
59  Workspace& ws,
60  Matrix& iy,
61  ArrayOfTensor4& iy_aux,
62  Ppath& ppath,
63  ArrayOfTensor3& diy_dx,
64  const Index& stokes_dim,
65  const Vector& f_grid,
66  const Index& atmosphere_dim,
67  const Vector& p_grid,
68  const Tensor3& z_field,
69  const Tensor3& t_field,
70  const Tensor4& vmr_field,
71  const ArrayOfArrayOfSpeciesTag& abs_species,
72  const Tensor3& wind_u_field,
73  const Tensor3& wind_v_field,
74  const Tensor3& wind_w_field,
75  const Tensor3& mag_u_field,
76  const Tensor3& mag_v_field,
77  const Tensor3& mag_w_field,
78  const Index& cloudbox_on,
79  const ArrayOfIndex& cloudbox_limits,
80  const Tensor4& pnd_field,
81  const Index& use_mean_scat_data,
82  const ArrayOfSingleScatteringData& scat_data_array,
83  const Matrix& particle_masses,
84  const String& iy_unit,
85  const ArrayOfString& iy_aux_vars,
86  const Index& jacobian_do,
87  const Agenda& ppath_agenda,
88  const Agenda& blackbody_radiation_agenda,
89  const Agenda& propmat_clearsky_agenda,
90  const Agenda& iy_main_agenda,
91  const Agenda& iy_space_agenda,
92  const Agenda& iy_surface_agenda,
93  const Index& iy_agenda_call1,
94  const Tensor3& iy_transmission,
95  const Vector& rte_pos,
96  const Vector& rte_los,
97  const Vector& rte_pos2,
98  const Numeric& rte_alonglos_v,
99  const Numeric& ppath_lraytrace,
100  const Matrix& fos_scatint_angles,
101  const Vector& fos_iyin_za_angles,
102  const Index& fos_za_interporder,
103  const Index& fos_n,
104  const Index& fos_i,
105  const Verbosity& verbosity )
106 {
107  // A temporary restriction
108  if( atmosphere_dim > 1 )
109  throw runtime_error( "FOS is so far only handling 1D atmospheres." );
110 
111  assert( fos_i >= 0 && fos_i <= fos_n );
112 
113 
114  // Determine propagation path
115  //
116  ppath_agendaExecute( ws, ppath, ppath_lraytrace, rte_pos, rte_los, rte_pos2,
117  0, 0, t_field, z_field, vmr_field, f_grid,
118  ppath_agenda );
119 
120  // Some basic sizes
121  //
122  const Index nf = f_grid.nelem();
123  const Index ns = stokes_dim;
124  const Index np = ppath.np;
125 
126  // The below copied from iyEmission. Activate for-loop when jacobians
127  // introduced.
128 
129  // Set up variable with index of species where we want abs_per_species.
130  // This variable can below be extended in iy_aux part.
131  //
132  ArrayOfIndex iaps(0);
133  //
134  //for( Index i=0; i<jac_species_i.nelem(); i++ )
135  // {
136  // if( jac_species_i[i] >= 0 )
137  // { iaps.push_back( jac_species_i[i] ); }
138  // }
139 
140  //=== iy_aux part ===========================================================
141  Index auxPressure = -1,
142  auxTemperature = -1,
143  auxAbsSum = -1,
144  auxBackground = -1,
145  auxIy = -1,
146  auxOptDepth = -1;
147  ArrayOfIndex auxAbsSpecies(0), auxAbsIsp(0);
148  ArrayOfIndex auxVmrSpecies(0), auxVmrIsp(0);
149  ArrayOfIndex auxPartCont(0), auxPartContI(0);
150  ArrayOfIndex auxPartField(0), auxPartFieldI(0);
151  //
152  if( !iy_agenda_call1 )
153  { iy_aux.resize( 0 ); }
154  else
155  {
156  const Index naux = iy_aux_vars.nelem();
157  iy_aux.resize( naux );
158  //
159  for( Index i=0; i<naux; i++ )
160  {
161  if( iy_aux_vars[i] == "Pressure" )
162  { auxPressure = i; iy_aux[i].resize( 1, 1, 1, np ); }
163  else if( iy_aux_vars[i] == "Temperature" )
164  { auxTemperature = i; iy_aux[i].resize( 1, 1, 1, np ); }
165  else if( iy_aux_vars[i].substr(0,13) == "VMR, species " )
166  {
167  Index ispecies;
168  istringstream is(iy_aux_vars[i].substr(13,2));
169  is >> ispecies;
170  if( ispecies < 0 || ispecies>=abs_species.nelem() )
171  {
172  ostringstream os;
173  os << "You have selected VMR of species with index "
174  << ispecies << ".\nThis species does not exist!";
175  throw runtime_error( os.str() );
176  }
177  auxVmrSpecies.push_back(i);
178  auxVmrIsp.push_back(ispecies);
179  iy_aux[i].resize( 1, 1, 1, np );
180  }
181  else if( iy_aux_vars[i] == "Absorption, summed" )
182  { auxAbsSum = i; iy_aux[i].resize( nf, ns, ns, np ); }
183  else if( iy_aux_vars[i].substr(0,20) == "Absorption, species " )
184  {
185  Index ispecies;
186  istringstream is(iy_aux_vars[i].substr(20,2));
187  is >> ispecies;
188  if( ispecies < 0 || ispecies>=abs_species.nelem() )
189  {
190  ostringstream os;
191  os << "You have selected absorption species with index "
192  << ispecies << ".\nThis species does not exist!";
193  throw runtime_error( os.str() );
194  }
195  auxAbsSpecies.push_back(i);
196  const Index ihit = find_first( iaps, ispecies );
197  if( ihit >= 0 )
198  { auxAbsIsp.push_back( ihit ); }
199  else
200  {
201  iaps.push_back(ispecies);
202  auxAbsIsp.push_back( iaps.nelem()-1 );
203  }
204  iy_aux[i].resize( nf, ns, ns, np );
205  }
206  else if( iy_aux_vars[i] == "Radiative background" )
207  { auxBackground = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
208  else if( iy_aux_vars[i] == "iy" && auxIy < 0 )
209  { auxIy = i; iy_aux[i].resize( nf, ns, 1, np ); }
210  else if( iy_aux_vars[i] == "Optical depth" )
211  { auxOptDepth = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
212  else if( iy_aux_vars[i].substr(0,14) == "Mass content, " )
213  {
214  Index icont;
215  istringstream is(iy_aux_vars[i].substr(14,2));
216  is >> icont;
217  if( icont < 0 || icont>=particle_masses.ncols() )
218  {
219  ostringstream os;
220  os << "You have selected particle mass content category with "
221  << "index " << icont << ".\nThis category is not defined!";
222  throw runtime_error( os.str() );
223  }
224  auxPartCont.push_back(i);
225  auxPartContI.push_back(icont);
226  iy_aux[i].resize( 1, 1, 1, np );
227  }
228  else if( iy_aux_vars[i].substr(0,10) == "PND, type " )
229  {
230  Index ip;
231  istringstream is(iy_aux_vars[i].substr(10,2));
232  is >> ip;
233  if( ip < 0 || ip>=pnd_field.nbooks() )
234  {
235  ostringstream os;
236  os << "You have selected particle number density field with "
237  << "index " << ip << ".\nThis field is not defined!";
238  throw runtime_error( os.str() );
239  }
240  auxPartField.push_back(i);
241  auxPartFieldI.push_back(ip);
242  iy_aux[i].resize( 1, 1, 1, np );
243  }
244  else
245  {
246  ostringstream os;
247  os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
248  << "\"\nThis choice is not recognised.";
249  throw runtime_error( os.str() );
250  }
251  }
252  }
253  //===========================================================================
254 
255  // Get atmospheric and RT quantities for each ppath point/step
256  //
257  Vector ppath_p, ppath_t;
258  Matrix ppath_vmr, ppath_pnd, ppath_wind, ppath_mag, ppath_f;
259  Matrix ppath_blackrad;
260  Tensor5 abs_per_species;
261  Tensor4 ppath_abs, trans_partial, trans_cumulat, pnd_ext_mat;
262  Tensor3 pnd_abs_vec;
263  Vector scalar_tau;
264  ArrayOfIndex clear2cloudbox;
265  //
267  ArrayOfArrayOfIndex extmat_case;
268  //
269  if( np > 1 )
270  {
271  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
272  ppath_wind, ppath_mag,
273  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
274  wind_u_field, wind_v_field, wind_w_field,
275  mag_u_field, mag_v_field, mag_w_field );
276  get_ppath_f( ppath_f, ppath, f_grid, atmosphere_dim,
277  rte_alonglos_v, ppath_wind );
278  get_ppath_abs( ws, ppath_abs, abs_per_species,
279  propmat_clearsky_agenda, ppath,
280  ppath_p, ppath_t, ppath_vmr, ppath_f,
281  ppath_mag, f_grid, stokes_dim, iaps );
282  get_ppath_blackrad( ws, ppath_blackrad, blackbody_radiation_agenda,
283  ppath, ppath_t, ppath_f );
284  if( !cloudbox_on )
285  {
286  get_ppath_trans( trans_partial, extmat_case, trans_cumulat,
287  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim );
288  }
289  else
290  {
291  get_ppath_ext( clear2cloudbox, pnd_abs_vec, pnd_ext_mat, scat_data,
292  ppath_pnd, ppath, ppath_t, stokes_dim, ppath_f,
293  atmosphere_dim, cloudbox_limits, pnd_field,
294  use_mean_scat_data, scat_data_array, verbosity );
295  get_ppath_trans2( trans_partial, extmat_case, trans_cumulat,
296  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim,
297  clear2cloudbox, pnd_ext_mat );
298  }
299  }
300  else // For cases totally outside the atmosphere,
301  { // set zero optical thickness and unit transmission
302  scalar_tau.resize( nf );
303  scalar_tau = 0;
304  trans_cumulat.resize( nf, ns, ns, np );
305  for( Index iv=0; iv<nf; iv++ )
306  { id_mat( trans_cumulat(iv,joker,joker,np-1) ); }
307  }
308 
309  // iy_transmission
310  //
311  Tensor3 iy_trans_new;
312  //
313  if( iy_agenda_call1 )
314  { iy_trans_new = trans_cumulat(joker,joker,joker,np-1); }
315  else
316  {
317  iy_transmission_mult( iy_trans_new, iy_transmission,
318  trans_cumulat(joker,joker,joker,np-1) ); }
319 
320  // Radiative background
321  //
322  {
323  Agenda iy_cbox_agenda;
324  get_iy_of_background( ws, iy, diy_dx,
325  iy_trans_new, jacobian_do, ppath, rte_pos2,
326  atmosphere_dim, t_field, z_field, vmr_field,
327  cloudbox_on, stokes_dim, f_grid, iy_main_agenda,
328  iy_space_agenda, iy_surface_agenda, iy_cbox_agenda,
329  verbosity );
330  }
331 
332  //=== iy_aux part ===========================================================
333  // Fill parts of iy_aux that are defined even for np=1.
334  // Radiative background
335  if( auxBackground >= 0 )
336  { iy_aux[auxBackground](0,0,0,0) = (Numeric)min( (Index)2,
337  ppath_what_background(ppath)-1); }
338  // Radiance
339  if( auxIy >= 0 )
340  { iy_aux[auxIy](joker,joker,0,np-1) = iy; }
341  // Transmission, total
342  if( auxOptDepth >= 0 )
343  { iy_aux[auxOptDepth](joker,0,0,0) = scalar_tau; }
344  //===========================================================================
345 
346 
347  // Do RT calculations
348  //
349  if( np > 1 )
350  {
351  //=== iy_aux part =======================================================
352  // iy_aux for point np-1:
353  // Pressure
354  if( auxPressure >= 0 )
355  { iy_aux[auxPressure](0,0,0,np-1) = ppath_p[np-1]; }
356  // Temperature
357  if( auxTemperature >= 0 )
358  { iy_aux[auxTemperature](0,0,0,np-1) = ppath_t[np-1]; }
359  // VMR
360  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
361  { iy_aux[auxVmrSpecies[j]](0,0,0,np-1) = ppath_vmr(auxVmrIsp[j],np-1); }
362  // Absorption
363  if( auxAbsSum >= 0 )
364  { for( Index iv=0; iv<nf; iv++ ) {
365  for( Index is1=0; is1<ns; is1++ ){
366  for( Index is2=0; is2<ns; is2++ ){
367  iy_aux[auxAbsSum](iv,is1,is2,np-1) =
368  ppath_abs(iv,is1,is2,np-1); } } } }
369  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
370  { for( Index iv=0; iv<nf; iv++ ) {
371  for( Index is1=0; is1<stokes_dim; is1++ ){
372  for( Index is2=0; is2<stokes_dim; is2++ ){
373  iy_aux[auxAbsSpecies[j]](iv,is1,is2,np-1) =
374  abs_per_species(auxAbsIsp[j],iv,is1,is2,np-1); } } } }
375  // Particle properties
376  if( cloudbox_on )
377  {
378  // Particle mass content
379  for( Index j=0; j<auxPartCont.nelem(); j++ )
380  { iy_aux[auxPartCont[j]](0,0,0,np-1) = ppath_pnd(joker,np-1) *
381  particle_masses(joker,auxPartContI[j]); }
382  // Particle field
383  for( Index j=0; j<auxPartField.nelem(); j++ )
384  { iy_aux[auxPartField[j]](0,0,0,np-1) =
385  ppath_pnd(auxPartFieldI[j],np-1); }
386  }
387  // Radiance for this point is handled above
388  //=======================================================================
389 
390 
391  // Scattering source term at ip (0) and ip+1 (1)
392  // (If any particles at ip=np-1, this is ignored. Could happen if
393  // particles at surface level.)
394  Matrix ssource0(nf,ns,0), ssource1(nf,ns);
395 
396  // Help variables for handling of *use_mean_scat_data*
397  Index nfs=nf, ivf=1;
398  if( use_mean_scat_data )
399  { nfs = 1; ivf = 0; }
400 
401  // Loop ppath steps
402  for( Index ip=np-2; ip>=0; ip-- )
403  {
404  // Path step average of emission source function: Bbar
405  Vector bbar(nf);
406  for( Index iv=0; iv<nf; iv++ )
407  { bbar[iv] = 0.5 * ( ppath_blackrad(iv,ip) +
408  ppath_blackrad(iv,ip+1) ); }
409 
410  // Check if any particles to consider
411  bool any_particles = clear2cloudbox[ip] >= 0 ||
412  clear2cloudbox[ip+1] >= 0;
413 
414  // -----------------------------------------------------------------
415  // i = N (only absorption/emission)
416  // -----------------------------------------------------------------
417  if( fos_i == fos_n )
418  {
419  // No particle absorption to consider
420  if( !any_particles )
421  {
422  emission_rtstep( iy, stokes_dim, bbar, extmat_case[ip],
423  trans_partial(joker,joker,joker,ip) );
424  }
425 
426  else // We want to include particle absorption, but not
427  { // extinction. trans_partial is then not valid.
428  Tensor3 t(nf,ns,ns);
429  ArrayOfIndex extmat_cas2(nf);
430  //
431  for( Index iv=0; iv<nf; iv++ )
432  {
433  // Particle absorption
434  //
435  Matrix pabs_mat(ns,ns,0);
436  //
437  if( clear2cloudbox[ip] >= 0 )
438  { ext_matFromabs_vec( pabs_mat, pnd_abs_vec(iv,joker,
439  clear2cloudbox[ip]), stokes_dim ); }
440  if( clear2cloudbox[ip+1] >= 0 )
441  { ext_matFromabs_vec( pabs_mat, pnd_abs_vec(iv,joker,
442  clear2cloudbox[ip+1]), stokes_dim ); }
443 
444  // Transmission of step
445  Matrix ext_mat(stokes_dim,stokes_dim);
446  for( Index is1=0; is1<stokes_dim; is1++ ) {
447  for( Index is2=0; is2<stokes_dim; is2++ ) {
448  ext_mat(is1,is2) = 0.5 * ( pabs_mat(is1,is2) +
449  ppath_abs(iv,is1,is2,ip) +
450  ppath_abs(iv,is1,is2,ip+1) );
451  } }
452  //
453  extmat_cas2[iv] = 0;
454  ext2trans( t(iv,joker,joker), extmat_cas2[iv],
455  ext_mat, ppath.lstep[ip] );
456  }
457 
458  // Perform RT
459  emission_rtstep( iy, stokes_dim, bbar, extmat_cas2, t );
460  }
461  }
462 
463 
464  // -----------------------------------------------------------------
465  // i < N
466  // -----------------------------------------------------------------
467  else
468  {
469  // Shift scattering source term (new 1 is old 0)
470  ssource1 = ssource0;
471 
472  // Clear-sky path step
473  if( !any_particles )
474  {
475  // Perform RT
476  emission_rtstep( iy, stokes_dim, bbar, extmat_case[ip],
477  trans_partial(joker,joker,joker,ip) );
478 
479  // Scattering source term at ip is zero:
480  ssource0 = 0;
481  }
482 
483  // Include scattering
484  else
485  {
486  // Determine scattering source term at ip
487  if( clear2cloudbox[ip] < 0 )
488  { ssource0 = 0; }
489  else
490  {
491  // Present position
492  // (Note that the Ppath positions (ppath.pos) for 1D have
493  // one column more than expected by most functions. Only
494  // the first atmosphere_dim values shall be copied.)
495  Vector pos = ppath.pos(ip,Range(0,atmosphere_dim));
496 
497  // Determine incoming radiation
498  //
499  const Index nin = fos_scatint_angles.nrows();
500  Tensor3 Y(nin,nf,ns);
501  {
502  // Do RT
503  const Index nza = fos_iyin_za_angles.nelem();
504  Tensor3 Y1(nza,nf,ns);
505  //
506  for( Index i=0; i<nza; i++ )
507  {
508  // LOS
509  Vector los( 1, fos_iyin_za_angles[i] );
510 
511  // Call recursively, with fos_i increased
512  //
513  Matrix iyl;
514  ArrayOfTensor4 iy_auxl;
515  Ppath ppathl;
516  ArrayOfTensor3 diy_dxl;
517  //
518  fos( ws, iyl, iy_auxl, ppathl, diy_dxl,
519  stokes_dim, f_grid, atmosphere_dim,
520  p_grid, z_field, t_field, vmr_field,
521  abs_species, wind_u_field, wind_v_field,
522  wind_w_field, mag_u_field, mag_v_field,
523  mag_w_field, cloudbox_on, cloudbox_limits,
524  pnd_field, use_mean_scat_data, scat_data_array,
525  particle_masses, iy_unit, iy_aux_vars,
526  jacobian_do, ppath_agenda,
527  blackbody_radiation_agenda,
528  propmat_clearsky_agenda, iy_main_agenda,
529  iy_space_agenda, iy_surface_agenda, 0,
530  iy_trans_new, pos, los, rte_pos2,
531  rte_alonglos_v, ppath_lraytrace,
532  fos_scatint_angles, fos_iyin_za_angles,
533  fos_za_interporder, fos_n, fos_i+1,
534  verbosity );
535 
536  Y1(i,joker,joker) = iyl;
537  }
538 
539  // Zenith angle interpolation of Y
540  ArrayOfGridPosPoly gp( nin );
541  gridpos_poly( gp, fos_iyin_za_angles,
542  fos_scatint_angles(joker,0),
543  fos_za_interporder );
544  Matrix itw( nin, fos_za_interporder+1 );
545  interpweights( itw, gp );
546  //
547  for( Index iv=0; iv<nf; iv++ )
548  {
549  for( Index is1=0; is1<stokes_dim; is1++ )
550  {
551  interp( Y(joker,iv,is1), itw,
552  Y1(joker,iv,is1), gp );
553  }
554  }
555  }
556 
557  // Direction of outgoing scattered radiation (which is
558  // reversed to LOS). Note that this outlos is only used
559  // for extracting scattering properties.
560  Vector outlos;
561  mirror_los( outlos, ppath.los(ip,joker), atmosphere_dim );
562 
563  // Determine phase matrix
564  Tensor4 P( nin, nfs, stokes_dim, stokes_dim );
565  Matrix P1( stokes_dim, stokes_dim );
566  //
567  for( Index ii=0; ii<nin; ii++ )
568  {
569  for( Index iv=0; iv<nfs; iv++ )
570  {
571  pha_mat_singleCalc( P1, outlos[0], outlos[1],
572  fos_scatint_angles(ii,0),
573  fos_scatint_angles(ii,1),
574  scat_data[iv], stokes_dim,
575  ppath_pnd(joker,ip),
576  ppath_t[ip], verbosity );
577  P(ii,iv,joker,joker) = P1;
578  }
579  }
580 
581 
582  // Scattering source term
583  ssource0 = 0.0;
584  for( Index iv=0; iv<nf; iv++ )
585  {
586  Vector sp(stokes_dim);
587  for( Index ii=0; ii<nin; ii++ )
588  {
589  mult( sp, P(ii,iv*ivf,joker,joker),
590  Y(ii,iv,joker) );
591  ssource0(iv,joker) += sp;
592  }
593  }
594  ssource0 *= 4*PI/(Numeric)nin;
595  }
596 
597  // RT of ppath step
598  for( Index iv=0; iv<nf; iv++ )
599  {
600  // Calculate average of absorption, extinction etc.
601  Matrix ext_mat( stokes_dim, stokes_dim );
602  Vector abs_vec( stokes_dim );
603  Vector sbar( stokes_dim, 0 );
604  //
605  // Contribution from abs_species
606  for( Index is1=0; is1<stokes_dim; is1++ )
607  {
608  abs_vec[is1] = 0.5 * (
609  ppath_abs(iv,is1,0,ip) +
610  ppath_abs(iv,is1,0,ip+1) );
611  for( Index is2=0; is2<stokes_dim; is2++ )
612  {
613  ext_mat(is1,is2) = 0.5 * (
614  ppath_abs(iv,is1,is2,ip) +
615  ppath_abs(iv,is1,is2,ip+1) );
616  }
617  }
618  // Particle contribution
619  if( clear2cloudbox[ip] >= 0 )
620  {
621  for( Index is1=0; is1<stokes_dim; is1++ )
622  {
623  sbar[is1] += 0.5 * ssource0(iv,is1);
624  abs_vec[is1] += 0.5 * (
625  pnd_abs_vec(iv,is1,clear2cloudbox[ip]) );
626  for( Index is2=0; is2<stokes_dim; is2++ )
627  {
628  ext_mat(is1,is2) += 0.5 * (
629  pnd_ext_mat(iv,is1,is2,clear2cloudbox[ip]) );
630  }
631  }
632  }
633  if( clear2cloudbox[ip+1] >= 0 )
634  {
635  for( Index is1=0; is1<stokes_dim; is1++ )
636  {
637  sbar[is1] += 0.5 * ssource1(iv,is1);
638  abs_vec[is1] += 0.5 * (
639  pnd_abs_vec(iv,is1,clear2cloudbox[ip+1]) );
640  for( Index is2=0; is2<stokes_dim; is2++ )
641  {
642  ext_mat(is1,is2) += 0.5 * (
643  pnd_ext_mat(iv,is1,is2,clear2cloudbox[ip+1]));
644  }
645  }
646  }
647 
648  // Perform RT
649  //
650  Matrix trans_mat = trans_partial(iv,joker,joker,ip);
651  rte_step_doit( iy(iv,joker), trans_mat, ext_mat, abs_vec,
652  sbar, ppath.lstep[ip], bbar[iv], true );
653  }
654  }
655  }
656 
657  //=== iy_aux part ===================================================
658  // Pressure
659  if( auxPressure >= 0 )
660  { iy_aux[auxPressure](0,0,0,ip) = ppath_p[ip]; }
661  // Temperature
662  if( auxTemperature >= 0 )
663  { iy_aux[auxTemperature](0,0,0,ip) = ppath_t[ip]; }
664  // VMR
665  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
666  { iy_aux[auxVmrSpecies[j]](0,0,0,ip) = ppath_vmr(auxVmrIsp[j],ip);}
667  // Absorption
668  if( auxAbsSum >= 0 )
669  { for( Index iv=0; iv<nf; iv++ ) {
670  for( Index is1=0; is1<ns; is1++ ){
671  for( Index is2=0; is2<ns; is2++ ){
672  iy_aux[auxAbsSum](iv,is1,is2,ip) =
673  ppath_abs(iv,is1,is2,ip); } } } }
674  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
675  { for( Index iv=0; iv<nf; iv++ ) {
676  for( Index is1=0; is1<stokes_dim; is1++ ){
677  for( Index is2=0; is2<stokes_dim; is2++ ){
678  iy_aux[auxAbsSpecies[j]](iv,is1,is2,ip) =
679  abs_per_species(auxAbsIsp[j],iv,is1,is2,ip); } } } }
680  // Particle properties
681  if( cloudbox_on )
682  {
683  // Particle mass content
684  for( Index j=0; j<auxPartCont.nelem(); j++ )
685  { iy_aux[auxPartCont[j]](0,0,0,ip) = ppath_pnd(joker,ip) *
686  particle_masses(joker,auxPartContI[j]); }
687  // Particle field
688  for( Index j=0; j<auxPartField.nelem(); j++ )
689  { iy_aux[auxPartField[j]](0,0,0,ip) =
690  ppath_pnd(auxPartFieldI[j],ip); }
691  }
692  // Radiance
693  if( auxIy >= 0 )
694  { iy_aux[auxIy](joker,joker,0,ip) = iy; }
695  //===================================================================
696  }
697  } // if np>1
698 
699 
700  // Unit conversions
701  if( iy_agenda_call1 )
702  {
703  // If any conversion, check that standard form of Planck used
704  if( !chk_if_std_blackbody_agenda( ws, blackbody_radiation_agenda ) )
705  {
706  ostringstream os;
707  os << "When any unit conversion is applied, "
708  << "*blackbody_radiation_agenda\nmust use "
709  << "*blackbody_radiationPlanck* or a corresponding WSM.\nA test "
710  << "call of the agenda indicates that this is not the case.";
711  throw runtime_error( os.str() );
712  }
713 
714  // Determine refractive index to use for the n2 radiance law
715  Numeric n = 1.0; // First guess is that sensor is in space
716  //
717  if( ppath.end_lstep == 0 ) // If true, sensor inside the atmosphere
718  { n = ppath.nreal[np-1]; }
719 
720  // Polarisation index variable
721  ArrayOfIndex i_pol(stokes_dim);
722  for( Index is=0; is<stokes_dim; is++ )
723  { i_pol[is] = is + 1; }
724 
725  // iy
726  apply_iy_unit( iy, iy_unit, f_grid, n, i_pol );
727 
728  // iy_aux
729  for( Index q=0; q<iy_aux.nelem(); q++ )
730  {
731  if( iy_aux_vars[q] == "iy")
732  {
733  for( Index ip=0; ip<ppath.np; ip++ )
734  { apply_iy_unit( iy_aux[q](joker,joker,0,ip), iy_unit, f_grid,
735  ppath.nreal[ip], i_pol ); }
736  }
737  }
738  }
739 }
740 
741 
742 
743 
744 /*===========================================================================
745  === The functions (in alphabetical order)
746  ===========================================================================*/
747 
748 
749 /* Workspace method: Doxygen documentation will be auto-generated */
750 void iyFOS(
751  Workspace& ws,
752  Matrix& iy,
753  ArrayOfTensor4& iy_aux,
754  Ppath& ppath,
755  ArrayOfTensor3& diy_dx,
756  const Index& stokes_dim,
757  const Vector& f_grid,
758  const Index& atmosphere_dim,
759  const Vector& p_grid,
760  const Tensor3& z_field,
761  const Tensor3& t_field,
762  const Tensor4& vmr_field,
763  const ArrayOfArrayOfSpeciesTag& abs_species,
764  const Tensor3& wind_u_field,
765  const Tensor3& wind_v_field,
766  const Tensor3& wind_w_field,
767  const Tensor3& mag_u_field,
768  const Tensor3& mag_v_field,
769  const Tensor3& mag_w_field,
770  const Index& cloudbox_on,
771  const ArrayOfIndex& cloudbox_limits,
772  const Tensor4& pnd_field,
773  const Index& use_mean_scat_data,
774  const ArrayOfSingleScatteringData& scat_data_array,
775  const Matrix& particle_masses,
776  const String& iy_unit,
777  const ArrayOfString& iy_aux_vars,
778  const Index& jacobian_do,
779  const Agenda& ppath_agenda,
780  const Agenda& blackbody_radiation_agenda,
781  const Agenda& propmat_clearsky_agenda,
782  const Agenda& iy_main_agenda,
783  const Agenda& iy_space_agenda,
784  const Agenda& iy_surface_agenda,
785  const Index& iy_agenda_call1,
786  const Tensor3& iy_transmission,
787  const Vector& rte_pos,
788  const Vector& rte_los,
789  const Vector& rte_pos2,
790  const Numeric& rte_alonglos_v,
791  const Numeric& ppath_lraytrace,
792  const Matrix& fos_scatint_angles,
793  const Vector& fos_iyin_za_angles,
794  const Index& fos_za_interporder,
795  const Index& fos_n,
796  const Verbosity& verbosity )
797 {
798  // Input checks
799  if( jacobian_do )
800  throw runtime_error(
801  "This method does not yet provide any jacobians (jacobian_do must be 0)" );
802  if( fos_scatint_angles.ncols() != 2 )
803  throw runtime_error( "The WSV *fos_scatint_angles* must have two columns." );
804  if( min(fos_scatint_angles(joker,0))<0 ||
805  max(fos_scatint_angles(joker,0))>180 )
806  throw runtime_error(
807  "The zenith angles in *fos_scatint_angles* shall be inside [0,180]." );
808  if( min(fos_scatint_angles(joker,1))<-180 ||
809  max(fos_scatint_angles(joker,1))>180 )
810  throw runtime_error(
811  "The azimuth angles in *fos_scatint_angles* shall be inside [-180,180]." );
812  if( min(fos_iyin_za_angles)<0 || max(fos_iyin_za_angles)>180 )
813  throw runtime_error(
814  "The zenith angles in *fos_iyin_za_angles* shall be inside [0,180]." );
815  if( fos_iyin_za_angles[0] != 0 )
816  throw runtime_error( "The first value in *fos_iyin_za_angles* must be 0." );
817  if( last(fos_iyin_za_angles) != 180 )
818  throw runtime_error( "The last value in *fos_iyin_za_angles* must be 180." );
819  if( fos_za_interporder < 1 )
820  throw runtime_error( "The argument *fos_za_interporder* must be >= 1." );
821  if( fos_iyin_za_angles.nelem() <= fos_za_interporder )
822  throw runtime_error( "The length of *fos_iyin_za_angles* must at least "
823  "be *fos_za_interporder* + 1." );
824  if( fos_n < 0 )
825  throw runtime_error( "The argument *fos_n* must be >= 0." );
826 
827  // Switch to order 0 if not primary call
828  // (This happens after surface reflection. If fos_n used (and >=1), new
829  // surface relections are created ..., and recursion never ends.)
830  Index n = fos_n;
831  if( !iy_agenda_call1 )
832  { n = 0; }
833 
834  fos( ws, iy, iy_aux, ppath, diy_dx, stokes_dim, f_grid, atmosphere_dim,
835  p_grid, z_field, t_field, vmr_field, abs_species, wind_u_field,
836  wind_v_field, wind_w_field, mag_u_field, mag_v_field, mag_w_field,
837  cloudbox_on, cloudbox_limits, pnd_field, use_mean_scat_data,
838  scat_data_array, particle_masses, iy_unit, iy_aux_vars, jacobian_do,
839  ppath_agenda, blackbody_radiation_agenda, propmat_clearsky_agenda,
840  iy_main_agenda, iy_space_agenda, iy_surface_agenda, iy_agenda_call1,
841  iy_transmission, rte_pos, rte_los, rte_pos2, rte_alonglos_v,
842  ppath_lraytrace, fos_scatint_angles, fos_iyin_za_angles,
843  fos_za_interporder, n, 0, verbosity );
844 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
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
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:281
bool chk_if_std_blackbody_agenda(Workspace &ws, const Agenda &blackbody_radiation_agenda)
Checks if blackbody_radiation_agenda returns frequency based radiance.
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
void iyFOS(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 ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Index &use_mean_scat_data, const ArrayOfSingleScatteringData &scat_data_array, const Matrix &particle_masses, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, 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 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 Matrix &fos_scatint_angles, const Vector &fos_iyin_za_angles, const Index &fos_za_interporder, const Index &fos_n, const Verbosity &verbosity)
WORKSPACE METHOD: iyFOS.
Definition: m_fos.cc:750
The Agenda class.
Definition: agenda_class.h:44
void iy_transmission_mult(Tensor3 &iy_trans_total, ConstTensor3View iy_trans_old, ConstTensor3View iy_trans_new)
iy_transmission_mult
Definition: rte.cc:2309
Index nelem() const
Number of elements.
Definition: array.h:176
#define q
Definition: continua.cc:21469
Matrix los
Definition: ppath.h:68
The Vector class.
Definition: matpackI.h:556
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
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
The Tensor4 class.
Definition: matpackIV.h:383
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:183
The range class.
Definition: matpackI.h:148
Vector lstep
Definition: ppath.h:70
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
Matrix pos
Definition: ppath.h:67
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
const Numeric DEG2RAD
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
Numeric end_lstep
Definition: ppath.h:73
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
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
void fos(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 ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Index &use_mean_scat_data, const ArrayOfSingleScatteringData &scat_data_array, const Matrix &particle_masses, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, 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 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 Matrix &fos_scatint_angles, const Vector &fos_iyin_za_angles, const Index &fos_za_interporder, const Index &fos_n, const Index &fos_i, const Verbosity &verbosity)
Definition: m_fos.cc:58
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
#define max(a, b)
Definition: continua.cc:20461
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Vector nreal
Definition: ppath.h:74
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
const Numeric RAD2DEG
const Numeric PI
Radiative transfer in cloudbox.
void emission_rtstep(Matrix &iy, const Index &stokes_dim, ConstVectorView bbar, ArrayOfIndex &extmat_case, ConstTensor3View t)
emission_rtstep
Definition: rte.cc:822
This can be used to make arrays out of anything.
Definition: array.h:40
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
#define ns
Definition: continua.cc:21931
#define min(a, b)
Definition: continua.cc:20460
Index np
Definition: ppath.h:61
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
void get_ppath_ext(ArrayOfIndex &clear2cloudbox, Tensor3 &pnd_abs_vec, Tensor4 &pnd_ext_mat, Array< ArrayOfSingleScatteringData > &scat_data, Matrix &ppath_pnd, const Ppath &ppath, ConstVectorView ppath_t, const Index &stokes_dim, ConstMatrixView ppath_f, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Index &use_mean_scat_data, const ArrayOfSingleScatteringData &scat_data_array, const Verbosity &verbosity)
get_ppath_ext
Definition: rte.cc:1546
void mirror_los(Vector &los_mirrored, ConstVectorView los, const Index &atmosphere_dim)
mirror_los
Definition: rte.cc:2389
void get_ppath_trans2(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, const ArrayOfIndex &clear2cloudbox, ConstTensor4View pnd_ext_mat)
get_ppath_trans2
Definition: rte.cc:1856
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
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
Header file for helper functions for OpenMP.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
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
void ext_matFromabs_vec(MatrixView ext_mat, ConstVectorView abs_vec, const Index &stokes_dim)
Derive extinction matrix from absorption vector.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1403
void rte_step_doit(VectorView stokes_vec, MatrixView trans_mat, ConstMatrixView ext_mat_av, ConstVectorView abs_vec_av, ConstVectorView sca_vec_av, const Numeric &lstep, const Numeric &rtp_planck_value, const bool &trans_is_precalc)
rte_step_doit
Definition: doit.cc:108
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Definition: rte.cc:903
void pha_mat_singleCalc(MatrixView Z, const Numeric za_sca, const Numeric aa_sca, const Numeric za_inc, const Numeric aa_inc, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index stokes_dim, ConstVectorView pnd_vec, const Numeric rtp_temperature, const Verbosity &verbosity)
pha_mat_singleCalc
Definition: montecarlo.cc:928
The Tensor5 class.
Definition: matpackV.h:451
Declaration of functions in rte.cc.