ARTS  2.2.66
m_transmitter.cc
Go to the documentation of this file.
1 /* Copyright (C) 2012
2  Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3  Stefan Buehler <sbuehler(at)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 
40 /*===========================================================================
41  === External declarations
42  ===========================================================================*/
43 
44 #include <cmath>
45 #include <stdexcept>
46 #include "arts.h"
47 #include "auto_md.h"
48 #include "complex.h"
49 #include "geodetic.h"
50 #include "jacobian.h"
51 #include "lin_alg.h"
52 #include "logic.h"
53 #include "math_funcs.h"
54 #include "messages.h"
55 #include "rte.h"
56 #include "sensor.h"
57 
58 extern const Numeric DEG2RAD;
59 extern const Numeric PI;
60 extern const Numeric RAD2DEG;
61 extern const Numeric SPEED_OF_LIGHT;
62 
63 
64 
65 /* Workspace method: Doxygen documentation will be auto-generated */
66 void iyRadioLink(
67  Workspace& ws,
68  Matrix& iy,
69  ArrayOfTensor4& iy_aux,
70  Ppath& ppath,
71  ArrayOfTensor3& diy_dx,
72  const Index& stokes_dim,
73  const Vector& f_grid,
74  const Index& atmosphere_dim,
75  const Vector& p_grid,
76  const Vector& lat_grid,
77  const Vector& lon_grid,
78  const Tensor3& z_field,
79  const Tensor3& t_field,
80  const Tensor4& vmr_field,
81  const ArrayOfArrayOfSpeciesTag& abs_species,
82  const Tensor3& wind_u_field,
83  const Tensor3& wind_v_field,
84  const Tensor3& wind_w_field,
85  const Tensor3& mag_u_field,
86  const Tensor3& mag_v_field,
87  const Tensor3& mag_w_field,
88  const Vector& refellipsoid,
89  const Matrix& z_surface,
90  const Index& cloudbox_on,
91  const ArrayOfIndex& cloudbox_limits,
92  const Tensor4& pnd_field,
93  const Index& use_mean_scat_data,
94  const ArrayOfSingleScatteringData& scat_data_array,
95  const Matrix& particle_masses,
96  const ArrayOfString& iy_aux_vars,
97  const Index& jacobian_do,
98  const Agenda& ppath_agenda,
99  const Agenda& ppath_step_agenda,
100  const Agenda& propmat_clearsky_agenda,
101  const Agenda& iy_transmitter_agenda,
102  const Index& iy_agenda_call1,
103  const Tensor3& iy_transmission,
104  const Vector& rte_pos,
105  const Vector& rte_los,
106  const Vector& rte_pos2,
107  const Numeric& rte_alonglos_v,
108  const Numeric& ppath_lraytrace,
109  const Index& defocus_method,
110  const Numeric& defocus_shift,
111  const Verbosity& verbosity )
112 {
113  // Throw error if unsupported features are requested
114  if( !iy_agenda_call1 )
115  throw runtime_error(
116  "Recursive usage not possible (iy_agenda_call1 must be 1)" );
117  if( iy_transmission.ncols() )
118  throw runtime_error( "*iy_transmission* must be empty" );
119  if( jacobian_do )
120  throw runtime_error( "This method does not provide any jacobians and "
121  "*jacobian_do* must be 0." );
122  if( defocus_method < 1 || defocus_method > 2 )
123  throw runtime_error( "Allowed choices for *defocus_method* is 1 and 2." );
124  diy_dx.resize(0);
125 
126 
127  //- Determine propagation path
128  ppath_agendaExecute( ws, ppath, ppath_lraytrace, rte_pos, rte_los,
129  rte_pos2, cloudbox_on, 0, t_field, z_field, vmr_field,
130  f_grid, ppath_agenda );
131 
132  //- Set np to zero if ground intersection
133  const Index radback = ppath_what_background(ppath);
134  // np should already be 1 fon non-OK cases, but for extra safety ...
135  if( radback == 0 || radback == 2 )
136  { ppath.np = 1; }
137 
138  // Some basic sizes
139  //
140  const Index nf = f_grid.nelem();
141  const Index ns = stokes_dim;
142  const Index np = ppath.np;
143 
144 
145  // The below copied from iyEmission. Activate for-loop when jacobians
146  // introduced.
147 
148  // Set up variable with index of species where we want abs_per_species.
149  // This variable can below be extended in iy_aux part.
150  //
151  ArrayOfIndex iaps(0);
152  //
153  //for( Index i=0; i<jac_species_i.nelem(); i++ )
154  // {
155  // if( jac_species_i[i] >= 0 )
156  // { iaps.push_back( jac_species_i[i] ); }
157  // }
158 
159  //=== iy_aux part ===========================================================
160  Index auxPressure = -1,
161  auxTemperature = -1,
162  auxAbsSum = -1,
163  auxPartExt = -1,
164  auxImpactParam = -1,
165  auxFreeSpaceLoss = -1,
166  auxFreeSpaceAtte = -1,
167  auxAtmosphericLoss = -1,
168  auxDefocusingLoss = -1,
169  auxFarRotTotal = -1,
170  auxFarRotSpeed = -1,
171  auxExtraPathDelay = -1,
172  auxBendingAngle = -1;
173  ArrayOfIndex auxAbsSpecies(0), auxAbsIsp(0);
174  ArrayOfIndex auxVmrSpecies(0), auxVmrIsp(0);
175  ArrayOfIndex auxPartCont(0), auxPartContI(0);
176  ArrayOfIndex auxPartField(0), auxPartFieldI(0);
177  //
178  const Index naux = iy_aux_vars.nelem();
179  iy_aux.resize( naux );
180  //
181  for( Index i=0; i<naux; i++ )
182  {
183  if( iy_aux_vars[i] == "Pressure" )
184  { auxPressure = i; iy_aux[i].resize( 1, 1, 1, np ); }
185  else if( iy_aux_vars[i] == "Temperature" )
186  { auxTemperature = i; iy_aux[i].resize( 1, 1, 1, np ); }
187  else if( iy_aux_vars[i].substr(0,13) == "VMR, species " )
188  {
189  Index ispecies;
190  istringstream is(iy_aux_vars[i].substr(13,2));
191  is >> ispecies;
192  if( ispecies < 0 || ispecies>=abs_species.nelem() )
193  {
194  ostringstream os;
195  os << "You have selected VMR of species with index "
196  << ispecies << ".\nThis species does not exist!";
197  throw runtime_error( os.str() );
198  }
199  auxVmrSpecies.push_back(i);
200  auxVmrIsp.push_back(ispecies);
201  iy_aux[i].resize( 1, 1, 1, np );
202  }
203  else if( iy_aux_vars[i] == "Absorption, summed" )
204  { auxAbsSum = i; iy_aux[i].resize( nf, ns, ns, np ); }
205  else if( iy_aux_vars[i] == "Particle extinction, summed" )
206  {
207  auxPartExt = i;
208  iy_aux[i].resize( nf, ns, ns, np );
209  iy_aux[i] = 0;
210  }
211  else if( iy_aux_vars[i].substr(0,20) == "Absorption, species " )
212  {
213  Index ispecies;
214  istringstream is(iy_aux_vars[i].substr(20,2));
215  is >> ispecies;
216  if( ispecies < 0 || ispecies>=abs_species.nelem() )
217  {
218  ostringstream os;
219  os << "You have selected absorption species with index "
220  << ispecies << ".\nThis species does not exist!";
221  throw runtime_error( os.str() );
222  }
223  auxAbsSpecies.push_back(i);
224  const Index ihit = find_first( iaps, ispecies );
225  if( ihit >= 0 )
226  { auxAbsIsp.push_back( ihit ); }
227  else
228  {
229  iaps.push_back(ispecies);
230  auxAbsIsp.push_back( iaps.nelem()-1 );
231  }
232  iy_aux[i].resize( nf, ns, ns, np );
233  }
234  else if( iy_aux_vars[i].substr(0,14) == "Mass content, " )
235  {
236  Index icont;
237  istringstream is(iy_aux_vars[i].substr(14,2));
238  is >> icont;
239  if( icont < 0 || icont>=particle_masses.ncols() )
240  {
241  ostringstream os;
242  os << "You have selected particle mass content category with "
243  << "index " << icont << ".\nThis category is not defined!";
244  throw runtime_error( os.str() );
245  }
246  auxPartCont.push_back(i);
247  auxPartContI.push_back(icont);
248  iy_aux[i].resize( 1, 1, 1, np );
249  }
250  else if( iy_aux_vars[i].substr(0,10) == "PND, type " )
251  {
252  Index ip;
253  istringstream is(iy_aux_vars[i].substr(10,2));
254  is >> ip;
255  if( ip < 0 || ip>=pnd_field.nbooks() )
256  {
257  ostringstream os;
258  os << "You have selected particle number density field with "
259  << "index " << ip << ".\nThis field is not defined!";
260  throw runtime_error( os.str() );
261  }
262  auxPartField.push_back(i);
263  auxPartFieldI.push_back(ip);
264  iy_aux[i].resize( 1, 1, 1, np );
265  }
266  else if( iy_aux_vars[i] == "Impact parameter" )
267  { auxImpactParam = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
268  else if( iy_aux_vars[i] == "Free space loss" )
269  { auxFreeSpaceLoss = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
270  else if( iy_aux_vars[i] == "Free space attenuation" )
271  { auxFreeSpaceAtte = i; iy_aux[i].resize( 1, 1, 1, np ); }
272  else if( iy_aux_vars[i] == "Atmospheric loss" )
273  { auxAtmosphericLoss = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
274  else if( iy_aux_vars[i] == "Defocusing loss" )
275  { auxDefocusingLoss = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
276  else if( iy_aux_vars[i] == "Faraday rotation" )
277  { auxFarRotTotal = i; iy_aux[i].resize( nf, 1, 1, 1 ); iy_aux[i] = 0; }
278  else if( iy_aux_vars[i] == "Faraday speed" )
279  { auxFarRotSpeed = i; iy_aux[i].resize( nf, 1, 1, np ); iy_aux[i] = 0; }
280  else if( iy_aux_vars[i] == "Extra path delay" )
281  { auxExtraPathDelay = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
282  else if( iy_aux_vars[i] == "Bending angle" )
283  { auxBendingAngle = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
284  else
285  {
286  ostringstream os;
287  os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
288  << "\"\nThis choice is not recognised.";
289  throw runtime_error( os.str() );
290  }
291  }
292 
293  // Special stuff to handle Faraday rotation
294  Index ife = -1; // When ready, ife refers abs_per_species
295  if( auxFarRotTotal>=0 || auxFarRotSpeed>=0 )
296  {
297  if( stokes_dim < 3 )
298  throw runtime_error(
299  "To include Faraday rotation, stokes_dim >= 3 is required." );
300 
301  // Determine species index of free electrons
302  for( Index sp = 0; sp < abs_species.nelem() && ife < 0; sp++ )
303  {
304  if (abs_species[sp][0].Type() == SpeciesTag::TYPE_FREE_ELECTRONS)
305  { ife = sp; }
306  }
307  // If not found, then aux values already set to zero
308  if( ife < 0 )
309  {
310  auxFarRotTotal = -1;
311  auxFarRotSpeed = -1;
312  }
313  else
314  {
315  const Index ihit = find_first( iaps, ife );
316  if( ihit >= 0 )
317  { ife = ihit; }
318  else
319  {
320  iaps.push_back(ife);
321  ife = iaps.nelem() - 1;
322  }
323  }
324  }
325  //===========================================================================
326 
327 
328  // Handle cases whn no link was establsihed
329  if( radback == 0 || radback == 2 )
330  {
331  Numeric fillvalue = 0;
332  if( radback == 0 )
333  { fillvalue = NAN; }
334  //
335  iy.resize( nf, stokes_dim );
336  iy = fillvalue;
337  //
338  for( Index i=0; i<naux; i++ )
339  { iy_aux[i] = fillvalue; }
340  //
341  return;
342  }
343 
344 
345  // Transmitted signal
346  //
347  iy_transmitter_agendaExecute( ws, iy, f_grid,
348  ppath.pos(np-1,Range(0,atmosphere_dim)),
349  ppath.los(np-1,joker), iy_transmitter_agenda );
350  if( iy.ncols() != stokes_dim || iy.nrows() != nf )
351  { throw runtime_error( "The size of *iy* returned from "
352  "*iy_transmitter_agenda* is not correct." ); }
353 
354 
355  // Get atmospheric and attenuation quantities for each ppath point/step
356  //
357  Vector ppath_p, ppath_t;
358  Matrix ppath_vmr, ppath_pnd, ppath_mag, ppath_wind, ppath_f;
359  Tensor5 abs_per_species;
360  Tensor4 ppath_abs, trans_partial, trans_cumulat, pnd_ext_mat;
361  Vector scalar_tau;
362  ArrayOfIndex clear2cloudbox;
363  //
364  if( np > 1 )
365  {
366  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
367  ppath_wind, ppath_mag,
368  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
369  wind_u_field, wind_v_field, wind_w_field,
370  mag_u_field, mag_v_field, mag_w_field );
371  get_ppath_f( ppath_f, ppath, f_grid, atmosphere_dim,
372  rte_alonglos_v, ppath_wind );
373  get_ppath_abs( ws, ppath_abs, abs_per_species,
374  propmat_clearsky_agenda, ppath,
375  ppath_p, ppath_t, ppath_vmr, ppath_f,
376  ppath_mag, f_grid, stokes_dim, iaps );
377  if( !cloudbox_on )
378  {
379  ArrayOfArrayOfIndex extmat_case;
380  get_ppath_trans( trans_partial, extmat_case, trans_cumulat,
381  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim );
382  }
383  else
384  {
385  Array<ArrayOfSingleScatteringData> scat_data;
386  ArrayOfArrayOfIndex extmat_case;
387  Tensor3 pnd_abs_vec;
388  //
389  get_ppath_ext( clear2cloudbox, pnd_abs_vec, pnd_ext_mat,
390  scat_data, ppath_pnd, ppath, ppath_t, stokes_dim,
391  ppath_f, atmosphere_dim, cloudbox_limits, pnd_field,
392  use_mean_scat_data, scat_data_array, verbosity );
393  get_ppath_trans2( trans_partial, extmat_case, trans_cumulat,
394  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim,
395  clear2cloudbox, pnd_ext_mat );
396  }
397  }
398 
399  // Ppath length variables
400  //
401  Numeric lbg; // Bent geometrical length of ray path
402  Numeric lba; // Bent apparent length of ray path
403  //
404  lbg = ppath.end_lstep;
405  lba = lbg;
406 
407  // Do RT calculations
408  //
409  if( np > 1 )
410  {
411  //=== iy_aux part =======================================================
412  // iy_aux for point np-1:
413  // Pressure
414  if( auxPressure >= 0 )
415  { iy_aux[auxPressure](0,0,0,np-1) = ppath_p[np-1]; }
416  // Temperature
417  if( auxTemperature >= 0 )
418  { iy_aux[auxTemperature](0,0,0,np-1) = ppath_t[np-1]; }
419  // VMR
420  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
421  { iy_aux[auxVmrSpecies[j]](0,0,0,np-1) = ppath_vmr(auxVmrIsp[j],np-1);}
422  // Absorption
423  if( auxAbsSum >= 0 )
424  { for( Index iv=0; iv<nf; iv++ ) {
425  for( Index is1=0; is1<ns; is1++ ){
426  for( Index is2=0; is2<ns; is2++ ){
427  iy_aux[auxAbsSum](iv,is1,is2,np-1) =
428  ppath_abs(iv,is1,is2,np-1); } } } }
429  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
430  { for( Index iv=0; iv<nf; iv++ ) {
431  for( Index is1=0; is1<stokes_dim; is1++ ){
432  for( Index is2=0; is2<stokes_dim; is2++ ){
433  iy_aux[auxAbsSpecies[j]](iv,is1,is2,np-1) =
434  abs_per_species(auxAbsIsp[j],iv,is1,is2,np-1); } } } }
435  // Particle properties
436  if( cloudbox_on )
437  {
438  // Extinction
439  if( auxPartExt >= 0 && clear2cloudbox[np-1] >= 0 )
440  {
441  const Index ic = clear2cloudbox[np-1];
442  for( Index iv=0; iv<nf; iv++ ) {
443  for( Index is1=0; is1<ns; is1++ ){
444  for( Index is2=0; is2<ns; is2++ ){
445  iy_aux[auxPartExt](iv,is1,is2,np-1) =
446  pnd_ext_mat(iv,is1,is2,ic); } } }
447  }
448  // Particle mass content
449  for( Index j=0; j<auxPartCont.nelem(); j++ )
450  { iy_aux[auxPartCont[j]](0,0,0,np-1) = ppath_pnd(joker,np-1) *
451  particle_masses(joker,auxPartContI[j]); }
452  // Particle field
453  for( Index j=0; j<auxPartField.nelem(); j++ )
454  { iy_aux[auxPartField[j]](0,0,0,np-1) =
455  ppath_pnd(auxPartFieldI[j],np-1); }
456  }
457  // Free space
458  if( auxFreeSpaceAtte >= 0 )
459  { iy_aux[auxFreeSpaceAtte](joker,0,0,np-1) = 2/lbg; }
460  // Faraday speed
461  if( auxFarRotSpeed >= 0 )
462  { for( Index iv=0; iv<nf; iv++ ) {
463  iy_aux[auxFarRotSpeed](iv,0,0,np-1) = 0.5 *
464  abs_per_species(ife,iv,1,2,np-1); } }
465  //=======================================================================
466 
467  // Loop ppath steps
468  for( Index ip=np-2; ip>=0; ip-- )
469  {
470  // Lengths
471  lbg += ppath.lstep[ip];
472  lba += ppath.lstep[ip] * (ppath.ngroup[ip]+ppath.ngroup[ip+1]) / 2.0;
473 
474  // Atmospheric loss of path step + Faraday rotation
475  if( stokes_dim == 1 )
476  {
477  for( Index iv=0; iv<nf; iv++ )
478  { iy(iv,0) = iy(iv,0) * trans_partial(iv,0,0,ip); }
479  }
480  else
481  {
482  for( Index iv=0; iv<nf; iv++ )
483  {
484  // Unpolarised:
485  if( is_diagonal( trans_partial(iv,joker,joker,ip) ) )
486  {
487  for( Index is=0; is<ns; is++ )
488  { iy(iv,is) = iy(iv,is) * trans_partial(iv,is,is,ip); }
489  }
490  // The general case:
491  else
492  {
493  Vector t1(ns);
494  mult( t1, trans_partial(iv,joker,joker,ip), iy(iv,joker));
495  iy(iv,joker) = t1;
496  }
497  }
498  }
499 
500  //=== iy_aux part ===================================================
501  // Pressure
502  if( auxPressure >= 0 )
503  { iy_aux[auxPressure](0,0,0,ip) = ppath_p[ip]; }
504  // Temperature
505  if( auxTemperature >= 0 )
506  { iy_aux[auxTemperature](0,0,0,ip) = ppath_t[ip]; }
507  // VMR
508  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
509  { iy_aux[auxVmrSpecies[j]](0,0,0,ip) = ppath_vmr(auxVmrIsp[j],ip);}
510  // Absorption
511  if( auxAbsSum >= 0 )
512  { for( Index iv=0; iv<nf; iv++ ) {
513  for( Index is1=0; is1<ns; is1++ ){
514  for( Index is2=0; is2<ns; is2++ ){
515  iy_aux[auxAbsSum](iv,is1,is2,ip) =
516  ppath_abs(iv,is1,is2,ip); } } } }
517  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
518  { for( Index iv=0; iv<nf; iv++ ) {
519  for( Index is1=0; is1<stokes_dim; is1++ ){
520  for( Index is2=0; is2<stokes_dim; is2++ ){
521  iy_aux[auxAbsSpecies[j]](iv,is1,is2,ip) =
522  abs_per_species(auxAbsIsp[j],iv,is1,is2,ip); } } } }
523  // Particle properties
524  if( cloudbox_on )
525  {
526  // Extinction
527  if( auxPartExt >= 0 && clear2cloudbox[ip] >= 0 )
528  {
529  const Index ic = clear2cloudbox[ip];
530  for( Index iv=0; iv<nf; iv++ ) {
531  for( Index is1=0; is1<ns; is1++ ){
532  for( Index is2=0; is2<ns; is2++ ){
533  iy_aux[auxPartExt](iv,is1,is2,ip) =
534  pnd_ext_mat(iv,is1,is2,ic); } } }
535  }
536  // Particle mass content
537  for( Index j=0; j<auxPartCont.nelem(); j++ )
538  { iy_aux[auxPartCont[j]](0,0,0,ip) = ppath_pnd(joker,ip) *
539  particle_masses(joker,auxPartContI[j]); }
540  // Particle field
541  for( Index j=0; j<auxPartField.nelem(); j++ )
542  { iy_aux[auxPartField[j]](0,0,0,ip) =
543  ppath_pnd(auxPartFieldI[j],ip); }
544  }
545  // Free space loss
546  if( auxFreeSpaceAtte >= 0 )
547  { iy_aux[auxFreeSpaceAtte](joker,0,0,ip) = 2/lbg; }
548  // Faraday rotation, total
549  if( auxFarRotTotal >= 0 )
550  { for( Index iv=0; iv<nf; iv++ ) {
551  iy_aux[auxFarRotTotal](iv,0,0,0) += RAD2DEG * ppath.lstep[ip] *
552  0.25 * ( abs_per_species(ife,iv,1,2,ip)+
553  abs_per_species(ife,iv,1,2,ip+1)); } }
554  // Faraday speed
555  if( auxFarRotSpeed >= 0 )
556  { for( Index iv=0; iv<nf; iv++ ) {
557  iy_aux[auxFarRotSpeed](iv,0,0,ip) = 0.5 *
558  abs_per_species(ife,iv,1,2,ip); } }
559  //===================================================================
560  }
561 
562 
563  //=== iy_aux part =======================================================
564  if( auxAtmosphericLoss >= 0 )
565  { iy_aux[auxAtmosphericLoss](joker,0,0,0) = iy(joker,0); }
566  if( auxImpactParam >= 0 )
567  {
568  assert( ppath.constant >= 0 );
569  iy_aux[auxImpactParam](joker,0,0,0) = ppath.constant;
570  }
571  //=======================================================================
572 
573 
574  // Remaing length of ppath
575  lbg += ppath.start_lstep;
576  lba += ppath.start_lstep;
577 
578 
579  // Determine total free space loss
580  Numeric fspl = 1 / ( 4 * PI * lbg*lbg );
581  //
582  if( auxFreeSpaceLoss >= 0 )
583  { iy_aux[auxFreeSpaceLoss] = fspl; }
584 
585 
586  // Determine defocusing loss
587  Numeric dfl = 1;
588  if( defocus_method == 1 )
589  {
590  defocusing_general( ws, dfl, ppath_step_agenda, atmosphere_dim,
591  p_grid, lat_grid, lon_grid, t_field, z_field,
592  vmr_field, f_grid, refellipsoid,
593  z_surface, ppath, ppath_lraytrace,
594  defocus_shift, verbosity );
595  }
596  else if( defocus_method == 2 )
597  { defocusing_sat2sat( ws, dfl, ppath_step_agenda, atmosphere_dim,
598  p_grid, lat_grid, lon_grid, t_field, z_field,
599  vmr_field, f_grid, refellipsoid,
600  z_surface, ppath, ppath_lraytrace,
601  defocus_shift, verbosity );
602  }
603  if( auxDefocusingLoss >= 0 )
604  { iy_aux[auxDefocusingLoss] = dfl; }
605 
606 
607 
608  // Include free space and defocusing losses
609  iy *= fspl*dfl;
610 
611 
612  // Extra path delay
613  if( auxExtraPathDelay >= 0 )
614  {
615  // Radius of rte_pos and rte_pos2
616  const Numeric r1 = ppath.end_pos[0] +
617  pos2refell_r( atmosphere_dim, refellipsoid,
618  lat_grid, lon_grid, ppath.end_pos );
619  const Numeric r2 = ppath.start_pos[0] +
620  pos2refell_r( atmosphere_dim, refellipsoid,
621  lat_grid, lon_grid, ppath.start_pos );
622 
623  // Geometrical distance between start and end point
624  Numeric lgd ;
625  if( atmosphere_dim <= 2 )
626  { distance2D( lgd, r1, ppath.end_pos[1], r2, ppath.start_pos[1] ); }
627  else
628  { distance3D( lgd, r1, ppath.end_pos[1], ppath.end_pos[2],
629  r2, ppath.start_pos[1], ppath.start_pos[2] ); }
630  //
631  iy_aux[auxExtraPathDelay] = ( lba - lgd ) / SPEED_OF_LIGHT;
632  }
633 
634 
635  // Bending angle
636  if( auxBendingAngle >= 0 )
637  {
638  Numeric ba = -999;
639  bending_angle1d( ba, ppath );
640  //
641  iy_aux[auxBendingAngle] = ba;
642  }
643  }
644 }
645 
646 
647 
648 
649 
650 /* Workspace method: Doxygen documentation will be auto-generated */
651 void iyTransmissionStandard(
652  Workspace& ws,
653  Matrix& iy,
654  ArrayOfTensor4& iy_aux,
655  Ppath& ppath,
656  ArrayOfTensor3& diy_dx,
657  const Index& stokes_dim,
658  const Vector& f_grid,
659  const Index& atmosphere_dim,
660  const Vector& p_grid,
661  const Tensor3& z_field,
662  const Tensor3& t_field,
663  const Tensor4& vmr_field,
664  const ArrayOfArrayOfSpeciesTag& abs_species,
665  const Tensor3& wind_u_field,
666  const Tensor3& wind_v_field,
667  const Tensor3& wind_w_field,
668  const Tensor3& mag_u_field,
669  const Tensor3& mag_v_field,
670  const Tensor3& mag_w_field,
671  const Index& cloudbox_on,
672  const ArrayOfIndex& cloudbox_limits,
673  const Tensor4& pnd_field,
674  const Index& use_mean_scat_data,
675  const ArrayOfSingleScatteringData& scat_data_array,
676  const Matrix& particle_masses,
677  const ArrayOfString& iy_aux_vars,
678  const Index& jacobian_do,
679  const ArrayOfRetrievalQuantity& jacobian_quantities,
680  const ArrayOfArrayOfIndex& jacobian_indices,
681  const Agenda& ppath_agenda,
682  const Agenda& propmat_clearsky_agenda,
683  const Agenda& iy_transmitter_agenda,
684  const Index& iy_agenda_call1,
685  const Tensor3& iy_transmission,
686  const Vector& rte_pos,
687  const Vector& rte_los,
688  const Vector& rte_pos2,
689  const Numeric& rte_alonglos_v,
690  const Numeric& ppath_lraytrace,
691  const Verbosity& verbosity )
692 {
693  // Throw error if unsupported features are requested
694  if( !iy_agenda_call1 )
695  throw runtime_error(
696  "Recursive usage not possible (iy_agenda_call1 must be 1)" );
697  if( iy_transmission.ncols() )
698  throw runtime_error( "*iy_transmission* must be empty" );
699 
700 
701  // Determine propagation path
702  //
703  ppath_agendaExecute( ws, ppath, ppath_lraytrace, rte_pos, rte_los, rte_pos2,
704  0, 0, t_field, z_field, vmr_field, f_grid,
705  ppath_agenda );
706 
707  // Some basic sizes
708  //
709  const Index nf = f_grid.nelem();
710  const Index ns = stokes_dim;
711  const Index np = ppath.np;
712  const Index nq = jacobian_quantities.nelem();
713 
714  // Transmitted signal
715  //
716  iy_transmitter_agendaExecute( ws, iy, f_grid,
717  ppath.pos(np-1,Range(0,atmosphere_dim)),
718  ppath.los(np-1,joker), iy_transmitter_agenda );
719  if( iy.ncols() != stokes_dim || iy.nrows() != nf )
720  {
721  ostringstream os;
722  os << "The size of *iy* returned from *iy_transmitter_agenda* is\n"
723  << "not correct:\n"
724  << " expected size = [" << nf << "," << stokes_dim << "]\n"
725  << " size of iy = [" << iy.nrows() << "," << iy.ncols()<< "]\n";
726  throw runtime_error( os.str() );
727  }
728 
729 
730  //### jacobian part #########################################################
731  // Initialise analytical jacobians (diy_dx and help variables)
732  //
733  Index j_analytical_do = 0;
734  ArrayOfTensor3 diy_dpath;
735  ArrayOfIndex jac_species_i(0), jac_is_t(0), jac_wind_i(0);
736  //
737  if( jacobian_do ) { FOR_ANALYTICAL_JACOBIANS_DO( j_analytical_do = 1; ) }
738  //
739  if( !j_analytical_do )
740  { diy_dx.resize( 0 ); }
741  else
742  {
743  diy_dpath.resize( nq );
744  jac_species_i.resize( nq );
745  jac_is_t.resize( nq );
746  jac_wind_i.resize( nq );
747  //
748  FOR_ANALYTICAL_JACOBIANS_DO(
749  diy_dpath[iq].resize( np, nf, ns );
750  diy_dpath[iq] = 0.0;
751  )
752  get_pointers_for_analytical_jacobians( jac_species_i, jac_is_t,
753  jac_wind_i, jacobian_quantities, abs_species );
754  if( iy_agenda_call1 )
755  {
756  diy_dx.resize( nq );
757  //
758  FOR_ANALYTICAL_JACOBIANS_DO( diy_dx[iq].resize(
759  jacobian_indices[iq][1]-jacobian_indices[iq][0]+1, nf, ns );
760  diy_dx[iq] = 0.0;
761  )
762  }
763  }
764  //###########################################################################
765 
766 
767  // Set up variable with index of species where we want abs_per_species.
768  // This variable can below be extended in iy_aux part.
769  //
770  ArrayOfIndex iaps(0);
771  //
772  for( Index i=0; i<jac_species_i.nelem(); i++ )
773  {
774  if( jac_species_i[i] >= 0 )
775  { iaps.push_back( jac_species_i[i] ); }
776  }
777 
778 
779  //=== iy_aux part ===========================================================
780  Index auxPressure = -1,
781  auxTemperature = -1,
782  auxAbsSum = -1,
783  auxPartExt = -1,
784  auxIy = -1,
785  auxTrans = -1,
786  auxOptDepth = -1,
787  auxFarRotTotal = -1,
788  auxFarRotSpeed = -1;
789  ArrayOfIndex auxAbsSpecies(0), auxAbsIsp(0);
790  ArrayOfIndex auxVmrSpecies(0), auxVmrIsp(0);
791  ArrayOfIndex auxPartCont(0), auxPartContI(0);
792  ArrayOfIndex auxPartField(0), auxPartFieldI(0);
793  //
794  const Index naux = iy_aux_vars.nelem();
795  iy_aux.resize( naux );
796  //
797  for( Index i=0; i<naux; i++ )
798  {
799  if( iy_aux_vars[i] == "Pressure" )
800  { auxPressure = i; iy_aux[i].resize( 1, 1, 1, np ); }
801  else if( iy_aux_vars[i] == "Temperature" )
802  { auxTemperature = i; iy_aux[i].resize( 1, 1, 1, np ); }
803  else if( iy_aux_vars[i].substr(0,13) == "VMR, species " )
804  {
805  Index ispecies;
806  istringstream is(iy_aux_vars[i].substr(13,2));
807  is >> ispecies;
808  if( ispecies < 0 || ispecies>=abs_species.nelem() )
809  {
810  ostringstream os;
811  os << "You have selected VMR of species with index "
812  << ispecies << ".\nThis species does not exist!";
813  throw runtime_error( os.str() );
814  }
815  auxVmrSpecies.push_back(i);
816  auxVmrIsp.push_back(ispecies);
817  iy_aux[i].resize( 1, 1, 1, np );
818  }
819  else if( iy_aux_vars[i] == "Absorption, summed" )
820  { auxAbsSum = i; iy_aux[i].resize( nf, ns, ns, np ); }
821  else if( iy_aux_vars[i] == "Particle extinction, summed" )
822  {
823  auxPartExt = i;
824  iy_aux[i].resize( nf, ns, ns, np );
825  iy_aux[i] = 0;
826  }
827  else if( iy_aux_vars[i].substr(0,20) == "Absorption, species " )
828  {
829  Index ispecies;
830  istringstream is(iy_aux_vars[i].substr(20,2));
831  is >> ispecies;
832  if( ispecies < 0 || ispecies>=abs_species.nelem() )
833  {
834  ostringstream os;
835  os << "You have selected absorption species with index "
836  << ispecies << ".\nThis species does not exist!";
837  throw runtime_error( os.str() );
838  }
839  auxAbsSpecies.push_back(i);
840  const Index ihit = find_first( iaps, ispecies );
841  if( ihit >= 0 )
842  { auxAbsIsp.push_back( ihit ); }
843  else
844  {
845  iaps.push_back(ispecies);
846  auxAbsIsp.push_back( iaps.nelem()-1 );
847  }
848  iy_aux[i].resize( nf, ns, ns, np );
849  }
850  else if( iy_aux_vars[i].substr(0,14) == "Mass content, " )
851  {
852  Index icont;
853  istringstream is(iy_aux_vars[i].substr(14,2));
854  is >> icont;
855  if( icont < 0 || icont>=particle_masses.ncols() )
856  {
857  ostringstream os;
858  os << "You have selected particle mass content category with "
859  << "index " << icont << ".\nThis category is not defined!";
860  throw runtime_error( os.str() );
861  }
862  auxPartCont.push_back(i);
863  auxPartContI.push_back(icont);
864  iy_aux[i].resize( 1, 1, 1, np );
865  }
866  else if( iy_aux_vars[i].substr(0,10) == "PND, type " )
867  {
868  Index ip;
869  istringstream is(iy_aux_vars[i].substr(10,2));
870  is >> ip;
871  if( ip < 0 || ip>=pnd_field.nbooks() )
872  {
873  ostringstream os;
874  os << "You have selected particle number density field with "
875  << "index " << ip << ".\nThis field is not defined!";
876  throw runtime_error( os.str() );
877  }
878  auxPartField.push_back(i);
879  auxPartFieldI.push_back(ip);
880  iy_aux[i].resize( 1, 1, 1, np );
881  }
882  else if( iy_aux_vars[i] == "iy" && auxIy < 0 )
883  { auxIy = i; iy_aux[i].resize( nf, ns, 1, np ); }
884  else if( iy_aux_vars[i] == "Transmission" && auxTrans < 0 )
885  { auxTrans = i; iy_aux[i].resize( nf, ns, ns, np ); }
886  else if( iy_aux_vars[i] == "Optical depth" )
887  { auxOptDepth = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
888  else if( iy_aux_vars[i] == "Faraday rotation" )
889  { auxFarRotTotal = i; iy_aux[i].resize( nf, 1, 1, 1 ); iy_aux[i] = 0; }
890  else if( iy_aux_vars[i] == "Faraday speed" )
891  { auxFarRotSpeed = i; iy_aux[i].resize( nf, 1, 1, np ); iy_aux[i] = 0; }
892  else
893  {
894  ostringstream os;
895  os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
896  << "\"\nThis choice is not recognised.";
897  throw runtime_error( os.str() );
898  }
899  }
900 
901  // Special stuff to handle Faraday rotation
902  Index ife = -1; // When ready, ife refers abs_per_species
903  if( auxFarRotTotal>=0 || auxFarRotSpeed>=0 )
904  {
905  if( stokes_dim < 3 )
906  throw runtime_error(
907  "To include Faraday rotation, stokes_dim >= 3 is required." );
908 
909  // Determine species index of free electrons
910  for( Index sp = 0; sp < abs_species.nelem() && ife < 0; sp++ )
911  {
912  if (abs_species[sp][0].Type() == SpeciesTag::TYPE_FREE_ELECTRONS)
913  { ife = sp; }
914  }
915  // If not found, then aux values already set to zero
916  if( ife < 0 )
917  {
918  auxFarRotTotal = -1;
919  auxFarRotSpeed = -1;
920  }
921  else
922  {
923  const Index ihit = find_first( iaps, ife );
924  if( ihit >= 0 )
925  { ife = ihit; }
926  else
927  {
928  iaps.push_back(ife);
929  ife = iaps.nelem() - 1;
930  }
931  }
932  }
933  //===========================================================================
934 
935 
936  // Get atmospheric and RT quantities for each ppath point/step
937  //
938  Vector ppath_p, ppath_t;
939  Matrix ppath_vmr, ppath_pnd, ppath_wind, ppath_mag, ppath_f;
940  Tensor5 abs_per_species;
941  Tensor4 ppath_abs, trans_partial, trans_cumulat, pnd_ext_mat;
942  Vector scalar_tau;
943  ArrayOfIndex clear2cloudbox;
944  ArrayOfArrayOfIndex extmat_case;
945  //
946  if( np > 1 )
947  {
948  get_ppath_atmvars( ppath_p, ppath_t, ppath_vmr,
949  ppath_wind, ppath_mag,
950  ppath, atmosphere_dim, p_grid, t_field, vmr_field,
951  wind_u_field, wind_v_field, wind_w_field,
952  mag_u_field, mag_v_field, mag_w_field );
953  get_ppath_f( ppath_f, ppath, f_grid, atmosphere_dim,
954  rte_alonglos_v, ppath_wind );
955  get_ppath_abs( ws, ppath_abs, abs_per_species,
956  propmat_clearsky_agenda, ppath,
957  ppath_p, ppath_t, ppath_vmr, ppath_f,
958  ppath_mag, f_grid, stokes_dim, iaps );
959  if( !cloudbox_on )
960  {
961  get_ppath_trans( trans_partial, extmat_case, trans_cumulat,
962  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim );
963  }
964  else
965  {
966  Array<ArrayOfSingleScatteringData> scat_data;
967  Tensor3 pnd_abs_vec;
968  //
969  get_ppath_ext( clear2cloudbox, pnd_abs_vec, pnd_ext_mat, scat_data,
970  ppath_pnd, ppath, ppath_t, stokes_dim, ppath_f,
971  atmosphere_dim, cloudbox_limits, pnd_field,
972  use_mean_scat_data, scat_data_array, verbosity );
973  get_ppath_trans2( trans_partial, extmat_case, trans_cumulat,
974  scalar_tau, ppath, ppath_abs, f_grid, stokes_dim,
975  clear2cloudbox, pnd_ext_mat );
976  }
977  }
978 
979 
980  //=== iy_aux part ===========================================================
981  // Fill parts of iy_aux that are defined even for np=1.
982  // Radiance
983  if( auxIy >= 0 )
984  { iy_aux[auxIy](joker,joker,0,np-1) = iy; }
985  if( auxOptDepth >= 0 )
986  {
987  if( np == 1 )
988  { iy_aux[auxOptDepth] = 0; }
989  else
990  { iy_aux[auxOptDepth](joker,0,0,0) = scalar_tau; }
991  }
992  if( auxTrans >= 0 ) // Complete tensor filled!
993  {
994  if( np == 1 )
995  { for( Index iv=0; iv<nf; iv++ ) {
996  id_mat( iy_aux[auxTrans](iv,joker,joker,0) ); } }
997  else
998  { iy_aux[auxTrans] = trans_cumulat; }
999  }
1000  // Faraday rotation, total
1001  if( auxFarRotTotal >= 0 )
1002  { for( Index iv=0; iv<nf; iv++ ) {
1003  iy_aux[auxFarRotTotal](iv,0,0,0) = 0; } }
1004  //===========================================================================
1005 
1006 
1007  // Do RT calculations
1008  //
1009  if( np > 1 )
1010  {
1011  //### jacobian part #####################################################
1012  // Create container for the derivatives with respect to changes
1013  // at each ppath point. And find matching abs_species-index and
1014  // "temperature flag" (for analytical jacobians).
1015  //
1016  const Numeric dt = 0.1;
1017  const Numeric dw = 5;
1018  Tensor4 ppath_at2, ppath_awu, ppath_awv, ppath_aww;
1019  //
1020  if( j_analytical_do )
1021  {
1022  // Determine if temperature is among the analytical jac. quantities.
1023  // If yes, get emission and absorption for disturbed temperature
1024  // Same for wind, but disturb only absorption
1025  for( Index iq=0; iq<jac_is_t.nelem(); iq++ )
1026  {
1027  if( jac_is_t[iq] )
1028  {
1029  Tensor5 dummy_abs_per_species;
1030  Vector t2 = ppath_t; t2 += dt;
1031  get_ppath_abs( ws, ppath_at2, dummy_abs_per_species,
1032  propmat_clearsky_agenda, ppath, ppath_p,
1033  t2, ppath_vmr, ppath_f, ppath_mag, f_grid,
1034  stokes_dim, ArrayOfIndex(0) );
1035  }
1036  else if( jac_wind_i[iq] )
1037  {
1038  if( jac_wind_i[iq] == 1 )
1039  {
1040  Tensor5 dummy_abs_per_species;
1041  Matrix f2, w2 = ppath_wind; w2(0,joker) += dw;
1042  get_ppath_f( f2, ppath, f_grid, atmosphere_dim,
1043  rte_alonglos_v, w2 );
1044  get_ppath_abs( ws, ppath_awu, dummy_abs_per_species,
1045  propmat_clearsky_agenda, ppath, ppath_p,
1046  ppath_t, ppath_vmr, f2, ppath_mag, f_grid,
1047  stokes_dim, ArrayOfIndex(0) );
1048  }
1049  else if( jac_wind_i[iq] == 2 )
1050  {
1051  Tensor5 dummy_abs_per_species;
1052  Matrix f2, w2 = ppath_wind; w2(1,joker) += dw;
1053  get_ppath_f( f2, ppath, f_grid, atmosphere_dim,
1054  rte_alonglos_v, w2 );
1055  get_ppath_abs( ws, ppath_awv, dummy_abs_per_species,
1056  propmat_clearsky_agenda, ppath, ppath_p,
1057  ppath_t, ppath_vmr, f2, ppath_mag, f_grid,
1058  stokes_dim, ArrayOfIndex(0) );
1059  }
1060  else if( jac_wind_i[iq] == 3 )
1061  {
1062  Tensor5 dummy_abs_per_species;
1063  Matrix f2, w2 = ppath_wind; w2(2,joker) += dw;
1064  get_ppath_f( f2, ppath, f_grid, atmosphere_dim,
1065  rte_alonglos_v, w2 );
1066  get_ppath_abs( ws, ppath_aww, dummy_abs_per_species,
1067  propmat_clearsky_agenda, ppath, ppath_p,
1068  ppath_t, ppath_vmr, f2, ppath_mag, f_grid,
1069  stokes_dim, ArrayOfIndex(0) );
1070  }
1071  }
1072  }
1073  }
1074  //#######################################################################
1075 
1076  //=== iy_aux part =======================================================
1077  // iy_aux for point np-1:
1078  // Pressure
1079  if( auxPressure >= 0 )
1080  { iy_aux[auxPressure](0,0,0,np-1) = ppath_p[np-1]; }
1081  // Temperature
1082  if( auxTemperature >= 0 )
1083  { iy_aux[auxTemperature](0,0,0,np-1) = ppath_t[np-1]; }
1084  // VMR
1085  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
1086  { iy_aux[auxVmrSpecies[j]](0,0,0,np-1) = ppath_vmr(auxVmrIsp[j],np-1); }
1087  // Absorption
1088  if( auxAbsSum >= 0 )
1089  { for( Index iv=0; iv<nf; iv++ ) {
1090  for( Index is1=0; is1<ns; is1++ ){
1091  for( Index is2=0; is2<ns; is2++ ){
1092  iy_aux[auxAbsSum](iv,is1,is2,np-1) =
1093  ppath_abs(iv,is1,is2,np-1); } } } }
1094  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
1095  { for( Index iv=0; iv<nf; iv++ ) {
1096  for( Index is1=0; is1<stokes_dim; is1++ ){
1097  for( Index is2=0; is2<stokes_dim; is2++ ){
1098  iy_aux[auxAbsSpecies[j]](iv,is1,is2,np-1) =
1099  abs_per_species(auxAbsIsp[j],iv,is1,is2,np-1); } } } }
1100  // Particle properties
1101  if( cloudbox_on )
1102  {
1103  // Extinction
1104  if( auxPartExt >= 0 && clear2cloudbox[np-1] >= 0 )
1105  {
1106  const Index ic = clear2cloudbox[np-1];
1107  for( Index iv=0; iv<nf; iv++ ) {
1108  for( Index is1=0; is1<ns; is1++ ){
1109  for( Index is2=0; is2<ns; is2++ ){
1110  iy_aux[auxPartExt](iv,is1,is2,np-1) =
1111  pnd_ext_mat(iv,is1,is2,ic); } } }
1112  }
1113  // Particle mass content
1114  for( Index j=0; j<auxPartCont.nelem(); j++ )
1115  { iy_aux[auxPartCont[j]](0,0,0,np-1) = ppath_pnd(joker,np-1) *
1116  particle_masses(joker,auxPartContI[j]); }
1117  // Particle field
1118  for( Index j=0; j<auxPartField.nelem(); j++ )
1119  { iy_aux[auxPartField[j]](0,0,0,np-1) =
1120  ppath_pnd(auxPartFieldI[j],np-1); }
1121  }
1122  // Radiance for this point is handled above
1123  // Faraday speed
1124  if( auxFarRotSpeed >= 0 )
1125  { for( Index iv=0; iv<nf; iv++ ) {
1126  iy_aux[auxFarRotSpeed](iv,0,0,np-1) = 0.5 *
1127  abs_per_species(ife,iv,1,2,np-1); } }
1128  //=======================================================================
1129 
1130 
1131  // Loop ppath steps
1132  for( Index ip=np-2; ip>=0; ip-- )
1133  {
1134  //### jacobian part #################################################
1135  if( j_analytical_do )
1136  {
1137  // Loop quantities
1138  //
1139  Index iiaps = -1; // Index with respect to abs_per_species.
1140  // Jacobian species stored first and in order
1141  for( Index iq=0; iq<nq; iq++ )
1142  {
1143  if( jacobian_quantities[iq].Analytical() )
1144  {
1145  //- Absorbing species -----------------------------------
1146  if( jac_species_i[iq] >= 0 )
1147  {
1148  // Index with respect to abs_per_species and ppath_vmr
1149  iiaps += 1;
1150  const Index isp = jac_species_i[iq];
1151 
1152  // Scaling factors to handle retrieval unit
1153  Numeric unitscf1, unitscf2;
1154  vmrunitscf( unitscf1,
1155  jacobian_quantities[iq].Mode(),
1156  ppath_vmr(isp,ip), ppath_p[ip],
1157  ppath_t[ip] );
1158  vmrunitscf( unitscf2,
1159  jacobian_quantities[iq].Mode(),
1160  ppath_vmr(isp,ip+1), ppath_p[ip+1],
1161  ppath_t[ip+1] );
1162 
1163  for( Index iv=0; iv<nf; iv++ )
1164  {
1165  // Diagonal transmission matrix
1166  if( extmat_case[ip][iv] == 1 )
1167  {
1168  const Numeric x = -0.5 * ppath.lstep[ip] *
1169  trans_cumulat(iv,0,0,ip+1);
1170  for( Index is=0; is<ns; is++ )
1171  {
1172  const Numeric z = x * iy(iv,is);
1173  diy_dpath[iq](ip ,iv,is) += z * unitscf1*
1174  abs_per_species(iiaps,iv,0,0,ip );
1175  diy_dpath[iq](ip+1,iv,is) += z * unitscf2*
1176  abs_per_species(iiaps,iv,0,0,ip+1);
1177  }
1178  }
1179 
1180  // General case
1181  else
1182  {
1183  // Size of disturbance, a relative number
1184  const Numeric dd = 1e-3;
1185  // Disturb for ip
1186  Matrix ext_mat(ns,ns), dtdx(ns,ns);
1187  //
1188  for( Index is1=0; is1<ns; is1++ ) {
1189  for( Index is2=0; is2<ns; is2++ ) {
1190  ext_mat(is1,is2) = 0.5 * ( dd *
1191  abs_per_species(iiaps,iv,is1,is2,ip ) +
1192  ppath_abs(iv,is1,is2,ip) +
1193  ppath_abs(iv,is1,is2,ip+1) );
1194  } }
1195  ext2trans( dtdx, extmat_case[ip][iv],
1196  ext_mat, ppath.lstep[ip] );
1197  for( Index is1=0; is1<ns; is1++ ) {
1198  for( Index is2=0; is2<ns; is2++ ) {
1199  dtdx(is1,is2) = (unitscf1/dd) *
1200  ( dtdx(is1,is2) -
1201  trans_partial(iv,is1,is2,ip) );
1202  } }
1203  Vector x(ns), y(ns);
1204  mult( x, dtdx, iy(iv,joker) );
1205  mult( y, trans_cumulat(iv,joker,joker,ip),
1206  x );
1207  diy_dpath[iq](ip,iv,joker) += y;
1208  //
1209  // Disturb for ip+1
1210  for( Index is1=0; is1<ns; is1++ ) {
1211  for( Index is2=0; is2<ns; is2++ ) {
1212  ext_mat(is1,is2) = 0.5 * ( dd *
1213  abs_per_species(iiaps,iv,is1,is2,ip+1) +
1214  ppath_abs(iv,is1,is2,ip) +
1215  ppath_abs(iv,is1,is2,ip+1) );
1216  } }
1217  ext2trans( dtdx, extmat_case[ip][iv],
1218  ext_mat, ppath.lstep[ip] );
1219  for( Index is1=0; is1<ns; is1++ ) {
1220  for( Index is2=0; is2<ns; is2++ ) {
1221  dtdx(is1,is2) = (unitscf2/dd) *
1222  ( dtdx(is1,is2) -
1223  trans_partial(iv,is1,is2,ip) );
1224  } }
1225  mult( x, dtdx, iy(iv,joker) );
1226  mult( y, trans_cumulat(iv,joker,joker,ip),
1227  x );
1228  diy_dpath[iq](ip+1,iv,joker) += y;
1229  }
1230  }
1231  }
1232 
1233  //- Winds -----------------------------------------------
1234  else if( jac_wind_i[iq] )
1235  {
1236  for( Index iv=0; iv<nf; iv++ )
1237  {
1238  // Pick out disturbed absorption to use.
1239  // V-component is first guess.
1240  Tensor4* ppath_awx = &ppath_awv;
1241  if( jac_wind_i[iq] == 1 )
1242  { ppath_awx = &ppath_awu; }
1243  else if( jac_wind_i[iq] == 3 )
1244  { ppath_awx = &ppath_aww; }
1245 
1246  // Diagonal transmission matrix
1247  if( extmat_case[ip][iv] == 1 )
1248  {
1249  const Numeric dkdx1 = (1/dw) * (
1250  (*ppath_awx)(iv,0,0,ip ) -
1251  ppath_abs(iv,0,0,ip ) );
1252  const Numeric dkdx2 = (1/dw) * (
1253  (*ppath_awx)(iv,0,0,ip+1) -
1254  ppath_abs(iv,0,0,ip+1) );
1255  const Numeric x = -0.5 * ppath.lstep[ip] *
1256  trans_cumulat(iv,0,0,ip+1);
1257  for( Index is=0; is<ns; is++ )
1258  {
1259  const Numeric z = x * iy(iv,is);
1260  diy_dpath[iq](ip ,iv,is) += z * dkdx1;
1261  diy_dpath[iq](ip+1,iv,is) += z * dkdx2;
1262  }
1263  }
1264 
1265  // General case
1266  else
1267  {
1268  // Disturb for ip
1269  Matrix ext_mat(ns,ns), dtdx(ns,ns);
1270  //
1271  for( Index is1=0; is1<ns; is1++ ) {
1272  for( Index is2=0; is2<ns; is2++ ) {
1273  ext_mat(is1,is2) = 0.5 * (
1274  (*ppath_awx)(iv,is1,is2,ip ) +
1275  ppath_abs(iv,is1,is2,ip+1) );
1276  } }
1277  ext2trans( dtdx, extmat_case[ip][iv],
1278  ext_mat, ppath.lstep[ip] );
1279  for( Index is1=0; is1<ns; is1++ ) {
1280  for( Index is2=0; is2<ns; is2++ ) {
1281  dtdx(is1,is2) = (1/dw) *
1282  ( dtdx(is1,is2) -
1283  trans_partial(iv,is1,is2,ip) );
1284  } }
1285  Vector x(ns), y(ns);
1286  mult( x, dtdx, iy(iv,joker) );
1287  mult( y, trans_cumulat(iv,joker,joker,ip),
1288  x );
1289  diy_dpath[iq](ip,iv,joker) += y;
1290  //
1291  // Disturb for ip+1
1292  for( Index is1=0; is1<ns; is1++ ) {
1293  for( Index is2=0; is2<ns; is2++ ) {
1294  ext_mat(is1,is2) = 0.5 * (
1295  ppath_abs(iv,is1,is2,ip ) +
1296  (*ppath_awx)(iv,is1,is2,ip+1) );
1297  } }
1298  ext2trans( dtdx, extmat_case[ip][iv],
1299  ext_mat, ppath.lstep[ip] );
1300  for( Index is1=0; is1<ns; is1++ ) {
1301  for( Index is2=0; is2<ns; is2++ ) {
1302  dtdx(is1,is2) = (1/dw) *
1303  ( dtdx(is1,is2) -
1304  trans_partial(iv,is1,is2,ip) );
1305  } }
1306  mult( x, dtdx, iy(iv,joker) );
1307  mult( y, trans_cumulat(iv,joker,joker,ip),
1308  x );
1309  diy_dpath[iq](ip+1,iv,joker) += y;
1310  }
1311  }
1312  }
1313 
1314  //- Temperature -----------------------------------------
1315  else if( jac_is_t[iq] )
1316  {
1317  for( Index iv=0; iv<nf; iv++ )
1318  {
1319  // Diagonal transmission matrix
1320  if( extmat_case[ip][iv] == 1 )
1321  {
1322  const Numeric dkdt1 = 1/dt * (
1323  ppath_at2(iv,0,0,ip ) -
1324  ppath_abs(iv,0,0,ip ) );
1325  const Numeric dkdt2 = 1/dt * (
1326  ppath_at2(iv,0,0,ip+1) -
1327  ppath_abs(iv,0,0,ip+1) );
1328  const Numeric x = -0.5 * ppath.lstep[ip] *
1329  trans_cumulat(iv,0,0,ip+1);
1330  for( Index is=0; is<ns; is++ )
1331  {
1332  const Numeric z = x * iy(iv,is);
1333  diy_dpath[iq](ip ,iv,is) += z * dkdt1;
1334  diy_dpath[iq](ip+1,iv,is) += z * dkdt2;
1335  }
1336  //
1337  // The terms associated with Delta-s:
1338  if( jacobian_quantities[iq].Subtag() ==
1339  "HSE on" )
1340  {
1341  const Numeric kbar = 0.5 * (
1342  ppath_abs(iv,0,0,ip ) +
1343  ppath_abs(iv,0,0,ip+1) );
1344  for( Index is=0; is<ns; is++ )
1345  {
1346  const Numeric z = x * iy(iv,is);
1347  diy_dpath[iq](ip ,iv,is) +=
1348  z * kbar / ppath_t[ip];
1349  diy_dpath[iq](ip+1,iv,is) +=
1350  z * kbar / ppath_t[ip+1];
1351  }
1352  } //hse
1353  }
1354  // General case
1355  else
1356  {
1357  // Disturb for ip
1358  Matrix ext_mat(ns,ns), dtdx(ns,ns);
1359  //
1360  for( Index is1=0; is1<ns; is1++ ) {
1361  for( Index is2=0; is2<ns; is2++ ) {
1362  ext_mat(is1,is2) = 0.5 * (
1363  ppath_at2(iv,is1,is2,ip ) +
1364  ppath_abs(iv,is1,is2,ip+1) );
1365  } }
1366  ext2trans( dtdx, extmat_case[ip][iv],
1367  ext_mat, ppath.lstep[ip] );
1368  for( Index is1=0; is1<ns; is1++ ) {
1369  for( Index is2=0; is2<ns; is2++ ) {
1370  dtdx(is1,is2) = (1/dt) *
1371  ( dtdx(is1,is2) -
1372  trans_partial(iv,is1,is2,ip) );
1373  } }
1374  Vector x(ns), y(ns);
1375  mult( x, dtdx, iy(iv,joker) );
1376  mult( y, trans_cumulat(iv,joker,joker,ip),
1377  x );
1378  diy_dpath[iq](ip,iv,joker) += y;
1379  //
1380  // Disturb for ip+1
1381  for( Index is1=0; is1<ns; is1++ ) {
1382  for( Index is2=0; is2<ns; is2++ ) {
1383  ext_mat(is1,is2) = 0.5 * (
1384  ppath_abs(iv,is1,is2,ip ) +
1385  ppath_at2(iv,is1,is2,ip+1) );
1386  } }
1387  ext2trans( dtdx, extmat_case[ip][iv],
1388  ext_mat, ppath.lstep[ip] );
1389  for( Index is1=0; is1<ns; is1++ ) {
1390  for( Index is2=0; is2<ns; is2++ ) {
1391  dtdx(is1,is2) = (1/dt) *
1392  ( dtdx(is1,is2) -
1393  trans_partial(iv,is1,is2,ip) );
1394  } }
1395  mult( x, dtdx, iy(iv,joker) );
1396  mult( y, trans_cumulat(iv,joker,joker,ip),
1397  x );
1398  diy_dpath[iq](ip+1,iv,joker) += y;
1399  //
1400  // The terms associated with Delta-s:
1401  if( jacobian_quantities[iq].Subtag() ==
1402  "HSE on" )
1403  {
1404  for( Index is1=0; is1<ns; is1++ ) {
1405  for( Index is2=0; is2<ns; is2++ ) {
1406  ext_mat(is1,is2) = 0.5 * (
1407  ppath_abs(iv,is1,is2,ip ) +
1408  ppath_abs(iv,is1,is2,ip+1) );
1409  } }
1410  // dl for disturbed tbar
1411  const Numeric tbar = 0.5 * (
1412  ppath_t[ip] + ppath_t[ip+1] );
1413  const Numeric dl = ppath.lstep[ip] *
1414  ( 1 + dt/tbar );
1415  ext2trans( dtdx, extmat_case[ip][iv],
1416  ext_mat, dl );
1417  for( Index is1=0; is1<ns; is1++ ) {
1418  for( Index is2=0; is2<ns; is2++ ) {
1419  dtdx(is1,is2) = (1/dt) *
1420  ( dtdx(is1,is2) -
1421  trans_partial(iv,is1,is2,ip) );
1422  } }
1423  mult( x, dtdx, iy(iv,joker) );
1424  mult( y, trans_cumulat(iv,joker,joker,ip),
1425  x );
1426  // Contribution shared between the two
1427  // points and is proportional to 1/t
1428  // See also AUG.
1429  for( Index is=0; is<ns; is++ )
1430  {
1431  diy_dpath[iq](ip ,iv,is) += y[is] *
1432  0.5 * tbar / ppath_t[ip];
1433  diy_dpath[iq](ip+1,iv,is) += y[is] *
1434  0.5 * tbar / ppath_t[ip+1];
1435  }
1436  } // HSE
1437  } // General case
1438  } // Frequency loop
1439  } // Temperature
1440  } // if this analytical
1441  } // for iq
1442  } // if any analytical
1443  //###################################################################
1444 
1445 
1446  // Spectrum at end of ppath step
1447  if( stokes_dim == 1 )
1448  {
1449  for( Index iv=0; iv<nf; iv++ )
1450  { iy(iv,0) = iy(iv,0) * trans_partial(iv,0,0,ip); }
1451  }
1452  else
1453  {
1454  for( Index iv=0; iv<nf; iv++ )
1455  {
1456  // Unpolarised:
1457  if( is_diagonal( trans_partial(iv,joker,joker,ip) ) )
1458  {
1459  for( Index is=0; is<ns; is++ )
1460  { iy(iv,is) = iy(iv,is) * trans_partial(iv,is,is,ip); }
1461  }
1462  // The general case:
1463  else
1464  {
1465  Vector t1(ns);
1466  mult( t1, trans_partial(iv,joker,joker,ip), iy(iv,joker));
1467  iy(iv,joker) = t1;
1468  }
1469  }
1470  }
1471 
1472 
1473  //=== iy_aux part ===================================================
1474  // Pressure
1475  if( auxPressure >= 0 )
1476  { iy_aux[auxPressure](0,0,0,ip) = ppath_p[ip]; }
1477  // Temperature
1478  if( auxTemperature >= 0 )
1479  { iy_aux[auxTemperature](0,0,0,ip) = ppath_t[ip]; }
1480  // VMR
1481  for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
1482  { iy_aux[auxVmrSpecies[j]](0,0,0,ip) = ppath_vmr(auxVmrIsp[j],ip);}
1483  // Absorption
1484  if( auxAbsSum >= 0 )
1485  { for( Index iv=0; iv<nf; iv++ ) {
1486  for( Index is1=0; is1<ns; is1++ ){
1487  for( Index is2=0; is2<ns; is2++ ){
1488  iy_aux[auxAbsSum](iv,is1,is2,ip) =
1489  ppath_abs(iv,is1,is2,ip); } } } }
1490  for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
1491  { for( Index iv=0; iv<nf; iv++ ) {
1492  for( Index is1=0; is1<stokes_dim; is1++ ){
1493  for( Index is2=0; is2<stokes_dim; is2++ ){
1494  iy_aux[auxAbsSpecies[j]](iv,is1,is2,ip) =
1495  abs_per_species(auxAbsIsp[j],iv,is1,is2,ip); } } } }
1496  // Particle properties
1497  if( cloudbox_on )
1498  {
1499  // Extinction
1500  if( auxPartExt >= 0 && clear2cloudbox[ip] >= 0 )
1501  {
1502  const Index ic = clear2cloudbox[ip];
1503  for( Index iv=0; iv<nf; iv++ ) {
1504  for( Index is1=0; is1<ns; is1++ ){
1505  for( Index is2=0; is2<ns; is2++ ){
1506  iy_aux[auxPartExt](iv,is1,is2,ip) =
1507  pnd_ext_mat(iv,is1,is2,ic); } } }
1508  }
1509  // Particle mass content
1510  for( Index j=0; j<auxPartCont.nelem(); j++ )
1511  { iy_aux[auxPartCont[j]](0,0,0,ip) = ppath_pnd(joker,ip) *
1512  particle_masses(joker,auxPartContI[j]); }
1513  // Particle field
1514  for( Index j=0; j<auxPartField.nelem(); j++ )
1515  { iy_aux[auxPartField[j]](0,0,0,ip) =
1516  ppath_pnd(auxPartFieldI[j],ip); }
1517  }
1518  // Radiance
1519  if( auxIy >= 0 )
1520  { iy_aux[auxIy](joker,joker,0,ip) = iy; }
1521  // Faraday rotation, total
1522  if( auxFarRotTotal >= 0 )
1523  { for( Index iv=0; iv<nf; iv++ ) {
1524  iy_aux[auxFarRotTotal](iv,0,0,0) += RAD2DEG * ppath.lstep[ip] *
1525  0.25 * ( abs_per_species(ife,iv,1,2,ip ) +
1526  abs_per_species(ife,iv,1,2,ip+1)); } }
1527  // Faraday speed
1528  if( auxFarRotSpeed >= 0 )
1529  { for( Index iv=0; iv<nf; iv++ ) {
1530  iy_aux[auxFarRotSpeed](iv,0,0,ip) = 0.5 *
1531  abs_per_species(ife,iv,1,2,ip); } }
1532  //===================================================================
1533  }
1534 
1535  //### jacobian part #####################################################
1536  // Map jacobians from ppath to retrieval grids
1537  // (this operation corresponds to the term Dx_i/Dx)
1538  if( j_analytical_do )
1539  {
1540  FOR_ANALYTICAL_JACOBIANS_DO(
1541  diy_from_path_to_rgrids( diy_dx[iq], jacobian_quantities[iq],
1542  diy_dpath[iq], atmosphere_dim, ppath, ppath_p );
1543  )
1544  }
1545  //#######################################################################
1546  } // if np>1
1547 }
1548 
1549 
1550 
1551 
1552 
1553 /* Workspace method: Doxygen documentation will be auto-generated */
1554 void iy_transmitterMultiplePol(
1555  Matrix& iy,
1556  const Index& stokes_dim,
1557  const Vector& f_grid,
1558  const ArrayOfIndex& sensor_pol,
1559  const Verbosity& )
1560 {
1561  const Index nf = f_grid.nelem();
1562 
1563  if( sensor_pol.nelem() != nf )
1564  throw runtime_error( "The length of *f_grid* and the number of elements "
1565  "in *sensor_pol* must be equal." );
1566 
1567  iy.resize( nf, stokes_dim );
1568  iy = 0;
1569 
1570  ArrayOfVector s2p;
1571  stokes2pol( s2p, 1 );
1572 
1573  for( Index i=0; i<nf; i++ )
1574  {
1575  for( Index j=0; j<s2p[sensor_pol[i]-1].nelem(); j++ )
1576  {
1577  iy(i,j) = s2p[sensor_pol[i]-1][j];
1578  }
1579  }
1580 }
1581 
1582 
1583 
1584 /* Workspace method: Doxygen documentation will be auto-generated */
1585 void iy_transmitterSinglePol(
1586  Matrix& iy,
1587  const Index& stokes_dim,
1588  const Vector& f_grid,
1589  const ArrayOfIndex& sensor_pol,
1590  const Verbosity& )
1591 {
1592  const Index nf = f_grid.nelem();
1593 
1594  if( sensor_pol.nelem() != 1 )
1595  throw runtime_error( "The number of elements in *sensor_pol* must be 1." );
1596 
1597  iy.resize( nf, stokes_dim );
1598  iy = 0;
1599 
1600  ArrayOfVector s2p;
1601  stokes2pol( s2p, 1 );
1602 
1603  for( Index j=0; j<s2p[sensor_pol[0]-1].nelem(); j++ )
1604  {
1605  iy(0,j) = s2p[sensor_pol[0]-1][j];
1606  for( Index i=1; i<nf; i++ )
1607  {
1608  iy(i,j) = iy(0,j);
1609  }
1610  }
1611 }
1612 
1613 
1614 
Matrix
The Matrix class.
Definition: matpackI.h:788
defocusing_general
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
get_ppath_abs
void get_ppath_abs(Workspace &ws, Tensor4 &ppath_abs, Tensor5 &abs_per_species, const Agenda &propmat_clearsky_agenda, const Ppath &ppath, ConstVectorView ppath_p, ConstVectorView ppath_t, ConstMatrixView ppath_vmr, ConstMatrixView ppath_f, ConstMatrixView ppath_mag, ConstVectorView f_grid, const Index &stokes_dim, const ArrayOfIndex &ispecies)
get_ppath_abs
Definition: rte.cc:1351
get_ppath_trans2
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
bending_angle1d
void bending_angle1d(Numeric &alpha, const Ppath &ppath)
bending_angle1d
Definition: rte.cc:353
auto_md.h
get_pointers_for_analytical_jacobians
void get_pointers_for_analytical_jacobians(ArrayOfIndex &abs_species_i, ArrayOfIndex &is_t, ArrayOfIndex &wind_i, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:292
Tensor3
The Tensor3 class.
Definition: matpackIII.h:348
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
stokes2pol
void stokes2pol(ArrayOfVector &s2p, const Numeric &nv)
stokes2pol
Definition: sensor.cc:1196
is_diagonal
bool is_diagonal(ConstMatrixView A)
Checks if a square matrix is diagonal.
Definition: logic.cc:411
joker
const Joker joker
RAD2DEG
const Numeric RAD2DEG
PI
const Numeric PI
defocusing_sat2sat
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
vmrunitscf
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
vmrunitscf
Definition: jacobian.cc:853
Ppath::end_lstep
Numeric end_lstep
Definition: ppath.h:73
diy_from_path_to_rgrids
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Definition: jacobian.cc:110
iy_transmitterSinglePol
void iy_transmitterSinglePol(Matrix &iy, const Index &stokes_dim, const Vector &f_grid, const ArrayOfIndex &sensor_pol, const Verbosity &)
WORKSPACE METHOD: iy_transmitterSinglePol.
Definition: m_transmitter.cc:1585
Ppath::start_lstep
Numeric start_lstep
Definition: ppath.h:66
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:114
SpeciesTag::TYPE_FREE_ELECTRONS
@ TYPE_FREE_ELECTRONS
Definition: abs_species_tags.h:121
Tensor4
The Tensor4 class.
Definition: matpackIV.h:383
ppath_agendaExecute
void ppath_agendaExecute(Workspace &ws, Ppath &ppath, const Numeric ppath_lraytrace, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index cloudbox_on, const Index ppath_inside_cloudbox_do, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:14625
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
iyRadioLink
void iyRadioLink(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 Vector &lat_grid, const Vector &lon_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 Vector &refellipsoid, const Matrix &z_surface, 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 ArrayOfString &iy_aux_vars, const Index &jacobian_do, const Agenda &ppath_agenda, const Agenda &ppath_step_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &iy_transmitter_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 Index &defocus_method, const Numeric &defocus_shift, const Verbosity &verbosity)
WORKSPACE METHOD: iyRadioLink.
Definition: m_transmitter.cc:66
Agenda
The Agenda class.
Definition: agenda_class.h:44
Ppath::end_pos
Vector end_pos
Definition: ppath.h:71
iyTransmissionStandard
void iyTransmissionStandard(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 ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Agenda &ppath_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &iy_transmitter_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Numeric &rte_alonglos_v, const Numeric &ppath_lraytrace, const Verbosity &verbosity)
WORKSPACE METHOD: iyTransmissionStandard.
Definition: m_transmitter.cc:651
complex.h
A class implementing complex numbers for ARTS.
Array
This can be used to make arrays out of anything.
Definition: array.h:107
iy_transmitterMultiplePol
void iy_transmitterMultiplePol(Matrix &iy, const Index &stokes_dim, const Vector &f_grid, const ArrayOfIndex &sensor_pol, const Verbosity &)
WORKSPACE METHOD: iy_transmitterMultiplePol.
Definition: m_transmitter.cc:1554
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
sensor.h
Sensor modelling functions.
mult
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
messages.h
Declarations having to do with the four output streams.
Ppath::ngroup
Vector ngroup
Definition: ppath.h:75
pos2refell_r
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition: geodetic.cc:1139
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
Ppath::los
Matrix los
Definition: ppath.h:68
Ppath::start_pos
Vector start_pos
Definition: ppath.h:64
get_ppath_ext
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
Ppath::np
Index np
Definition: ppath.h:61
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
ArrayOfIndex
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:40
ns
#define ns
Definition: continua.cc:21931
jacobian.h
Declarations required for the calculation of jacobians.
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
lin_alg.h
Linear algebra functions.
Verbosity
Definition: messages.h:50
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
math_funcs.h
Ppath::lstep
Vector lstep
Definition: ppath.h:70
geodetic.h
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
iy_transmitter_agendaExecute
void iy_transmitter_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:14171
Tensor5
The Tensor5 class.
Definition: matpackV.h:451
Range
The range class.
Definition: matpackI.h:148
DEG2RAD
const Numeric DEG2RAD
get_ppath_atmvars
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstTensor3View wind_u_field, ConstTensor3View wind_v_field, ConstTensor3View wind_w_field, ConstTensor3View mag_u_field, ConstTensor3View mag_v_field, ConstTensor3View mag_w_field)
get_ppath_atmvars
Definition: rte.cc:1233
Ppath::constant
Numeric constant
Definition: ppath.h:62
get_ppath_trans
void get_ppath_trans(Tensor4 &trans_partial, ArrayOfArrayOfIndex &extmat_case, Tensor4 &trans_cumulat, Vector &scalar_tau, const Ppath &ppath, ConstTensor4View &ppath_abs, ConstVectorView f_grid, const Index &stokes_dim)
get_ppath_trans
Definition: rte.cc:1770
logic.h
Header file for logic.cc.
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
rte.h
Declaration of functions in rte.cc.
ext2trans
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Definition: rte.cc:903
Workspace
Workspace class.
Definition: workspace_ng.h:47
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
find_first
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:281
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
Vector
The Vector class.
Definition: matpackI.h:556
distance3D
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
Ppath::pos
Matrix pos
Definition: ppath.h:67
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:176
get_ppath_f
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, ConstVectorView f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, ConstMatrixView ppath_wind)
get_ppath_f
Definition: rte.cc:1690
distance2D
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:201
arts.h
The global header file for ARTS.