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