ARTS  2.2.66
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 
36 /*===========================================================================
37  === External declarations
38  ===========================================================================*/
39 
40 #include <cmath>
41 #include <stdexcept>
42 #include "auto_md.h"
43 #include "check_input.h"
44 #include "geodetic.h"
45 #include "logic.h"
46 #include "math_funcs.h"
47 #include "montecarlo.h"
48 #include "physics_funcs.h"
49 #include "ppath.h"
50 #include "rte.h"
51 #include "special_interp.h"
52 #include "lin_alg.h"
53 
54 extern const Numeric SPEED_OF_LIGHT;
55 
56 
57 
58 /*===========================================================================
59  === The functions in alphabetical order
60  ===========================================================================*/
61 
62 
64 
81 void adjust_los(
82  VectorView los,
83  const Index & atmosphere_dim )
84 {
85  if( atmosphere_dim == 1 )
86  {
87  if( los[0] < 0 ) { los[0] = -los[0]; }
88  else if( los[0] > 180 ) { los[0] = 360-los[0]; }
89  }
90  else if( atmosphere_dim == 2 )
91  {
92  if( los[0] < -180 ) { los[0] = los[0] + 360; }
93  else if( los[0] > 180 ) { los[0] = los[0] - 360; }
94  }
95  else
96  {
97  // If any of the angles out-of-bounds, use cart2zaaa to resolve
98  if( abs(los[0]-90) > 90 || abs(los[1]) > 180 )
99  {
100  Numeric dx, dy, dz;
101  zaaa2cart( dx, dy, dz, los[0], los[1] );
102  cart2zaaa( los[0], los[1], dx, dy, dz );
103  }
104  }
105 }
106 
107 
108 
110 
128  MatrixView iy,
129  const String& y_unit,
130  ConstVectorView f_grid,
131  const Numeric& n,
132  const ArrayOfIndex& i_pol )
133 {
134  // The code is largely identical between the two apply_iy_unit functions.
135  // If any change here, remember to update the other function.
136 
137  const Index nf = iy.nrows();
138  const Index ns = iy.ncols();
139 
140  assert( f_grid.nelem() == nf );
141  assert( i_pol.nelem() == ns );
142 
143  if( y_unit == "1" )
144  {
145  if( n != 1 )
146  { iy *= (n*n); }
147  }
148 
149  else if( y_unit == "RJBT" )
150  {
151  for( Index iv=0; iv<nf; iv++ )
152  {
153  const Numeric scfac = invrayjean( 1, f_grid[iv] );
154  for( Index is=0; is<ns; is++ )
155  {
156  if( i_pol[is] < 5 ) // Stokes components
157  { iy(iv,is) *= scfac; }
158  else // Measuement single pols
159  { iy(iv,is) *= 2*scfac; }
160  }
161  }
162  }
163 
164  else if( y_unit == "PlanckBT" )
165  {
166  for( Index iv=0; iv<nf; iv++ )
167  {
168  for( Index is=ns-1; is>=0; is-- ) // Order must here be reversed
169  {
170  if( i_pol[is] == 1 )
171  { iy(iv,is) = invplanck( iy(iv,is), f_grid[iv] ); }
172  else if( i_pol[is] < 5 )
173  {
174  assert( i_pol[0] == 1 );
175  iy(iv,is) =
176  invplanck( 0.5*(iy(iv,0)+iy(iv,is)), f_grid[iv] ) -
177  invplanck( 0.5*(iy(iv,0)-iy(iv,is)), f_grid[iv] );
178  }
179  else
180  { iy(iv,is) = invplanck( 2*iy(iv,is), f_grid[iv] ); }
181  }
182  }
183  }
184 
185  else if ( y_unit == "W/(m^2 m sr)" )
186  {
187  for( Index iv=0; iv<nf; iv++ )
188  {
189  const Numeric scfac = n*n * f_grid[iv] * (f_grid[iv]/SPEED_OF_LIGHT);
190  for( Index is=0; is<ns; is++ )
191  { iy(iv,is) *= scfac; }
192  }
193  }
194 
195  else if ( y_unit == "W/(m^2 m-1 sr)" )
196  {
197  iy *= ( n * n * SPEED_OF_LIGHT );
198  }
199 
200  else
201  {
202  ostringstream os;
203  os << "Unknown option: y_unit = \"" << y_unit << "\"\n"
204  << "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
205  << "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\"";
206 
207  throw runtime_error( os.str() );
208  }
209 }
210 
211 
212 
214 
233  Tensor3View J,
234  ConstMatrixView iy,
235  const String& y_unit,
236  ConstVectorView f_grid,
237  const Numeric& n,
238  const ArrayOfIndex& i_pol )
239 {
240  // The code is largely identical between the two apply_iy_unit functions.
241  // If any change here, remember to update the other function.
242 
243  const Index nf = iy.nrows();
244  const Index ns = iy.ncols();
245  const Index np = J.npages();
246 
247  assert( J.nrows() == nf );
248  assert( J.ncols() == ns );
249  assert( f_grid.nelem() == nf );
250  assert( i_pol.nelem() == ns );
251 
252  if( y_unit == "1" )
253  {
254  if( n != 1 )
255  { J *= (n*n); }
256  }
257 
258  else if( y_unit == "RJBT" )
259  {
260  for( Index iv=0; iv<nf; iv++ )
261  {
262  const Numeric scfac = invrayjean( 1, f_grid[iv] );
263  for( Index is=0; is<ns; is++ )
264  {
265  if( i_pol[is] < 5 ) // Stokes componenets
266  {
267  for( Index ip=0; ip<np; ip++ )
268  { J(ip,iv,is) *= scfac; }
269  }
270  else // Measuement single pols
271  {
272  for( Index ip=0; ip<np; ip++ )
273  { J(ip,iv,is) *= 2*scfac; }
274  }
275  }
276  }
277  }
278 
279  else if( y_unit == "PlanckBT" )
280  {
281  for( Index iv=0; iv<f_grid.nelem(); iv++ )
282  {
283  for( Index is=ns-1; is>=0; is-- )
284  {
285  Numeric scfac = 1;
286  if( i_pol[is] == 1 )
287  { scfac = dinvplanckdI( iy(iv,is), f_grid[iv] ); }
288  else if( i_pol[is] < 5 )
289  {
290  assert( i_pol[0] == 1 );
291  scfac =
292  dinvplanckdI( 0.5*(iy(iv,0)+iy(iv,is)), f_grid[iv] ) +
293  dinvplanckdI( 0.5*(iy(iv,0)-iy(iv,is)), f_grid[iv] );
294  }
295  else
296  { scfac = dinvplanckdI( 2*iy(iv,is), f_grid[iv] ); }
297  //
298  for( Index ip=0; ip<np; ip++ )
299  { J(ip,iv,is) *= scfac; }
300  }
301  }
302  }
303 
304  else if ( y_unit == "W/(m^2 m sr)" )
305  {
306  for( Index iv=0; iv<nf; iv++ )
307  {
308  const Numeric scfac = n*n * f_grid[iv] * (f_grid[iv]/SPEED_OF_LIGHT);
309  for( Index ip=0; ip<np; ip++ )
310  {
311  for( Index is=0; is<ns; is++ )
312  { J(ip,iv,is) *= scfac; }
313  }
314  }
315  }
316 
317  else if ( y_unit == "W/(m^2 m-1 sr)" )
318  {
319  J *= ( n *n * SPEED_OF_LIGHT );
320  }
321 
322  else
323  {
324  ostringstream os;
325  os << "Unknown option: y_unit = \"" << y_unit << "\"\n"
326  << "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
327  << "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\"";
328 
329  throw runtime_error( os.str() );
330  }
331 }
332 
333 
334 
336 
354  Numeric& alpha,
355  const Ppath& ppath )
356 {
357  Numeric theta;
358  if( ppath.dim < 3 )
359  { theta = abs( ppath.start_pos[1] - ppath.end_pos[1] ); }
360  else
361  { theta = sphdist( ppath.start_pos[1], ppath.start_pos[2],
362  ppath.end_pos[1], ppath.end_pos[2] ); }
363 
364  // Eq 17 in Kursinski et al., TAO, 2000:
365  alpha = ppath.start_los[0] - ppath.end_los[0] + theta;
366 
367  // This as
368  // phi_r = 180 - ppath.end_los[0]
369  // phi_t = ppath.start_los[0]
370 }
371 
372 
373 
375 
405  Workspace& ws,
406  Vector& pos,
407  Vector& rte_los,
408  Index& background,
409  ConstVectorView rte_pos,
410  const Numeric& lo0,
411  const Agenda& ppath_step_agenda,
412  const Numeric& ppath_lraytrace,
413  const Index& atmosphere_dim,
414  ConstVectorView p_grid,
415  ConstVectorView lat_grid,
416  ConstVectorView lon_grid,
417  ConstTensor3View t_field,
418  ConstTensor3View z_field,
419  ConstTensor4View vmr_field,
420  ConstVectorView f_grid,
421  ConstVectorView refellipsoid,
422  ConstMatrixView z_surface,
423  const Verbosity& verbosity )
424 {
425  // Special treatment of 1D around zenith/nadir
426  // (zenith angles outside [0,180] are changed by *adjust_los*)
427  bool invert_lat = false;
428  if( atmosphere_dim == 1 && ( rte_los[0] < 0 || rte_los[0] > 180 ) )
429  { invert_lat = true; }
430 
431  // Handle cases where angles have moved out-of-bounds due to disturbance
432  adjust_los( rte_los, atmosphere_dim );
433 
434  // Calculate the ppath for disturbed rte_los
435  Ppath ppx;
436  //
437  ppath_calc( ws, ppx, ppath_step_agenda, atmosphere_dim, p_grid, lat_grid,
438  lon_grid, t_field, z_field, vmr_field,
439  f_grid, refellipsoid, z_surface, 0, ArrayOfIndex(0),
440  rte_pos, rte_los, ppath_lraytrace, 0, verbosity );
441  //
442  background = ppath_what_background( ppx );
443 
444  // Calcualte cumulative optical path for ppx
445  Vector lox( ppx.np );
446  Index ilast = ppx.np-1;
447  lox[0] = ppx.end_lstep;
448  for( Index i=1; i<=ilast; i++ )
449  { lox[i] = lox[i-1] + ppx.lstep[i-1] * ( ppx.nreal[i-1] +
450  ppx.nreal[i] ) / 2.0; }
451 
452  pos.resize( max( Index(2), atmosphere_dim ) );
453 
454  // Reciever at a longer distance (most likely out in space):
455  if( lox[ilast] < lo0 )
456  {
457  const Numeric dl = lo0 - lox[ilast];
458  if( atmosphere_dim < 3 )
459  {
460  Numeric x, z, dx, dz;
461  poslos2cart( x, z, dx, dz, ppx.r[ilast], ppx.pos(ilast,1),
462  ppx.los(ilast,0) );
463  cart2pol( pos[0], pos[1], x+dl*dx, z+dl*dz, ppx.pos(ilast,1),
464  ppx.los(ilast,0) );
465  }
466  else
467  {
468  Numeric x, y, z, dx, dy, dz;
469  poslos2cart( x, y, z, dx, dy, dz, ppx.r[ilast], ppx.pos(ilast,1),
470  ppx.pos(ilast,2), ppx.los(ilast,0), ppx.los(ilast,1) );
471  cart2sph( pos[0], pos[1], pos[2], x+dl*dx, y+dl*dy, z+dl*dz,
472  ppx.pos(ilast,1), ppx.pos(ilast,2),
473  ppx.los(ilast,0), ppx.los(ilast,1) );
474  }
475  }
476 
477  // Interpolate to lo0
478  else
479  {
480  GridPos gp;
481  Vector itw(2);
482  gridpos( gp, lox, lo0 );
483  interpweights( itw, gp );
484  //
485  pos[0] = interp( itw, ppx.r, gp );
486  pos[1] = interp( itw, ppx.pos(joker,1), gp );
487  if( atmosphere_dim == 3 )
488  { pos[2] = interp( itw, ppx.pos(joker,2), gp ); }
489  }
490 
491  if( invert_lat )
492  { pos[1] = -pos[1]; }
493 }
494 
495 
497 
531  Workspace& ws,
532  Numeric& dlf,
533  const Agenda& ppath_step_agenda,
534  const Index& atmosphere_dim,
535  ConstVectorView p_grid,
536  ConstVectorView lat_grid,
537  ConstVectorView lon_grid,
538  ConstTensor3View t_field,
539  ConstTensor3View z_field,
540  ConstTensor4View vmr_field,
541  ConstVectorView f_grid,
542  ConstVectorView refellipsoid,
543  ConstMatrixView z_surface,
544  const Ppath& ppath,
545  const Numeric& ppath_lraytrace,
546  const Numeric& dza,
547  const Verbosity& verbosity )
548 {
549  // Optical and physical path between transmitter and reciver
550  Numeric lo = ppath.start_lstep + ppath.end_lstep;
551  Numeric lp = lo;
552  for( Index i=0; i<=ppath.np-2; i++ )
553  { lp += ppath.lstep[i];
554  lo += ppath.lstep[i] * ( ppath.nreal[i] + ppath.nreal[i+1] ) / 2.0;
555  }
556  // Extract rte_pos and rte_los
557  const Vector rte_pos = ppath.start_pos[Range(0,atmosphere_dim)];
558  //
559  Vector rte_los0(max(Index(1),atmosphere_dim-1)), rte_los;
560  mirror_los( rte_los, ppath.start_los, atmosphere_dim );
561  rte_los0 = rte_los[Range(0,max(Index(1),atmosphere_dim-1))];
562 
563  // A new ppath with positive zenith angle off-set
564  //
565  Vector pos1;
566  Index backg1;
567  //
568  rte_los = rte_los0;
569  rte_los[0] += dza;
570  //
571  defocusing_general_sub( ws, pos1, rte_los, backg1, rte_pos, lo,
572  ppath_step_agenda, ppath_lraytrace, atmosphere_dim,
573  p_grid, lat_grid, lon_grid, t_field, z_field,
574  vmr_field, f_grid, refellipsoid,
575  z_surface, verbosity );
576 
577  // Same thing with negative zenit angle off-set
578  Vector pos2;
579  Index backg2;
580  //
581  rte_los = rte_los0; // Use rte_los0 as rte_los can have been "adjusted"
582  rte_los[0] -= dza;
583  //
584  defocusing_general_sub( ws, pos2, rte_los, backg2, rte_pos, lo,
585  ppath_step_agenda, ppath_lraytrace, atmosphere_dim,
586  p_grid, lat_grid, lon_grid, t_field, z_field,
587  vmr_field, f_grid, refellipsoid,
588  z_surface, verbosity );
589 
590  // Calculate distance between pos1 and 2, and derive the loss factor
591  // All appears OK:
592  if( backg1 == backg2 )
593  {
594  Numeric l12;
595  if( atmosphere_dim < 3 )
596  { distance2D( l12, pos1[0], pos1[1], pos2[0], pos2[1] ); }
597  else
598  { distance3D( l12, pos1[0], pos1[1], pos1[2],
599  pos2[0], pos2[1], pos2[2] ); }
600  //
601  dlf = lp*2*DEG2RAD*dza / l12;
602  }
603  // If different backgrounds, then only use the second calculation
604  else
605  {
606  Numeric l12;
607  if( atmosphere_dim == 1 )
608  {
609  const Numeric r = refellipsoid[0];
610  distance2D( l12, r+ppath.end_pos[0], 0, pos2[0], pos2[1] );
611  }
612  else if( atmosphere_dim == 2 )
613  {
614  const Numeric r = refell2r( refellipsoid, ppath.end_pos[1] );
615  distance2D( l12, r+ppath.end_pos[0], ppath.end_pos[1],
616  pos2[0], pos2[1] );
617  }
618  else
619  {
620  const Numeric r = refell2r( refellipsoid, ppath.end_pos[1] );
621  distance3D( l12, r+ppath.end_pos[0], ppath.end_pos[1],
622  ppath.end_pos[2], pos2[0], pos2[1], pos2[2] );
623  }
624  //
625  dlf = lp*DEG2RAD*dza / l12;
626  }
627 }
628 
629 
630 
632 
664  Workspace& ws,
665  Numeric& dlf,
666  const Agenda& ppath_step_agenda,
667  const Index& atmosphere_dim,
668  ConstVectorView p_grid,
669  ConstVectorView lat_grid,
670  ConstVectorView lon_grid,
671  ConstTensor3View t_field,
672  ConstTensor3View z_field,
673  ConstTensor4View vmr_field,
674  ConstVectorView f_grid,
675  ConstVectorView refellipsoid,
676  ConstMatrixView z_surface,
677  const Ppath& ppath,
678  const Numeric& ppath_lraytrace,
679  const Numeric& dza,
680  const Verbosity& verbosity )
681 {
682  if( ppath.end_los[0] < 90 || ppath.start_los[0] > 90 )
683  throw runtime_error( "The function *defocusing_sat2sat* can only be used "
684  "for limb sounding geometry." );
685 
686  // Index of tangent point
687  Index it;
688  find_tanpoint( it, ppath );
689  assert( it >= 0 );
690 
691  // Length between tangent point and transmitter/reciver
692  Numeric lt = ppath.start_lstep, lr = ppath.end_lstep;
693  for( Index i=it; i<=ppath.np-2; i++ )
694  { lt += ppath.lstep[i]; }
695  for( Index i=0; i<it; i++ )
696  { lr += ppath.lstep[i]; }
697 
698  // Bending angle and impact parameter for centre ray
699  Numeric alpha0, a0;
700  bending_angle1d( alpha0, ppath );
701  alpha0 *= DEG2RAD;
702  a0 = ppath.constant;
703 
704  // Azimuth loss term (Eq 18.5 in Kursinski et al.)
705  const Numeric lf = lr*lt / (lr+lt);
706  const Numeric alt = 1 / ( 1 - alpha0*lf / refellipsoid[0] );
707 
708  // Calculate two new ppaths to get dalpha/da
709  Numeric alpha1, a1, alpha2, a2, dada;
710  Ppath ppt;
711  Vector rte_pos = ppath.end_pos[Range(0,atmosphere_dim)];
712  Vector rte_los = ppath.end_los;
713  //
714  rte_los[0] -= dza;
715  adjust_los( rte_los, atmosphere_dim );
716  ppath_calc( ws, ppt, ppath_step_agenda, atmosphere_dim, p_grid, lat_grid,
717  lon_grid, t_field, z_field, vmr_field,
718  f_grid, refellipsoid, z_surface, 0, ArrayOfIndex(0),
719  rte_pos, rte_los, ppath_lraytrace, 0, verbosity );
720  bending_angle1d( alpha2, ppt );
721  alpha2 *= DEG2RAD;
722  a2 = ppt.constant;
723  //
724  rte_los[0] += 2*dza;
725  adjust_los( rte_los, atmosphere_dim );
726  ppath_calc( ws, ppt, ppath_step_agenda, atmosphere_dim, p_grid, lat_grid,
727  lon_grid, t_field, z_field, vmr_field,
728  f_grid, refellipsoid, z_surface, 0, ArrayOfIndex(0),
729  rte_pos, rte_los, ppath_lraytrace, 0, verbosity );
730  // This path can hit the surface. And we need to check if ppt is OK.
731  // (remember this function only deals with sat-to-sat links and OK
732  // background here is be space)
733  // Otherwise use the centre ray as the second one.
734  if( ppath_what_background(ppt) == 1 )
735  {
736  bending_angle1d( alpha1, ppt );
737  alpha1 *= DEG2RAD;
738  a1 = ppt.constant;
739  dada = (alpha2-alpha1) / (a2-a1);
740  }
741  else
742  {
743  dada = (alpha2-alpha0) / (a2-a0);
744  }
745 
746  // Zenith loss term (Eq 18 in Kursinski et al.)
747  const Numeric zlt = 1 / ( 1 - dada*lf );
748 
749  // Total defocusing loss
750  dlf = zlt * alt;
751 }
752 
753 
754 
756 
777  ConstVectorView los,
778  const Numeric& u,
779  const Numeric& v,
780  const Numeric& w,
781  const Index& atmosphere_dim )
782 {
783  // Strength of field
784  const Numeric f = sqrt( u*u + v*v + w*w );
785 
786  // Zenith and azimuth angle for field (in radians)
787  const Numeric za_f = acos( w/f );
788  const Numeric aa_f = atan2( u, v );
789 
790  // Zenith and azimuth angle for photon direction (in radians)
791  Vector los_p;
792  mirror_los( los_p, los, atmosphere_dim );
793  const Numeric za_p = DEG2RAD * los_p[0];
794  const Numeric aa_p = DEG2RAD * los_p[1];
795 
796  return f * ( cos(za_f) * cos(za_p) +
797  sin(za_f) * sin(za_p) * cos(aa_f-aa_p) );
798 }
799 
800 
801 
803 
823  Matrix& iy,
824  const Index& stokes_dim,
825  ConstVectorView bbar,
826  ArrayOfIndex& extmat_case,
827  ConstTensor3View t )
828 {
829  const Index nf = bbar.nelem();
830 
831  assert( t.ncols() == stokes_dim && t.nrows() == stokes_dim );
832  assert( t.npages() == nf );
833  assert( extmat_case.nelem() == nf );
834 
835  // Spectrum at end of ppath step
836  if( stokes_dim == 1 )
837  {
838  for( Index iv=0; iv<nf; iv++ )
839  { iy(iv,0) = iy(iv,0) * t(iv,0,0) + bbar[iv] * ( 1 - t(iv,0,0) ); }
840  }
841 
842  else
843  {
844 #pragma omp parallel for \
845  if (!arts_omp_in_parallel() \
846  && nf >= arts_omp_get_max_threads())
847  for( Index iv=0; iv<nf; iv++ )
848  {
849  assert( extmat_case[iv]>=1 && extmat_case[iv]<=3 );
850  // Unpolarised absorption:
851  if( extmat_case[iv] == 1 )
852  {
853  iy(iv,0) = iy(iv,0) * t(iv,0,0) + bbar[iv] * ( 1 - t(iv,0,0) );
854  for( Index is=1; is<stokes_dim; is++ )
855  { iy(iv,is) = iy(iv,is) * t(iv,is,is); }
856  }
857  // The general case:
858  else
859  {
860  // Transmitted term
861  Vector tt(stokes_dim);
862  mult( tt, t(iv,joker,joker), iy(iv,joker));
863  // Add emission, first Stokes element
864  iy(iv,0) = tt[0] + bbar[iv] * ( 1 - t(iv,0,0) );
865  // Remaining Stokes elements
866  for( Index i=1; i<stokes_dim; i++ )
867  { iy(iv,i) = tt[i] - bbar[iv] * t(iv,i,0); }
868 
869  }
870  }
871  }
872 }
873 
874 
875 
877 
904  MatrixView trans_mat,
905  Index& icase,
906  ConstMatrixView ext_mat,
907  const Numeric& lstep )
908 {
909  const Index stokes_dim = ext_mat.ncols();
910 
911  assert( ext_mat.nrows()==stokes_dim );
912  assert( trans_mat.nrows()==stokes_dim && trans_mat.ncols()==stokes_dim );
913 
914  // Theoretically ext_mat(0,0) >= 0, but to demand this can cause problems for
915  // iterative retrievals, and the assert is skipped. Negative should be a
916  // result of negative vmr, and this issue is checked in basics_checkedCalc.
917  //assert( ext_mat(0,0) >= 0 );
918 
919  assert( icase>=0 && icase<=3 );
920  assert( !is_singular( ext_mat ) );
921  assert( lstep >= 0 );
922 
923  // Analyse ext_mat?
924  if( icase == 0 )
925  {
926  icase = 1; // Start guess is diagonal
927 
928  //--- Scalar case ----------------------------------------------------------
929  if( stokes_dim == 1 )
930  {}
931 
932  //--- Vector RT ------------------------------------------------------------
933  else
934  {
935  // Check symmetries and analyse structure of exp_mat:
936  assert( ext_mat(1,1) == ext_mat(0,0) );
937  assert( ext_mat(1,0) == ext_mat(0,1) );
938 
939  if( ext_mat(1,0) != 0 )
940  { icase = 2; }
941 
942  if( stokes_dim >= 3 )
943  {
944  assert( ext_mat(2,2) == ext_mat(0,0) );
945  assert( ext_mat(2,1) == -ext_mat(1,2) );
946  assert( ext_mat(2,0) == ext_mat(0,2) );
947 
948  if( ext_mat(2,0) != 0 || ext_mat(2,1) != 0 )
949  { icase = 3; }
950 
951  if( stokes_dim > 3 )
952  {
953  assert( ext_mat(3,3) == ext_mat(0,0) );
954  assert( ext_mat(3,2) == -ext_mat(2,3) );
955  assert( ext_mat(3,1) == -ext_mat(1,3) );
956  assert( ext_mat(3,0) == ext_mat(0,3) );
957 
958  if( icase < 3 ) // if icase==3, already at most complex case
959  {
960  if( ext_mat(3,0) != 0 || ext_mat(3,1) != 0 )
961  { icase = 3; }
962  else if( ext_mat(3,2) != 0 )
963  { icase = 2; }
964  }
965  }
966  }
967  }
968  }
969 
970 
971  // Calculation options:
972  if( icase == 1 )
973  {
974  trans_mat = 0;
975  trans_mat(0,0) = exp( -ext_mat(0,0) * lstep );
976  for( Index i=1; i<stokes_dim; i++ )
977  { trans_mat(i,i) = trans_mat(0,0); }
978  }
979 
980  else if( icase == 2 )
981  {
982  // Expressions below are found in "Polarization in Spectral Lines" by
983  // Landi Degl'Innocenti and Landolfi (2004).
984  const Numeric tI = exp( -ext_mat(0,0) * lstep );
985  const Numeric HQ = ext_mat(0,1) * lstep;
986  trans_mat(0,0) = tI * cosh( HQ );
987  trans_mat(1,1) = trans_mat(0,0);
988  trans_mat(1,0) = -tI * sinh( HQ );
989  trans_mat(0,1) = trans_mat(1,0);
990  if( stokes_dim >= 3 )
991  {
992  trans_mat(2,0) = 0;
993  trans_mat(2,1) = 0;
994  trans_mat(0,2) = 0;
995  trans_mat(1,2) = 0;
996  const Numeric RQ = ext_mat(2,3) * lstep;
997  trans_mat(2,2) = tI * cos( RQ );
998  if( stokes_dim > 3 )
999  {
1000  trans_mat(3,0) = 0;
1001  trans_mat(3,1) = 0;
1002  trans_mat(0,3) = 0;
1003  trans_mat(1,3) = 0;
1004  trans_mat(3,3) = trans_mat(2,2);
1005  trans_mat(3,2) = tI * sin( RQ );
1006  trans_mat(2,3) = -trans_mat(3,2);
1007  }
1008  }
1009  }
1010  else
1011  {
1012  Matrix ext_mat_ds = ext_mat;
1013  ext_mat_ds *= -lstep;
1014  //
1015  Index q = 10; // index for the precision of the matrix exp function
1016  //
1017  matrix_exp( trans_mat, ext_mat_ds, q );
1018  }
1019 }
1020 
1021 
1022 
1023 
1024 
1026 
1046 void get_iy(
1047  Workspace& ws,
1048  Matrix& iy,
1049  ConstTensor3View t_field,
1050  ConstTensor3View z_field,
1051  ConstTensor4View vmr_field,
1052  const Index& cloudbox_on,
1053  ConstVectorView f_grid,
1054  ConstVectorView rte_pos,
1055  ConstVectorView rte_los,
1056  ConstVectorView rte_pos2,
1057  const Agenda& iy_main_agenda )
1058 {
1059  ArrayOfTensor3 diy_dx;
1060  ArrayOfTensor4 iy_aux;
1061  Ppath ppath;
1062  Tensor3 iy_transmission(0,0,0);
1063 
1064  iy_main_agendaExecute( ws, iy, iy_aux, ppath, diy_dx, 1, iy_transmission,
1065  ArrayOfString(0), cloudbox_on, 0, t_field, z_field,
1066  vmr_field, f_grid, rte_pos, rte_los, rte_pos2,
1067  iy_main_agenda );
1068 }
1069 
1070 
1071 
1072 
1074 
1107  Workspace& ws,
1108  Matrix& iy,
1109  ArrayOfTensor3& diy_dx,
1110  ConstTensor3View iy_transmission,
1111  const Index& jacobian_do,
1112  const Ppath& ppath,
1113  ConstVectorView rte_pos2,
1114  const Index& atmosphere_dim,
1115  ConstTensor3View t_field,
1116  ConstTensor3View z_field,
1117  ConstTensor4View vmr_field,
1118  const Index& cloudbox_on,
1119  const Index& stokes_dim,
1120  ConstVectorView f_grid,
1121  const Agenda& iy_main_agenda,
1122  const Agenda& iy_space_agenda,
1123  const Agenda& iy_surface_agenda,
1124  const Agenda& iy_cloudbox_agenda,
1125  const Verbosity& verbosity)
1126 {
1127  CREATE_OUT3;
1128 
1129  // Some sizes
1130  const Index nf = f_grid.nelem();
1131  const Index np = ppath.np;
1132 
1133  // Set rtp_pos and rtp_los to match the last point in ppath.
1134  //
1135  // Note that the Ppath positions (ppath.pos) for 1D have one column more
1136  // than expected by most functions. Only the first atmosphere_dim values
1137  // shall be copied.
1138  //
1139  Vector rtp_pos, rtp_los;
1140  rtp_pos.resize( atmosphere_dim );
1141  rtp_pos = ppath.pos(np-1,Range(0,atmosphere_dim));
1142  rtp_los.resize( ppath.los.ncols() );
1143  rtp_los = ppath.los(np-1,joker);
1144 
1145  out3 << "Radiative background: " << ppath.background << "\n";
1146 
1147 
1148  // Handle the different background cases
1149  //
1150  String agenda_name;
1151  //
1152  switch( ppath_what_background( ppath ) )
1153  {
1154 
1155  case 1: //--- Space ----------------------------------------------------
1156  {
1157  agenda_name = "iy_space_agenda";
1158  chk_not_empty( agenda_name, iy_space_agenda );
1159  iy_space_agendaExecute( ws, iy, f_grid, rtp_pos, rtp_los,
1160  iy_space_agenda );
1161  }
1162  break;
1163 
1164  case 2: //--- The surface -----------------------------------------------
1165  {
1166  agenda_name = "iy_surface_agenda";
1167  chk_not_empty( agenda_name, iy_surface_agenda );
1168  iy_surface_agendaExecute( ws, iy, diy_dx, iy_transmission, cloudbox_on,
1169  jacobian_do, t_field, z_field, vmr_field,
1170  f_grid, iy_main_agenda, rtp_pos, rtp_los,
1171  rte_pos2, iy_surface_agenda );
1172  }
1173  break;
1174 
1175  case 3: //--- Cloudbox boundary or interior ------------------------------
1176  case 4:
1177  {
1178  agenda_name = "iy_cloudbox_agenda";
1179  chk_not_empty( agenda_name, iy_cloudbox_agenda );
1180  iy_cloudbox_agendaExecute( ws, iy, f_grid, rtp_pos, rtp_los,
1181  iy_cloudbox_agenda );
1182  }
1183  break;
1184 
1185  default: //--- ????? ----------------------------------------------------
1186  // Are we here, the coding is wrong somewhere
1187  assert( false );
1188  }
1189 
1190  if( iy.ncols() != stokes_dim || iy.nrows() != nf )
1191  {
1192  ostringstream os;
1193  os << "The size of *iy* returned from *" << agenda_name << "* is\n"
1194  << "not correct:\n"
1195  << " expected size = [" << nf << "," << stokes_dim << "]\n"
1196  << " size of iy = [" << iy.nrows() << "," << iy.ncols()<< "]\n";
1197  throw runtime_error( os.str() );
1198  }
1199 }
1200 
1201 
1202 
1204 
1234  Vector& ppath_p,
1235  Vector& ppath_t,
1236  Matrix& ppath_vmr,
1237  Matrix& ppath_wind,
1238  Matrix& ppath_mag,
1239  const Ppath& ppath,
1240  const Index& atmosphere_dim,
1241  ConstVectorView p_grid,
1242  ConstTensor3View t_field,
1243  ConstTensor4View vmr_field,
1244  ConstTensor3View wind_u_field,
1245  ConstTensor3View wind_v_field,
1246  ConstTensor3View wind_w_field,
1247  ConstTensor3View mag_u_field,
1248  ConstTensor3View mag_v_field,
1249  ConstTensor3View mag_w_field )
1250 {
1251  const Index np = ppath.np;
1252  // Pressure:
1253  ppath_p.resize(np);
1254  Matrix itw_p(np,2);
1255  interpweights( itw_p, ppath.gp_p );
1256  itw2p( ppath_p, p_grid, ppath.gp_p, itw_p );
1257 
1258  // Temperature:
1259  ppath_t.resize(np);
1260  Matrix itw_field;
1261  interp_atmfield_gp2itw( itw_field, atmosphere_dim,
1262  ppath.gp_p, ppath.gp_lat, ppath.gp_lon );
1263  interp_atmfield_by_itw( ppath_t, atmosphere_dim, t_field,
1264  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1265 
1266  // VMR fields:
1267  const Index ns = vmr_field.nbooks();
1268  ppath_vmr.resize(ns,np);
1269  for( Index is=0; is<ns; is++ )
1270  {
1271  interp_atmfield_by_itw( ppath_vmr(is, joker), atmosphere_dim,
1272  vmr_field( is, joker, joker, joker ),
1273  ppath.gp_p, ppath.gp_lat, ppath.gp_lon,
1274  itw_field );
1275  }
1276 
1277  // Winds:
1278  ppath_wind.resize(3,np);
1279  ppath_wind = 0;
1280  //
1281  if( wind_u_field.npages() > 0 )
1282  {
1283  interp_atmfield_by_itw( ppath_wind(0,joker), atmosphere_dim, wind_u_field,
1284  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1285  }
1286  if( wind_v_field.npages() > 0 )
1287  {
1288  interp_atmfield_by_itw( ppath_wind(1,joker), atmosphere_dim, wind_v_field,
1289  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1290  }
1291  if( wind_w_field.npages() > 0 )
1292  {
1293  interp_atmfield_by_itw( ppath_wind(2,joker), atmosphere_dim, wind_w_field,
1294  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1295  }
1296 
1297  // Magnetic field:
1298  ppath_mag.resize(3,np);
1299  ppath_mag = 0;
1300  //
1301  if( mag_u_field.npages() > 0 )
1302  {
1303  interp_atmfield_by_itw( ppath_mag(0,joker), atmosphere_dim, mag_u_field,
1304  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1305  }
1306  if( mag_v_field.npages() > 0 )
1307  {
1308  interp_atmfield_by_itw( ppath_mag(1,joker), atmosphere_dim, mag_v_field,
1309  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1310  }
1311  if( mag_w_field.npages() > 0 )
1312  {
1313  interp_atmfield_by_itw( ppath_mag(2,joker), atmosphere_dim, mag_w_field,
1314  ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
1315  }
1316 }
1317 
1318 
1319 
1321 
1352  Workspace& ws,
1353  Tensor4& ppath_abs,
1354  Tensor5& abs_per_species,
1355  const Agenda& propmat_clearsky_agenda,
1356  const Ppath& ppath,
1357  ConstVectorView ppath_p,
1358  ConstVectorView ppath_t,
1359  ConstMatrixView ppath_vmr,
1360  ConstMatrixView ppath_f,
1361  ConstMatrixView ppath_mag,
1362  ConstVectorView f_grid,
1363  const Index& stokes_dim,
1364  const ArrayOfIndex& ispecies )
1365 {
1366  // Sizes
1367  const Index nf = f_grid.nelem();
1368  const Index np = ppath.np;
1369  const Index nabs = ppath_vmr.nrows();
1370  const Index nisp = ispecies.nelem();
1371 
1372  DEBUG_ONLY(
1373  for( Index i=0; i<nisp; i++ )
1374  {
1375  assert( ispecies[i] >= 0 );
1376  assert( ispecies[i] < nabs );
1377  }
1378  )
1379 
1380  // Size variable
1381  try
1382  {
1383  ppath_abs.resize( nf, stokes_dim, stokes_dim, np );
1384  abs_per_species.resize( nisp, nf, stokes_dim, stokes_dim, np );
1385  }
1386  catch (std::bad_alloc x)
1387  {
1388  ostringstream os;
1389  os << "Run-time error in function: get_ppath_abs" << endl
1390  << "Memory allocation failed for ppath_abs("
1391  << nabs << ", " << nf << ", " << stokes_dim << ", "
1392  << stokes_dim << ", " << np << ")" << endl;
1393  throw runtime_error(os.str());
1394  }
1395 
1396  String fail_msg;
1397  bool failed = false;
1398 
1399  // Loop ppath points
1400  //
1401  Workspace l_ws (ws);
1402  Agenda l_propmat_clearsky_agenda (propmat_clearsky_agenda);
1403  //
1404  if (np)
1405 #pragma omp parallel for \
1406  if (!arts_omp_in_parallel() \
1407  && np >= arts_omp_get_max_threads()) \
1408  firstprivate(l_ws, l_propmat_clearsky_agenda)
1409  for( Index ip=0; ip<np; ip++ )
1410  {
1411  if (failed) continue;
1412 
1413  // Call agenda
1414  //
1415  Tensor4 propmat_clearsky;
1416  //
1417  try {
1418  Vector rtp_vmr(0);
1419  if( nabs )
1420  {
1421  propmat_clearsky_agendaExecute( l_ws, propmat_clearsky,
1422  ppath_f(joker,ip), ppath_mag(joker,ip), ppath.los(ip,joker),
1423  ppath_p[ip], ppath_t[ip], ppath_vmr(joker,ip),
1424  l_propmat_clearsky_agenda );
1425  }
1426  else
1427  {
1428  propmat_clearsky_agendaExecute( l_ws, propmat_clearsky,
1429  ppath_f(joker,ip), ppath_mag(joker,ip), ppath.los(ip,joker),
1430  ppath_p[ip], ppath_t[ip], Vector(0), l_propmat_clearsky_agenda );
1431  }
1432  } catch (runtime_error e) {
1433 #pragma omp critical (get_ppath_abs_fail)
1434  { failed = true; fail_msg = e.what();}
1435  }
1436 
1437  // Copy to output argument
1438  if( !failed )
1439  {
1440  assert( propmat_clearsky.ncols() == stokes_dim );
1441  assert( propmat_clearsky.nrows() == stokes_dim );
1442  assert( propmat_clearsky.npages() == nf );
1443  assert( propmat_clearsky.nbooks() == max(nabs,Index(1)) );
1444 
1445  for( Index i1=0; i1<nf; i1++ )
1446  {
1447  for( Index i2=0; i2<stokes_dim; i2++ )
1448  {
1449  for( Index i3=0; i3<stokes_dim; i3++ )
1450  {
1451  ppath_abs(i1,i2,i3,ip) = propmat_clearsky(joker,i1,i2,i3).sum();
1452 
1453  for( Index ia=0; ia<nisp; ia++ )
1454  {
1455  abs_per_species(ia,i1,i2,i3,ip) =
1456  propmat_clearsky(ispecies[ia],i1,i2,i3);
1457  }
1458 
1459  }
1460  }
1461  }
1462  }
1463  }
1464 
1465  if (failed)
1466  throw runtime_error(fail_msg);
1467 }
1468 
1469 
1470 
1472 
1488  Workspace& ws,
1489  Matrix& ppath_blackrad,
1490  const Agenda& blackbody_radiation_agenda,
1491  const Ppath& ppath,
1492  ConstVectorView ppath_t,
1493  ConstMatrixView ppath_f )
1494 {
1495  // Sizes
1496  const Index nf = ppath_f.nrows();
1497  const Index np = ppath.np;
1498 
1499  // Loop path and call agenda
1500  //
1501  ppath_blackrad.resize( nf, np );
1502  //
1503  for( Index ip=0; ip<np; ip++ )
1504  {
1505  Vector bvector;
1506 
1507  blackbody_radiation_agendaExecute( ws, bvector, ppath_t[ip],
1508  ppath_f(joker,ip),
1509  blackbody_radiation_agenda );
1510  ppath_blackrad(joker,ip) = bvector;
1511  }
1512 }
1513 
1514 
1515 
1517 
1547  ArrayOfIndex& clear2cloudbox,
1548  Tensor3& pnd_abs_vec,
1549  Tensor4& pnd_ext_mat,
1551  Matrix& ppath_pnd,
1552  const Ppath& ppath,
1553  ConstVectorView ppath_t,
1554  const Index& stokes_dim,
1555  ConstMatrixView ppath_f,
1556  const Index& atmosphere_dim,
1557  const ArrayOfIndex& cloudbox_limits,
1558  const Tensor4& pnd_field,
1559  const Index& use_mean_scat_data,
1560  const ArrayOfSingleScatteringData& scat_data_array,
1561  const Verbosity& verbosity )
1562 {
1563  const Index nf = ppath_f.nrows();
1564  const Index np = ppath.np;
1565 
1566  // Pnd along the ppath
1567  ppath_pnd.resize( pnd_field.nbooks(), np );
1568  ppath_pnd = 0;
1569 
1570  // A variable that maps from total ppath to extension data index.
1571  // If outside cloudbox or all pnd=0, this variable holds -1.
1572  // Otherwise it gives the index in pnd_ext_mat etc.
1573  clear2cloudbox.resize( np );
1574 
1575  // Determine ppath_pnd
1576  Index nin = 0;
1577  for( Index ip=0; ip<np; ip++ )
1578  {
1579  Matrix itw( 1, Index(pow(2.0,Numeric(atmosphere_dim))) );
1580 
1581  ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
1582  GridPos gp_lat, gp_lon;
1583  if( atmosphere_dim >= 2 ) { gridpos_copy( gp_lat, ppath.gp_lat[ip] ); }
1584  if( atmosphere_dim == 3 ) { gridpos_copy( gp_lon, ppath.gp_lon[ip] ); }
1585  if( is_gp_inside_cloudbox( ppath.gp_p[ip], gp_lat, gp_lon,
1586  cloudbox_limits, true, atmosphere_dim ) )
1587  {
1589  gpc_p[0], gpc_lat[0], gpc_lon[0],
1590  ppath.gp_p[ip], gp_lat, gp_lon,
1591  atmosphere_dim, cloudbox_limits );
1592  for( Index i=0; i<pnd_field.nbooks(); i++ )
1593  {
1594  interp_atmfield_by_itw( ppath_pnd(i,ip), atmosphere_dim,
1595  pnd_field(i,joker,joker,joker),
1596  gpc_p, gpc_lat, gpc_lon, itw );
1597  }
1598  if( max(ppath_pnd(joker,ip)) > 0 )
1599  { clear2cloudbox[ip] = nin; nin++; }
1600  else
1601  { clear2cloudbox[ip] = -1; }
1602  }
1603  else
1604  { clear2cloudbox[ip] = -1; }
1605  }
1606 
1607  // Particle single scattering properties (are independent of position)
1608  //
1609  if( use_mean_scat_data )
1610  {
1611  const Numeric f = (mean(ppath_f(0,joker))+mean(ppath_f(nf-1,joker)))/2.0;
1612  scat_data.resize( 1 );
1613  scat_data_array_monoCalc( scat_data[0], scat_data_array, Vector(1,f), 0,
1614  verbosity );
1615  }
1616  else
1617  {
1618  scat_data.resize( nf );
1619  for( Index iv=0; iv<nf; iv++ )
1620  {
1621  const Numeric f = mean(ppath_f(iv,joker));
1622  scat_data_array_monoCalc( scat_data[iv], scat_data_array, Vector(1,f), 0,
1623  verbosity );
1624  }
1625  }
1626 
1627  // Resize absorption and extension tensors
1628  pnd_abs_vec.resize( nf, stokes_dim, nin );
1629  pnd_ext_mat.resize( nf, stokes_dim, stokes_dim, nin );
1630 
1631  // Loop ppath points
1632  //
1633  for( Index ip=0; ip<np; ip++ )
1634  {
1635  const Index i = clear2cloudbox[ip];
1636  if( i>=0 )
1637  {
1638  // Direction of outgoing scattered radiation (which is reversed to
1639  // LOS). Note that rtp_los2 is only used for extracting scattering
1640  // properties.
1641  Vector rtp_los2;
1642  mirror_los( rtp_los2, ppath.los(ip,joker), atmosphere_dim );
1643 
1644  // Extinction and absorption
1645  if( use_mean_scat_data )
1646  {
1647  Vector abs_vec( stokes_dim );
1648  Matrix ext_mat( stokes_dim, stokes_dim );
1649  opt_propCalc( ext_mat, abs_vec, rtp_los2[0], rtp_los2[1],
1650  scat_data[0], stokes_dim, ppath_pnd(joker,ip),
1651  ppath_t[ip], verbosity);
1652  for( Index iv=0; iv<nf; iv++ )
1653  {
1654  pnd_ext_mat(iv,joker,joker,i) = ext_mat;
1655  pnd_abs_vec(iv,joker,i) = abs_vec;
1656  }
1657  }
1658  else
1659  {
1660  for( Index iv=0; iv<nf; iv++ )
1661  {
1662  opt_propCalc( pnd_ext_mat(iv,joker,joker,i),
1663  pnd_abs_vec(iv,joker,i), rtp_los2[0],
1664  rtp_los2[1], scat_data[iv], stokes_dim,
1665  ppath_pnd(joker,ip), ppath_t[ip], verbosity );
1666  }
1667  }
1668  }
1669  }
1670 }
1671 
1672 
1673 
1675 
1691  Matrix& ppath_f,
1692  const Ppath& ppath,
1693  ConstVectorView f_grid,
1694  const Index& atmosphere_dim,
1695  const Numeric& rte_alonglos_v,
1696  ConstMatrixView ppath_wind )
1697 {
1698  // Sizes
1699  const Index nf = f_grid.nelem();
1700  const Index np = ppath.np;
1701 
1702  ppath_f.resize(nf,np);
1703 
1704  // Doppler relevant velocity
1705  //
1706  for( Index ip=0; ip<np; ip++ )
1707  {
1708  // Start by adding rte_alonglos_v (most likely sensor effects)
1709  Numeric v_doppler = rte_alonglos_v;
1710 
1711  // Include wind
1712  if( ppath_wind(1,ip) != 0 || ppath_wind(0,ip) != 0 ||
1713  ppath_wind(2,ip) != 0 )
1714  {
1715  // The dot product below is valid for the photon direction. Winds
1716  // along this direction gives a positive contribution.
1717  v_doppler += dotprod_with_los( ppath.los(ip,joker), ppath_wind(0,ip),
1718  ppath_wind(1,ip), ppath_wind(2,ip), atmosphere_dim );
1719  }
1720 
1721  // Determine frequency grid
1722  if( v_doppler == 0 )
1723  { ppath_f(joker,ip) = f_grid; }
1724  else
1725  {
1726  // Positive v_doppler means that sensor measures lower rest
1727  // frequencies
1728  const Numeric a = 1 - v_doppler / SPEED_OF_LIGHT;
1729  for( Index iv=0; iv<nf; iv++ )
1730  { ppath_f(iv,ip) = a * f_grid[iv]; }
1731  }
1732  }
1733 }
1734 
1735 
1736 
1738 
1771  Tensor4& trans_partial,
1772  ArrayOfArrayOfIndex& extmat_case,
1773  Tensor4& trans_cumulat,
1774  Vector& scalar_tau,
1775  const Ppath& ppath,
1776  ConstTensor4View& ppath_abs,
1777  ConstVectorView f_grid,
1778  const Index& stokes_dim )
1779 {
1780  // Sizes
1781  const Index nf = f_grid.nelem();
1782  const Index np = ppath.np;
1783 
1784  // Init variables
1785  //
1786  trans_partial.resize( nf, stokes_dim, stokes_dim, np-1 );
1787  trans_cumulat.resize( nf, stokes_dim, stokes_dim, np );
1788  //
1789  extmat_case.resize(np-1);
1790  for( Index i=0; i<np-1; i++ )
1791  { extmat_case[i].resize(nf); }
1792  //
1793  scalar_tau.resize( nf );
1794  scalar_tau = 0;
1795 
1796  // Loop ppath points (in the anti-direction of photons)
1797  //
1798  for( Index ip=0; ip<np; ip++ )
1799  {
1800  // If first point, calculate sum of absorption and set transmission
1801  // to identity matrix.
1802  if( ip == 0 )
1803  {
1804  for( Index iv=0; iv<nf; iv++ )
1805  { id_mat( trans_cumulat(iv,joker,joker,ip) ); }
1806  }
1807 
1808  else
1809  {
1810  for( Index iv=0; iv<nf; iv++ )
1811  {
1812  // Transmission due to absorption
1813  Matrix ext_mat(stokes_dim,stokes_dim);
1814  for( Index is1=0; is1<stokes_dim; is1++ ) {
1815  for( Index is2=0; is2<stokes_dim; is2++ ) {
1816  ext_mat(is1,is2) = 0.5 * ( ppath_abs(iv,is1,is2,ip-1) +
1817  ppath_abs(iv,is1,is2,ip ) );
1818  } }
1819  scalar_tau[iv] += ppath.lstep[ip-1] * ext_mat(0,0);
1820  extmat_case[ip-1][iv] = 0;
1821  ext2trans( trans_partial(iv,joker,joker,ip-1),
1822  extmat_case[ip-1][iv], ext_mat, ppath.lstep[ip-1] );
1823 
1824  // Cumulative transmission
1825  // (note that multiplication below depends on ppath loop order)
1826  mult( trans_cumulat(iv,joker,joker,ip),
1827  trans_cumulat(iv,joker,joker,ip-1),
1828  trans_partial(iv,joker,joker,ip-1) );
1829  }
1830  }
1831  }
1832 }
1833 
1834 
1835 
1837 
1857  Tensor4& trans_partial,
1858  ArrayOfArrayOfIndex& extmat_case,
1859  Tensor4& trans_cumulat,
1860  Vector& scalar_tau,
1861  const Ppath& ppath,
1862  ConstTensor4View& ppath_abs,
1863  ConstVectorView f_grid,
1864  const Index& stokes_dim,
1865  const ArrayOfIndex& clear2cloudbox,
1866  ConstTensor4View pnd_ext_mat )
1867 {
1868  // Sizes
1869  const Index nf = f_grid.nelem();
1870  const Index np = ppath.np;
1871 
1872  // Init variables
1873  //
1874  trans_partial.resize( nf, stokes_dim, stokes_dim, np-1 );
1875  trans_cumulat.resize( nf, stokes_dim, stokes_dim, np );
1876  //
1877  extmat_case.resize(np-1);
1878  for( Index i=0; i<np-1; i++ )
1879  { extmat_case[i].resize(nf); }
1880  //
1881  scalar_tau.resize( nf );
1882  scalar_tau = 0;
1883 
1884  // Loop ppath points (in the anti-direction of photons)
1885  //
1886  Tensor3 extsum_old, extsum_this( nf, stokes_dim, stokes_dim );
1887  //
1888  for( Index ip=0; ip<np; ip++ )
1889  {
1890  // If first point, calculate sum of absorption and set transmission
1891  // to identity matrix.
1892  if( ip == 0 )
1893  {
1894  for( Index iv=0; iv<nf; iv++ )
1895  {
1896  for( Index is1=0; is1<stokes_dim; is1++ ) {
1897  for( Index is2=0; is2<stokes_dim; is2++ ) {
1898  extsum_this(iv,is1,is2) = ppath_abs(iv,is1,is2,ip);
1899  } }
1900  id_mat( trans_cumulat(iv,joker,joker,ip) );
1901  }
1902  // First point should not be "cloudy", but just in case:
1903  if( clear2cloudbox[ip] >= 0 )
1904  {
1905  const Index ic = clear2cloudbox[ip];
1906  for( Index iv=0; iv<nf; iv++ ) {
1907  for( Index is1=0; is1<stokes_dim; is1++ ) {
1908  for( Index is2=0; is2<stokes_dim; is2++ ) {
1909  extsum_this(iv,is1,is2) += pnd_ext_mat(iv,is1,is2,ic);
1910  } } }
1911  }
1912  }
1913 
1914  else
1915  {
1916  const Index ic = clear2cloudbox[ip];
1917  //
1918  for( Index iv=0; iv<nf; iv++ )
1919  {
1920  // Transmission due to absorption and scattering
1921  Matrix ext_mat(stokes_dim,stokes_dim); // -1*tau
1922  for( Index is1=0; is1<stokes_dim; is1++ ) {
1923  for( Index is2=0; is2<stokes_dim; is2++ ) {
1924  extsum_this(iv,is1,is2) = ppath_abs(iv,is1,is2,ip);
1925  if( ic >= 0 )
1926  { extsum_this(iv,is1,is2) += pnd_ext_mat(iv,is1,is2,ic); }
1927 
1928  ext_mat(is1,is2) = 0.5 * ( extsum_old(iv,is1,is2) +
1929  extsum_this(iv,is1,is2) );
1930  } }
1931  scalar_tau[iv] += ppath.lstep[ip-1] * ext_mat(0,0);
1932  extmat_case[ip-1][iv] = 0;
1933  ext2trans( trans_partial(iv,joker,joker,ip-1),
1934  extmat_case[ip-1][iv], ext_mat, ppath.lstep[ip-1] );
1935 
1936  // Note that multiplication below depends on ppath loop order
1937  mult( trans_cumulat(iv,joker,joker,ip),
1938  trans_cumulat(iv,joker,joker,ip-1),
1939  trans_partial(iv,joker,joker,ip-1) );
1940  }
1941  }
1942 
1943  extsum_old = extsum_this;
1944  }
1945 }
1946 
1947 
1948 
1950 
1961  const Sparse& sensor_response,
1962  const Index& mblock_index )
1963 {
1964  const Index n1y = sensor_response.nrows();
1965  return Range( n1y*mblock_index, n1y );
1966 }
1967 
1968 
1970  bool& failed,
1971  String& fail_msg,
1972  ArrayOfArrayOfTensor4& iy_aux_array,
1973  Workspace& ws,
1974  Vector& iyb,
1975  ArrayOfMatrix& diyb_dx,
1976  const Index& mblock_index,
1977  const Index& atmosphere_dim,
1978  ConstTensor3View t_field,
1979  ConstTensor3View z_field,
1980  ConstTensor4View vmr_field,
1981  const Index& cloudbox_on,
1982  const Index& stokes_dim,
1983  ConstVectorView f_grid,
1984  ConstMatrixView sensor_pos,
1985  ConstMatrixView sensor_los,
1986  ConstMatrixView transmitter_pos,
1987  ConstVectorView mblock_za_grid,
1988  ConstVectorView mblock_aa_grid,
1989  const Index& antenna_dim,
1990  const Agenda& iy_main_agenda,
1991  const Index& j_analytical_do,
1992  const ArrayOfRetrievalQuantity& jacobian_quantities,
1993  const ArrayOfArrayOfIndex& jacobian_indices,
1994  const ArrayOfString& iy_aux_vars,
1995  const Index& naa,
1996  const Index& iza,
1997  const Index& nf)
1998 {
1999  // The try block here is necessary to correctly handle
2000  // exceptions inside the parallel region.
2001  try
2002  {
2003  for( Index iaa=0; iaa<naa; iaa++ )
2004  {
2005  //--- LOS of interest
2006  //
2007  Vector los( sensor_los.ncols() );
2008  //
2009  los = sensor_los( mblock_index, joker );
2010  los[0] += mblock_za_grid[iza];
2011  //
2012  if( antenna_dim == 2 ) // map_daa handles also "adjustment"
2013  { map_daa( los[0], los[1], los[0], los[1],
2014  mblock_aa_grid[iaa] ); }
2015  else
2016  { adjust_los( los, atmosphere_dim ); }
2017 
2018  //--- rtp_pos 1 and 2
2019  //
2020  Vector rtp_pos, rtp_pos2(0);
2021  //
2022  rtp_pos = sensor_pos( mblock_index, joker );
2023  if( transmitter_pos.nrows() )
2024  { rtp_pos2 = transmitter_pos( mblock_index, joker ); }
2025 
2026  // Calculate iy and associated variables
2027  //
2028  Matrix iy;
2029  ArrayOfTensor3 diy_dx;
2030  Ppath ppath;
2031  Tensor3 iy_transmission(0,0,0);
2032  Index iang = iza*naa + iaa;
2033  //
2034  iy_main_agendaExecute(ws, iy, iy_aux_array[iang], ppath,
2035  diy_dx, 1, iy_transmission, iy_aux_vars,
2036  cloudbox_on, j_analytical_do, t_field,
2037  z_field, vmr_field, f_grid, rtp_pos, los,
2038  rtp_pos2, iy_main_agenda );
2039 
2040  // Check that aux data can be handled and has correct size
2041  for( Index q=0; q<iy_aux_array[iang].nelem(); q++ )
2042  {
2043  if( iy_aux_array[iang][q].ncols() != 1 ||
2044  iy_aux_array[iang][q].nrows() != 1 )
2045  {
2046  throw runtime_error( "For calculations using yCalc, "
2047  "*iy_aux_vars* can not include\nvariables of "
2048  "along-the-path or extinction matrix type.");
2049  }
2050  assert( iy_aux_array[iang][q].npages() == 1 ||
2051  iy_aux_array[iang][q].npages() == stokes_dim );
2052  assert( iy_aux_array[iang][q].nbooks() == 1 ||
2053  iy_aux_array[iang][q].nbooks() == nf );
2054  }
2055 
2056  // Start row in iyb etc. for present LOS
2057  //
2058  const Index row0 = iang * nf * stokes_dim;
2059 
2060  // Jacobian part
2061  //
2062  if( j_analytical_do )
2063  {
2065  for( Index ip=0; ip<jacobian_indices[iq][1] -
2066  jacobian_indices[iq][0]+1; ip++ )
2067  {
2068  for( Index is=0; is<stokes_dim; is++ )
2069  {
2070  diyb_dx[iq](Range(row0+is,nf,stokes_dim),ip)=
2071  diy_dx[iq](ip,joker,is);
2072  }
2073  }
2074  )
2075  }
2076 
2077  // iy : copy to iyb
2078  for( Index is=0; is<stokes_dim; is++ )
2079  { iyb[Range(row0+is,nf,stokes_dim)] = iy(joker,is); }
2080 
2081  } // End aa loop
2082  } // End try
2083 
2084  catch (runtime_error e)
2085  {
2086 #pragma omp critical (iyb_calc_fail)
2087  { fail_msg = e.what(); failed = true; }
2088  }
2089 }
2090 
2091 
2093 
2102  Workspace& ws,
2103  Vector& iyb,
2104  ArrayOfVector& iyb_aux,
2105  ArrayOfMatrix& diyb_dx,
2106  const Index& mblock_index,
2107  const Index& atmosphere_dim,
2108  ConstTensor3View t_field,
2109  ConstTensor3View z_field,
2110  ConstTensor4View vmr_field,
2111  const Index& cloudbox_on,
2112  const Index& stokes_dim,
2113  ConstVectorView f_grid,
2114  ConstMatrixView sensor_pos,
2115  ConstMatrixView sensor_los,
2116  ConstMatrixView transmitter_pos,
2117  ConstVectorView mblock_za_grid,
2118  ConstVectorView mblock_aa_grid,
2119  const Index& antenna_dim,
2120  const Agenda& iy_main_agenda,
2121  const Index& j_analytical_do,
2122  const ArrayOfRetrievalQuantity& jacobian_quantities,
2123  const ArrayOfArrayOfIndex& jacobian_indices,
2124  const ArrayOfString& iy_aux_vars,
2125  const Verbosity& verbosity)
2126 {
2127  CREATE_OUT3;
2128 
2129  // Sizes
2130  const Index nf = f_grid.nelem();
2131  const Index nza = mblock_za_grid.nelem();
2132  Index naa = mblock_aa_grid.nelem();
2133  if( antenna_dim == 1 )
2134  { naa = 1; }
2135  const Index niyb = nf * nza * naa * stokes_dim;
2136  // Set up size of containers for data of 1 measurement block.
2137  // (can not be made below due to parallalisation)
2138  iyb.resize( niyb );
2139  //
2140  if( j_analytical_do )
2141  {
2142  diyb_dx.resize( jacobian_indices.nelem() );
2144  diyb_dx[iq].resize( niyb, jacobian_indices[iq][1] -
2145  jacobian_indices[iq][0] + 1 );
2146  )
2147  }
2148  else
2149  { diyb_dx.resize( 0 ); }
2150 
2151  // For iy_aux we don't know the number of quantities, and we have to store
2152  // all outout
2153  ArrayOfArrayOfTensor4 iy_aux_array( nza*naa );
2154 
2155  // We have to make a local copy of the Workspace and the agendas because
2156  // only non-reference types can be declared firstprivate in OpenMP
2157  Workspace l_ws (ws);
2158  Agenda l_iy_main_agenda (iy_main_agenda);
2159 
2160  String fail_msg;
2161  bool failed = false;
2162  if (nza >= arts_omp_get_max_threads() || nza*10 >= nf)
2163  {
2164  out3 << " Parallelizing za loop (" << nza << " iterations, "
2165  << nf << " frequencies)\n";
2166 
2167  // Start of actual calculations
2168 #pragma omp parallel for \
2169 if (!arts_omp_in_parallel()) \
2170 firstprivate(l_ws, l_iy_main_agenda)
2171  for( Index iza=0; iza<nza; iza++ )
2172  {
2173  // Skip remaining iterations if an error occurred
2174  if (failed) continue;
2175 
2176  iyb_calc_za_loop_body(failed,
2177  fail_msg,
2178  iy_aux_array,
2179  l_ws,
2180  iyb,
2181  diyb_dx,
2182  mblock_index,
2183  atmosphere_dim,
2184  t_field,
2185  z_field,
2186  vmr_field,
2187  cloudbox_on,
2188  stokes_dim,
2189  f_grid,
2190  sensor_pos,
2191  sensor_los,
2192  transmitter_pos,
2193  mblock_za_grid,
2194  mblock_aa_grid,
2195  antenna_dim,
2196  l_iy_main_agenda,
2197  j_analytical_do,
2198  jacobian_quantities,
2199  jacobian_indices,
2200  iy_aux_vars,
2201  naa,
2202  iza,
2203  nf);
2204 
2205  } // End za loop
2206  }
2207  else
2208  {
2209  out3 << " Not parallelizing za loop (" << nza << " iterations, "
2210  << nf << " frequencies)\n";
2211 
2212  for( Index iza=0; iza<nza; iza++ )
2213  {
2214  // Skip remaining iterations if an error occurred
2215  if (failed) continue;
2216 
2217  iyb_calc_za_loop_body(failed,
2218  fail_msg,
2219  iy_aux_array,
2220  ws,
2221  iyb,
2222  diyb_dx,
2223  mblock_index,
2224  atmosphere_dim,
2225  t_field,
2226  z_field,
2227  vmr_field,
2228  cloudbox_on,
2229  stokes_dim,
2230  f_grid,
2231  sensor_pos,
2232  sensor_los,
2233  transmitter_pos,
2234  mblock_za_grid,
2235  mblock_aa_grid,
2236  antenna_dim,
2237  iy_main_agenda,
2238  j_analytical_do,
2239  jacobian_quantities,
2240  jacobian_indices,
2241  iy_aux_vars,
2242  naa,
2243  iza,
2244  nf);
2245 
2246  } // End za loop
2247  }
2248 
2249  if( failed )
2250  throw runtime_error("Run-time error in function: iyb_calc\n" + fail_msg);
2251 
2252  // Compile iyb_aux
2253  //
2254  const Index nq = iy_aux_array[0].nelem();
2255  iyb_aux.resize( nq );
2256  //
2257  for( Index q=0; q<nq; q++ )
2258  {
2259  iyb_aux[q].resize( niyb );
2260  //
2261  for( Index iza=0; iza<nza; iza++ )
2262  {
2263  for( Index iaa=0; iaa<naa; iaa++ )
2264  {
2265  const Index iang = iza*naa + iaa;
2266  const Index row0 = iang * nf * stokes_dim;
2267  for( Index iv=0; iv<nf; iv++ )
2268  {
2269  const Index row1 = row0 + iv*stokes_dim;
2270  const Index i1 = min( iv, iy_aux_array[iang][q].nbooks()-1 );
2271  for( Index is=0; is<stokes_dim; is++ )
2272  {
2273  Index i2 = min( is, iy_aux_array[iang][q].npages()-1 );
2274  iyb_aux[q][row1+is] = iy_aux_array[iang][q](i1,i2,0,0);
2275  }
2276  }
2277  }
2278  }
2279  }
2280 }
2281 
2282 
2283 
2285 
2310  Tensor3& iy_trans_total,
2311  ConstTensor3View iy_trans_old,
2312  ConstTensor3View iy_trans_new )
2313 {
2314  const Index nf = iy_trans_old.npages();
2315  const Index ns = iy_trans_old.ncols();
2316 
2317  assert( ns == iy_trans_old.nrows() );
2318  assert( nf == iy_trans_new.npages() );
2319  assert( ns == iy_trans_new.nrows() );
2320  assert( ns == iy_trans_new.ncols() );
2321 
2322  iy_trans_total.resize( nf, ns, ns );
2323 
2324  for( Index iv=0; iv<nf; iv++ )
2325  {
2326  mult( iy_trans_total(iv,joker,joker), iy_trans_old(iv,joker,joker),
2327  iy_trans_new(iv,joker,joker) );
2328  }
2329 }
2330 
2331 
2332 
2335 
2348 void los3d(
2349  Vector& los3d,
2350  ConstVectorView los,
2351  const Index& atmosphere_dim )
2352 {
2353  los3d.resize(2);
2354  //
2355  los3d[0] = abs( los[0] );
2356  //
2357  if( atmosphere_dim == 1 )
2358  { los3d[1] = 0; }
2359  else if( atmosphere_dim == 2 )
2360  {
2361  if( los[0] >= 0 )
2362  { los3d[1] = 0; }
2363  else
2364  { los3d[1] = 180; }
2365  }
2366  else if( atmosphere_dim == 3 )
2367  { los3d[1] = los[1]; }
2368 }
2369 
2370 
2371 
2373 
2390  Vector& los_mirrored,
2391  ConstVectorView los,
2392  const Index& atmosphere_dim )
2393 {
2394  los_mirrored.resize(2);
2395  //
2396  if( atmosphere_dim == 1 )
2397  {
2398  los_mirrored[0] = 180 - los[0];
2399  los_mirrored[1] = 180;
2400  }
2401  else if( atmosphere_dim == 2 )
2402  {
2403  los_mirrored[0] = 180 - fabs( los[0] );
2404  if( los[0] >= 0 )
2405  { los_mirrored[1] = 180; }
2406  else
2407  { los_mirrored[1] = 0; }
2408  }
2409  else if( atmosphere_dim == 3 )
2410  {
2411  los_mirrored[0] = 180 - los[0];
2412  los_mirrored[1] = los[1] + 180;
2413  if( los_mirrored[1] > 180 )
2414  { los_mirrored[1] -= 360; }
2415  }
2416 }
2417 
2418 
2419 
2421 
2439  Numeric& lat,
2440  Numeric& lon,
2441  const Index& atmosphere_dim,
2442  ConstVectorView lat_grid,
2443  ConstVectorView lat_true,
2444  ConstVectorView lon_true,
2445  ConstVectorView pos )
2446 {
2447  assert( pos.nelem() == atmosphere_dim );
2448 
2449  if( atmosphere_dim == 1 )
2450  {
2451  assert( lat_true.nelem() == 1 );
2452  assert( lon_true.nelem() == 1 );
2453  //
2454  lat = lat_true[0];
2455  lon = lon_true[0];
2456  }
2457 
2458  else if( atmosphere_dim == 2 )
2459  {
2460  assert( lat_true.nelem() == lat_grid.nelem() );
2461  assert( lon_true.nelem() == lat_grid.nelem() );
2462  GridPos gp;
2463  Vector itw(2);
2464  gridpos( gp, lat_grid, pos[1] );
2465  interpweights( itw, gp );
2466  lat = interp( itw, lat_true, gp );
2467  lon = interp( itw, lon_true, gp );
2468  }
2469 
2470  else
2471  {
2472  lat = pos[1];
2473  lon = pos[2];
2474  }
2475 }
2476 
2477 
2479 
2497  Matrix& iy,
2498  ConstTensor3View I,
2499  ConstMatrixView surface_los,
2500  ConstTensor4View surface_rmatrix,
2501  ConstMatrixView surface_emission )
2502 {
2503  // Some sizes
2504  const Index nf = I.nrows();
2505  const Index stokes_dim = I.ncols();
2506  const Index nlos = surface_los.nrows();
2507 
2508  iy = surface_emission;
2509 
2510  // Loop *surface_los*-es. If no such LOS, we are ready.
2511  if( nlos > 0 )
2512  {
2513  for( Index ilos=0; ilos<nlos; ilos++ )
2514  {
2515  Vector rtmp(stokes_dim); // Reflected Stokes vector for 1 frequency
2516 
2517  for( Index iv=0; iv<nf; iv++ )
2518  {
2519  mult( rtmp, surface_rmatrix(ilos,iv,joker,joker), I(ilos,iv,joker) );
2520  iy(iv,joker) += rtmp;
2521  }
2522  }
2523  }
2524 }
2525 
2526 
2527 
2529 
2543  Numeric& l,
2544  Vector& los,
2545  const Numeric& u,
2546  const Numeric& v,
2547  const Numeric& w )
2548 {
2549  l= sqrt( u*u + v*v + w*w );
2550  //
2551  los.resize(2);
2552  //
2553  los[0] = acos( w / l );
2554  los[1] = atan2( u, v );
2555 }
2556 
2557 
2558 
Numeric dotprod_with_los(ConstVectorView los, const Numeric &u, const Numeric &v, const Numeric &w, const Index &atmosphere_dim)
dotprod_with_los
Definition: rte.cc:776
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
ArrayOfGridPos gp_lat
Definition: ppath.h:77
The VectorView class.
Definition: matpackI.h:372
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:201
Numeric constant
Definition: ppath.h:62
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
The Agenda class.
Definition: agenda_class.h:44
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:194
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
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:2430
void surface_calc(Matrix &iy, ConstTensor3View I, ConstMatrixView surface_los, ConstTensor4View surface_rmatrix, ConstMatrixView surface_emission)
surface_calc
Definition: rte.cc:2496
#define q
Definition: continua.cc:21469
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
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
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
zaaa2cart
Definition: ppath.cc:484
Matrix los
Definition: ppath.h:68
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Vector end_pos
Definition: ppath.h:71
Index dim
Definition: ppath.h:60
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 Sparse class.
Definition: matpackII.h:55
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
Vector lstep
Definition: ppath.h:70
void iyb_calc_za_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfTensor4 &iy_aux_array, Workspace &ws, Vector &iyb, 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 Index &naa, const Index &iza, const Index &nf)
Definition: rte.cc:1969
Linear algebra functions.
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:47
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
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
void cart2zaaa(Numeric &za, Numeric &aa, const Numeric &dx, const Numeric &dy, const Numeric &dz)
cart2zaaa
Definition: ppath.cc:443
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
Matrix pos
Definition: ppath.h:67
void defocusing_general(Workspace &ws, Numeric &dlf, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Ppath &ppath, const Numeric &ppath_lraytrace, const Numeric &dza, const Verbosity &verbosity)
defocusing_general
Definition: rte.cc:530
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
Numeric dinvplanckdI(const Numeric &i, const Numeric &f)
dinvplanckdI
void defocusing_general_sub(Workspace &ws, Vector &pos, Vector &rte_los, Index &background, ConstVectorView rte_pos, const Numeric &lo0, const Agenda &ppath_step_agenda, const Numeric &ppath_lraytrace, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Verbosity &verbosity)
defocusing_general_sub
Definition: rte.cc:404
Vector start_pos
Definition: ppath.h:64
Vector r
Definition: ppath.h:69
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:146
void bending_angle1d(Numeric &alpha, const Ppath &ppath)
bending_angle1d
Definition: rte.cc:353
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:75
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:114
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
const Numeric SPEED_OF_LIGHT
Structure to store a grid position.
Definition: interpolation.h:74
A constant view of a Tensor4.
Definition: matpackIV.h:141
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:40
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
Numeric end_lstep
Definition: ppath.h:73
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:768
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
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 los3d(Vector &los3d, ConstVectorView los, const Index &atmosphere_dim)
iy_transmission_mult_scalar_tau los3d
Definition: rte.cc:2348
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
void map_daa(Numeric &za, Numeric &aa, const Numeric &za0, const Numeric &aa0, const Numeric &aa_grid)
Definition: ppath.cc:576
Vector end_los
Definition: ppath.h:72
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
String background
Definition: ppath.h:63
#define DEBUG_ONLY(...)
Definition: debug.h:37
void pos2true_latlon(Numeric &lat, Numeric &lon, const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true, ConstVectorView pos)
pos2true_latlon
Definition: rte.cc:2438
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.
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1276
The Tensor3View class.
Definition: matpackIII.h:232
#define max(a, b)
Definition: continua.cc:20461
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
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:56
Numeric start_lstep
Definition: ppath.h:66
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void propmat_clearsky_agendaExecute(Workspace &ws, Tensor4 &propmat_clearsky, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:13039
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:13886
#define abs(x)
Definition: continua.cc:20458
const Joker joker
void find_tanpoint(Index &it, const Ppath ppath)
find_tanpoint
Definition: ppath.cc:705
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
void blackbody_radiation_agendaExecute(Workspace &ws, Vector &blackbody_radiation, const Numeric rtp_temperature, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:13205
#define dx
Definition: continua.cc:21928
void opt_propCalc(MatrixView ext_mat_mono, VectorView abs_vec_mono, const Numeric za, const Numeric aa, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index stokes_dim, ConstVectorView pnd_vec, const Numeric rtp_temperature, const Verbosity &verbosity)
opt_propCalc
Definition: montecarlo.cc:684
The Matrix class.
Definition: matpackI.h:788
Vector nreal
Definition: ppath.h:74
void distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
Definition: geodetic.cc:658
Array< String > ArrayOfString
An array of Strings.
Definition: mystring.h:321
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1199
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
interp_atmfield_gp2itw
Header file for special_interp.cc.
Propagation path structure and functions.
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
Header file for logic.cc.
void emission_rtstep(Matrix &iy, const Index &stokes_dim, ConstVectorView bbar, ArrayOfIndex &extmat_case, ConstTensor3View t)
emission_rtstep
Definition: rte.cc:822
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
This can be used to make arrays out of anything.
Definition: array.h:40
bool is_singular(ConstMatrixView A)
Checks if a square matrix is singular.
Definition: logic.cc:377
void adjust_los(VectorView los, const Index &atmosphere_dim)
adjust_los
Definition: rte.cc:81
void vectorfield2los(Numeric &l, Vector &los, const Numeric &u, const Numeric &v, const Numeric &w)
vectorfield2los
Definition: rte.cc:2542
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:82
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
A constant view of a Tensor3.
Definition: matpackIII.h:139
A constant view of a Vector.
Definition: matpackI.h:292
Index np
Definition: ppath.h:61
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
ArrayOfGridPos gp_lon
Definition: ppath.h:78
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
get_rowindex_for_mblock
Definition: rte.cc:1960
Numeric mean(const ConstVectorView &x)
Mean function, vector version.
Definition: matpackI.cc:1972
A constant view of a Matrix.
Definition: matpackI.h:596
Vector start_los
Definition: ppath.h:65
Numeric invplanck(const Numeric &i, const Numeric &f)
invplanck
Workspace class.
Definition: workspace_ng.h:47
void cart2sph(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
cart2sph
Definition: geodetic.cc:596
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
void iy_cloudbox_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:13702
void iy_surface_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmission, 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 Agenda &iy_main_agenda, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition: auto_md.cc:14070
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
#define CREATE_OUT3
Definition: messages.h:214
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 ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
ppath_calc
Definition: ppath.cc:5344
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
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
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
void get_iy(Workspace &ws, Matrix &iy, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, const Index &cloudbox_on, ConstVectorView f_grid, ConstVectorView rte_pos, ConstVectorView rte_los, ConstVectorView rte_pos2, const Agenda &iy_main_agenda)
get_iy
Definition: rte.cc:1046
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 poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:387
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
ArrayOfGridPos gp_p
Definition: ppath.h:76
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
void defocusing_sat2sat(Workspace &ws, Numeric &dlf, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Ppath &ppath, const Numeric &ppath_lraytrace, const Numeric &dza, const Verbosity &verbosity)
defocusing_sat2sat
Definition: rte.cc:663
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1403
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Definition: rte.cc:903
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
itw2p
The Tensor5 class.
Definition: matpackV.h:451
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580