ARTS 2.5.0 (git: 9ee3ac6c)
m_fos.cc
Go to the documentation of this file.
1/* Copyright (C) 2013
2 Patrick Eriksson <patrick.eriksson@chalmers.se>
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
19/*===========================================================================
20 === File description
21 ===========================================================================*/
22
34/*===========================================================================
35 === External declarations
36 ===========================================================================*/
37
38#include "arts.h"
39#include "arts_omp.h"
40#include "auto_md.h"
41#include "doit.h"
42#include "m_xml.h"
43#include "math_funcs.h"
44#include "montecarlo.h"
45#include "rte.h"
46#include <cmath>
47#include <stdexcept>
48
49extern const Numeric DEG2RAD;
50extern const Numeric RAD2DEG;
51extern const Numeric PI;
52
53// FOS implemented as an internal function, to allow an recursive algorithm
54/*
55void fos(
56 Workspace& ws,
57 Matrix& iy,
58 ArrayOfTensor4& iy_aux,
59 Ppath& ppath,
60 ArrayOfTensor3& diy_dx,
61 const Index& stokes_dim,
62 const Vector& f_grid,
63 const Index& atmosphere_dim,
64 const Vector& p_grid,
65 const Tensor3& z_field,
66 const Tensor3& t_field,
67 const Tensor4& vmr_field,
68 const ArrayOfArrayOfSpeciesTag& abs_species,
69 const Tensor3& wind_u_field,
70 const Tensor3& wind_v_field,
71 const Tensor3& wind_w_field,
72 const Tensor3& mag_u_field,
73 const Tensor3& mag_v_field,
74 const Tensor3& mag_w_field,
75 const Index& cloudbox_on,
76 const ArrayOfIndex& cloudbox_limits,
77 const Tensor4& pnd_field,
78 const ArrayOfArrayOfSingleScatteringData& scat_data,
79 const Matrix& particle_masses,
80 const String& iy_unit,
81 const ArrayOfString& iy_aux_vars,
82 const Index& jacobian_do,
83 const Agenda& ppath_agenda,
84 const Agenda& propmat_clearsky_agenda,
85 const Agenda& iy_main_agenda,
86 const Agenda& iy_space_agenda,
87 const Agenda& iy_surface_agenda,
88 const Index& iy_agenda_call1,
89 const Tensor3& iy_transmittance,
90 const Vector& rte_pos,
91 const Vector& rte_los,
92 const Vector& rte_pos2,
93 const Numeric& rte_alonglos_v,
94 const Numeric& ppath_lmax,
95 const Numeric& ppath_lraytrace,
96 const Matrix& fos_scatint_angles,
97 const Vector& fos_iyin_za_angles,
98 const Index& fos_za_interporder,
99 const Index& fos_n,
100 const Index& fos_i,
101 const Verbosity& verbosity )
102{
103 // A temporary restriction
104 if( atmosphere_dim > 1 )
105 throw runtime_error( "FOS is so far only handling 1D atmospheres." );
106
107 ARTS_ASSERT( fos_i >= 0 && fos_i <= fos_n );
108
109
110 // Determine propagation path
111 //
112 ppath_agendaExecute( ws, ppath, ppath_lmax, ppath_lraytrace,
113 rte_pos, rte_los, rte_pos2,
114 0, 0, t_field, z_field, vmr_field, f_grid,
115 ppath_agenda );
116
117 // Some basic sizes
118 //
119 const Index nf = f_grid.nelem();
120 const Index ns = stokes_dim;
121 const Index np = ppath.np;
122
123 // The below copied from iyEmission. Activate for-loop when jacobians
124 // introduced.
125
126 // Set up variable with index of species where we want abs_per_species.
127 // This variable can below be extended in iy_aux part.
128 //
129 ArrayOfIndex iaps(0);
130 //
131 //for( Index i=0; i<jac_species_i.nelem(); i++ )
132 // {
133 // if( jac_species_i[i] >= 0 )
134 // { iaps.push_back( jac_species_i[i] ); }
135 // }
136
137 //=== iy_aux part ===========================================================
138 Index auxPressure = -1,
139 auxTemperature = -1,
140 auxAbsSum = -1,
141 auxBackground = -1,
142 auxIy = -1,
143 auxOptDepth = -1;
144 ArrayOfIndex auxAbsSpecies(0), auxAbsIsp(0);
145 ArrayOfIndex auxVmrSpecies(0), auxVmrIsp(0);
146 ArrayOfIndex auxPartCont(0), auxPartContI(0);
147 ArrayOfIndex auxPartField(0), auxPartFieldI(0);
148 //
149 if( !iy_agenda_call1 )
150 { iy_aux.resize( 0 ); }
151 else
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].substr(0,20) == "Absorption, species " )
181 {
182 Index ispecies;
183 istringstream is(iy_aux_vars[i].substr(20,2));
184 is >> ispecies;
185 if( ispecies < 0 || ispecies>=abs_species.nelem() )
186 {
187 ostringstream os;
188 os << "You have selected absorption species with index "
189 << ispecies << ".\nThis species does not exist!";
190 throw runtime_error( os.str() );
191 }
192 auxAbsSpecies.push_back(i);
193 const Index ihit = find_first( iaps, ispecies );
194 if( ihit >= 0 )
195 { auxAbsIsp.push_back( ihit ); }
196 else
197 {
198 iaps.push_back(ispecies);
199 auxAbsIsp.push_back( iaps.nelem()-1 );
200 }
201 iy_aux[i].resize( nf, ns, ns, np );
202 }
203 else if( iy_aux_vars[i] == "Radiative background" )
204 { auxBackground = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
205 else if( iy_aux_vars[i] == "iy" && auxIy < 0 )
206 { auxIy = i; iy_aux[i].resize( nf, ns, 1, np ); }
207 else if( iy_aux_vars[i] == "Optical depth" )
208 { auxOptDepth = i; iy_aux[i].resize( nf, 1, 1, 1 ); }
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
242 {
243 ostringstream os;
244 os << "In *iy_aux_vars* you have included: \"" << iy_aux_vars[i]
245 << "\"\nThis choice is not recognised.";
246 throw runtime_error( os.str() );
247 }
248 }
249 }
250 //===========================================================================
251
252 // Get atmospheric and RT quantities for each ppath point/step
253 //
254 Vector ppath_p, ppath_t;
255 Matrix ppath_vmr, ppath_pnd, ppath_wind, ppath_mag, ppath_f, ppath_nlte;
256 Matrix ppath_blackrad;
257 ArrayOfArrayOfPropagationMatrix abs_per_species, dummy_dppath_ext_dx;
258 ArrayOfPropagationMatrix ppath_ext, pnd_ext_mat;
259 Tensor4 trans_partial, trans_cumulat;
260 ArrayOfArrayOfStokesVector dummy_dppath_nlte_dx;
261 Tensor3 pnd_abs_vec;
262 ArrayOfStokesVector ppath_nlte_source;
263 Vector scalar_tau;
264 ArrayOfIndex clear2cloudy, lte;
265 ArrayOfMatrix dummy_ppath_dpnd_dx;
266 ArrayOfTensor4 dummy_dpnd_field_dx;
267 const Tensor4 nlte_field_empty(0,0,0,0);
268 //
269 Array<ArrayOfArrayOfSingleScatteringData> scat_data_single;
270 ArrayOfArrayOfIndex extmat_case;
271 //
272 if( np > 1 )
273 {
274 get_ppath_atmvars( ppath_p, ppath_t, ppath_nlte, ppath_vmr,
275 ppath_wind, ppath_mag,
276 ppath, atmosphere_dim, p_grid, t_field,
277 nlte_field_empty, vmr_field,
278 wind_u_field, wind_v_field, wind_w_field,
279 mag_u_field, mag_v_field, mag_w_field );
280
281 get_ppath_f( ppath_f, ppath, f_grid, atmosphere_dim,
282 rte_alonglos_v, ppath_wind );
283
284 get_ppath_pmat( ws, ppath_ext, ppath_nlte_source, lte, abs_per_species,
285 dummy_dppath_ext_dx, dummy_dppath_nlte_dx,
286 propmat_clearsky_agenda, ArrayOfRetrievalQuantity(0), ppath,
287 ppath_p, ppath_t, ppath_nlte, ppath_vmr, ppath_f,
288 ppath_mag, f_grid, stokes_dim, iaps );
289
290 for( Index i=0; i<lte.nelem(); i++ )
291 {
292 if( lte[i] == 0 )
293 throw runtime_error( "FOS can so far only handle LTE conditions." );
294 }
295
296 get_ppath_blackrad( ppath_blackrad, ppath, ppath_t, ppath_f );
297
298 if( !cloudbox_on )
299 {
300 get_ppath_trans( trans_partial, extmat_case, trans_cumulat,
301 scalar_tau, ppath, ppath_ext, f_grid, stokes_dim );
302 }
303 else
304 {
305 get_ppath_cloudvars( clear2cloudy, ppath_pnd, dummy_ppath_dpnd_dx,
306 ppath, atmosphere_dim, cloudbox_limits,
307 pnd_field, dummy_dpnd_field_dx );
308 get_ppath_partopt( pnd_abs_vec, pnd_ext_mat, scat_data_single,
309 clear2cloudy, ppath_pnd, ppath, ppath_t,
310 stokes_dim, ppath_f, atmosphere_dim, scat_data,
311 verbosity );
312 get_ppath_trans2( trans_partial, extmat_case, trans_cumulat,
313 scalar_tau, ppath, ppath_ext, f_grid, stokes_dim,
314 clear2cloudy, pnd_ext_mat );
315 }
316 }
317 else // For cases totally outside the atmosphere,
318 { // set zero optical thickness and unit transmission
319 scalar_tau.resize( nf );
320 scalar_tau = 0;
321 trans_cumulat.resize( nf, ns, ns, np );
322 for( Index iv=0; iv<nf; iv++ )
323 { id_mat( trans_cumulat(iv,joker,joker,np-1) ); }
324 }
325
326 // iy_transmittance
327 //
328 Tensor3 iy_trans_new;
329 //
330 if( iy_agenda_call1 )
331 { iy_trans_new = trans_cumulat(joker,joker,joker,np-1); }
332 else
333 {
334 iy_transmittance_mult( iy_trans_new, iy_transmittance,
335 trans_cumulat(joker,joker,joker,np-1) ); }
336
337 // Radiative background
338 //
339 {
340 Agenda iy_cbox_agenda;
341 const Index iy_id = 0;
342 get_iy_of_background( ws, iy, diy_dx,
343 iy_trans_new, iy_id, jacobian_do, ppath, rte_pos2,
344 atmosphere_dim, t_field, z_field, vmr_field,
345 cloudbox_on, stokes_dim, f_grid, iy_unit,
346 iy_main_agenda, iy_space_agenda, iy_surface_agenda,
347 iy_cbox_agenda, verbosity );
348 }
349
350 //=== iy_aux part ===========================================================
351 // Fill parts of iy_aux that are defined even for np=1.
352 // Radiative background
353 if( auxBackground >= 0 )
354 { iy_aux[auxBackground](0,0,0,0) = (Numeric)min( (Index)2,
355 ppath_what_background(ppath)-1); }
356 // Radiance
357 if( auxIy >= 0 )
358 { iy_aux[auxIy](joker,joker,0,np-1) = iy; }
359 // Transmission, total
360 if( auxOptDepth >= 0 )
361 { iy_aux[auxOptDepth](joker,0,0,0) = scalar_tau; }
362 //===========================================================================
363
364
365 // Do RT calculations
366 //
367 if( np > 1 )
368 {
369 //=== iy_aux part =======================================================
370 // iy_aux for point np-1:
371 // Pressure
372 if( auxPressure >= 0 )
373 { iy_aux[auxPressure](0,0,0,np-1) = ppath_p[np-1]; }
374 // Temperature
375 if( auxTemperature >= 0 )
376 { iy_aux[auxTemperature](0,0,0,np-1) = ppath_t[np-1]; }
377 // VMR
378 for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
379 { iy_aux[auxVmrSpecies[j]](0,0,0,np-1) = ppath_vmr(auxVmrIsp[j],np-1); }
380 // Absorption
381 if( auxAbsSum >= 0 )
382 { for( Index iv=0; iv<nf; iv++ ) {
383 for( Index is1=0; is1<ns; is1++ ){
384 for( Index is2=0; is2<ns; is2++ ){
385 iy_aux[auxAbsSum](iv,is1,is2,np-1) =
386 ppath_ext[np-1](iv,is1,is2); } } } }
387 for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
388 { for( Index iv=0; iv<nf; iv++ ) {
389 for( Index is1=0; is1<stokes_dim; is1++ ){
390 for( Index is2=0; is2<stokes_dim; is2++ ){
391 iy_aux[auxAbsSpecies[j]](iv,is1,is2,np-1) =
392 abs_per_species[np-1][auxAbsIsp[j]](iv,is1,is2); } } } }
393 // Particle properties
394 if( cloudbox_on )
395 {
396 // Particle mass content
397 for( Index j=0; j<auxPartCont.nelem(); j++ )
398 { iy_aux[auxPartCont[j]](0,0,0,np-1) = ppath_pnd(joker,np-1) *
399 particle_masses(joker,auxPartContI[j]); }
400 // Particle number density
401 for( Index j=0; j<auxPartField.nelem(); j++ )
402 { iy_aux[auxPartField[j]](0,0,0,np-1) =
403 ppath_pnd(auxPartFieldI[j],np-1); }
404 }
405 // Radiance for this point is handled above
406 //=======================================================================
407
408
409 // Scattering source term at ip (0) and ip+1 (1)
410 // (If any particles at ip=np-1, this is ignored. Could happen if
411 // particles at surface level.)
412 Matrix ssource0(nf,ns,0), ssource1(nf,ns);
413
414 // Dummy variables for non-LTE
415 const bool nonlte = false;
416 const Matrix sourcebar_dummy(0,0);
417 const Tensor3 extbar_dummy(0,0,0);
418
419 // Loop ppath steps
420 for( Index ip=np-2; ip>=0; ip-- )
421 {
422 // Path step average of emission source function: Bbar
423 Vector bbar(nf);
424 for( Index iv=0; iv<nf; iv++ )
425 { bbar[iv] = 0.5 * ( ppath_blackrad(iv,ip) +
426 ppath_blackrad(iv,ip+1) ); }
427
428 // Check if any particles to consider
429 bool any_particles = clear2cloudy[ip] >= 0 ||
430 clear2cloudy[ip+1] >= 0;
431
432 // -----------------------------------------------------------------
433 // i = N (only absorption/emission)
434 // -----------------------------------------------------------------
435 if( fos_i == fos_n )
436 {
437 // No particle absorption to consider
438 if( !any_particles )
439 {
440 emission_rtstep( iy, stokes_dim, bbar, extmat_case[ip],
441 trans_partial(joker,joker,joker,ip),
442 nonlte, extbar_dummy, sourcebar_dummy );
443 }
444
445 else // We want to include particle absorption, but not
446 { // extinction. trans_partial is then not valid.
447 Tensor3 t(nf,ns,ns);
448 ArrayOfIndex extmat_cas2(nf);
449 //
450 for( Index iv=0; iv<nf; iv++ )
451 {
452 // Particle absorption
453 //
454 Matrix pabs_mat(ns,ns,0);
455 //
456 if( clear2cloudy[ip] >= 0 )
457 { ext_matFromabs_vec( pabs_mat, pnd_abs_vec(iv,joker,
458 clear2cloudy[ip]), stokes_dim ); }
459 if( clear2cloudy[ip+1] >= 0 )
460 { ext_matFromabs_vec( pabs_mat, pnd_abs_vec(iv,joker,
461 clear2cloudy[ip+1]), stokes_dim ); }
462
463 // Transmission of step
464 Matrix ext_mat(stokes_dim,stokes_dim);
465 for( Index is1=0; is1<stokes_dim; is1++ ) {
466 for( Index is2=0; is2<stokes_dim; is2++ ) {
467 ext_mat(is1,is2) = 0.5 * ( pabs_mat(is1,is2) +
468 ppath_ext[ip](iv,is1,is2) +
469 ppath_ext[ip+1](iv,is1,is2) );
470 } }
471 //
472 extmat_cas2[iv] = 0;
473 ext2trans( t(iv,joker,joker), extmat_cas2[iv],
474 ext_mat, ppath.lstep[ip] );
475 }
476
477 // Perform RT
478 emission_rtstep( iy, stokes_dim, bbar, extmat_cas2, t,
479 nonlte, extbar_dummy, sourcebar_dummy );
480 }
481 }
482
483
484 // -----------------------------------------------------------------
485 // i < N
486 // -----------------------------------------------------------------
487 else
488 {
489 // Shift scattering source term (new 1 is old 0)
490 ssource1 = ssource0;
491
492 // Clear-sky path step
493 if( !any_particles )
494 {
495 // Perform RT
496 emission_rtstep( iy, stokes_dim, bbar, extmat_case[ip],
497 trans_partial(joker,joker,joker,ip),
498 nonlte, extbar_dummy, sourcebar_dummy );
499
500 // Scattering source term at ip is zero:
501 ssource0 = 0;
502 }
503
504 // Include scattering
505 else
506 {
507 // Determine scattering source term at ip
508 if( clear2cloudy[ip] < 0 )
509 { ssource0 = 0; }
510 else
511 {
512 // Present position
513 // (Note that the Ppath positions (ppath.pos) for 1D have
514 // one column more than expected by most functions. Only
515 // the first atmosphere_dim values shall be copied.)
516 Vector pos = ppath.pos(ip,Range(0,atmosphere_dim));
517
518 // Determine incoming radiation
519 //
520 const Index nin = fos_scatint_angles.nrows();
521 Tensor3 Y(nin,nf,ns);
522 {
523 // Do RT
524 const Index nza = fos_iyin_za_angles.nelem();
525 Tensor3 Y1(nza,nf,ns);
526 //
527 for( Index i=0; i<nza; i++ )
528 {
529 // LOS
530 Vector los( 1, fos_iyin_za_angles[i] );
531
532 // Call recursively, with fos_i increased
533 //
534 Matrix iyl;
535 ArrayOfTensor4 iy_auxl;
536 Ppath ppathl;
537 ArrayOfTensor3 diy_dxl;
538 //
539 fos( ws, iyl, iy_auxl, ppathl, diy_dxl,
540 stokes_dim, f_grid, atmosphere_dim,
541 p_grid, z_field, t_field, vmr_field,
542 abs_species, wind_u_field, wind_v_field,
543 wind_w_field, mag_u_field, mag_v_field,
544 mag_w_field, cloudbox_on, cloudbox_limits,
545 pnd_field,
546 scat_data,
547 particle_masses, iy_unit, iy_aux_vars,
548 jacobian_do, ppath_agenda,
549 propmat_clearsky_agenda, iy_main_agenda,
550 iy_space_agenda, iy_surface_agenda, 0,
551 iy_trans_new, pos, los, rte_pos2,
552 rte_alonglos_v, ppath_lmax, ppath_lraytrace,
553 fos_scatint_angles, fos_iyin_za_angles,
554 fos_za_interporder, fos_n, fos_i+1,
555 verbosity );
556
557 Y1(i,joker,joker) = iyl;
558 }
559
560 // Zenith angle interpolation of Y
561 //ArrayOfGridPosPoly gp( nin ); FIXME: REMOVE WHEN UNCOMMENTING THIS CODE
562 //gridpos_poly( gp, fos_iyin_za_angles,
563 // fos_scatint_angles(joker,0),
564 // fos_za_interporder );
565 const auto lag = Interpolation::LagrangeVector(fos_scatint_angles(joker, 0), fos_iyin_za_angles, 0.5, false, Interpolation::LagrangeType::Linear);
566 //Matrix itw( nin, fos_za_interporder+1 );
567 //interpweights( itw, gp );
568 const auto itw = interpweights( lag )
569 //
570 for( Index iv=0; iv<nf; iv++ )
571 {
572 for( Index is1=0; is1<stokes_dim; is1++ )
573 {
574 //interp( Y(joker,iv,is1), itw,
575 // Y1(joker,iv,is1), gp );
576 reinterp( Y(joker,iv,is1), Y1(joker,iv,is1), itw, lag );
577 }
578 }
579 }
580
581 // Direction of outgoing scattered radiation (which is
582 // reversed to LOS). Note that this outlos is only used
583 // for extracting scattering properties.
584 Vector outlos;
585 mirror_los( outlos, ppath.los(ip,joker), atmosphere_dim );
586
587 // Determine phase matrix
588 Tensor4 P( nin, nf, stokes_dim, stokes_dim );
589 Matrix P1( stokes_dim, stokes_dim );
590 //
591 for( Index ii=0; ii<nin; ii++ )
592 {
593 for( Index iv=0; iv<nf; iv++ )
594 {
595 pha_mat_singleCalc( P1, outlos[0], outlos[1],
596 fos_scatint_angles(ii,0),
597 fos_scatint_angles(ii,1),
598 scat_data_single[iv], stokes_dim,
599 ppath_pnd(joker,ip),
600 ppath_t[ip], verbosity );
601 P(ii,iv,joker,joker) = P1;
602 }
603 }
604
605
606 // Scattering source term
607 ssource0 = 0.0;
608 for( Index iv=0; iv<nf; iv++ )
609 {
610 Vector sp(stokes_dim);
611 for( Index ii=0; ii<nin; ii++ )
612 {
613 mult( sp, P(ii,iv,joker,joker),
614 Y(ii,iv,joker) );
615 ssource0(iv,joker) += sp;
616 }
617 }
618 ssource0 *= 4*PI/(Numeric)nin;
619 }
620
621 // RT of ppath step
622 for( Index iv=0; iv<nf; iv++ )
623 {
624 // Calculate average of absorption, extinction etc.
625 Matrix ext_mat( stokes_dim, stokes_dim );
626 Vector abs_vec( stokes_dim );
627 Vector sbar( stokes_dim, 0 );
628 //
629 // Contribution from abs_species
630 for( Index is1=0; is1<stokes_dim; is1++ )
631 {
632 abs_vec[is1] = 0.5 * (
633 ppath_ext[ip](iv,is1,0) +
634 ppath_ext[ip+1](iv,is1,0) );
635 for( Index is2=0; is2<stokes_dim; is2++ )
636 {
637 ext_mat(is1,is2) = 0.5 * (
638 ppath_ext[ip](iv,is1,is2) +
639 ppath_ext[ip+1](iv,is1,is2) );
640 }
641 }
642 // Particle contribution
643 if( clear2cloudy[ip] >= 0 )
644 {
645 for( Index is1=0; is1<stokes_dim; is1++ )
646 {
647 sbar[is1] += 0.5 * ssource0(iv,is1);
648 abs_vec[is1] += 0.5 * (
649 pnd_abs_vec(iv,is1,clear2cloudy[ip]) );
650 for( Index is2=0; is2<stokes_dim; is2++ )
651 {
652 ext_mat(is1,is2) += 0.5 * (
653 pnd_ext_mat[clear2cloudy[ip]](iv,is1,is2) );
654 }
655 }
656 }
657 if( clear2cloudy[ip+1] >= 0 )
658 {
659 for( Index is1=0; is1<stokes_dim; is1++ )
660 {
661 sbar[is1] += 0.5 * ssource1(iv,is1);
662 abs_vec[is1] += 0.5 * (
663 pnd_abs_vec(iv,is1,clear2cloudy[ip+1]) );
664 for( Index is2=0; is2<stokes_dim; is2++ )
665 {
666 ext_mat(is1,is2) += 0.5 * (
667 pnd_ext_mat[clear2cloudy[ip+1]](iv,is1,is2));
668 }
669 }
670 }
671
672 // Perform RT
673 //
674 Matrix trans_mat = trans_partial(iv,joker,joker,ip);
675 rte_step_doit( iy(iv,joker), trans_mat, ext_mat, abs_vec,
676 sbar, ppath.lstep[ip], bbar[iv], true );
677 }
678 }
679 }
680
681 //=== iy_aux part ===================================================
682 // Pressure
683 if( auxPressure >= 0 )
684 { iy_aux[auxPressure](0,0,0,ip) = ppath_p[ip]; }
685 // Temperature
686 if( auxTemperature >= 0 )
687 { iy_aux[auxTemperature](0,0,0,ip) = ppath_t[ip]; }
688 // VMR
689 for( Index j=0; j<auxVmrSpecies.nelem(); j++ )
690 { iy_aux[auxVmrSpecies[j]](0,0,0,ip) = ppath_vmr(auxVmrIsp[j],ip);}
691 // Absorption
692 if( auxAbsSum >= 0 )
693 { for( Index iv=0; iv<nf; iv++ ) {
694 for( Index is1=0; is1<ns; is1++ ){
695 for( Index is2=0; is2<ns; is2++ ){
696 iy_aux[auxAbsSum](iv,is1,is2,ip) =
697 ppath_ext[ip](iv,is1,is2); } } } }
698 for( Index j=0; j<auxAbsSpecies.nelem(); j++ )
699 { for( Index iv=0; iv<nf; iv++ ) {
700 for( Index is1=0; is1<stokes_dim; is1++ ){
701 for( Index is2=0; is2<stokes_dim; is2++ ){
702 iy_aux[auxAbsSpecies[j]](iv,is1,is2,ip) =
703 abs_per_species[ip][auxAbsIsp[j]](iv,is1,is2); } } } }
704 // Particle properties
705 if( cloudbox_on )
706 {
707 // Particle mass content
708 for( Index j=0; j<auxPartCont.nelem(); j++ )
709 { iy_aux[auxPartCont[j]](0,0,0,ip) = ppath_pnd(joker,ip) *
710 particle_masses(joker,auxPartContI[j]); }
711 // Particle number density
712 for( Index j=0; j<auxPartField.nelem(); j++ )
713 { iy_aux[auxPartField[j]](0,0,0,ip) =
714 ppath_pnd(auxPartFieldI[j],ip); }
715 }
716 // Radiance
717 if( auxIy >= 0 )
718 { iy_aux[auxIy](joker,joker,0,ip) = iy; }
719 //===================================================================
720 }
721 } // if np>1
722
723
724 // Unit conversions
725 if( iy_agenda_call1 )
726 {
727 // Determine refractive index to use for the n2 radiance law
728 Numeric n = 1.0; // First guess is that sensor is in space
729 //
730 if( ppath.end_lstep == 0 ) // If true, sensor inside the atmosphere
731 { n = ppath.nreal[np-1]; }
732
733 // Polarisation index variable
734 ArrayOfIndex i_pol(stokes_dim);
735 for( Index is=0; is<stokes_dim; is++ )
736 { i_pol[is] = is + 1; }
737
738 // iy
739 apply_iy_unit( iy, iy_unit, f_grid, n, i_pol );
740
741 // iy_aux
742 for( Index q=0; q<iy_aux.nelem(); q++ )
743 {
744 if( iy_aux_vars[q] == "iy")
745 {
746 for( Index ip=0; ip<ppath.np; ip++ )
747 { apply_iy_unit( iy_aux[q](joker,joker,0,ip), iy_unit, f_grid,
748 ppath.nreal[ip], i_pol ); }
749 }
750 }
751 }
752}
753*/
754
755/*===========================================================================
756 === The functions (in alphabetical order)
757 ===========================================================================*/
758
759/* Workspace method: Doxygen documentation will be auto-generated */
760/*
761void iyFOS(
762 Workspace& ws,
763 Matrix& iy,
764 ArrayOfTensor4& iy_aux,
765 Ppath& ppath,
766 ArrayOfTensor3& diy_dx,
767 const Index& stokes_dim,
768 const Vector& f_grid,
769 const Index& atmosphere_dim,
770 const Vector& p_grid,
771 const Tensor3& z_field,
772 const Tensor3& t_field,
773 const Tensor4& vmr_field,
774 const ArrayOfArrayOfSpeciesTag& abs_species,
775 const Tensor3& wind_u_field,
776 const Tensor3& wind_v_field,
777 const Tensor3& wind_w_field,
778 const Tensor3& mag_u_field,
779 const Tensor3& mag_v_field,
780 const Tensor3& mag_w_field,
781 const Index& cloudbox_on,
782 const ArrayOfIndex& cloudbox_limits,
783 const Tensor4& pnd_field,
784 const ArrayOfArrayOfSingleScatteringData& scat_data,
785 const Matrix& particle_masses,
786 const String& iy_unit,
787 const ArrayOfString& iy_aux_vars,
788 const Index& jacobian_do,
789 const Agenda& ppath_agenda,
790 const Agenda& propmat_clearsky_agenda,
791 const Agenda& iy_main_agenda,
792 const Agenda& iy_space_agenda,
793 const Agenda& iy_surface_agenda,
794 const Index& iy_agenda_call1,
795 const Tensor3& iy_transmittance,
796 const Vector& rte_pos,
797 const Vector& rte_los,
798 const Vector& rte_pos2,
799 const Numeric& rte_alonglos_v,
800 const Numeric& ppath_lmax,
801 const Numeric& ppath_lraytrace,
802 const Matrix& fos_scatint_angles,
803 const Vector& fos_iyin_za_angles,
804 const Index& fos_za_interporder,
805 const Index& fos_n,
806 const Verbosity& verbosity )
807{
808 // Input checks
809 if( jacobian_do )
810 throw runtime_error(
811 "This method does not yet provide any jacobians (jacobian_do must be 0)" );
812 if( fos_scatint_angles.ncols() != 2 )
813 throw runtime_error( "The WSV *fos_scatint_angles* must have two columns." );
814 if( min(fos_scatint_angles(joker,0))<0 ||
815 max(fos_scatint_angles(joker,0))>180 )
816 throw runtime_error(
817 "The zenith angles in *fos_scatint_angles* shall be inside [0,180]." );
818 if( min(fos_scatint_angles(joker,1))<-180 ||
819 max(fos_scatint_angles(joker,1))>180 )
820 throw runtime_error(
821 "The azimuth angles in *fos_scatint_angles* shall be inside [-180,180]." );
822 if( min(fos_iyin_za_angles)<0 || max(fos_iyin_za_angles)>180 )
823 throw runtime_error(
824 "The zenith angles in *fos_iyin_za_angles* shall be inside [0,180]." );
825 if( fos_iyin_za_angles[0] != 0 )
826 throw runtime_error( "The first value in *fos_iyin_za_angles* must be 0." );
827 if( last(fos_iyin_za_angles) != 180 )
828 throw runtime_error( "The last value in *fos_iyin_za_angles* must be 180." );
829 if( fos_za_interporder < 1 )
830 throw runtime_error( "The argument *fos_za_interporder* must be >= 1." );
831 if( fos_iyin_za_angles.nelem() <= fos_za_interporder )
832 throw runtime_error( "The length of *fos_iyin_za_angles* must at least "
833 "be *fos_za_interporder* + 1." );
834 if( fos_n < 0 )
835 throw runtime_error( "The argument *fos_n* must be >= 0." );
836
837 // Switch to order 0 if not primary call
838 // (This happens after surface reflection. If fos_n used (and >=1), new
839 // surface relections are created ..., and recursion never ends.)
840 Index n = fos_n;
841 if( !iy_agenda_call1 )
842 { n = 0; }
843
844 fos( ws, iy, iy_aux, ppath, diy_dx, stokes_dim, f_grid, atmosphere_dim,
845 p_grid, z_field, t_field, vmr_field, abs_species, wind_u_field,
846 wind_v_field, wind_w_field, mag_u_field, mag_v_field, mag_w_field,
847 cloudbox_on, cloudbox_limits, pnd_field,
848 scat_data, particle_masses, iy_unit, iy_aux_vars, jacobian_do,
849 ppath_agenda, propmat_clearsky_agenda,
850 iy_main_agenda, iy_space_agenda, iy_surface_agenda, iy_agenda_call1,
851 iy_transmittance, rte_pos, rte_los, rte_pos2, rte_alonglos_v,
852 ppath_lmax, ppath_lraytrace, fos_scatint_angles, fos_iyin_za_angles,
853 fos_za_interporder, n, 0, verbosity );
854}
855*/
856
857/* Workspace method: Doxygen documentation will be auto-generated */
859 Matrix& iy,
860 ArrayOfMatrix& iy_aux,
861 ArrayOfTensor3& diy_dx,
862 Vector& ppvar_p,
863 Vector& ppvar_t,
864 EnergyLevelMap& ppvar_nlte,
865 Matrix& ppvar_vmr,
866 Matrix& ppvar_wind,
867 Matrix& ppvar_mag,
868 Matrix& ppvar_pnd,
869 Matrix& ppvar_f,
870 Tensor3& ppvar_iy,
871 Tensor4& ppvar_trans_cumulat,
872 Tensor4& ppvar_trans_partial,
873 const Index& iy_id,
874 const Index& stokes_dim,
875 const Vector& f_grid,
876 const Index& atmosphere_dim,
877 const Vector& p_grid,
878 const Tensor3& t_field,
879 const EnergyLevelMap& nlte_field,
880 const Tensor4& vmr_field,
881 const ArrayOfArrayOfSpeciesTag& abs_species,
882 const Tensor3& wind_u_field,
883 const Tensor3& wind_v_field,
884 const Tensor3& wind_w_field,
885 const Tensor3& mag_u_field,
886 const Tensor3& mag_v_field,
887 const Tensor3& mag_w_field,
888 const Index& cloudbox_on,
889 const ArrayOfIndex& cloudbox_limits,
890 const Tensor4& pnd_field,
891 const ArrayOfTensor4& dpnd_field_dx,
892 const ArrayOfString& scat_species,
893 const ArrayOfArrayOfSingleScatteringData& scat_data,
894 const String& iy_unit,
895 const ArrayOfString& iy_aux_vars,
896 const Index& jacobian_do,
897 const ArrayOfRetrievalQuantity& jacobian_quantities,
898 const Agenda& propmat_clearsky_agenda,
899 const Agenda& water_p_eq_agenda,
900 const String& rt_integration_option,
901 const Agenda& iy_main_agenda,
902 const Agenda& iy_space_agenda,
903 const Agenda& iy_surface_agenda,
904 const Agenda& iy_cloudbox_agenda,
905 const Index& iy_agenda_call1,
906 const Tensor3& iy_transmittance,
907 const Ppath& ppath,
908 const Vector& rte_pos2,
909 const Numeric& rte_alonglos_v,
910 const Tensor3& surface_props_data,
911 const Tensor7& cloudbox_field,
912 const Vector& za_grid,
913 const Index& Naa,
914 const Index& t_interp_order,
915 const Verbosity& verbosity) {
916 // If cloudbox off, switch to use clearsky method
917 if (!cloudbox_on) {
918 Tensor4 dummy;
920 iy,
921 iy_aux,
922 diy_dx,
923 ppvar_p,
924 ppvar_t,
925 ppvar_nlte,
926 ppvar_vmr,
927 ppvar_wind,
928 ppvar_mag,
929 ppvar_f,
930 ppvar_iy,
931 ppvar_trans_cumulat,
932 ppvar_trans_partial,
933 iy_id,
934 stokes_dim,
935 f_grid,
936 atmosphere_dim,
937 p_grid,
938 t_field,
939 nlte_field,
940 vmr_field,
941 abs_species,
942 wind_u_field,
943 wind_v_field,
944 wind_w_field,
945 mag_u_field,
946 mag_v_field,
947 mag_w_field,
948 cloudbox_on,
949 iy_unit,
950 iy_aux_vars,
951 jacobian_do,
952 jacobian_quantities,
953 ppath,
954 rte_pos2,
955 propmat_clearsky_agenda,
956 water_p_eq_agenda,
957 rt_integration_option,
958 iy_main_agenda,
959 iy_space_agenda,
960 iy_surface_agenda,
961 iy_cloudbox_agenda,
962 iy_agenda_call1,
963 iy_transmittance,
964 rte_alonglos_v,
965 surface_props_data,
966 verbosity);
967 return;
968 }
969 // Init Jacobian quantities?
970 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<1>(jacobian_quantities) : 0;
971
972 // Some basic sizes
973 const Index nf = f_grid.nelem();
974 const Index ns = stokes_dim;
975 const Index np = ppath.np;
976 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
977
978 // Radiative background index
979 const Index rbi = ppath_what_background(ppath);
980
981 // Throw error if unsupported features are requested
982 if (atmosphere_dim != 1)
983 throw runtime_error(
984 "With cloudbox on, this method handles only 1D calculations.");
985 if (Naa < 3) throw runtime_error("Naa must be > 2.");
986 if (jacobian_do)
987 if (dpnd_field_dx.nelem() != jacobian_quantities.nelem())
988 throw runtime_error(
989 "*dpnd_field_dx* not properly initialized:\n"
990 "Number of elements in dpnd_field_dx must be equal number of jacobian"
991 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
992 " is calculated/set.");
993 if (rbi < 1 || rbi > 9)
994 throw runtime_error(
995 "ppath.background is invalid. Check your "
996 "calculation of *ppath*?");
997 if (rbi == 3 || rbi == 4)
998 throw runtime_error(
999 "The propagation path ends inside or at boundary of "
1000 "the cloudbox.\nFor this method, *ppath* must be "
1001 "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
1002 // iy_aux_vars checked below
1003 // Checks of i_field
1004 if (cloudbox_field.ncols() != stokes_dim)
1005 throw runtime_error(
1006 "Obtained *cloudbox_field* number of Stokes elements inconsistent with "
1007 "*stokes_dim*.");
1008 if (cloudbox_field.nrows() != 1)
1009 throw runtime_error(
1010 "Obtained *cloudbox_field* has wrong number of azimuth angles.");
1011 if (cloudbox_field.npages() != za_grid.nelem())
1012 throw runtime_error(
1013 "Obtained *cloudbox_field* number of zenith angles inconsistent with "
1014 "*za_grid*.");
1015 if (cloudbox_field.nbooks() != 1)
1016 throw runtime_error(
1017 "Obtained *cloudbox_field* has wrong number of longitude points.");
1018 if (cloudbox_field.nshelves() != 1)
1019 throw runtime_error(
1020 "Obtained *cloudbox_field* has wrong number of latitude points.");
1021 if (cloudbox_field.nvitrines() != cloudbox_limits[1] - cloudbox_limits[0] + 1)
1022 throw runtime_error(
1023 "Obtained *cloudbox_field* number of pressure points inconsistent with "
1024 "*cloudbox_limits*.");
1025 if (cloudbox_field.nlibraries() != nf)
1026 throw runtime_error(
1027 "Obtained *cloudbox_field* number of frequency points inconsistent with "
1028 "*f_grid*.");
1029
1030 // Set diy_dpath if we are doing are doing jacobian calculations
1031 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
1032
1033 // Set the species pointers if we are doing jacobian
1034 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
1035
1036 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
1037 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
1038
1039 // Checks that the scattering species are treated correctly if their derivatives are needed (we can here discard the Array)
1040 if (j_analytical_do and iy_agenda_call1) get_pointers_for_scat_species(jacobian_quantities, scat_species, cloudbox_on);
1041
1042 // Init iy_aux and fill where possible
1043 const Index naux = iy_aux_vars.nelem();
1044 iy_aux.resize(naux);
1045 //
1046 for (Index i = 0; i < naux; i++) {
1047 iy_aux[i].resize(nf, ns);
1048
1049 if (iy_aux_vars[i] == "Optical depth") { /*pass*/
1050 } // Filled below
1051 else if (iy_aux_vars[i] == "Radiative background")
1052 iy_aux[i] = (Numeric)min((Index)2, rbi - 1);
1053 else {
1054 ostringstream os;
1055 os << "The only allowed strings in *iy_aux_vars* are:\n"
1056 << " \"Radiative background\"\n"
1057 << " \"Optical depth\"\n"
1058 << "but you have selected: \"" << iy_aux_vars[i] << "\"";
1059 throw runtime_error(os.str());
1060 }
1061 }
1062
1063 // Get atmospheric and radiative variables along the propagation path
1064 ppvar_trans_cumulat.resize(np, nf, ns, ns);
1065 ppvar_trans_partial.resize(np, nf, ns, ns);
1066 ppvar_iy.resize(nf, ns, np);
1067
1069 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
1072 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
1075
1076 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
1078 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
1080
1081 ArrayOfIndex clear2cloudy;
1082 //
1083 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
1084 ppvar_p.resize(0);
1085 ppvar_t.resize(0);
1086 ppvar_vmr.resize(0, 0);
1087 ppvar_wind.resize(0, 0);
1088 ppvar_mag.resize(0, 0);
1089 ppvar_f.resize(0, 0);
1090 ppvar_trans_cumulat = 0;
1091 ppvar_trans_partial = 0;
1092 for (Index iv = 0; iv < nf; iv++) {
1093 for (Index is = 0; is < ns; is++) {
1094 ppvar_trans_cumulat(0,iv,is,is) = 1;
1095 ppvar_trans_partial(0,iv,is,is) = 1;
1096 }
1097 }
1098 } else {
1099 // Basic atmospheric variables
1100 get_ppath_atmvars(ppvar_p,
1101 ppvar_t,
1102 ppvar_nlte,
1103 ppvar_vmr,
1104 ppvar_wind,
1105 ppvar_mag,
1106 ppath,
1107 atmosphere_dim,
1108 p_grid,
1109 t_field,
1110 nlte_field,
1111 vmr_field,
1112 wind_u_field,
1113 wind_v_field,
1114 wind_w_field,
1115 mag_u_field,
1116 mag_v_field,
1117 mag_w_field);
1118
1120 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
1121
1122 // here, the cloudbox is on, ie we don't need to check and branch this here
1123 // anymore.
1124 ArrayOfMatrix ppvar_dpnd_dx;
1125 //
1126 get_ppath_cloudvars(clear2cloudy,
1127 ppvar_pnd,
1128 ppvar_dpnd_dx,
1129 ppath,
1130 atmosphere_dim,
1131 cloudbox_limits,
1132 pnd_field,
1133 dpnd_field_dx);
1134
1135 // Size radiative variables always used
1136 Vector B(nf);
1137 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
1138 StokesVector a(nf, ns), S(nf, ns), Sp(nf, ns);
1139 ArrayOfIndex lte(np);
1140
1141 // Init variables only used if analytical jacobians done
1142 Vector dB_dT(0);
1143 ArrayOfPropagationMatrix dK_this_dx(nq), dK_past_dx(nq), dKp_dx(nq);
1144 ArrayOfStokesVector da_dx(nq), dS_dx(nq), dSp_dx(nq);
1145
1146 // HSE variables
1147 Index temperature_derivative_position = -1;
1148 bool do_hse = false;
1149
1150 if (j_analytical_do) {
1151 dB_dT.resize(nf);
1153 dK_this_dx[iq] = PropagationMatrix(nf, ns);
1154 dK_past_dx[iq] = PropagationMatrix(nf, ns);
1155 dKp_dx[iq] = PropagationMatrix(nf, ns);
1156 da_dx[iq] = StokesVector(nf, ns);
1157 dS_dx[iq] = StokesVector(nf, ns);
1158 dSp_dx[iq] = StokesVector(nf, ns);
1159 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
1160 temperature_derivative_position = iq;
1161 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
1162 })
1163 }
1164 const bool temperature_jacobian =
1165 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
1166
1167 // Loop ppath points and determine radiative properties
1168 for (Index ip = 0; ip < np; ip++) {
1170 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
1171
1173 K_this,
1174 S,
1175 lte[ip],
1176 dK_this_dx,
1177 dS_dx,
1178 propmat_clearsky_agenda,
1179 jacobian_quantities,
1180 ppvar_f(joker, ip),
1181 ppvar_mag(joker, ip),
1182 ppath.los(ip, joker),
1183 ppvar_nlte[ip],
1184 ppvar_vmr(joker, ip),
1185 ppvar_t[ip],
1186 ppvar_p[ip],
1187 j_analytical_do);
1188
1189 if (j_analytical_do)
1191 dS_dx,
1192 jacobian_quantities,
1193 ppvar_f(joker, ip),
1194 ppath.los(ip, joker),
1195 lte[ip],
1196 atmosphere_dim,
1197 j_analytical_do);
1198
1199 if (clear2cloudy[ip] + 1) {
1201 Kp,
1202 da_dx,
1203 dKp_dx,
1204 jacobian_quantities,
1205 ppvar_pnd(joker, Range(ip, 1)),
1206 ppvar_dpnd_dx,
1207 ip,
1208 scat_data,
1209 ppath.los(ip, joker),
1210 ppvar_t[Range(ip, 1)],
1211 atmosphere_dim,
1212 jacobian_do);
1213 a += K_this;
1214 K_this += Kp;
1215
1216 if (j_analytical_do)
1217 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] += dK_this_dx[iq];
1218 dK_this_dx[iq] += dKp_dx[iq];)
1219
1220 Vector aa_grid;
1221 nlinspace(aa_grid, 0, 360, Naa);
1222 //
1224 dSp_dx,
1225 jacobian_quantities,
1226 ppvar_pnd(joker, ip),
1227 ppvar_dpnd_dx,
1228 ip,
1229 scat_data,
1230 cloudbox_field,
1231 za_grid,
1232 aa_grid,
1233 ppath.los(Range(ip, 1), joker),
1234 ppath.gp_p[ip],
1235 ppvar_t[Range(ip, 1)],
1236 atmosphere_dim,
1237 jacobian_do,
1238 t_interp_order);
1239 S += Sp;
1240
1241 if (j_analytical_do)
1242 FOR_ANALYTICAL_JACOBIANS_DO(dS_dx[iq] += dSp_dx[iq];)
1243 } else { // no particles present at this level
1244 a = K_this;
1245 if (j_analytical_do)
1246 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_this_dx[iq];)
1247 }
1248
1249 if (ip not_eq 0) {
1250 const Numeric dr_dT_past =
1251 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
1252 const Numeric dr_dT_this =
1253 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
1254 stepwise_transmission(lyr_tra[ip],
1255 dlyr_tra_above[ip],
1256 dlyr_tra_below[ip],
1257 K_past,
1258 K_this,
1259 dK_past_dx,
1260 dK_this_dx,
1261 ppath.lstep[ip - 1],
1262 dr_dT_past,
1263 dr_dT_this,
1264 temperature_derivative_position);
1265 }
1266
1267 stepwise_source(src_rad[ip],
1268 dsrc_rad[ip],
1269 K_this,
1270 a,
1271 S,
1272 dK_this_dx,
1273 da_dx,
1274 dS_dx,
1275 B,
1276 dB_dT,
1277 jacobian_quantities,
1278 jacobian_do);
1279
1280 swap(K_past, K_this);
1281 swap(dK_past_dx, dK_this_dx);
1282 }
1283 }
1284
1285 const ArrayOfTransmissionMatrix tot_tra =
1287
1288 // iy_transmittance
1289 Tensor3 iy_trans_new;
1290 if (iy_agenda_call1)
1291 iy_trans_new = tot_tra[np - 1];
1292 else
1293 iy_transmittance_mult(iy_trans_new, iy_transmittance, tot_tra[np - 1]);
1294
1295 // Copy transmission to iy_aux
1296 for (Index i = 0; i < naux; i++)
1297 if (iy_aux_vars[i] == "Optical depth")
1298 for (Index iv = 0; iv < nf; iv++)
1299 iy_aux[i](iv, joker) = -log(ppvar_trans_cumulat(np - 1, iv, 0, 0));
1300
1301 // Radiative background
1303 iy,
1304 diy_dx,
1305 iy_trans_new,
1306 iy_id,
1307 jacobian_do,
1308 jacobian_quantities,
1309 ppath,
1310 rte_pos2,
1311 atmosphere_dim,
1312 nlte_field,
1313 cloudbox_on,
1314 stokes_dim,
1315 f_grid,
1316 iy_unit,
1317 surface_props_data,
1318 iy_main_agenda,
1319 iy_space_agenda,
1320 iy_surface_agenda,
1321 iy_cloudbox_agenda,
1322 iy_agenda_call1,
1323 verbosity);
1324
1325 lvl_rad[np - 1] = iy;
1326
1327 // Radiative transfer calculations
1328 for (Index ip = np - 2; ip >= 0; ip--) {
1329 lvl_rad[ip] = lvl_rad[ip + 1];
1330 update_radiation_vector(lvl_rad[ip],
1331 dlvl_rad[ip],
1332 dlvl_rad[ip + 1],
1333 src_rad[ip],
1334 src_rad[ip + 1],
1335 dsrc_rad[ip],
1336 dsrc_rad[ip + 1],
1337 lyr_tra[ip + 1],
1338 tot_tra[ip],
1339 dlyr_tra_above[ip + 1],
1340 dlyr_tra_below[ip + 1],
1345 Numeric(),
1346 Vector(),
1347 Vector(),
1348 0,
1349 0,
1351 }
1352
1353 // Copy back to ARTS external style
1354 iy = lvl_rad[0];
1355 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
1356 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
1357 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
1358 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
1359 if (j_analytical_do)
1360 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
1361 dlvl_rad[ip][iq];);
1362 }
1363
1364 // Finalize analytical Jacobians
1365 if (j_analytical_do)
1367 diy_dx,
1368 diy_dpath,
1369 ns,
1370 nf,
1371 np,
1372 atmosphere_dim,
1373 ppath,
1374 ppvar_p,
1375 ppvar_t,
1376 ppvar_vmr,
1377 iy_agenda_call1,
1378 iy_transmittance,
1379 water_p_eq_agenda,
1380 jacobian_quantities,
1381 jac_species_i);
1382
1383 // Unit conversions
1384 if (iy_agenda_call1)
1386 diy_dx,
1387 ppvar_iy,
1388 ns,
1389 np,
1390 f_grid,
1391 ppath,
1392 jacobian_quantities,
1393 j_analytical_do,
1394 iy_unit);
1395}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:39
Array< Tensor3 > ArrayOfTensor3
An array of Tensor3.
Definition: array.h:73
The global header file for ARTS.
Header file for helper functions for OpenMP.
void iyEmissionStandard(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_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, 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 String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandard.
Definition: m_rte.cc:160
The Agenda class.
Definition: agenda_class.h:44
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:48
Index nrows() const
Returns the number of rows.
Definition: matpackVII.cc:57
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVII.cc:45
Index npages() const
Returns the number of pages.
Definition: matpackVII.cc:54
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:60
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:42
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:51
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
Radiation Vector for Stokes dimension 1-4.
The range class.
Definition: matpackI.h:165
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:339
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:658
The Tensor4 class.
Definition: matpackIV.h:421
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1058
The Tensor7 class.
Definition: matpackVII.h:2382
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
Radiative transfer in cloudbox.
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:564
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:605
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1147
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:596
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:618
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:563
#define ns
#define min(a, b)
void iyHybrid(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 &iy_id, 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 Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Ppath &ppath, const Vector &rte_pos2, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Tensor7 &cloudbox_field, const Vector &za_grid, const Index &Naa, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyHybrid.
Definition: m_fos.cc:858
const Numeric RAD2DEG
const Numeric PI
const Numeric DEG2RAD
Workspace methods and template functions for supergeneric XML IO.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:225
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
const Joker joker
Temperature
Definition: jacobian.h:57
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:1476
#define a
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:1316
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:1131
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const ConstTensor3View &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const ConstVectorView &rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const ConstVectorView &f_grid, const String &iy_unit, const ConstTensor3View &surface_props_data, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Verbosity &verbosity)
Determines iy of the "background" of a propgation path.
Definition: rte.cc:728
void get_stepwise_scattersky_source(StokesVector &Sp, ArrayOfStokesVector &dSp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstTensor7View &cloudbox_field, const ConstVectorView &za_grid, const ConstVectorView &aa_grid, const ConstMatrixView &ppath_line_of_sight, const GridPos &ppath_pressure, const Vector &temperature, const Index &atmosphere_dim, const bool &jacobian_do, const Index &t_interp_order)
Calculates the stepwise scattering source terms.
Definition: rte.cc:1433
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:843
void iy_transmittance_mult(Tensor3 &iy_trans_total, const ConstTensor3View &iy_trans_old, const ConstTensor3View &iy_trans_new)
Multiplicates iy_transmittance with transmissions.
Definition: rte.cc:1973
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:1066
void rtmethods_unit_conversion(Matrix &iy, ArrayOfTensor3 &diy_dx, Tensor3 &ppvar_iy, const Index &ns, const Index &np, const Vector &f_grid, const Ppath &ppath, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &j_analytical_do, const String &iy_unit)
This function handles the unit conversion to be done at the end of some radiative transfer WSMs.
Definition: rte.cc:2206
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:2068
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:55
void get_stepwise_blackbody_radiation(VectorView B, VectorView dB_dT, const ConstVectorView &ppath_f_grid, const Numeric &ppath_temperature, const bool &do_temperature_derivative)
Get the blackbody radiation at propagation path point.
Definition: rte.cc:1116
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:962
Declaration of functions in rte.cc.
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Matrix los
Line-of-sight at each ppath point.
Definition: ppath.h:66
Index np
Number of points describing the ppath.
Definition: ppath.h:52
Vector lstep
The length between ppath points.
Definition: ppath.h:70
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath.h:82
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK, const ArrayOfStokesVector &da, const ArrayOfStokesVector &dS, const ConstVectorView &B, const ConstVectorView &dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Set the stepwise source.
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
Array< TransmissionMatrix > ArrayOfTransmissionMatrix