ARTS 2.5.11 (git: 725533f0)
m_transmitter.cc
Go to the documentation of this file.
1/*===========================================================================
2 === File description
3 ===========================================================================*/
4
14/*===========================================================================
15 === External declarations
16 ===========================================================================*/
17
18#include "arts.h"
19#include "arts_constants.h"
20#include "arts_conversions.h"
21#include "auto_md.h"
22#include "geodetic.h"
23#include "jacobian.h"
24#include "lin_alg.h"
25#include "logic.h"
26#include "math_funcs.h"
27#include "matpack_complex.h"
28#include "messages.h"
29#include "rte.h"
30#include "sensor.h"
31#include <cmath>
32#include <stdexcept>
33
34inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
35inline constexpr Numeric PI=Constant::pi;
36inline constexpr Numeric RAD2DEG=Conversion::rad2deg(1);
38
39/* Workspace method: Doxygen documentation will be auto-generated */
40/* Has not been updated since v2.2
41void iyRadioLink(
42 Workspace& ws,
43 Matrix& iy,
44 ArrayOfTensor4& iy_aux,
45 Ppath& ppath,
46 ArrayOfTensor3& diy_dx,
47 const Index& stokes_dim,
48 const Vector& f_grid,
49 const Index& atmosphere_dim,
50 const Vector& p_grid,
51 const Vector& lat_grid,
52 const Vector& lon_grid,
53 const Tensor3& z_field,
54 const Tensor3& t_field,
55 const Tensor4& vmr_field,
56 const ArrayOfArrayOfSpeciesTag& abs_species,
57 const Tensor3& wind_u_field,
58 const Tensor3& wind_v_field,
59 const Tensor3& wind_w_field,
60 const Tensor3& mag_u_field,
61 const Tensor3& mag_v_field,
62 const Tensor3& mag_w_field,
63 const Vector& refellipsoid,
64 const Matrix& z_surface,
65 const Index& cloudbox_on,
66 const ArrayOfIndex& cloudbox_limits,
67 const Tensor4& pnd_field,
68 const ArrayOfArrayOfSingleScatteringData& scat_data,
69 const Matrix& particle_masses,
70 const ArrayOfString& iy_aux_vars,
71 const Index& jacobian_do,
72 const Agenda& ppath_agenda,
73 const Agenda& ppath_step_agenda,
74 const Agenda& propmat_clearsky_agenda,
75 const Agenda& iy_transmitter_agenda,
76 const Index& iy_agenda_call1,
77 const Tensor3& iy_transmittance,
78 const Vector& rte_pos,
79 const Vector& rte_los,
80 const Vector& rte_pos2,
81 const Numeric& rte_alonglos_v,
82 const Numeric& ppath_lmax,
83 const Numeric& ppath_lraytrace,
84 const Index& defocus_method,
85 const Numeric& defocus_shift,
86 const Verbosity& verbosity )
87{
88 // Throw error if unsupported features are requested
89 if( !iy_agenda_call1 )
90 throw runtime_error(
91 "Recursive usage not possible (iy_agenda_call1 must be 1)" );
92 if( !iy_transmittance.empty() )
93 throw runtime_error( "*iy_transmittance* must be empty" );
94 if( jacobian_do )
95 throw runtime_error( "This method does not provide any jacobians and "
96 "*jacobian_do* must be 0." );
97 if( defocus_method < 1 || defocus_method > 2 )
98 throw runtime_error( "Allowed choices for *defocus_method* is 1 and 2." );
99 diy_dx.resize(0);
100
101
102 //- Determine propagation path
103 ppath_agendaExecute( ws, ppath, ppath_lmax, ppath_lraytrace, rte_pos, rte_los,
104 rte_pos2, cloudbox_on, 0, t_field, z_field, vmr_field,
105 f_grid, ppath_agenda );
106
107 //- Set np to zero if ground intersection
108 const Index radback = ppath_what_background(ppath);
109 // np should already be 1 fon non-OK cases, but for extra safety ...
110 if( radback == 0 || radback == 2 )
111 { ppath.np = 1; }
112
113 // Some basic sizes
114 //
115 const Index nf = f_grid.nelem();
116 const Index ns = stokes_dim;
117 const Index np = ppath.np;
118
119
120 // The below copied from iyEmission. Activate for-loop when jacobians
121 // introduced.
122
123 // Set up variable with index of species where we want abs_per_species.
124 // This variable can below be extended in iy_aux part.
125 //
126 ArrayOfIndex iaps(0);
127 //
128 //for( Index i=0; i<jac_species_i.nelem(); i++ )
129 // {
130 // if( jac_species_i[i] >= 0 )
131 // { iaps.push_back( jac_species_i[i] ); }
132 // }
133
134 //=== iy_aux part ===========================================================
135 Index auxPressure = -1,
136 auxTemperature = -1,
137 auxAbsSum = -1,
138 auxPartExt = -1,
139 auxImpactParam = -1,
140 auxFreeSpaceLoss = -1,
141 auxFreeSpaceAtte = -1,
142 auxAtmosphericLoss = -1,
143 auxDefocusingLoss = -1,
144 auxFarRotTotal = -1,
145 auxFarRotSpeed = -1,
146 auxExtraPathDelay = -1,
147 auxBendingAngle = -1;
148 ArrayOfIndex auxAbsSpecies(0), auxAbsIsp(0);
149 ArrayOfIndex auxVmrSpecies(0), auxVmrIsp(0);
150 ArrayOfIndex auxPartCont(0), auxPartContI(0);
151 ArrayOfIndex auxPartField(0), auxPartFieldI(0);
152 //
153 const Index naux = iy_aux_vars.nelem();
154 iy_aux.resize( naux );
155 //
156 for( Index i=0; i<naux; i++ )
157 {
158 if( iy_aux_vars[i] == "Pressure" )
159 { auxPressure = i; iy_aux[i].resize( 1, 1, 1, np ); }
160 else if( iy_aux_vars[i] == "Temperature" )
161 { auxTemperature = i; iy_aux[i].resize( 1, 1, 1, np ); }
162 else if( iy_aux_vars[i].substr(0,13) == "VMR, species " )
163 {
164 Index ispecies;
165 istringstream is(iy_aux_vars[i].substr(13,2));
166 is >> ispecies;
167 if( ispecies < 0 || ispecies>=abs_species.nelem() )
168 {
169 ostringstream os;
170 os << "You have selected VMR of species with index "
171 << ispecies << ".\nThis species does not exist!";
172 throw runtime_error( os.str() );
173 }
174 auxVmrSpecies.push_back(i);
175 auxVmrIsp.push_back(ispecies);
176 iy_aux[i].resize( 1, 1, 1, np );
177 }
178 else if( iy_aux_vars[i] == "Absorption, summed" )
179 { auxAbsSum = i; iy_aux[i].resize( nf, ns, ns, np ); }
180 else if( iy_aux_vars[i] == "Particle extinction, summed" )
181 {
182 auxPartExt = i;
183 iy_aux[i].resize( nf, ns, ns, np );
184 iy_aux[i] = 0;
185 }
186 else if( iy_aux_vars[i].substr(0,20) == "Absorption, species " )
187 {
188 Index ispecies;
189 istringstream is(iy_aux_vars[i].substr(20,2));
190 is >> ispecies;
191 if( ispecies < 0 || ispecies>=abs_species.nelem() )
192 {
193 ostringstream os;
194 os << "You have selected absorption species with index "
195 << ispecies << ".\nThis species does not exist!";
196 throw runtime_error( os.str() );
197 }
198 auxAbsSpecies.push_back(i);
199 const Index ihit = find_first( iaps, ispecies );
200 if( ihit >= 0 )
201 { auxAbsIsp.push_back( ihit ); }
202 else
203 {
204 iaps.push_back(ispecies);
205 auxAbsIsp.push_back( iaps.nelem()-1 );
206 }
207 iy_aux[i].resize( nf, ns, ns, np );
208 }
209 else if( iy_aux_vars[i].substr(0,14) == "Mass content, " )
210 {
211 Index icont;
212 istringstream is(iy_aux_vars[i].substr(14,2));
213 is >> icont;
214 if( icont < 0 || icont>=particle_masses.ncols() )
215 {
216 ostringstream os;
217 os << "You have selected particle mass content category with "
218 << "index " << icont << ".\nThis category is not defined!";
219 throw runtime_error( os.str() );
220 }
221 auxPartCont.push_back(i);
222 auxPartContI.push_back(icont);
223 iy_aux[i].resize( 1, 1, 1, np );
224 }
225 else if( iy_aux_vars[i].substr(0,10) == "PND, type " )
226 {
227 Index ip;
228 istringstream is(iy_aux_vars[i].substr(10,2));
229 is >> ip;
230 if( ip < 0 || ip>=pnd_field.nbooks() )
231 {
232 ostringstream os;
233 os << "You have selected particle number density field with "
234 << "index " << ip << ".\nThis field is not defined!";
235 throw runtime_error( os.str() );
236 }
237 auxPartField.push_back(i);
238 auxPartFieldI.push_back(ip);
239 iy_aux[i].resize( 1, 1, 1, np );
240 }
241 else if( iy_aux_vars[i] == "Impact parameter" )
242 { auxImpactParam = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
243 else if( iy_aux_vars[i] == "Free space loss" )
244 { auxFreeSpaceLoss = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
245 else if( iy_aux_vars[i] == "Free space attenuation" )
246 { auxFreeSpaceAtte = i; iy_aux[i].resize( 1, 1, 1, np ); }
247 else if( iy_aux_vars[i] == "Atmospheric loss" )
248 { auxAtmosphericLoss = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
249 else if( iy_aux_vars[i] == "Defocusing loss" )
250 { auxDefocusingLoss = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
251 else if( iy_aux_vars[i] == "Faraday rotation" )
252 { auxFarRotTotal = i; iy_aux[i].resize( nf, 1, 1, 1 ); iy_aux[i] = 0; }
253 else if( iy_aux_vars[i] == "Faraday speed" )
254 { auxFarRotSpeed = i; iy_aux[i].resize( nf, 1, 1, np ); iy_aux[i] = 0; }
255 else if( iy_aux_vars[i] == "Extra path delay" )
256 { auxExtraPathDelay = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
257 else if( iy_aux_vars[i] == "Bending angle" )
258 { auxBendingAngle = i; iy_aux[i].resize( 1, 1, 1, 1 ); }
259 else
260 {
261 ostringstream os;
262 os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
263 << "\"\nThis choice is not recognised.";
264 throw runtime_error( os.str() );
265 }
266 }
267
268 // Special stuff to handle Faraday rotation
269 Index ife = -1; // When ready, ife refers abs_per_species
270 if( auxFarRotTotal>=0 || auxFarRotSpeed>=0 )
271 {
272 if( stokes_dim < 3 )
273 throw runtime_error(
274 "To include Faraday rotation, stokes_dim >= 3 is required." );
275
276 // Determine species index of free electrons
277 for( Index sp = 0; sp < abs_species.nelem() && ife < 0; sp++ )
278 {
279 if (abs_species[sp][0].Type() == SpeciesTag::TYPE_FREE_ELECTRONS)
280 { ife = sp; }
281 }
282 // If not found, then aux values already set to zero
283 if( ife < 0 )
284 {
285 auxFarRotTotal = -1;
286 auxFarRotSpeed = -1;
287 }
288 else
289 {
290 const Index ihit = find_first( iaps, ife );
291 if( ihit >= 0 )
292 { ife = ihit; }
293 else
294 {
295 iaps.push_back(ife);
296 ife = iaps.nelem() - 1;
297 }
298 }
299 }
300 //===========================================================================
301
302
303 // Handle cases whn no link was establsihed
304 if( radback == 0 || radback == 2 )
305 {
306 Numeric fillvalue = 0;
307 if( radback == 0 )
308 { fillvalue = NAN; }
309 //
310 iy.resize( nf, stokes_dim );
311 iy = fillvalue;
312 //
313 for( Index i=0; i<naux; i++ )
314 { iy_aux[i] = fillvalue; }
315 //
316 return;
317 }
318
319
320 // Transmitted signal
321 //
322 iy_transmitter_agendaExecute( ws, iy, f_grid,
323 ppath.pos(np-1,Range(0,atmosphere_dim)),
324 ppath.los(np-1,joker), iy_transmitter_agenda );
325 if( iy.ncols() != stokes_dim || iy.nrows() != nf )
326 { throw runtime_error( "The size of *iy* returned from "
327 "*iy_transmitter_agenda* is not correct." ); }
328
329
330 // Get atmospheric and attenuation quantities for each ppath point/step
331 //
332 Vector ppath_p, ppath_t;
333 Matrix ppath_vmr, ppath_pnd, ppath_mag, ppath_wind, ppath_f, ppath_t_nlte;
334 ArrayOfArrayOfPropagationMatrix abs_per_species, dummy_dppath_ext_dx;
335 ArrayOfPropagationMatrix ppath_ext, pnd_ext_mat;
336 Tensor4 trans_partial, trans_cumulat;
337 ArrayOfArrayOfStokesVector dummy_dppath_nlte_dx;
338 ArrayOfStokesVector dummy_ppath_nlte_source;
339 Vector scalar_tau;
340 ArrayOfIndex clear2cloudy, dummy_lte;
341 ArrayOfMatrix dummy_ppath_dpnd_dx;
342 ArrayOfTensor4 dummy_dpnd_field_dx;
343 const Tensor4 t_nlte_field_empty(0,0,0,0);
344 //
345 if( np > 1 )
346 {
347 get_ppath_atmvars( ppath_p, ppath_t, ppath_t_nlte, ppath_vmr,
348 ppath_wind, ppath_mag,
349 ppath, atmosphere_dim, p_grid, t_field, t_nlte_field_empty,
350 vmr_field, wind_u_field, wind_v_field, wind_w_field,
351 mag_u_field, mag_v_field, mag_w_field );
352
353 get_ppath_f( ppath_f, ppath, f_grid, atmosphere_dim,
354 rte_alonglos_v, ppath_wind );
355
356 get_ppath_pmat( ws, ppath_ext, dummy_ppath_nlte_source, dummy_lte, abs_per_species,
357 dummy_dppath_ext_dx, dummy_dppath_nlte_dx,
358 propmat_clearsky_agenda, ArrayOfRetrievalQuantity(0), ppath,
359 ppath_p, ppath_t, ppath_t_nlte, ppath_vmr, ppath_f,
360 ppath_mag, f_grid, stokes_dim, iaps );
361
362 if( !cloudbox_on )
363 {
364 ArrayOfArrayOfIndex extmat_case;
365 get_ppath_trans( trans_partial, extmat_case, trans_cumulat,
366 scalar_tau, ppath, ppath_ext, f_grid, stokes_dim );
367 }
368 else
369 {
370 // Extract basic scattering data
371 Array<ArrayOfArrayOfSingleScatteringData> scat_data_single;
372 ArrayOfArrayOfIndex extmat_case;
373 Tensor3 pnd_abs_vec;
374
375 get_ppath_cloudvars( clear2cloudy, ppath_pnd, dummy_ppath_dpnd_dx,
376 ppath, atmosphere_dim, cloudbox_limits,
377 pnd_field, dummy_dpnd_field_dx );
378 get_ppath_partopt( pnd_abs_vec, pnd_ext_mat, scat_data_single,
379 clear2cloudy, ppath_pnd, ppath, ppath_t,
380 stokes_dim, ppath_f, atmosphere_dim, scat_data,
381 verbosity );
382
383 get_ppath_trans2( trans_partial, extmat_case, trans_cumulat,
384 scalar_tau, ppath, ppath_ext, f_grid, stokes_dim,
385 clear2cloudy, pnd_ext_mat );
386 }
387 }
388
389 // Ppath length variables
390 //
391 Numeric lbg; // Bent geometrical length of ray path
392 Numeric lba; // Bent apparent length of ray path
393 //
394 lbg = ppath.end_lstep;
395 lba = lbg;
396
397 // Do RT calculations
398 //
399 if( np > 1 )
400 {
401 //=== iy_aux part =======================================================
402 // iy_aux for point np-1:
403 // Pressure
404 if( auxPressure >= 0 )
405 { iy_aux[auxPressure](0,0,0,np-1) = ppath_p[np-1]; }
406 // Temperature
407 if( auxTemperature >= 0 )
408 { iy_aux[auxTemperature](0,0,0,np-1) = ppath_t[np-1]; }
409 // VMR
410 for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
411 { iy_aux[auxVmrSpecies[j]](0,0,0,np-1) = ppath_vmr(auxVmrIsp[j],np-1);}
412 // Absorption
413 if( auxAbsSum >= 0 )
414 { for( Index iv=0; iv<nf; iv++ ) {
415 for( Index is1=0; is1<ns; is1++ ){
416 for( Index is2=0; is2<ns; is2++ ){
417 iy_aux[auxAbsSum](iv,is1,is2,np-1) =
418 ppath_ext[np-1](iv,is1,is2); } } } }
419 for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
420 { for( Index iv=0; iv<nf; iv++ ) {
421 for( Index is1=0; is1<stokes_dim; is1++ ){
422 for( Index is2=0; is2<stokes_dim; is2++ ){
423 iy_aux[auxAbsSpecies[j]](iv,is1,is2,np-1) =
424 abs_per_species[np-1][auxAbsIsp[j]](iv,is1,is2); } } } }
425 // Particle properties
426 if( cloudbox_on )
427 {
428 // Extinction
429 if( auxPartExt >= 0 && clear2cloudy[np-1] >= 0 )
430 {
431 const Index ic = clear2cloudy[np-1];
432 for( Index iv=0; iv<nf; iv++ ) {
433 for( Index is1=0; is1<ns; is1++ ){
434 for( Index is2=0; is2<ns; is2++ ){
435 iy_aux[auxPartExt](iv,is1,is2,np-1) =
436 pnd_ext_mat[ic](iv,is1,is2); } } }
437 }
438 // Particle mass content
439 for( Index j=0; j<auxPartCont.nelem(); j++ )
440 { iy_aux[auxPartCont[j]](0,0,0,np-1) = ppath_pnd(joker,np-1) *
441 particle_masses(joker,auxPartContI[j]); }
442 // Particle number density
443 for( Index j=0; j<auxPartField.nelem(); j++ )
444 { iy_aux[auxPartField[j]](0,0,0,np-1) =
445 ppath_pnd(auxPartFieldI[j],np-1); }
446 }
447 // Free space
448 if( auxFreeSpaceAtte >= 0 )
449 { iy_aux[auxFreeSpaceAtte](joker,0,0,np-1) = 2/lbg; }
450 // Faraday speed
451 if( auxFarRotSpeed >= 0 )
452 { for( Index iv=0; iv<nf; iv++ ) {
453 iy_aux[auxFarRotSpeed](iv,0,0,np-1) = 0.5 *
454 abs_per_species[np-1][ife](iv,1,2); } }
455 //=======================================================================
456
457 // Loop ppath steps
458 for( Index ip=np-2; ip>=0; ip-- )
459 {
460 // Lengths
461 lbg += ppath.lstep[ip];
462 lba += ppath.lstep[ip] * (ppath.ngroup[ip]+ppath.ngroup[ip+1]) / 2.0;
463
464 // Atmospheric loss of path step + Faraday rotation
465 if( stokes_dim == 1 )
466 {
467 for( Index iv=0; iv<nf; iv++ )
468 { iy(iv,0) = iy(iv,0) * trans_partial(iv,0,0,ip); }
469 }
470 else
471 {
472 for( Index iv=0; iv<nf; iv++ )
473 {
474 // Unpolarised:
475 if( is_diagonal( trans_partial(iv,joker,joker,ip) ) )
476 {
477 for( Index is=0; is<ns; is++ )
478 { iy(iv,is) = iy(iv,is) * trans_partial(iv,is,is,ip); }
479 }
480 // The general case:
481 else
482 {
483 Vector t1(ns);
484 mult( t1, trans_partial(iv,joker,joker,ip), iy(iv,joker));
485 iy(iv,joker) = t1;
486 }
487 }
488 }
489
490 //=== iy_aux part ===================================================
491 // Pressure
492 if( auxPressure >= 0 )
493 { iy_aux[auxPressure](0,0,0,ip) = ppath_p[ip]; }
494 // Temperature
495 if( auxTemperature >= 0 )
496 { iy_aux[auxTemperature](0,0,0,ip) = ppath_t[ip]; }
497 // VMR
498 for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
499 { iy_aux[auxVmrSpecies[j]](0,0,0,ip) = ppath_vmr(auxVmrIsp[j],ip);}
500 // Absorption
501 if( auxAbsSum >= 0 )
502 { for( Index iv=0; iv<nf; iv++ ) {
503 for( Index is1=0; is1<ns; is1++ ){
504 for( Index is2=0; is2<ns; is2++ ){
505 iy_aux[auxAbsSum](iv,is1,is2,ip) =
506 ppath_ext[ip](iv,is1,is2); } } } }
507 for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
508 { for( Index iv=0; iv<nf; iv++ ) {
509 for( Index is1=0; is1<stokes_dim; is1++ ){
510 for( Index is2=0; is2<stokes_dim; is2++ ){
511 iy_aux[auxAbsSpecies[j]](iv,is1,is2,ip) =
512 abs_per_species[ip][auxAbsIsp[j]](iv,is1,is2); } } } }
513 // Particle properties
514 if( cloudbox_on )
515 {
516 // Extinction
517 if( auxPartExt >= 0 && clear2cloudy[ip] >= 0 )
518 {
519 const Index ic = clear2cloudy[ip];
520 for( Index iv=0; iv<nf; iv++ ) {
521 for( Index is1=0; is1<ns; is1++ ){
522 for( Index is2=0; is2<ns; is2++ ){
523 iy_aux[auxPartExt](iv,is1,is2,ip) =
524 pnd_ext_mat[ic](iv,is1,is2); } } }
525 }
526 // Particle mass content
527 for( Index j=0; j<auxPartCont.nelem(); j++ )
528 { iy_aux[auxPartCont[j]](0,0,0,ip) = ppath_pnd(joker,ip) *
529 particle_masses(joker,auxPartContI[j]); }
530 // Particle number density
531 for( Index j=0; j<auxPartField.nelem(); j++ )
532 { iy_aux[auxPartField[j]](0,0,0,ip) =
533 ppath_pnd(auxPartFieldI[j],ip); }
534 }
535 // Free space loss
536 if( auxFreeSpaceAtte >= 0 )
537 { iy_aux[auxFreeSpaceAtte](joker,0,0,ip) = 2/lbg; }
538 // Faraday rotation, total
539 if( auxFarRotTotal >= 0 )
540 { for( Index iv=0; iv<nf; iv++ ) {
541 iy_aux[auxFarRotTotal](iv,0,0,0) += RAD2DEG * ppath.lstep[ip] *
542 0.25 * ( abs_per_species[ip][ife](iv,1,2)+
543 abs_per_species[ip+1][ife](iv,1,2)); } }
544 // Faraday speed
545 if( auxFarRotSpeed >= 0 )
546 { for( Index iv=0; iv<nf; iv++ ) {
547 iy_aux[auxFarRotSpeed](iv,0,0,ip) = 0.5 *
548 abs_per_species[ip][ife](iv,1,2); } }
549 //===================================================================
550 }
551
552
553 //=== iy_aux part =======================================================
554 if( auxAtmosphericLoss >= 0 )
555 { iy_aux[auxAtmosphericLoss](joker,0,0,0) = iy(joker,0); }
556 if( auxImpactParam >= 0 )
557 {
558 ARTS_ASSERT( ppath.constant >= 0 );
559 iy_aux[auxImpactParam](joker,0,0,0) = ppath.constant;
560 }
561 //=======================================================================
562
563
564 // Remaing length of ppath
565 lbg += ppath.start_lstep;
566 lba += ppath.start_lstep;
567
568
569 // Determine total free space loss
570 Numeric fspl = 1 / ( 4 * PI * lbg*lbg );
571 //
572 if( auxFreeSpaceLoss >= 0 )
573 { iy_aux[auxFreeSpaceLoss] = fspl; }
574
575
576 // Determine defocusing loss
577 Numeric dfl = 1;
578 if( defocus_method == 1 )
579 {
580 defocusing_general( ws, dfl, ppath_step_agenda, atmosphere_dim,
581 p_grid, lat_grid, lon_grid, t_field, z_field,
582 vmr_field, f_grid, refellipsoid,
583 z_surface, ppath, ppath_lmax, ppath_lraytrace,
584 defocus_shift, verbosity );
585 }
586 else if( defocus_method == 2 )
587 { defocusing_sat2sat( ws, dfl, ppath_step_agenda, atmosphere_dim,
588 p_grid, lat_grid, lon_grid, t_field, z_field,
589 vmr_field, f_grid, refellipsoid,
590 z_surface, ppath, ppath_lmax, ppath_lraytrace,
591 defocus_shift, verbosity );
592 }
593 if( auxDefocusingLoss >= 0 )
594 { iy_aux[auxDefocusingLoss] = dfl; }
595
596
597
598 // Include free space and defocusing losses
599 iy *= fspl*dfl;
600
601
602 // Extra path delay
603 if( auxExtraPathDelay >= 0 )
604 {
605 // Radius of rte_pos and rte_pos2
606 const Numeric r1 = ppath.end_pos[0] +
607 pos2refell_r( atmosphere_dim, refellipsoid,
608 lat_grid, lon_grid, ppath.end_pos );
609 const Numeric r2 = ppath.start_pos[0] +
610 pos2refell_r( atmosphere_dim, refellipsoid,
611 lat_grid, lon_grid, ppath.start_pos );
612
613 // Geometrical distance between start and end point
614 Numeric lgd ;
615 if( atmosphere_dim <= 2 )
616 { distance2D( lgd, r1, ppath.end_pos[1], r2, ppath.start_pos[1] ); }
617 else
618 { distance3D( lgd, r1, ppath.end_pos[1], ppath.end_pos[2],
619 r2, ppath.start_pos[1], ppath.start_pos[2] ); }
620 //
621 iy_aux[auxExtraPathDelay] = ( lba - lgd ) / SPEED_OF_LIGHT;
622 }
623
624
625 // Bending angle
626 if( auxBendingAngle >= 0 )
627 {
628 Numeric ba = -999;
629 bending_angle1d( ba, ppath );
630 //
631 iy_aux[auxBendingAngle] = ba;
632 }
633 }
634}
635*/
636
637/* Workspace method: Doxygen documentation will be auto-generated */
639 Matrix& iy,
640 ArrayOfMatrix& iy_aux,
641 ArrayOfTensor3& diy_dx,
642 Vector& ppvar_p,
643 Vector& ppvar_t,
644 EnergyLevelMap& ppvar_nlte,
645 Matrix& ppvar_vmr,
646 Matrix& ppvar_wind,
647 Matrix& ppvar_mag,
648 Matrix& ppvar_pnd,
649 Matrix& ppvar_f,
650 Tensor3& ppvar_iy,
651 Tensor4& ppvar_trans_cumulat,
652 Tensor4& ppvar_trans_partial,
653 const Index& stokes_dim,
654 const Vector& f_grid,
655 const Index& atmosphere_dim,
656 const Vector& p_grid,
657 const Tensor3& t_field,
658 const EnergyLevelMap& nlte_field,
659 const Tensor4& vmr_field,
660 const ArrayOfArrayOfSpeciesTag& abs_species,
661 const Tensor3& wind_u_field,
662 const Tensor3& wind_v_field,
663 const Tensor3& wind_w_field,
664 const Tensor3& mag_u_field,
665 const Tensor3& mag_v_field,
666 const Tensor3& mag_w_field,
667 const Index& cloudbox_on,
668 const ArrayOfIndex& cloudbox_limits,
669 const Index& gas_scattering_do,
670 const Tensor4& pnd_field,
671 const ArrayOfTensor4& dpnd_field_dx,
672 const ArrayOfString& scat_species,
673 const ArrayOfArrayOfSingleScatteringData& scat_data,
674 const ArrayOfString& iy_aux_vars,
675 const Index& jacobian_do,
676 const ArrayOfRetrievalQuantity& jacobian_quantities,
677 const Ppath& ppath,
678 const Matrix& iy_transmitter,
679 const Agenda& propmat_clearsky_agenda,
680 const Agenda& water_p_eq_agenda,
681 const Agenda& gas_scattering_agenda,
682 const Index& iy_agenda_call1,
683 const Tensor3& iy_transmittance,
684 const Numeric& rte_alonglos_v,
685 const Verbosity&) {
686 // Init Jacobian quantities?
687 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<1>(jacobian_quantities) : 0;
688
689 // Some basic sizes
690 const Index nf = f_grid.nelem();
691 const Index ns = stokes_dim;
692 const Index np = ppath.np;
693 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
694
695 // Radiative background index
696 const Index rbi = ppath_what_background(ppath);
697
698 // Checks of input
699 // Throw error if unsupported features are requested
700 if (!iy_agenda_call1)
701 throw runtime_error(
702 "Recursive usage not possible (iy_agenda_call1 must be 1)");
703 if (!iy_transmittance.empty())
704 throw runtime_error("*iy_transmittance* must be empty");
705 if (rbi < 1 || rbi > 9)
706 throw runtime_error(
707 "ppath.background is invalid. Check your "
708 "calculation of *ppath*?");
709 if (jacobian_do) {
710 if (dpnd_field_dx.nelem() != jacobian_quantities.nelem())
711 throw runtime_error(
712 "*dpnd_field_dx* not properly initialized:\n"
713 "Number of elements in dpnd_field_dx must be equal number of jacobian"
714 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
715 " is calculated/set.");
716 }
717
718 ARTS_USER_ERROR_IF(jacobian_quantities.nelem() && gas_scattering_do, R"--(
719Jacobian calculation are not supported when gas scattering or suns are included.
720This feature will be added in a future version.
721)--");
722
723 // iy_aux_vars checked below
724
725 // Transmitted signal
726// iy_transmitter=iy_transmitter;
727 if (iy_transmitter.ncols() != ns || iy_transmitter.nrows() != nf) {
728 ostringstream os;
729 os << "The size of *iy_transmitter* returned from *iy_transmitter_agenda* is\n"
730 << "not correct:\n"
731 << " expected size = [" << nf << "," << stokes_dim << "]\n"
732 << " size of iy_transmitter = [" << iy_transmitter.nrows() << "," << iy_transmitter.ncols() << "]\n";
733 throw runtime_error(os.str());
734 }
735
736 // Set diy_dpath if we are doing are doing jacobian calculations
737 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
738
739 // Set the species pointers if we are doing jacobian
740 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
741
742 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
743 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
744
745 // Checks that the scattering species are treated correctly if their derivatives are needed (we can here discard the Array)
746 if (j_analytical_do and iy_agenda_call1) get_pointers_for_scat_species(jacobian_quantities, scat_species, cloudbox_on);
747
748 // Init iy_aux and fill where possible
749 const Index naux = iy_aux_vars.nelem();
750 iy_aux.resize(naux);
751 //
752 Index auxOptDepth = -1;
753 //
754 for (Index i = 0; i < naux; i++) {
755 iy_aux[i].resize(nf, ns);
756 iy_aux[i] = 0;
757
758 if (iy_aux_vars[i] == "Radiative background")
759 iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
760 else if (iy_aux_vars[i] == "Optical depth")
761 auxOptDepth = i;
762 else {
763 ostringstream os;
764 os << "The only allowed strings in *iy_aux_vars* are:\n"
765 << " \"Radiative background\"\n"
766 << " \"Optical depth\"\n"
767 << "but you have selected: \"" << iy_aux_vars[i] << "\"";
768 throw runtime_error(os.str());
769 }
770 }
771
772 // Get atmospheric and radiative variables along the propagation path
773 //
774 ppvar_trans_cumulat.resize(np, nf, ns, ns);
775 ppvar_trans_partial.resize(np, nf, ns, ns);
776
777 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
779 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
780
786
787 ArrayOfIndex clear2cloudy;
788 //
789 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
790 ppvar_p.resize(0);
791 ppvar_t.resize(0);
792 ppvar_vmr.resize(0, 0);
793 ppvar_wind.resize(0, 0);
794 ppvar_mag.resize(0, 0);
795 ppvar_pnd.resize(0, 0);
796 ppvar_f.resize(0, 0);
797 ppvar_iy.resize(0, 0, 0);
798 ppvar_trans_cumulat = 0;
799 ppvar_trans_partial = 0;
800 for (Index iv = 0; iv < nf; iv++) {
801 for (Index is = 0; is < ns; is++) {
802 ppvar_trans_cumulat(0,iv,is,is) = 1;
803 ppvar_trans_partial(0,iv,is,is) = 1;
804 }
805 }
806
807 } else {
808 // ppvar_iy
809 ppvar_iy.resize(nf, ns, np);
810 ppvar_iy(joker, joker, np - 1) = iy_transmitter;
811
812 // Basic atmospheric variables
813 get_ppath_atmvars(ppvar_p,
814 ppvar_t,
815 ppvar_nlte,
816 ppvar_vmr,
817 ppvar_wind,
818 ppvar_mag,
819 ppath,
820 atmosphere_dim,
821 p_grid,
822 t_field,
823 nlte_field,
824 vmr_field,
825 wind_u_field,
826 wind_v_field,
827 wind_w_field,
828 mag_u_field,
829 mag_v_field,
830 mag_w_field);
831
833 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
834
835 // pnd_field
836 ArrayOfMatrix ppvar_dpnd_dx(0);
837 //
838 if (cloudbox_on)
839 get_ppath_cloudvars(clear2cloudy,
840 ppvar_pnd,
841 ppvar_dpnd_dx,
842 ppath,
843 atmosphere_dim,
844 cloudbox_limits,
845 pnd_field,
846 dpnd_field_dx);
847 else {
848 clear2cloudy.resize(np);
849 for (Index ip = 0; ip < np; ip++) clear2cloudy[ip] = -1;
850 }
851
852
853
854
855 // Size radiative variables always used
856 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
857 StokesVector a(nf, ns), S(nf, ns), Sp(nf, ns);
858 ArrayOfIndex lte(np);
859
860 // size gas scattering variables
861 TransmissionMatrix gas_scattering_mat;
862 Vector sca_fct_dummy;
863 PropagationMatrix K_sca;
864 if (gas_scattering_do) {
865 K_sca = PropagationMatrix(nf, ns);
866 }
867
868 Vector gas_scattering_los_in, gas_scattering_los_out;
869
870 // Init variables only used if analytical jacobians done
871 Vector dB_dT(0);
872 ArrayOfPropagationMatrix dK_this_dx(nq), dK_past_dx(nq), dKp_dx(nq);
873 ArrayOfStokesVector da_dx(nq), dS_dx(nq), dSp_dx(nq);
874 //
875 Index temperature_derivative_position = -1;
876 bool do_hse = false;
877
878 if (j_analytical_do) {
879 dB_dT.resize(nf);
880 FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] = PropagationMatrix(nf, ns);
881 dK_past_dx[iq] = PropagationMatrix(nf, ns);
882 dKp_dx[iq] = PropagationMatrix(nf, ns);
883 da_dx[iq] = StokesVector(nf, ns);
884 dS_dx[iq] = StokesVector(nf, ns);
885 dSp_dx[iq] = StokesVector(nf, ns);
886 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
887 temperature_derivative_position = iq;
888 do_hse = jacobian_quantities[iq].Subtag() ==
889 "HSE on";
890 })
891 }
892
893 // Loop ppath points and determine radiative properties
894 for (Index ip = 0; ip < np; ip++) {
896 K_this,
897 S,
898 lte[ip],
899 dK_this_dx,
900 dS_dx,
901 propmat_clearsky_agenda,
902 jacobian_quantities,
903 Vector{ppvar_f(joker, ip)},
904 Vector{ppvar_mag(joker, ip)},
905 Vector{ppath.los(ip, joker)},
906 ppvar_nlte[ip],
907 Vector{ppvar_vmr(joker, ip)},
908 ppvar_t[ip],
909 ppvar_p[ip],
910 j_analytical_do);
911
912 // get Extinction from gas scattering
913 if (gas_scattering_do) {
915 K_sca,
916 gas_scattering_mat,
917 sca_fct_dummy,
918 f_grid,
919 ppvar_p[ip],
920 ppvar_t[ip],
921 Vector{ppvar_vmr(joker, ip)},
922 gas_scattering_los_in,
923 gas_scattering_los_out,
924 0,
925 gas_scattering_agenda);
926
927 K_this += K_sca;
928 }
929
930 if (j_analytical_do) {
932 dS_dx,
933 jacobian_quantities,
934 ppvar_f(joker, ip),
935 ppath.los(ip, joker),
936 lte[ip],
937 atmosphere_dim,
938 j_analytical_do);
939 }
940
941 if (clear2cloudy[ip] + 1) {
943 Kp,
944 da_dx,
945 dKp_dx,
946 jacobian_quantities,
947 ppvar_pnd(joker, Range(ip, 1)),
948 ppvar_dpnd_dx,
949 ip,
950 scat_data,
951 ppath.los(ip, joker),
952 ppvar_t[Range(ip, 1)],
953 atmosphere_dim,
954 jacobian_do);
955 K_this += Kp;
956
957 if (j_analytical_do) {
958 FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] += dKp_dx[iq];)
959 }
960 }
961
962 // Transmission
963 if (ip not_eq 0) {
964 const Numeric dr_dT_past =
965 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
966 const Numeric dr_dT_this =
967 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
968 stepwise_transmission(lyr_tra[ip],
969 dlyr_tra_above[ip],
970 dlyr_tra_below[ip],
971 K_past,
972 K_this,
973 dK_past_dx,
974 dK_this_dx,
975 ppath.lstep[ip - 1],
976 dr_dT_past,
977 dr_dT_this,
978 temperature_derivative_position);
979 }
980
981 swap(K_past, K_this);
982 swap(dK_past_dx, dK_this_dx);
983 }
984 }
985
986 const ArrayOfTransmissionMatrix tot_tra =
988
989 // iy_aux: Optical depth
990 if (auxOptDepth >= 0) {
991 for (Index iv = 0; iv < nf; iv++)
992 iy_aux[auxOptDepth](iv, 0) = -std::log(tot_tra[np - 1](iv, 0, 0));
993 }
994
995 lvl_rad[np - 1] = iy_transmitter;
996
997 // Radiative transfer calculations
998 for (Index ip = np - 2; ip >= 0; ip--) {
999 lvl_rad[ip] = lvl_rad[ip + 1];
1000 update_radiation_vector(lvl_rad[ip],
1001 dlvl_rad[ip],
1002 dlvl_rad[ip + 1],
1007 lyr_tra[ip + 1],
1008 tot_tra[ip],
1009 dlyr_tra_above[ip + 1],
1010 dlyr_tra_below[ip + 1],
1011 PropagationMatrix(),
1012 PropagationMatrix(),
1013 ArrayOfPropagationMatrix(),
1014 ArrayOfPropagationMatrix(),
1015 Numeric(),
1016 Vector(),
1017 Vector(),
1018 0,
1019 0,
1021 }
1022
1023 // Copy back to ARTS external style
1024 iy = lvl_rad[0];
1025 //
1026 if (np > 1) {
1027 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
1028 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
1029 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
1030 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
1031 if (j_analytical_do)
1032 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
1033 dlvl_rad[ip][iq];);
1034 }
1035 }
1036
1037 // Finalize analytical Jacobians
1038 if (j_analytical_do) {
1040 diy_dx,
1041 diy_dpath,
1042 ns,
1043 nf,
1044 np,
1045 atmosphere_dim,
1046 ppath,
1047 ppvar_p,
1048 ppvar_t,
1049 ppvar_vmr,
1050 iy_agenda_call1,
1051 iy_transmittance,
1052 water_p_eq_agenda,
1053 jacobian_quantities,
1054 jac_species_i);
1055 }
1056}
1057
1058/* Workspace method: Doxygen documentation will be auto-generated */
1059void iy_transmitterMultiplePol(Matrix& iy_transmitter,
1060 const Index& stokes_dim,
1061 const Vector& f_grid,
1062 const ArrayOfIndex& instrument_pol,
1063 const Verbosity&) {
1064 const Index nf = f_grid.nelem();
1065
1066 if (instrument_pol.nelem() != nf)
1067 throw runtime_error(
1068 "The length of *f_grid* and the number of elements "
1069 "in *instrument_pol* must be equal.");
1070
1071 iy_transmitter.resize(nf, stokes_dim);
1072
1073 for (Index i = 0; i < nf; i++) {
1074 stokes2pol(iy_transmitter(i, joker), stokes_dim, instrument_pol[i], 1);
1075 }
1076}
1077
1078/* Workspace method: Doxygen documentation will be auto-generated */
1079void iy_transmitterSinglePol(Matrix& iy_transmitter,
1080 const Index& stokes_dim,
1081 const Vector& f_grid,
1082 const ArrayOfIndex& instrument_pol,
1083 const Verbosity&) {
1084 const Index nf = f_grid.nelem();
1085
1086 if (instrument_pol.nelem() != 1)
1087 throw runtime_error(
1088 "The number of elements in *instrument_pol* must be 1.");
1089
1090 iy_transmitter.resize(nf, stokes_dim);
1091
1092 stokes2pol(iy_transmitter(0, joker), stokes_dim, instrument_pol[0], 1);
1093
1094 for (Index i = 1; i < nf; i++) {
1095 iy_transmitter(i, joker) = iy_transmitter(0, joker);
1096 }
1097}
Array< Index > ArrayOfIndex
An array of Index.
Definition array.h:119
base min(const Array< base > &x)
Min function.
Definition array.h:144
The global header file for ARTS.
Constants of physical expressions as constexpr.
Common ARTS conversions.
void gas_scattering_agendaExecute(Workspace &ws, PropagationMatrix &gas_scattering_coef, TransmissionMatrix &gas_scattering_mat, Vector &gas_scattering_fct_legendre, const Vector &f_grid, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Vector &gas_scattering_los_in, const Vector &gas_scattering_los_out, const Index gas_scattering_output_type, const Agenda &input_agenda)
Definition auto_md.cc:26834
The Agenda class.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Workspace class.
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition jacobian.cc:548
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition jacobian.cc:589
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition jacobian.cc:580
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
Definition jacobian.cc:602
Routines for setting up the jacobian.
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition jacobian.h:532
void iy_transmitterMultiplePol(Matrix &iy_transmitter, const Index &stokes_dim, const Vector &f_grid, const ArrayOfIndex &instrument_pol, const Verbosity &)
WORKSPACE METHOD: iy_transmitterMultiplePol.
constexpr Numeric SPEED_OF_LIGHT
constexpr Numeric DEG2RAD
constexpr Numeric RAD2DEG
void iy_transmitterSinglePol(Matrix &iy_transmitter, const Index &stokes_dim, const Vector &f_grid, const ArrayOfIndex &instrument_pol, const Verbosity &)
WORKSPACE METHOD: iy_transmitterSinglePol.
void iyTransmissionStandard(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_pnd, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_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 Index &gas_scattering_do, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Matrix &iy_transmitter, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Verbosity &)
WORKSPACE METHOD: iyTransmissionStandard.
constexpr Numeric PI
Declarations having to do with the four output streams.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition ppath.cc:1460
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstMatrixView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstVectorView &ppath_line_of_sight, const ConstVectorView &ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
Definition rte.cc:1188
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
Definition rte.cc:828
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, const ConstVectorView &f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, const ConstMatrixView &ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
Definition rte.cc:1051
void rtmethods_jacobian_finalisation(Workspace &ws, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &atmosphere_dim, const Ppath &ppath, const Vector &ppvar_p, const Vector &ppvar_t, const Matrix &ppvar_vmr, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &jac_species_i)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
Definition rte.cc:1948
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &ppath_f_grid, const Vector &ppath_magnetic_field, const Vector &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const Vector &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition rte.cc:1116
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_line_of_sight, const Index &lte, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
Definition rte.cc:38
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Definition rte.cc:947
Declaration of functions in rte.cc.
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Definition sensor.cc:978
Sensor modelling functions.
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
Index np
Number of points describing the ppath.
Vector lstep
The length between ppath points.
Radiation Vector for Stokes dimension 1-4.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void update_radiation_vector(RadiationVector &I, ArrayOfRadiationVector &dI1, ArrayOfRadiationVector &dI2, const RadiationVector &J1, const RadiationVector &J2, const ArrayOfRadiationVector &dJ1, const ArrayOfRadiationVector &dJ2, const TransmissionMatrix &T, const TransmissionMatrix &PiT, const ArrayOfTransmissionMatrix &dT1, const ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the Radiation Vector.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
Array< RadiationVector > ArrayOfRadiationVector
#define a
Array< TransmissionMatrix > ArrayOfTransmissionMatrix