ARTS 2.5.11 (git: 6827797f)
rte.cc
Go to the documentation of this file.
1
9/*===========================================================================
10 === External declarations
11 ===========================================================================*/
12
13#include "rte.h"
14#include "arts_constexpr_math.h"
15#include "arts_constants.h"
16#include "arts_conversions.h"
17#include "auto_md.h"
18#include "check_input.h"
19#include "geodetic.h"
20#include "lin_alg.h"
21#include "logic.h"
22#include "math_funcs.h"
23#include "montecarlo.h"
24#include "physics_funcs.h"
25#include "ppath.h"
26#include "refraction.h"
27#include "special_interp.h"
28#include <cmath>
29#include <stdexcept>
30
32inline constexpr Numeric TEMP_0_C=Constant::temperature_at_0c;
33
34/*===========================================================================
35 === The functions in alphabetical order
36 ===========================================================================*/
37
39 ArrayOfPropagationMatrix& dK_dx,
40 ArrayOfStokesVector& dS_dx,
41 const ArrayOfRetrievalQuantity& jacobian_quantities,
42 const ConstVectorView& ppath_f_grid,
43 const ConstVectorView& ppath_line_of_sight,
44 const Index& lte,
45 const Index& atmosphere_dim,
46 const bool& jacobian_do) {
47 if (not jacobian_do) return;
48
49 // All relevant quantities are extracted first
50 const Index nq = jacobian_quantities.nelem();
51
52 // Computational temporary vector
53 Vector a;
54
55 for (Index i = 0; i < nq; i++) {
56 if (jacobian_quantities[i] == Jacobian::Type::Sensor or jacobian_quantities[i] == Jacobian::Special::SurfaceString) continue;
57
58 if (jacobian_quantities[i].is_wind()) {
59 const auto scale = get_stepwise_f_partials(ppath_line_of_sight, ppath_f_grid, jacobian_quantities[i].Target().atm, atmosphere_dim);
60 dK_dx[i] *= scale;
61 if (not lte) dS_dx[i] *= scale;
62 }
63 }
64}
65
66void adjust_los(VectorView los, const Index& atmosphere_dim) {
67 if (atmosphere_dim == 1) {
68 if (los[0] < 0) {
69 los[0] = -los[0];
70 } else if (los[0] > 180) {
71 los[0] = 360 - los[0];
72 }
73 } else if (atmosphere_dim == 2) {
74 if (los[0] < -180) {
75 los[0] = los[0] + 360;
76 } else if (los[0] > 180) {
77 los[0] = los[0] - 360;
78 }
79 } else {
80 // If any of the angles out-of-bounds, use cart2zaaa to resolve
81 if (abs(los[0] - 90) > 90 || abs(los[1]) > 180) {
82 Numeric dx, dy, dz;
83 zaaa2cart(dx, dy, dz, los[0], los[1]);
84 cart2zaaa(los[0], los[1], dx, dy, dz);
85 }
86 }
87}
88
89void apply_iy_unit(MatrixView iy,
90 const String& iy_unit,
91 const ConstVectorView& f_grid,
92 const Numeric& n,
93 const ArrayOfIndex& i_pol) {
94 // The code is largely identical between the two apply_iy_unit functions.
95 // If any change here, remember to update the other function.
96
97 const Index nf = iy.nrows();
98 const Index ns = iy.ncols();
99
100 ARTS_ASSERT(f_grid.nelem() == nf);
101 ARTS_ASSERT(i_pol.nelem() == ns);
102
103 if (iy_unit == "1") {
104 if (n != 1) {
105 iy *= (n * n);
106 }
107 }
108
109 else if (iy_unit == "RJBT") {
110 for (Index iv = 0; iv < nf; iv++) {
111 const Numeric scfac = invrayjean(1, f_grid[iv]);
112 for (Index is = 0; is < ns; is++) {
113 if (i_pol[is] < 5) // Stokes components
114 {
115 iy(iv, is) *= scfac;
116 } else // Measuement single pols
117 {
118 iy(iv, is) *= 2 * scfac;
119 }
120 }
121 }
122 }
123
124 else if (iy_unit == "PlanckBT") {
125 for (Index iv = 0; iv < nf; iv++) {
126 for (Index is = ns - 1; is >= 0; is--) // Order must here be reversed
127 {
128 if (i_pol[is] == 1) {
129 iy(iv, is) = invplanck(iy(iv, is), f_grid[iv]);
130 } else if (i_pol[is] < 5) {
131 ARTS_ASSERT(i_pol[0] == 1);
132 iy(iv, is) = invplanck(0.5 * (iy(iv, 0) + iy(iv, is)), f_grid[iv]) -
133 invplanck(0.5 * (iy(iv, 0) - iy(iv, is)), f_grid[iv]);
134 } else {
135 iy(iv, is) = invplanck(2 * iy(iv, is), f_grid[iv]);
136 }
137 }
138 }
139 }
140
141 else if (iy_unit == "W/(m^2 m sr)") {
142 for (Index iv = 0; iv < nf; iv++) {
143 const Numeric scfac = n * n * f_grid[iv] * (f_grid[iv] / SPEED_OF_LIGHT);
144 for (Index is = 0; is < ns; is++) {
145 iy(iv, is) *= scfac;
146 }
147 }
148 }
149
150 else if (iy_unit == "W/(m^2 m-1 sr)") {
151 iy *= (n * n * SPEED_OF_LIGHT);
152 }
153
154 else {
156 "Unknown option: iy_unit = \"", iy_unit, "\"\n"
157 "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
158 "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\"")
159 }
160}
161
162void apply_iy_unit2(Tensor3View J,
163 const ConstMatrixView& iy,
164 const String& iy_unit,
165 const ConstVectorView& f_grid,
166 const Numeric& n,
167 const ArrayOfIndex& i_pol) {
168 // The code is largely identical between the two apply_iy_unit functions.
169 // If any change here, remember to update the other function.
170
171 const Index nf = iy.nrows();
172 const Index ns = iy.ncols();
173 const Index np = J.npages();
174
175 ARTS_ASSERT(J.nrows() == nf);
176 ARTS_ASSERT(J.ncols() == ns);
177 ARTS_ASSERT(f_grid.nelem() == nf);
178 ARTS_ASSERT(i_pol.nelem() == ns);
179
180 if (iy_unit == "1") {
181 if (n != 1) {
182 J *= (n * n);
183 }
184 }
185
186 else if (iy_unit == "RJBT") {
187 for (Index iv = 0; iv < nf; iv++) {
188 const Numeric scfac = invrayjean(1, f_grid[iv]);
189 for (Index is = 0; is < ns; is++) {
190 if (i_pol[is] < 5) // Stokes componenets
191 {
192 for (Index ip = 0; ip < np; ip++) {
193 J(ip, iv, is) *= scfac;
194 }
195 } else // Measuement single pols
196 {
197 for (Index ip = 0; ip < np; ip++) {
198 J(ip, iv, is) *= 2 * scfac;
199 }
200 }
201 }
202 }
203 }
204
205 else if (iy_unit == "PlanckBT") {
206 for (Index iv = 0; iv < f_grid.nelem(); iv++) {
207 for (Index is = ns - 1; is >= 0; is--) {
208 Numeric scfac = 1;
209 if (i_pol[is] == 1) {
210 scfac = dinvplanckdI(iy(iv, is), f_grid[iv]);
211 } else if (i_pol[is] < 5) {
212 ARTS_ASSERT(i_pol[0] == 1);
213 scfac = dinvplanckdI(0.5 * (iy(iv, 0) + iy(iv, is)), f_grid[iv]) +
214 dinvplanckdI(0.5 * (iy(iv, 0) - iy(iv, is)), f_grid[iv]);
215 } else {
216 scfac = dinvplanckdI(2 * iy(iv, is), f_grid[iv]);
217 }
218 //
219 for (Index ip = 0; ip < np; ip++) {
220 J(ip, iv, is) *= scfac;
221 }
222 }
223 }
224 }
225
226 else if (iy_unit == "W/(m^2 m sr)") {
227 for (Index iv = 0; iv < nf; iv++) {
228 const Numeric scfac = n * n * f_grid[iv] * (f_grid[iv] / SPEED_OF_LIGHT);
229 for (Index ip = 0; ip < np; ip++) {
230 for (Index is = 0; is < ns; is++) {
231 J(ip, iv, is) *= scfac;
232 }
233 }
234 }
235 }
236
237 else if (iy_unit == "W/(m^2 m-1 sr)") {
238 J *= (n * n * SPEED_OF_LIGHT);
239 }
240
241 else {
243 "Unknown option: iy_unit = \"", iy_unit, "\"\n"
244 "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
245 "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\"")
246 }
247}
248
249void bending_angle1d(Numeric& alpha, const Ppath& ppath) {
250 Numeric theta;
251 if (ppath.dim < 3) {
252 theta = abs(ppath.start_pos[1] - ppath.end_pos[1]);
253 } else {
254 theta = sphdist(ppath.start_pos[1],
255 ppath.start_pos[2],
256 ppath.end_pos[1],
257 ppath.end_pos[2]);
258 }
259
260 // Eq 17 in Kursinski et al., TAO, 2000:
261 alpha = ppath.start_los[0] - ppath.end_los[0] + theta;
262
263 // This as
264 // phi_r = 180 - ppath.end_los[0]
265 // phi_t = ppath.start_los[0]
266}
267
295 Vector& pos,
296 Vector& rte_los,
297 Index& background,
298 const Vector& rte_pos,
299 const Numeric& lo0,
300 const Agenda& ppath_step_agenda,
301 const Numeric& ppath_lmax,
302 const Numeric& ppath_lraytrace,
303 const Index& atmosphere_dim,
304 const Vector& p_grid,
305 const Vector& lat_grid,
306 const Vector& lon_grid,
307 const Tensor3& z_field,
308 const Vector& f_grid,
309 const Vector& refellipsoid,
310 const Matrix& z_surface,
311 const Verbosity& verbosity) {
312 // Special treatment of 1D around zenith/nadir
313 // (zenith angles outside [0,180] are changed by *adjust_los*)
314 bool invert_lat = false;
315 if (atmosphere_dim == 1 && (rte_los[0] < 0 || rte_los[0] > 180)) {
316 invert_lat = true;
317 }
318
319 // Handle cases where angles have moved out-of-bounds due to disturbance
320 adjust_los(rte_los, atmosphere_dim);
321
322 // Calculate the ppath for disturbed rte_los
323 Ppath ppx;
324 //
325 ppath_calc(ws,
326 ppx,
327 ppath_step_agenda,
328 atmosphere_dim,
329 p_grid,
330 lat_grid,
331 lon_grid,
332 z_field,
333 f_grid,
334 refellipsoid,
335 z_surface,
336 0,
337 ArrayOfIndex(0),
338 rte_pos,
339 rte_los,
340 ppath_lmax,
341 ppath_lraytrace,
342 false,
343 verbosity);
344 //
345 background = ppath_what_background(ppx);
346
347 // Calcualte cumulative optical path for ppx
348 Vector lox(ppx.np);
349 Index ilast = ppx.np - 1;
350 lox[0] = ppx.end_lstep;
351 for (Index i = 1; i <= ilast; i++) {
352 lox[i] =
353 lox[i - 1] + ppx.lstep[i - 1] * (ppx.nreal[i - 1] + ppx.nreal[i]) / 2.0;
354 }
355
356 pos.resize(max(Index(2), atmosphere_dim));
357
358 // Reciever at a longer distance (most likely out in space):
359 if (lox[ilast] < lo0) {
360 const Numeric dl = lo0 - lox[ilast];
361 if (atmosphere_dim < 3) {
362 Numeric x, z, dx, dz;
364 x, z, dx, dz, ppx.r[ilast], ppx.pos(ilast, 1), ppx.los(ilast, 0));
365 cart2pol(pos[0],
366 pos[1],
367 x + dl * dx,
368 z + dl * dz,
369 ppx.pos(ilast, 1),
370 ppx.los(ilast, 0));
371 } else {
372 Numeric x, y, z, dx, dy, dz;
373 poslos2cart(x,
374 y,
375 z,
376 dx,
377 dy,
378 dz,
379 ppx.r[ilast],
380 ppx.pos(ilast, 1),
381 ppx.pos(ilast, 2),
382 ppx.los(ilast, 0),
383 ppx.los(ilast, 1));
384 cart2sph(pos[0],
385 pos[1],
386 pos[2],
387 x + dl * dx,
388 y + dl * dy,
389 z + dl * dz,
390 ppx.pos(ilast, 1),
391 ppx.pos(ilast, 2),
392 ppx.los(ilast, 0),
393 ppx.los(ilast, 1));
394 }
395 }
396
397 // Interpolate to lo0
398 else {
399 GridPos gp;
400 Vector itw(2);
401 gridpos(gp, lox, lo0);
402 interpweights(itw, gp);
403 //
404 pos[0] = interp(itw, ppx.r, gp);
405 pos[1] = interp(itw, ppx.pos(joker, 1), gp);
406 if (atmosphere_dim == 3) {
407 pos[2] = interp(itw, ppx.pos(joker, 2), gp);
408 }
409 }
410
411 if (invert_lat) {
412 pos[1] = -pos[1];
413 }
414}
415
417 Numeric& dlf,
418 const Agenda& ppath_step_agenda,
419 const Index& atmosphere_dim,
420 const Vector& p_grid,
421 const Vector& lat_grid,
422 const Vector& lon_grid,
423 const Tensor3& z_field,
424 const Vector& f_grid,
425 const Vector& refellipsoid,
426 const Matrix& z_surface,
427 const Ppath& ppath,
428 const Numeric& ppath_lmax,
429 const Numeric& ppath_lraytrace,
430 const Numeric& dza,
431 const Verbosity& verbosity) {
432 // Optical and physical path between transmitter and reciver
433 Numeric lo = ppath.start_lstep + ppath.end_lstep;
434 Numeric lp = lo;
435 for (Index i = 0; i <= ppath.np - 2; i++) {
436 lp += ppath.lstep[i];
437 lo += ppath.lstep[i] * (ppath.nreal[i] + ppath.nreal[i + 1]) / 2.0;
438 }
439 // Extract rte_pos and rte_los
440 const Vector rte_pos{ppath.start_pos[Range(0, atmosphere_dim)]};
441 //
442 Vector rte_los0(max(Index(1), atmosphere_dim - 1)), rte_los;
443 mirror_los(rte_los, ppath.start_los, atmosphere_dim);
444 rte_los0 = rte_los[Range(0, max(Index(1), atmosphere_dim - 1))];
445
446 // A new ppath with positive zenith angle off-set
447 //
448 Vector pos1;
449 Index backg1;
450 //
451 rte_los = rte_los0;
452 rte_los[0] += dza;
453 //
455 pos1,
456 rte_los,
457 backg1,
458 rte_pos,
459 lo,
460 ppath_step_agenda,
461 ppath_lmax,
462 ppath_lraytrace,
463 atmosphere_dim,
464 p_grid,
465 lat_grid,
466 lon_grid,
467 z_field,
468 f_grid,
469 refellipsoid,
470 z_surface,
471 verbosity);
472
473 // Same thing with negative zenit angle off-set
474 Vector pos2;
475 Index backg2;
476 //
477 rte_los = rte_los0; // Use rte_los0 as rte_los can have been "adjusted"
478 rte_los[0] -= dza;
479 //
481 pos2,
482 rte_los,
483 backg2,
484 rte_pos,
485 lo,
486 ppath_step_agenda,
487 ppath_lmax,
488 ppath_lraytrace,
489 atmosphere_dim,
490 p_grid,
491 lat_grid,
492 lon_grid,
493 z_field,
494 f_grid,
495 refellipsoid,
496 z_surface,
497 verbosity);
498
499 // Calculate distance between pos1 and 2, and derive the loss factor
500 // All appears OK:
501 if (backg1 == backg2) {
502 Numeric l12;
503 if (atmosphere_dim < 3) {
504 distance2D(l12, pos1[0], pos1[1], pos2[0], pos2[1]);
505 } else {
506 distance3D(l12, pos1[0], pos1[1], pos1[2], pos2[0], pos2[1], pos2[2]);
507 }
508 //
509 dlf = lp * 2 * Conversion::deg2rad(1) * dza / l12;
510 }
511 // If different backgrounds, then only use the second calculation
512 else {
513 Numeric l12;
514 if (atmosphere_dim == 1) {
515 const Numeric r = refellipsoid[0];
516 distance2D(l12, r + ppath.end_pos[0], 0, pos2[0], pos2[1]);
517 } else if (atmosphere_dim == 2) {
518 const Numeric r = refell2r(refellipsoid, ppath.end_pos[1]);
519 distance2D(l12, r + ppath.end_pos[0], ppath.end_pos[1], pos2[0], pos2[1]);
520 } else {
521 const Numeric r = refell2r(refellipsoid, ppath.end_pos[1]);
522 distance3D(l12,
523 r + ppath.end_pos[0],
524 ppath.end_pos[1],
525 ppath.end_pos[2],
526 pos2[0],
527 pos2[1],
528 pos2[2]);
529 }
530 //
531 dlf = lp * Conversion::deg2rad(1) * dza / l12;
532 }
533}
534
536 Numeric& dlf,
537 const Agenda& ppath_step_agenda,
538 const Index& atmosphere_dim,
539 const Vector& p_grid,
540 const Vector& lat_grid,
541 const Vector& lon_grid,
542 const Tensor3& z_field,
543 const Vector& f_grid,
544 const Vector& refellipsoid,
545 const Matrix& z_surface,
546 const Ppath& ppath,
547 const Numeric& ppath_lmax,
548 const Numeric& ppath_lraytrace,
549 const Numeric& dza,
550 const Verbosity& verbosity) {
551 ARTS_USER_ERROR_IF (ppath.end_los[0] < 90 || ppath.start_los[0] > 90,
552 "The function *defocusing_sat2sat* can only be used "
553 "for limb sounding geometry.");
554
555 // Index of tangent point
556 Index it;
557 find_tanpoint(it, ppath);
558 ARTS_ASSERT(it >= 0);
559
560 // Length between tangent point and transmitter/reciver
561 Numeric lt = ppath.start_lstep, lr = ppath.end_lstep;
562 for (Index i = it; i <= ppath.np - 2; i++) {
563 lt += ppath.lstep[i];
564 }
565 for (Index i = 0; i < it; i++) {
566 lr += ppath.lstep[i];
567 }
568
569 // Bending angle and impact parameter for centre ray
570 Numeric alpha0, a0;
571 bending_angle1d(alpha0, ppath);
572 alpha0 *= Conversion::deg2rad(1);
573 a0 = ppath.constant;
574
575 // Azimuth loss term (Eq 18.5 in Kursinski et al.)
576 const Numeric lf = lr * lt / (lr + lt);
577 const Numeric alt = 1 / (1 - alpha0 * lf / refellipsoid[0]);
578
579 // Calculate two new ppaths to get dalpha/da
580 Numeric alpha1, a1, alpha2, a2, dada;
581 Ppath ppt;
582 Vector rte_pos{ppath.end_pos[Range(0, atmosphere_dim)]};
583 Vector rte_los{ppath.end_los};
584 //
585 rte_los[0] -= dza;
586 adjust_los(rte_los, atmosphere_dim);
587 ppath_calc(ws,
588 ppt,
589 ppath_step_agenda,
590 atmosphere_dim,
591 p_grid,
592 lat_grid,
593 lon_grid,
594 z_field,
595 f_grid,
596 refellipsoid,
597 z_surface,
598 0,
599 ArrayOfIndex(0),
600 rte_pos,
601 rte_los,
602 ppath_lmax,
603 ppath_lraytrace,
604 false,
605 verbosity);
606 bending_angle1d(alpha2, ppt);
607 alpha2 *= Conversion::deg2rad(1);
608 a2 = ppt.constant;
609 //
610 rte_los[0] += 2 * dza;
611 adjust_los(rte_los, atmosphere_dim);
612 ppath_calc(ws,
613 ppt,
614 ppath_step_agenda,
615 atmosphere_dim,
616 p_grid,
617 lat_grid,
618 lon_grid,
619 z_field,
620 f_grid,
621 refellipsoid,
622 z_surface,
623 0,
624 ArrayOfIndex(0),
625 rte_pos,
626 rte_los,
627 ppath_lmax,
628 ppath_lraytrace,
629 false,
630 verbosity);
631 // This path can hit the surface. And we need to check if ppt is OK.
632 // (remember this function only deals with sat-to-sat links and OK
633 // background here is be space)
634 // Otherwise use the centre ray as the second one.
635 if (ppath_what_background(ppt) == 1) {
636 bending_angle1d(alpha1, ppt);
637 alpha1 *= Conversion::deg2rad(1);
638 a1 = ppt.constant;
639 dada = (alpha2 - alpha1) / (a2 - a1);
640 } else {
641 dada = (alpha2 - alpha0) / (a2 - a0);
642 }
643
644 // Zenith loss term (Eq 18 in Kursinski et al.)
645 const Numeric zlt = 1 / (1 - dada * lf);
646
647 // Total defocusing loss
648 dlf = zlt * alt;
649}
650
651Numeric dotprod_with_los(const ConstVectorView& los,
652 const Numeric& u,
653 const Numeric& v,
654 const Numeric& w,
655 const Index& atmosphere_dim) {
656 // Strength of field
657 const Numeric f = sqrt(u * u + v * v + w * w);
658
659 // Zenith and azimuth angle for field (in radians)
660 const Numeric za_f = acos(w / f);
661 const Numeric aa_f = atan2(u, v);
662
663 // Zenith and azimuth angle for photon direction (in radians)
664 Vector los_p;
665 mirror_los(los_p, los, atmosphere_dim);
666 const Numeric za_p = Conversion::deg2rad(1) * los_p[0];
667 const Numeric aa_p = Conversion::deg2rad(1) * los_p[1];
668
669 return f * (cos(za_f) * cos(za_p) + sin(za_f) * sin(za_p) * cos(aa_f - aa_p));
670}
671
673 Matrix& iy,
674 const Index& cloudbox_on,
675 const Vector& f_grid,
676 const EnergyLevelMap& nlte_field,
677 const Vector& rte_pos,
678 const Vector& rte_los,
679 const Vector& rte_pos2,
680 const String& iy_unit,
681 const Agenda& iy_main_agenda) {
682 ArrayOfTensor3 diy_dx;
683 ArrayOfMatrix iy_aux;
684 Ppath ppath;
685 Vector geo_pos;
686 Tensor3 iy_transmittance(0, 0, 0);
687 const Index iy_agenda_call1 = 1;
688 const ArrayOfString iy_aux_vars(0);
689 const Index iy_id = 0;
690 const Index jacobian_do = 0;
691
693 iy,
694 iy_aux,
695 ppath,
696 diy_dx,
697 geo_pos,
698 iy_agenda_call1,
699 iy_transmittance,
700 iy_aux_vars,
701 iy_id,
702 iy_unit,
703 cloudbox_on,
704 jacobian_do,
705 f_grid,
706 nlte_field,
707 rte_pos,
708 rte_los,
709 rte_pos2,
710 iy_main_agenda);
711}
712
714 Matrix& iy,
715 ArrayOfTensor3& diy_dx,
716 const Tensor3& iy_transmittance,
717 const Index& iy_id,
718 const Index& jacobian_do,
719 const ArrayOfRetrievalQuantity& jacobian_quantities,
720 const Ppath& ppath,
721 const Vector& rte_pos2,
722 const Index& atmosphere_dim,
723 const EnergyLevelMap& nlte_field,
724 const Index& cloudbox_on,
725 const Index& stokes_dim,
726 const Vector& f_grid,
727 const String& iy_unit,
728 const Tensor3& surface_props_data,
729 const Agenda& iy_main_agenda,
730 const Agenda& iy_space_agenda,
731 const Agenda& iy_surface_agenda,
732 const Agenda& iy_cloudbox_agenda,
733 const Index& iy_agenda_call1,
734 const Verbosity& verbosity) {
736
737 // Some sizes
738 const Index nf = f_grid.nelem();
739 const Index np = ppath.np;
740
741 // Set rtp_pos and rtp_los to match the last point in ppath.
742 //
743 // Note that the Ppath positions (ppath.pos) for 1D have one column more
744 // than expected by most functions. Only the first atmosphere_dim values
745 // shall be copied.
746 //
747 Vector rtp_pos, rtp_los;
748 rtp_pos.resize(atmosphere_dim);
749 rtp_pos = ppath.pos(np - 1, Range(0, atmosphere_dim));
750 rtp_los.resize(ppath.los.ncols());
751 rtp_los = ppath.los(np - 1, joker);
752
753 out3 << "Radiative background: " << ppath.background << "\n";
754
755 // Handle the different background cases
756 //
757 String agenda_name;
758 //
759 switch (ppath_what_background(ppath)) {
760 case 1: //--- Space ----------------------------------------------------
761 {
762 agenda_name = "iy_space_agenda";
763 chk_not_empty(agenda_name, iy_space_agenda);
764 iy_space_agendaExecute(ws, iy, f_grid, rtp_pos, rtp_los, iy_space_agenda);
765 } break;
766
767 case 2: //--- The surface -----------------------------------------------
768 {
769 agenda_name = "iy_surface_agenda";
770 chk_not_empty(agenda_name, iy_surface_agenda);
771 //
772 const Index los_id = iy_id % (Index)1000;
773 Index iy_id_new = iy_id + (Index)9 * los_id;
774 //
775 // Surface jacobian stuff:
776 ArrayOfString dsurface_names(0);
777 if (jacobian_do && iy_agenda_call1) {
778 for (Index i = 0; i < jacobian_quantities.nelem(); i++) {
779 if (jacobian_quantities[i] == Jacobian::Special::SurfaceString) {
780 dsurface_names.push_back(jacobian_quantities[i].Subtag());
781 }
782 }
783 }
784 ArrayOfTensor4 dsurface_rmatrix_dx(dsurface_names.nelem());
785 ArrayOfMatrix dsurface_emission_dx(dsurface_names.nelem());
786 //
788 iy,
789 diy_dx,
790 dsurface_rmatrix_dx,
791 dsurface_emission_dx,
792 iy_unit,
793 iy_transmittance,
794 iy_id_new,
795 cloudbox_on,
796 jacobian_do,
797 iy_main_agenda,
798 f_grid,
799 nlte_field,
800 rtp_pos,
801 rtp_los,
802 rte_pos2,
803 surface_props_data,
804 dsurface_names,
805 iy_surface_agenda);
806 } break;
807
808 case 3: //--- Cloudbox boundary or interior ------------------------------
809 case 4: {
810 agenda_name = "iy_cloudbox_agenda";
811 chk_not_empty(agenda_name, iy_cloudbox_agenda);
813 ws, iy, f_grid, rtp_pos, rtp_los, iy_cloudbox_agenda);
814 } break;
815
816 default: //--- ????? ----------------------------------------------------
817 // Are we here, the coding is wrong somewhere
818 ARTS_ASSERT(false);
819 }
820
821 ARTS_USER_ERROR_IF (iy.ncols() != stokes_dim || iy.nrows() != nf,
822 "The size of *iy* returned from *", agenda_name, "* is\n"
823 "not correct:\n"
824 " expected size = [", nf, ",", stokes_dim, "]\n"
825 " size of iy = [", iy.nrows(), ",", iy.ncols(), "]\n")
826}
827
828void get_ppath_atmvars(Vector& ppath_p,
829 Vector& ppath_t,
830 EnergyLevelMap& ppath_nlte,
831 Matrix& ppath_vmr,
832 Matrix& ppath_wind,
833 Matrix& ppath_mag,
834 const Ppath& ppath,
835 const Index& atmosphere_dim,
836 const ConstVectorView& p_grid,
837 const ConstTensor3View& t_field,
838 const EnergyLevelMap& nlte_field,
839 const ConstTensor4View& vmr_field,
840 const ConstTensor3View& wind_u_field,
841 const ConstTensor3View& wind_v_field,
842 const ConstTensor3View& wind_w_field,
843 const ConstTensor3View& mag_u_field,
844 const ConstTensor3View& mag_v_field,
845 const ConstTensor3View& mag_w_field) {
846 const Index np = ppath.np;
847 // Pressure:
848 ppath_p.resize(np);
849 Matrix itw_p(np, 2);
850 interpweights(itw_p, ppath.gp_p);
851 itw2p(ppath_p, p_grid, ppath.gp_p, itw_p);
852
853 // Temperature:
854 ppath_t.resize(np);
855 Matrix itw_field;
857 itw_field, atmosphere_dim, ppath.gp_p, ppath.gp_lat, ppath.gp_lon);
859 atmosphere_dim,
860 t_field,
861 ppath.gp_p,
862 ppath.gp_lat,
863 ppath.gp_lon,
864 itw_field);
865
866 // VMR fields:
867 const Index ns = vmr_field.nbooks();
868 ppath_vmr.resize(ns, np);
869 for (Index is = 0; is < ns; is++) {
870 interp_atmfield_by_itw(ppath_vmr(is, joker),
871 atmosphere_dim,
872 vmr_field(is, joker, joker, joker),
873 ppath.gp_p,
874 ppath.gp_lat,
875 ppath.gp_lon,
876 itw_field);
877 }
878
879 // NLTE temperatures
880 ppath_nlte = nlte_field.InterpToGridPos(atmosphere_dim, ppath.gp_p, ppath.gp_lat, ppath.gp_lon);
881
882 // Winds:
883 ppath_wind.resize(3, np);
884 ppath_wind = 0;
885 //
886 if (wind_u_field.npages() > 0) {
887 interp_atmfield_by_itw(ppath_wind(0, joker),
888 atmosphere_dim,
889 wind_u_field,
890 ppath.gp_p,
891 ppath.gp_lat,
892 ppath.gp_lon,
893 itw_field);
894 }
895 if (wind_v_field.npages() > 0) {
896 interp_atmfield_by_itw(ppath_wind(1, joker),
897 atmosphere_dim,
898 wind_v_field,
899 ppath.gp_p,
900 ppath.gp_lat,
901 ppath.gp_lon,
902 itw_field);
903 }
904 if (wind_w_field.npages() > 0) {
905 interp_atmfield_by_itw(ppath_wind(2, joker),
906 atmosphere_dim,
907 wind_w_field,
908 ppath.gp_p,
909 ppath.gp_lat,
910 ppath.gp_lon,
911 itw_field);
912 }
913
914 // Magnetic field:
915 ppath_mag.resize(3, np);
916 ppath_mag = 0;
917 //
918 if (mag_u_field.npages() > 0) {
919 interp_atmfield_by_itw(ppath_mag(0, joker),
920 atmosphere_dim,
921 mag_u_field,
922 ppath.gp_p,
923 ppath.gp_lat,
924 ppath.gp_lon,
925 itw_field);
926 }
927 if (mag_v_field.npages() > 0) {
928 interp_atmfield_by_itw(ppath_mag(1, joker),
929 atmosphere_dim,
930 mag_v_field,
931 ppath.gp_p,
932 ppath.gp_lat,
933 ppath.gp_lon,
934 itw_field);
935 }
936 if (mag_w_field.npages() > 0) {
937 interp_atmfield_by_itw(ppath_mag(2, joker),
938 atmosphere_dim,
939 mag_w_field,
940 ppath.gp_p,
941 ppath.gp_lat,
942 ppath.gp_lon,
943 itw_field);
944 }
945}
946
948 Matrix& ppath_pnd,
949 ArrayOfMatrix& ppath_dpnd_dx,
950 const Ppath& ppath,
951 const Index& atmosphere_dim,
952 const ArrayOfIndex& cloudbox_limits,
953 const Tensor4& pnd_field,
954 const ArrayOfTensor4& dpnd_field_dx) {
955 const Index np = ppath.np;
956
957 // Pnd along the ppath
958 ppath_pnd.resize(pnd_field.nbooks(), np);
959 ppath_pnd = 0;
960 ppath_dpnd_dx.resize(dpnd_field_dx.nelem());
961
962 bool any_dpnd = false;
963 for (Index iq = 0; iq < dpnd_field_dx.nelem(); iq++) {
964 if (dpnd_field_dx[iq].empty()) {
965 ppath_dpnd_dx[iq].resize(0, 0);
966 } else {
967 any_dpnd = true;
968 ppath_dpnd_dx[iq].resize(pnd_field.nbooks(), np);
969 }
970 }
971
972 // A variable that can map from ppath to particle containers.
973 // If outside cloudbox or all (d)pnd=0, this variable holds -1.
974 clear2cloudy.resize(np);
975
976 // Determine ppath_pnd and ppath_dpnd_dx
977 Index nin = 0;
978 for (Index ip = 0; ip < np; ip++) // PPath point
979 {
980 Matrix itw(1, Index(pow(2.0, Numeric(atmosphere_dim))));
981
982 ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
983 GridPos gp_lat, gp_lon;
984 if (atmosphere_dim >= 2) {
985 gridpos_copy(gp_lat, ppath.gp_lat[ip]);
986 }
987 if (atmosphere_dim == 3) {
988 gridpos_copy(gp_lon, ppath.gp_lon[ip]);
989 }
990
991 if (is_gp_inside_cloudbox(ppath.gp_p[ip],
992 gp_lat,
993 gp_lon,
994 cloudbox_limits,
995 true,
996 atmosphere_dim)) {
997 interp_cloudfield_gp2itw(itw(0, joker),
998 gpc_p[0],
999 gpc_lat[0],
1000 gpc_lon[0],
1001 ppath.gp_p[ip],
1002 gp_lat,
1003 gp_lon,
1004 atmosphere_dim,
1005 cloudbox_limits);
1006 for (Index i = 0; i < pnd_field.nbooks(); i++) {
1007 interp_atmfield_by_itw(ExhaustiveVectorView{ppath_pnd(i, ip)},
1008 atmosphere_dim,
1009 pnd_field(i, joker, joker, joker),
1010 gpc_p,
1011 gpc_lat,
1012 gpc_lon,
1013 itw);
1014 }
1015 bool any_ppath_dpnd = false;
1016 if (any_dpnd) {
1017 for (Index iq = 0; iq < dpnd_field_dx.nelem();
1018 iq++) // Jacobian parameter
1019 {
1020 if (!dpnd_field_dx[iq].empty()) {
1021 for (Index i = 0; i < pnd_field.nbooks();
1022 i++) // Scattering element
1023 {
1024 interp_atmfield_by_itw(ExhaustiveVectorView{ppath_dpnd_dx[iq](i, ip)},
1025 atmosphere_dim,
1026 dpnd_field_dx[iq](i, joker, joker, joker),
1027 gpc_p,
1028 gpc_lat,
1029 gpc_lon,
1030 itw);
1031 }
1032 if (max(ppath_dpnd_dx[iq](joker, ip)) > 0. ||
1033 min(ppath_dpnd_dx[iq](joker, ip)) < 0.)
1034 any_ppath_dpnd = true;
1035 }
1036 }
1037 }
1038 if (max(ppath_pnd(joker, ip)) > 0. || min(ppath_pnd(joker, ip)) < 0. ||
1039 any_ppath_dpnd) {
1040 clear2cloudy[ip] = nin;
1041 nin++;
1042 } else {
1043 clear2cloudy[ip] = -1;
1044 }
1045 } else {
1046 clear2cloudy[ip] = -1;
1047 }
1048 }
1049}
1050
1051void get_ppath_f(Matrix& ppath_f,
1052 const Ppath& ppath,
1053 const ConstVectorView& f_grid,
1054 const Index& atmosphere_dim,
1055 const Numeric& rte_alonglos_v,
1056 const ConstMatrixView& ppath_wind) {
1057 // Sizes
1058 const Index nf = f_grid.nelem();
1059 const Index np = ppath.np;
1060
1061 ppath_f.resize(nf, np);
1062
1063 // Doppler relevant velocity
1064 //
1065 for (Index ip = 0; ip < np; ip++) {
1066 // Start by adding rte_alonglos_v (most likely sensor effects)
1067 Numeric v_doppler = rte_alonglos_v;
1068
1069 // Include wind
1070 if (ppath_wind(1, ip) != 0 || ppath_wind(0, ip) != 0 ||
1071 ppath_wind(2, ip) != 0) {
1072 // The dot product below is valid for the photon direction. Winds
1073 // along this direction gives a positive contribution.
1074 v_doppler += dotprod_with_los(ppath.los(ip, joker),
1075 ppath_wind(0, ip),
1076 ppath_wind(1, ip),
1077 ppath_wind(2, ip),
1078 atmosphere_dim);
1079 }
1080
1081 // Determine frequency grid
1082 if (v_doppler == 0) {
1083 ppath_f(joker, ip) = f_grid;
1084 } else {
1085 // Positive v_doppler means that sensor measures lower rest
1086 // frequencies
1087 const Numeric a = 1 - v_doppler / SPEED_OF_LIGHT;
1088 for (Index iv = 0; iv < nf; iv++) {
1089 ppath_f(iv, ip) = a * f_grid[iv];
1090 }
1091 }
1092 }
1093}
1094
1095Range get_rowindex_for_mblock(const Sparse& sensor_response,
1096 const Index& mblock_index) {
1097 const Index n1y = sensor_response.nrows();
1098 return Range(n1y * mblock_index, n1y);
1099}
1100
1102 VectorView dB_dT,
1103 const ConstVectorView& ppath_f_grid,
1104 const Numeric& ppath_temperature,
1105 const bool& do_temperature_derivative) {
1106 const Index nf = ppath_f_grid.nelem();
1107
1108 for (Index i = 0; i < nf; i++)
1109 B[i] = planck(ppath_f_grid[i], ppath_temperature);
1110
1111 if (do_temperature_derivative)
1112 for (Index i = 0; i < nf; i++)
1113 dB_dT[i] = dplanck_dt(ppath_f_grid[i], ppath_temperature);
1114}
1115
1117 Workspace& ws,
1118 PropagationMatrix& K,
1119 StokesVector& S,
1120 Index& lte,
1121 ArrayOfPropagationMatrix& dK_dx,
1122 ArrayOfStokesVector& dS_dx,
1123 const Agenda& propmat_clearsky_agenda,
1124 const ArrayOfRetrievalQuantity& jacobian_quantities,
1125 const Vector& ppath_f_grid,
1126 const Vector& ppath_magnetic_field,
1127 const Vector& ppath_line_of_sight,
1128 const EnergyLevelMap& ppath_nlte,
1129 const Vector& ppath_vmrs,
1130 const Numeric& ppath_temperature,
1131 const Numeric& ppath_pressure,
1132 const bool& jacobian_do) {
1133 // Perform the propagation matrix computations
1135 K,
1136 S,
1137 dK_dx,
1138 dS_dx,
1139 jacobian_do ? jacobian_quantities : ArrayOfRetrievalQuantity(0),
1140 {},
1141 ppath_f_grid,
1142 ppath_magnetic_field,
1143 ppath_line_of_sight,
1144 ppath_pressure,
1145 ppath_temperature,
1146 ppath_nlte,
1147 ppath_vmrs,
1148 propmat_clearsky_agenda);
1149
1150 // If there are no NLTE values, then set the LTE flag as true
1151 lte = S.allZeroes(); // FIXME: Should be nlte_do?
1152}
1153
1154Vector get_stepwise_f_partials(const ConstVectorView& line_of_sight,
1155 const ConstVectorView& f_grid,
1156 const Jacobian::Atm wind_type,
1157 const Index& atmosphere_dim) {
1158 // Doppler relevant velocity
1159 Numeric dv_doppler_dx = 0.0;
1160
1161 Vector deriv(f_grid);
1162
1163 switch (wind_type) {
1164 case Jacobian::Atm::WindMagnitude:
1165 dv_doppler_dx = 1.0;
1166 break;
1167 case Jacobian::Atm::WindU:
1168 dv_doppler_dx =
1169 (dotprod_with_los(line_of_sight, 1, 0, 0, atmosphere_dim));
1170 break;
1171 case Jacobian::Atm::WindV:
1172 dv_doppler_dx =
1173 (dotprod_with_los(line_of_sight, 0, 1, 0, atmosphere_dim));
1174 break;
1175 case Jacobian::Atm::WindW:
1176 dv_doppler_dx =
1177 (dotprod_with_los(line_of_sight, 0, 0, 1, atmosphere_dim));
1178 break;
1179 default:
1180 ARTS_ASSERT(false, "Not allowed to call this function without a wind parameter as wind_type");
1181 break;
1182 }
1183
1184 deriv *= - dv_doppler_dx / Constant::c;
1185 return deriv;
1186}
1187
1189 StokesVector& ap,
1190 PropagationMatrix& Kp,
1191 ArrayOfStokesVector& dap_dx,
1192 ArrayOfPropagationMatrix& dKp_dx,
1193 const ArrayOfRetrievalQuantity& jacobian_quantities,
1194 const ConstMatrixView& ppath_1p_pnd, // the ppath_pnd at this ppath point
1195 const ArrayOfMatrix&
1196 ppath_dpnd_dx, // the full ppath_dpnd_dx, ie all ppath points
1197 const Index ppath_1p_id,
1198 const ArrayOfArrayOfSingleScatteringData& scat_data,
1199 const ConstVectorView& ppath_line_of_sight,
1200 const ConstVectorView& ppath_temperature,
1201 const Index& atmosphere_dim,
1202 const bool& jacobian_do) {
1203 const Index nf = Kp.NumberOfFrequencies(), stokes_dim = Kp.StokesDimensions();
1204
1205 //StokesVector da_aux(nf, stokes_dim);
1206 //PropagationMatrix dK_aux(nf, stokes_dim);
1207
1209
1210 // Direction of outgoing scattered radiation (which is reversed to
1211 // LOS). Only used for extracting scattering properties.
1212 Vector dir;
1213 mirror_los(dir, ppath_line_of_sight, atmosphere_dim);
1214 Matrix dir_array(1, 2, 0.);
1215 dir_array(0, joker) = dir;
1216
1217 ArrayOfArrayOfTensor5 ext_mat_Nse;
1218 ArrayOfArrayOfTensor4 abs_vec_Nse;
1219 ArrayOfArrayOfIndex ptypes_Nse;
1220 Matrix t_ok;
1221 ArrayOfTensor5 ext_mat_ssbulk;
1222 ArrayOfTensor4 abs_vec_ssbulk;
1223 ArrayOfIndex ptype_ssbulk;
1224 Tensor5 ext_mat_bulk;
1225 Tensor4 abs_vec_bulk;
1226 Index ptype_bulk;
1227
1228 // get per-scat-elem data here. and fold with pnd.
1229 // keep per-scat-elem data to fold with dpnd_dx further down in
1230 // analyt-jac-loop.
1231 opt_prop_NScatElems(ext_mat_Nse,
1232 abs_vec_Nse,
1233 ptypes_Nse,
1234 t_ok,
1235 scat_data,
1236 stokes_dim,
1237 Vector{ppath_temperature},
1238 dir_array,
1239 -1);
1240
1241 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
1242 abs_vec_ssbulk,
1243 ptype_ssbulk,
1244 ext_mat_Nse,
1245 abs_vec_Nse,
1246 ptypes_Nse,
1247 ppath_1p_pnd,
1248 t_ok);
1249 opt_prop_Bulk(ext_mat_bulk,
1250 abs_vec_bulk,
1251 ptype_bulk,
1252 ext_mat_ssbulk,
1253 abs_vec_ssbulk,
1254 ptype_ssbulk);
1255
1256 const Index nf_ssd = abs_vec_bulk.nbooks(); // number of freqs in extracted
1257 // optprops. if 1, we need to
1258 // duplicate the ext/abs output.
1259
1260 for (Index iv = 0; iv < nf; iv++) {
1261 if (nf_ssd > 1) {
1262 ap.SetAtPosition(abs_vec_bulk(iv, 0, 0, joker), iv);
1263 Kp.SetAtPosition(ext_mat_bulk(iv, 0, 0, joker, joker), iv);
1264 } else {
1265 ap.SetAtPosition(abs_vec_bulk(0, 0, 0, joker), iv);
1266 Kp.SetAtPosition(ext_mat_bulk(0, 0, 0, joker, joker), iv);
1267 }
1268 }
1269
1270 if (jacobian_do)
1272 if (ppath_dpnd_dx[iq].empty()) {
1273 dap_dx[iq].SetZero();
1274 dKp_dx[iq].SetZero();
1275 } else {
1276 // check, whether we have any non-zero ppath_dpnd_dx in this
1277 // pnd-affecting x? might speed up things a little bit.
1278 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
1279 abs_vec_ssbulk,
1280 ptype_ssbulk,
1281 ext_mat_Nse,
1282 abs_vec_Nse,
1283 ptypes_Nse,
1284 ppath_dpnd_dx[iq](joker, Range(ppath_1p_id, 1)),
1285 t_ok);
1286 opt_prop_Bulk(ext_mat_bulk,
1287 abs_vec_bulk,
1288 ptype_bulk,
1289 ext_mat_ssbulk,
1290 abs_vec_ssbulk,
1291 ptype_ssbulk);
1292 for (Index iv = 0; iv < nf; iv++) {
1293 if (nf_ssd > 1) {
1294 dap_dx[iq].SetAtPosition(abs_vec_bulk(iv, 0, 0, joker), iv);
1295 dKp_dx[iq].SetAtPosition(ext_mat_bulk(iv, 0, 0, joker, joker),
1296 iv);
1297 } else {
1298 dap_dx[iq].SetAtPosition(abs_vec_bulk(0, 0, 0, joker), iv);
1299 dKp_dx[iq].SetAtPosition(ext_mat_bulk(0, 0, 0, joker, joker), iv);
1300 }
1301 }
1302 })
1303}
1304
1306 StokesVector& Sp,
1307 ArrayOfStokesVector& dSp_dx,
1308 const ArrayOfRetrievalQuantity& jacobian_quantities,
1309 const ConstVectorView& ppath_1p_pnd, // the ppath_pnd at this ppath point
1310 const ArrayOfMatrix&
1311 ppath_dpnd_dx, // the full ppath_dpnd_dx, ie all ppath points
1312 const Index ppath_1p_id,
1313 const ArrayOfArrayOfSingleScatteringData& scat_data,
1314 const ConstTensor7View& cloudbox_field,
1315 const ConstVectorView& za_grid,
1316 const ConstVectorView& aa_grid,
1317 const ConstMatrixView& ppath_line_of_sight,
1318 const GridPos& ppath_pressure,
1319 const Vector& temperature,
1320 const Index& atmosphere_dim,
1321 const bool& jacobian_do,
1322 const Index& t_interp_order) {
1323 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
1324 "This function handles so far only 1D atmospheres.");
1325
1326 const Index nf = Sp.NumberOfFrequencies();
1327 const Index stokes_dim = Sp.StokesDimensions();
1328 const Index ne = ppath_1p_pnd.nelem();
1329 ARTS_ASSERT(TotalNumberOfElements(scat_data) == ne);
1330 const Index nza = za_grid.nelem();
1331 const Index naa = aa_grid.nelem();
1332 const Index nq = jacobian_do ? jacobian_quantities.nelem() : 0;
1333
1334 // interpolate incident field to this ppath point (no need to do this
1335 // separately per scatelem)
1336 GridPos gp_p;
1337 gridpos_copy(gp_p, ppath_pressure);
1338 Vector itw_p(2);
1339 interpweights(itw_p, gp_p);
1340 Tensor3 inc_field(nf, nza, stokes_dim, 0.);
1341 for (Index iv = 0; iv < nf; iv++) {
1342 for (Index iza = 0; iza < nza; iza++) {
1343 for (Index i = 0; i < stokes_dim; i++) {
1344 inc_field(iv, iza, i) =
1345 interp(itw_p, cloudbox_field(iv, joker, 0, 0, iza, 0, i), gp_p);
1346 }
1347 }
1348 }
1349
1350 // create matrix of incident directions (flat representation of the
1351 // za_grid * aa_grid matrix)
1352 Matrix idir(nza * naa, 2);
1353 Index ia = 0;
1354 for (Index iza = 0; iza < nza; iza++) {
1355 for (Index iaa = 0; iaa < naa; iaa++) {
1356 idir(ia, 0) = za_grid[iza];
1357 idir(ia, 1) = aa_grid[iaa];
1358 ia++;
1359 }
1360 }
1361
1362 // setting prop (aka scattered) direction
1363 Matrix pdir(1, 2);
1364 //if( ppath_line_of_sight.ncols()==2 )
1365 // pdir(0,joker) = ppath_line_of_sight;
1366 //else // 1D case only (currently the only handled case). azimuth not defined.
1367 {
1368 pdir(0, 0) = ppath_line_of_sight(0, 0);
1369 pdir(0, 1) = 0.;
1370 }
1371
1372 // some further variables needed for pha_mat extraction
1373 Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
1374 Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
1375 Tensor6 pha_mat_1se(nf_ssd, 1, 1, nza * naa, stokes_dim, stokes_dim);
1376 Vector t_ok(1);
1377 Index ptype;
1378 Tensor3 scat_source_1se(ne, nf, stokes_dim, 0.);
1379
1380 Index ise_flat = 0;
1381 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++) {
1382 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
1383 // determine whether we have some valid pnd for this
1384 // scatelem (in pnd or dpnd)
1385 Index val_pnd = 0;
1386 if (ppath_1p_pnd[ise_flat] != 0) {
1387 val_pnd = 1;
1388 } else if (jacobian_do) {
1389 for (Index iq = 0; (!val_pnd) && (iq < nq); iq++) {
1390 if ((not(jacobian_quantities[iq] == Jacobian::Type::Sensor) and not(jacobian_quantities[iq] == Jacobian::Special::SurfaceString)) &&
1391 !ppath_dpnd_dx[iq].empty() &&
1392 ppath_dpnd_dx[iq](ise_flat, ppath_1p_id) != 0) {
1393 val_pnd = 1;
1394 }
1395 }
1396 }
1397
1398 if (val_pnd) {
1399 pha_mat_1ScatElem(pha_mat_1se,
1400 ptype,
1401 t_ok,
1402 scat_data[i_ss][i_se],
1403 temperature,
1404 pdir,
1405 idir,
1406 0,
1407 t_interp_order);
1408 ARTS_USER_ERROR_IF (t_ok[0] == 0,
1409 "Interpolation error for (flat-array) scattering "
1410 "element #", ise_flat, "\n"
1411 "at location/temperature point #", ppath_1p_id, "\n")
1412
1413 Index this_iv = 0;
1414 for (Index iv = 0; iv < nf; iv++) {
1415 if (!duplicate_freqs) {
1416 this_iv = iv;
1417 }
1418 Tensor3 product_fields(nza, naa, stokes_dim, 0.);
1419
1420 ia = 0;
1421 for (Index iza = 0; iza < nza; iza++) {
1422 for (Index iaa = 0; iaa < naa; iaa++) {
1423 for (Index i = 0; i < stokes_dim; i++) {
1424 for (Index j = 0; j < stokes_dim; j++) {
1425 product_fields(iza, iaa, i) +=
1426 pha_mat_1se(this_iv, 0, 0, ia, i, j) *
1427 inc_field(iv, iza, j);
1428 }
1429 }
1430 ia++;
1431 }
1432 }
1433
1434 for (Index i = 0; i < stokes_dim; i++) {
1435 scat_source_1se(ise_flat, iv, i) = AngIntegrate_trapezoid(
1436 product_fields(joker, joker, i), za_grid, aa_grid);
1437 }
1438 } // for iv
1439 } // if val_pnd
1440
1441 ise_flat++;
1442
1443 } // for i_se
1444 } // for i_ss
1445
1446 for (Index iv = 0; iv < nf; iv++) {
1447 Vector scat_source(stokes_dim, 0.);
1448 for (ise_flat = 0; ise_flat < ne; ise_flat++) {
1449 for (Index i = 0; i < stokes_dim; i++) {
1450 scat_source[i] +=
1451 scat_source_1se(ise_flat, iv, i) * ppath_1p_pnd[ise_flat];
1452 }
1453 }
1454
1455 Sp.SetAtPosition(scat_source, iv);
1456
1457 if (jacobian_do) {
1459 if (ppath_dpnd_dx[iq].empty()) { dSp_dx[iq].SetZero(); } else {
1460 scat_source = 0.;
1461 for (ise_flat = 0; ise_flat < ne; ise_flat++) {
1462 for (Index i = 0; i < stokes_dim; i++) {
1463 scat_source[i] += scat_source_1se(ise_flat, iv, i) *
1464 ppath_dpnd_dx[iq](ise_flat, ppath_1p_id);
1465 dSp_dx[iq].SetAtPosition(scat_source, iv);
1466 }
1467 }
1468 })
1469 }
1470 } // for iv
1471}
1472
1473void iyb_calc_body(bool& failed,
1474 String& fail_msg,
1475 ArrayOfArrayOfMatrix& iy_aux_array,
1476 Workspace& ws,
1477 Ppath& ppath,
1478 Vector& iyb,
1479 ArrayOfMatrix& diyb_dx,
1480 Vector& geo_pos,
1481 const Index& mblock_index,
1482 const Index& atmosphere_dim,
1483 const EnergyLevelMap& nlte_field,
1484 const Index& cloudbox_on,
1485 const Index& stokes_dim,
1486 const Matrix& sensor_pos,
1487 const Matrix& sensor_los,
1488 const Matrix& transmitter_pos,
1489 const Matrix& mblock_dlos,
1490 const String& iy_unit,
1491 const Agenda& iy_main_agenda,
1492 const Index& j_analytical_do,
1493 const ArrayOfRetrievalQuantity& jacobian_quantities,
1494 const ArrayOfArrayOfIndex& jacobian_indices,
1495 const Vector& f_grid,
1496 const ArrayOfString& iy_aux_vars,
1497 const Index& ilos,
1498 const Index& nf) {
1499 // The try block here is necessary to correctly handle
1500 // exceptions inside the parallel region.
1501 try {
1502 //--- LOS of interest
1503 //
1504 Vector los(sensor_los.ncols());
1505 //
1506 los = sensor_los(mblock_index, joker);
1507 if (mblock_dlos.ncols() == 1) {
1508 los[0] += mblock_dlos(ilos, 0);
1509 adjust_los(los, atmosphere_dim);
1510 } else {
1511 add_za_aa(los[0],
1512 los[1],
1513 los[0],
1514 los[1],
1515 mblock_dlos(ilos, 0),
1516 mblock_dlos(ilos, 1));
1517 }
1518
1519 //--- rtp_pos 1 and 2
1520 //
1521 Vector rtp_pos, rtp_pos2(0);
1522 //
1523 rtp_pos = sensor_pos(mblock_index, joker);
1524 if (!transmitter_pos.empty()) {
1525 rtp_pos2 = transmitter_pos(mblock_index, joker);
1526 }
1527
1528 // Calculate iy and associated variables
1529 //
1530 Matrix iy;
1531 ArrayOfTensor3 diy_dx;
1532 Tensor3 iy_transmittance(0, 0, 0);
1533 const Index iy_agenda_call1 = 1;
1534 const Index iy_id =
1535 (Index)1e6 * (mblock_index + 1) + (Index)1e3 * (ilos + 1);
1536 //
1538 iy,
1539 iy_aux_array[ilos],
1540 ppath,
1541 diy_dx,
1542 geo_pos,
1543 iy_agenda_call1,
1544 iy_transmittance,
1545 iy_aux_vars,
1546 iy_id,
1547 iy_unit,
1548 cloudbox_on,
1549 j_analytical_do,
1550 f_grid,
1551 nlte_field,
1552 rtp_pos,
1553 los,
1554 rtp_pos2,
1555 iy_main_agenda);
1556
1557 // Start row in iyb etc. for present LOS
1558 //
1559 const Index row0 = ilos * nf * stokes_dim;
1560
1561 // Jacobian part
1562 //
1563 if (j_analytical_do) {
1565 for (Index ip = 0;
1566 ip < jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1;
1567 ip++) {
1568 for (Index is = 0; is < stokes_dim; is++) {
1569 diyb_dx[iq](Range(row0 + is, nf, stokes_dim), ip) =
1570 diy_dx[iq](ip, joker, is);
1571 }
1572 })
1573 }
1574
1575 // iy : copy to iyb
1576 for (Index is = 0; is < stokes_dim; is++) {
1577 iyb[Range(row0 + is, nf, stokes_dim)] = iy(joker, is);
1578 }
1579
1580 } // End try
1581
1582 catch (const std::exception& e) {
1583#pragma omp critical(iyb_calc_fail)
1584 {
1585 fail_msg = e.what();
1586 failed = true;
1587 }
1588 }
1589}
1590
1592 Vector& iyb,
1593 ArrayOfVector& iyb_aux,
1594 ArrayOfMatrix& diyb_dx,
1595 Matrix& geo_pos_matrix,
1596 const Index& mblock_index,
1597 const Index& atmosphere_dim,
1598 const EnergyLevelMap& nlte_field,
1599 const Index& cloudbox_on,
1600 const Index& stokes_dim,
1601 const Vector& f_grid,
1602 const Matrix& sensor_pos,
1603 const Matrix& sensor_los,
1604 const Matrix& transmitter_pos,
1605 const Matrix& mblock_dlos,
1606 const String& iy_unit,
1607 const Agenda& iy_main_agenda,
1608 const Index& j_analytical_do,
1609 const ArrayOfRetrievalQuantity& jacobian_quantities,
1610 const ArrayOfArrayOfIndex& jacobian_indices,
1611 const ArrayOfString& iy_aux_vars,
1612 const Verbosity& verbosity) {
1614
1615 // Sizes
1616 const Index nf = f_grid.nelem();
1617 const Index nlos = mblock_dlos.nrows();
1618 const Index niyb = nf * nlos * stokes_dim;
1619 // Set up size of containers for data of 1 measurement block.
1620 // (can not be made below due to parallalisation)
1621 iyb.resize(niyb);
1622 //
1623 if (j_analytical_do) {
1624 diyb_dx.resize(jacobian_indices.nelem());
1625 FOR_ANALYTICAL_JACOBIANS_DO2(diyb_dx[iq].resize(
1626 niyb, jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1);)
1627 } else {
1628 diyb_dx.resize(0);
1629 }
1630 // Assume empty geo_pos.
1631 geo_pos_matrix.resize(nlos, 5);
1632 geo_pos_matrix = NAN;
1633
1634 // For iy_aux we don't know the number of quantities, and we have to store
1635 // all outout
1636 ArrayOfArrayOfMatrix iy_aux_array(nlos);
1637
1638 String fail_msg;
1639 bool failed = false;
1640 if (nlos >= arts_omp_get_max_threads() || nlos * 10 >= nf) {
1641 out3 << " Parallelizing los loop (" << nlos << " iterations, " << nf
1642 << " frequencies)\n";
1643
1645
1646 // Start of actual calculations
1647#pragma omp parallel for if (!arts_omp_in_parallel()) firstprivate(wss)
1648 for (Index ilos = 0; ilos < nlos; ilos++) {
1649 // Skip remaining iterations if an error occurred
1650 if (failed) continue;
1651
1652 Ppath ppath;
1653 Vector geo_pos;
1654 iyb_calc_body(failed,
1655 fail_msg,
1656 iy_aux_array,
1657 wss,
1658 ppath,
1659 iyb,
1660 diyb_dx,
1661 geo_pos,
1662 mblock_index,
1663 atmosphere_dim,
1664 nlte_field,
1665 cloudbox_on,
1666 stokes_dim,
1667 sensor_pos,
1668 sensor_los,
1669 transmitter_pos,
1670 mblock_dlos,
1671 iy_unit,
1672 iy_main_agenda,
1673 j_analytical_do,
1674 jacobian_quantities,
1675 jacobian_indices,
1676 f_grid,
1677 iy_aux_vars,
1678 ilos,
1679 nf);
1680
1681 if (geo_pos.nelem()) geo_pos_matrix(ilos, joker) = geo_pos;
1682
1683 // Skip remaining iterations if an error occurred
1684 if (failed) continue;
1685 }
1686 } else {
1687 out3 << " Not parallelizing los loop (" << nlos << " iterations, " << nf
1688 << " frequencies)\n";
1689
1690 for (Index ilos = 0; ilos < nlos; ilos++) {
1691 // Skip remaining iterations if an error occurred
1692 if (failed) continue;
1693
1694 Ppath ppath;
1695 Vector geo_pos;
1696 iyb_calc_body(failed,
1697 fail_msg,
1698 iy_aux_array,
1699 ws,
1700 ppath,
1701 iyb,
1702 diyb_dx,
1703 geo_pos,
1704 mblock_index,
1705 atmosphere_dim,
1706 nlte_field,
1707 cloudbox_on,
1708 stokes_dim,
1709 sensor_pos,
1710 sensor_los,
1711 transmitter_pos,
1712 mblock_dlos,
1713 iy_unit,
1714 iy_main_agenda,
1715 j_analytical_do,
1716 jacobian_quantities,
1717 jacobian_indices,
1718 f_grid,
1719 iy_aux_vars,
1720 ilos,
1721 nf);
1722
1723 if (geo_pos.nelem()) geo_pos_matrix(ilos, joker) = geo_pos;
1724
1725 // Skip remaining iterations if an error occurred
1726 if (failed) continue;
1727 }
1728 }
1729
1730 ARTS_USER_ERROR_IF (failed,
1731 "Run-time error in function: iyb_calc\n", fail_msg);
1732
1733 // Compile iyb_aux
1734 //
1735 const Index nq = iy_aux_array[0].nelem();
1736 iyb_aux.resize(nq);
1737 //
1738 for (Index q = 0; q < nq; q++) {
1739 iyb_aux[q].resize(niyb);
1740 //
1741 for (Index ilos = 0; ilos < nlos; ilos++) {
1742 const Index row0 = ilos * nf * stokes_dim;
1743 for (Index iv = 0; iv < nf; iv++) {
1744 const Index row1 = row0 + iv * stokes_dim;
1745 const Index i1 = min(iv, iy_aux_array[ilos][q].nrows() - 1);
1746 for (Index is = 0; is < stokes_dim; is++) {
1747 Index i2 = min(is, iy_aux_array[ilos][q].ncols() - 1);
1748 iyb_aux[q][row1 + is] = iy_aux_array[ilos][q](i1, i2);
1749 }
1750 }
1751 }
1752 }
1753}
1754
1755void iy_transmittance_mult(Tensor3& iy_trans_total,
1756 const ConstTensor3View& iy_trans_old,
1757 const ConstTensor3View& iy_trans_new) {
1758 const Index nf = iy_trans_old.npages();
1759 const Index ns = iy_trans_old.ncols();
1760
1761 ARTS_ASSERT(ns == iy_trans_old.nrows());
1762 ARTS_ASSERT(nf == iy_trans_new.npages());
1763 ARTS_ASSERT(ns == iy_trans_new.nrows());
1764 ARTS_ASSERT(ns == iy_trans_new.ncols());
1765
1766 iy_trans_total.resize(nf, ns, ns);
1767
1768 for (Index iv = 0; iv < nf; iv++) {
1769 mult(iy_trans_total(iv, joker, joker),
1770 iy_trans_old(iv, joker, joker),
1771 iy_trans_new(iv, joker, joker));
1772 }
1773}
1774
1775void iy_transmittance_mult(Matrix& iy_new,
1776 const ConstTensor3View& iy_trans,
1777 const ConstMatrixView& iy_old) {
1778 const Index nf = iy_trans.npages();
1779 const Index ns = iy_trans.ncols();
1780
1781 ARTS_ASSERT(ns == iy_trans.nrows());
1782 ARTS_ASSERT(nf == iy_old.nrows());
1783 ARTS_ASSERT(ns == iy_old.ncols());
1784
1785 iy_new.resize(nf, ns);
1786
1787 for (Index iv = 0; iv < nf; iv++) {
1788 mult(iy_new(iv, joker), iy_trans(iv, joker, joker), iy_old(iv, joker));
1789 }
1790}
1791
1792void mirror_los(Vector& los_mirrored,
1793 const ConstVectorView& los,
1794 const Index& atmosphere_dim) {
1795 los_mirrored.resize(2);
1796 //
1797 if (atmosphere_dim == 1) {
1798 los_mirrored[0] = 180 - los[0];
1799 los_mirrored[1] = 180;
1800 } else if (atmosphere_dim == 2) {
1801 los_mirrored[0] = 180 - fabs(los[0]);
1802 if (los[0] >= 0) {
1803 los_mirrored[1] = 180;
1804 } else {
1805 los_mirrored[1] = 0;
1806 }
1807 } else if (atmosphere_dim == 3) {
1808 los_mirrored[0] = 180 - los[0];
1809 los_mirrored[1] = los[1] + 180;
1810 if (los_mirrored[1] > 180) {
1811 los_mirrored[1] -= 360;
1812 }
1813 }
1814}
1815
1817 const Index& stokes_dim,
1818 const Numeric& rotangle) {
1819 ARTS_ASSERT(stokes_dim > 1);
1820 ARTS_ASSERT(stokes_dim <= 4);
1821 ARTS_ASSERT(H.nrows() == stokes_dim);
1822 ARTS_ASSERT(H.ncols() == stokes_dim);
1823 ARTS_ASSERT(H(0, 1) == 0);
1824 ARTS_ASSERT(H(1, 0) == 0);
1825 //
1826 H.rw(0, 0) = 1;
1827 const Numeric a = Conversion::cosd(2 * rotangle);
1828 H.rw(1, 1) = a;
1829 if (stokes_dim > 2) {
1830 ARTS_ASSERT(H(2, 0) == 0);
1831 ARTS_ASSERT(H(0, 2) == 0);
1832
1833 const Numeric b = Conversion::sind(2 * rotangle);
1834 H.rw(1, 2) = b;
1835 H.rw(2, 1) = -b;
1836 H.rw(2, 2) = a;
1837 if (stokes_dim > 3) {
1838 // More values should be checked, but to save time we just ARTS_ASSERT one
1839 ARTS_ASSERT(H(2, 3) == 0);
1840 H.rw(3, 3) = 1;
1841 }
1842 }
1843}
1844
1845void mueller_modif2stokes(Matrix& Cs,
1846 const Index& stokes_dim) {
1847 ARTS_ASSERT(stokes_dim >= 1);
1848 ARTS_ASSERT(stokes_dim <= 4);
1849 //
1850 Cs.resize(stokes_dim, stokes_dim);
1851 Cs(0,0) = 1;
1852 if (stokes_dim > 1 ) {
1853 Cs(0,1) = Cs(1,0) = 1;
1854 Cs(1,1) = -1;
1855 if (stokes_dim > 2 ) {
1856 Cs(0,2) = Cs(1,2) = Cs(2,0) = Cs(2,1) = 0;
1857 Cs(2,2) = 1;
1858 if (stokes_dim > 3 ) {
1859 Cs(0,3) = Cs(1,3) = Cs(2,3) = Cs(3,0) = Cs(3,1) = Cs(3,2) = 0;
1860 Cs(3,3) = 1;
1861 }
1862 }
1863 }
1864}
1865
1866void mueller_rotation(Matrix& L,
1867 const Index& stokes_dim,
1868 const Numeric& rotangle) {
1869 ARTS_ASSERT(stokes_dim >= 1);
1870 ARTS_ASSERT(stokes_dim <= 4);
1871 //
1872 L.resize(stokes_dim, stokes_dim);
1873 L(0, 0) = 1;
1874 if (stokes_dim > 1 ) {
1875 const Numeric alpha = 2 * Conversion::deg2rad(1) * rotangle;
1876 const Numeric c2 = cos(alpha);
1877 L(0,1) = L(1,0) = 0;
1878 L(1,1) = c2;
1879 if (stokes_dim > 2 ) {
1880 const Numeric s2 = sin(alpha);
1881 L(0,2) = L(2,0) = 0;
1882 L(1,2) = s2;
1883 L(2,1) = -s2;
1884 L(2,2) = c2;
1885 if (stokes_dim > 3 ) {
1886 L(0,3) = L(1,3) = L(2,3) = L(3,0) = L(3,1) = L(3,2) = 0;
1887 L(3,3) = 1;
1888 }
1889 }
1890 }
1891}
1892
1893void mueller_stokes2modif(Matrix& Cm,
1894 const Index& stokes_dim) {
1895 ARTS_ASSERT(stokes_dim >= 1);
1896 ARTS_ASSERT(stokes_dim <= 4);
1897 //
1898 Cm.resize(stokes_dim, stokes_dim);
1899 Cm(0,0) = 0.5;
1900 if (stokes_dim > 1 ) {
1901 Cm(0,1) = Cm(1,0) = 0.5;
1902 Cm(1,1) = -0.5;
1903 if (stokes_dim > 2 ) {
1904 Cm(0,2) = Cm(1,2) = Cm(2,0) = Cm(2,1) = 0;
1905 Cm(2,2) = 1;
1906 if (stokes_dim > 3 ) {
1907 Cm(0,3) = Cm(1,3) = Cm(2,3) = Cm(3,0) = Cm(3,1) = Cm(3,2) = 0;
1908 Cm(3,3) = 1;
1909 }
1910 }
1911 }
1912}
1913
1914void pos2true_latlon(Numeric& lat,
1915 Numeric& lon,
1916 const Index& atmosphere_dim,
1917 const ConstVectorView& lat_grid,
1918 const ConstVectorView& lat_true,
1919 const ConstVectorView& lon_true,
1920 const ConstVectorView& pos) {
1921 ARTS_ASSERT(pos.nelem() == atmosphere_dim);
1922
1923 if (atmosphere_dim == 1) {
1924 ARTS_ASSERT(lat_true.nelem() == 1);
1925 ARTS_ASSERT(lon_true.nelem() == 1);
1926 //
1927 lat = lat_true[0];
1928 lon = lon_true[0];
1929 }
1930
1931 else if (atmosphere_dim == 2) {
1932 ARTS_ASSERT(lat_true.nelem() == lat_grid.nelem());
1933 ARTS_ASSERT(lon_true.nelem() == lat_grid.nelem());
1934 GridPos gp;
1935 Vector itw(2);
1936 gridpos(gp, lat_grid, pos[1]);
1937 interpweights(itw, gp);
1938 lat = interp(itw, lat_true, gp);
1939 lon = interp(itw, lon_true, gp);
1940 }
1941
1942 else {
1943 lat = pos[1];
1944 lon = pos[2];
1945 }
1946}
1947
1949 Workspace& ws,
1950 ArrayOfTensor3& diy_dx,
1951 ArrayOfTensor3& diy_dpath,
1952 const Index& ns,
1953 const Index& nf,
1954 const Index& np,
1955 const Index& atmosphere_dim,
1956 const Ppath& ppath,
1957 const Vector& ppvar_p,
1958 const Vector& ppvar_t,
1959 const Matrix& ppvar_vmr,
1960 const Index& iy_agenda_call1,
1961 const Tensor3& iy_transmittance,
1962 const Agenda& water_p_eq_agenda,
1963 const ArrayOfRetrievalQuantity& jacobian_quantities,
1964 const ArrayOfIndex& jac_species_i) {
1965 // Weight with iy_transmittance
1966 if (!iy_agenda_call1) {
1967 Matrix X, Y;
1968 //
1970 Y.resize(ns, diy_dpath[iq].npages());
1971 for (Index iv = 0; iv < nf; iv++) {
1972 X = transpose(diy_dpath[iq](joker, iv, joker));
1973 mult(Y, iy_transmittance(iv, joker, joker), X);
1974 diy_dpath[iq](joker, iv, joker) = transpose(Y);
1975 })
1976 }
1977
1978 // Handle abs species retrieval units, both internally and impact on T-jacobian
1979 //
1980 Tensor3 water_p_eq(0, 0, 0);
1981 //
1982 // Conversion for abs species itself
1983 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
1984 // Let x be VMR, and z the selected retrieval unit.
1985 // We have then that diy/dz = diy/dx * dx/dz
1986 //
1987 if (not(jacobian_quantities[iq] == Jacobian::Type::Sensor) and
1988 not(jacobian_quantities[iq] == Jacobian::Special::SurfaceString) and
1989 jac_species_i[iq] >= 0) {
1990 if (jacobian_quantities[iq].Mode() == "vmr") {
1991 }
1992
1993 else if (jacobian_quantities[iq].Mode() == "rel") {
1994 // Here x = vmr*z
1995 for (Index ip = 0; ip < np; ip++) {
1996 diy_dpath[iq](ip, joker, joker) *= ppvar_vmr(jac_species_i[iq], ip);
1997 }
1998 }
1999
2000 else if (jacobian_quantities[iq].Mode() == "nd") {
2001 // Here x = z/nd_tot
2002 for (Index ip = 0; ip < np; ip++) {
2003 diy_dpath[iq](ip, joker, joker) /=
2004 number_density(ppvar_p[ip], ppvar_t[ip]);
2005 }
2006 }
2007
2008 else if (jacobian_quantities[iq].Mode() == "rh") {
2009 // Here x = (p_sat/p) * z
2010 Tensor3 t_data(ppvar_t.nelem(), 1, 1);
2011 t_data(joker, 0, 0) = ppvar_t;
2012 water_p_eq_agendaExecute(ws, water_p_eq, t_data, water_p_eq_agenda);
2013 for (Index ip = 0; ip < np; ip++) {
2014 diy_dpath[iq](ip, joker, joker) *= water_p_eq(ip, 0, 0) / ppvar_p[ip];
2015 }
2016 }
2017
2018 else if (jacobian_quantities[iq].Mode() == "q") {
2019 // Here we use the approximation of x = z/0.622
2020 diy_dpath[iq](joker, joker, joker) /= 0.622;
2021 }
2022
2023 else {
2024 ARTS_ASSERT(0);
2025 }
2026 }
2027 }
2028
2029 // Correction of temperature Jacobian
2030 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
2031 // Let a be unit for abs species, and iy = f(T,a(T))
2032 // We have then that diy/dT = df/dT + df/da*da/dT
2033 // diy_dpath holds already df/dT. Remains is to add
2034 // df/da*da/dT for which abs species having da/dT != 0
2035 // This is only true for "nd" and "rh"
2036 //
2037 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
2038 // Loop abs species, again
2039 for (Index ia = 0; ia < jacobian_quantities.nelem(); ia++) {
2040 if (jac_species_i[ia] >= 0) {
2041 if (jacobian_quantities[ia].Mode() == "nd") {
2042 for (Index ip = 0; ip < np; ip++) {
2043 Matrix ddterm{diy_dpath[ia](ip, joker, joker)};
2044 ddterm *= ppvar_vmr(jac_species_i[ia], ip) *
2045 (number_density(ppvar_p[ip], ppvar_t[ip] + 1) -
2046 number_density(ppvar_p[ip], ppvar_t[ip]));
2047 diy_dpath[iq](ip, joker, joker) += ddterm;
2048 }
2049 } else if (jacobian_quantities[ia].Mode() == "rh") {
2050 Tensor3 t_data(ppvar_t.nelem(), 1, 1);
2051 t_data(joker, 0, 0) = ppvar_t;
2052 // Calculate water sat. pressure if not already done
2053 if (water_p_eq.npages() == 0) {
2055 ws, water_p_eq, t_data, water_p_eq_agenda);
2056 }
2057 // Sat.pressure for +1K
2058 Tensor3 water_p_eq1K;
2059 t_data(joker, 0, 0) += 1;
2061 ws, water_p_eq1K, t_data, water_p_eq_agenda);
2062
2063 for (Index ip = 0; ip < np; ip++) {
2064 const Numeric p_eq = water_p_eq(ip, 0, 0);
2065 const Numeric p_eq1K = water_p_eq1K(ip, 0, 0);
2066 Matrix ddterm{diy_dpath[ia](ip, joker, joker)};
2067 ddterm *= ppvar_vmr(jac_species_i[ia], ip) *
2068 (ppvar_p[ip] / pow(p_eq, 2.0)) * (p_eq1K - p_eq);
2069 diy_dpath[iq](ip, joker, joker) += ddterm;
2070 }
2071 }
2072 }
2073 }
2074 }
2075 }
2076
2077 // Map to retrieval grids
2079 jacobian_quantities[iq],
2080 diy_dpath[iq],
2081 atmosphere_dim,
2082 ppath,
2083 ppvar_p);)
2084}
2085
2087 Matrix& iy,
2088 ArrayOfTensor3& diy_dx,
2089 Tensor3& ppvar_iy,
2090 const Index& ns,
2091 const Index& np,
2092 const Vector& f_grid,
2093 const Ppath& ppath,
2094 const ArrayOfRetrievalQuantity& jacobian_quantities,
2095 const Index& j_analytical_do,
2096 const String& iy_unit) {
2097 // Determine refractive index to use for the n2 radiance law
2098 Numeric n = 1.0; // First guess is that sensor is in space
2099 //
2100 if (ppath.end_lstep == 0) // If true, sensor inside the atmosphere
2101 {
2102 n = ppath.nreal[np - 1];
2103 }
2104
2105 // Polarisation index variable
2106 ArrayOfIndex i_pol(ns);
2107 for (Index is = 0; is < ns; is++) {
2108 i_pol[is] = is + 1;
2109 }
2110
2111 // Jacobian part (must be converted to Tb before iy for PlanckBT)
2112 //
2113 if (j_analytical_do) {
2115 apply_iy_unit2(diy_dx[iq], iy, iy_unit, f_grid, n, i_pol);)
2116 }
2117
2118 // iy
2119 apply_iy_unit(iy, iy_unit, f_grid, n, i_pol);
2120
2121 // ppvar_iy
2122 for (Index ip = 0; ip < ppath.np; ip++) {
2124 ppvar_iy(joker, joker, ip), iy_unit, f_grid, ppath.nreal[ip], i_pol);
2125 }
2126}
2127
2128void yCalc_mblock_loop_body(bool& failed,
2129 String& fail_msg,
2130 ArrayOfArrayOfVector& iyb_aux_array,
2131 Workspace& ws,
2132 Vector& y,
2133 Vector& y_f,
2134 ArrayOfIndex& y_pol,
2135 Matrix& y_pos,
2136 Matrix& y_los,
2137 Matrix& y_geo,
2138 Matrix& jacobian,
2139 const Index& atmosphere_dim,
2140 const EnergyLevelMap& nlte_field,
2141 const Index& cloudbox_on,
2142 const Index& stokes_dim,
2143 const Vector& f_grid,
2144 const Matrix& sensor_pos,
2145 const Matrix& sensor_los,
2146 const Matrix& transmitter_pos,
2147 const Matrix& mblock_dlos,
2148 const Sparse& sensor_response,
2149 const Vector& sensor_response_f,
2150 const ArrayOfIndex& sensor_response_pol,
2151 const Matrix& sensor_response_dlos,
2152 const String& iy_unit,
2153 const Agenda& iy_main_agenda,
2154 const Agenda& jacobian_agenda,
2155 const Index& jacobian_do,
2156 const ArrayOfRetrievalQuantity& jacobian_quantities,
2157 const ArrayOfArrayOfIndex& jacobian_indices,
2158 const ArrayOfString& iy_aux_vars,
2159 const Verbosity& verbosity,
2160 const Index& mblock_index,
2161 const Index& n1y,
2162 const Index& j_analytical_do) {
2163 try {
2164 // Calculate monochromatic pencil beam data for 1 measurement block
2165 //
2166 Vector iyb, iyb_error, yb(n1y);
2167 ArrayOfMatrix diyb_dx;
2168 Matrix geo_pos_matrix;
2169 //
2170 iyb_calc(ws,
2171 iyb,
2172 iyb_aux_array[mblock_index],
2173 diyb_dx,
2174 geo_pos_matrix,
2175 mblock_index,
2176 atmosphere_dim,
2177 nlte_field,
2178 cloudbox_on,
2179 stokes_dim,
2180 f_grid,
2181 sensor_pos,
2182 sensor_los,
2183 transmitter_pos,
2184 mblock_dlos,
2185 iy_unit,
2186 iy_main_agenda,
2187 j_analytical_do,
2188 jacobian_quantities,
2189 jacobian_indices,
2190 iy_aux_vars,
2191 verbosity);
2192
2193 // Apply sensor response matrix on iyb, and put into y
2194 //
2195 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
2196 const Index row0 = rowind.offset;
2197 //
2198 mult(yb, sensor_response, iyb);
2199 //
2200 y[rowind] = yb; // *yb* also used below, as input to jacobian_agenda
2201
2202 // Fill information variables. And search for NaNs in *y*.
2203 //
2204 for (Index i = 0; i < n1y; i++) {
2205 const Index ii = row0 + i;
2206 ARTS_USER_ERROR_IF (std::isnan(y[ii]),
2207 "One or several NaNs found in *y*.");
2208 y_f[ii] = sensor_response_f[i];
2209 y_pol[ii] = sensor_response_pol[i];
2210 y_pos(ii, joker) = sensor_pos(mblock_index, joker);
2211 y_los(ii, joker) = sensor_los(mblock_index, joker);
2212 y_los(ii, 0) += sensor_response_dlos(i, 0);
2213 if (sensor_response_dlos.ncols() > 1) {
2214 y_los(ii, 1) += sensor_response_dlos(i, 1);
2215 }
2216 }
2217
2218 // Apply sensor response matrix on diyb_dx, and put into jacobian
2219 // (that is, analytical jacobian part)
2220 //
2221 if (j_analytical_do) {
2223 mult(jacobian(rowind,
2224 Range(jacobian_indices[iq][0],
2225 jacobian_indices[iq][1] -
2226 jacobian_indices[iq][0] + 1)),
2227 sensor_response,
2228 diyb_dx[iq]);)
2229 }
2230
2231 // Calculate remaining parts of *jacobian*
2232 //
2233 if (jacobian_do) {
2235 ws, jacobian, mblock_index, iyb, yb, jacobian_agenda);
2236 }
2237
2238 // Handle geo-positioning
2239 if (!std::isnan(geo_pos_matrix(0, 0))) // No data are flagged as NaN
2240 {
2241 // We set geo_pos based on the max value in sensor_response
2242 const Index nfs = f_grid.nelem() * stokes_dim;
2243 for (Index i = 0; i < n1y; i++) {
2244 Index jmax = -1;
2245 Numeric rmax = -99e99;
2246 for (Index j = 0; j < sensor_response.ncols(); j++) {
2247 if (sensor_response(i, j) > rmax) {
2248 rmax = sensor_response(i, j);
2249 jmax = j;
2250 }
2251 }
2252 const auto jhit = Index(floor(jmax / nfs));
2253 y_geo(row0 + i, joker) = geo_pos_matrix(jhit, joker);
2254 }
2255 }
2256 }
2257
2258 catch (const std::exception& e) {
2259#pragma omp critical(yCalc_fail)
2260 {
2261 fail_msg = e.what();
2262 failed = true;
2263 }
2264 }
2265}
2266
2267void ze_cfac(Vector& fac,
2268 const Vector& f_grid,
2269 const Numeric& ze_tref,
2270 const Numeric& k2) {
2271 const Index nf = f_grid.nelem();
2272
2273 ARTS_ASSERT(fac.nelem() == nf);
2274
2275 // Refractive index for water (if needed)
2276 Matrix complex_n(0, 0);
2277 if (k2 <= 0) {
2278 complex_n_water_liebe93(complex_n, f_grid, ze_tref);
2279 }
2280
2281 // Common conversion factor
2282 static constexpr Numeric a = 4e18 / Math::pow4(Constant::pi);
2283
2284 for (Index iv = 0; iv < nf; iv++) {
2285 // Reference dielectric factor.
2286 Numeric K2;
2287 if (k2 >= 0) {
2288 K2 = k2;
2289 } else {
2290 Complex n(complex_n(iv, 0), complex_n(iv, 1));
2291 Complex n2 = n * n;
2292 Complex K = (n2 - Numeric(1.0)) / (n2 + Numeric(2.0));
2293 Numeric absK = abs(K);
2294 K2 = absK * absK;
2295 }
2296
2297 // Wavelength
2298 Numeric la = SPEED_OF_LIGHT / f_grid[iv];
2299
2300 fac[iv] = a * la * la * la * la / K2;
2301 }
2302}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:119
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:224
base max(const Array< base > &x)
Max function.
Definition: array.h:128
base min(const Array< base > &x)
Min function.
Definition: array.h:144
Constants of physical expressions as constexpr.
Simple constexpr math that we can make use of in Arts.
Common ARTS conversions.
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:29
void iy_cloudbox_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:25135
void iy_surface_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const ArrayOfTensor4 &dsurface_rmatrix_dx, const ArrayOfMatrix &dsurface_emission_dx, const String &iy_unit, const Tensor3 &iy_transmittance, const Index iy_id, const Index cloudbox_on, const Index jacobian_do, const Agenda &iy_main_agenda, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const Tensor3 &surface_props_data, const ArrayOfString &dsurface_names, const Agenda &input_agenda)
Definition: auto_md.cc:25457
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:25421
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, Vector &geo_pos, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const String &iy_unit, const Index cloudbox_on, const Index jacobian_do, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition: auto_md.cc:25309
void jacobian_agendaExecute(Workspace &ws, Matrix &jacobian, const Index mblock_index, const Vector &iyb, const Vector &yb, const Agenda &input_agenda)
Definition: auto_md.cc:25519
void propmat_clearsky_agendaExecute(Workspace &ws, PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfSpeciesTag &select_abs_species, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:25794
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:26129
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:608
The Agenda class.
Definition: agenda_class.h:52
This can be used to make arrays out of anything.
Definition: array.h:31
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
Workspace class.
Definition: workspace_ng.h:36
bool is_gp_inside_cloudbox(const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, const bool &include_boundaries, const Index &atmosphere_dim)
Definition: cloudbox.cc:443
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:84
#define ARTS_USER_ERROR(...)
Definition: debug.h:151
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:135
void distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
Definition: geodetic.cc:652
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:165
void cart2sph(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
cart2sph
Definition: geodetic.cc:585
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1248
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1318
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:338
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:57
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Maps jacobian data for points along the propagation path, to jacobian retrieval grid data.
Definition: jacobian.cc:296
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition: jacobian.h:540
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:532
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:529
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:46
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:318
#define CREATE_OUT3
Definition: messages.h:189
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 Numeric c
Speed of light convenience name [m/s].
constexpr Numeric temperature_at_0c
Global constant, Temperature in Celsius of 0 Kelvin.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
auto sind(auto x) noexcept
Returns the sine of the deg2rad of the input.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
constexpr auto pow4(auto x) noexcept
power of four
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
void pha_mat_1ScatElem(Tensor6View pha_mat, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_start, const Index &t_interp_order)
Preparing phase matrix from one scattering element.
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
Numeric dinvplanckdI(const Numeric &i, const Numeric &f)
dinvplanckdI
Numeric planck(const Numeric &f, const Numeric &t)
planck
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
Numeric invplanck(const Numeric &i, const Numeric &f)
invplanck
Numeric dplanck_dt(const Numeric &f, const Numeric &t)
dplanck_dt
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:51
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to a cartesian unit vector.
Definition: ppath.cc:296
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1460
void add_za_aa(Numeric &za, Numeric &aa, const Numeric &za0, const Numeric &aa0, const Numeric &dza, const Numeric &daa)
Adds up zenith and azimuth angles.
Definition: ppath.cc:375
void cart2zaaa(Numeric &za, Numeric &aa, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Converts a cartesian directional vector to zenith and azimuth.
Definition: ppath.cc:283
void find_tanpoint(Index &it, const Ppath &ppath)
Identifies the tangent point of a propagation path.
Definition: ppath.cc:494
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
This is the core for the WSM ppathStepByStep.
Definition: ppath.cc:5162
Propagation path structure and functions.
void complex_n_water_liebe93(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_water_liebe93
Definition: refraction.cc:57
Refraction functions.
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstMatrixView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstVectorView &ppath_line_of_sight, const ConstVectorView &ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
Definition: rte.cc:1188
void ze_cfac(Vector &fac, const Vector &f_grid, const Numeric &ze_tref, const Numeric &k2)
Calculates factor to convert back-scattering to Ze.
Definition: rte.cc:2267
void iyb_calc_body(bool &failed, String &fail_msg, ArrayOfArrayOfMatrix &iy_aux_array, Workspace &ws, Ppath &ppath, Vector &iyb, ArrayOfMatrix &diyb_dx, Vector &geo_pos, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Vector &f_grid, const ArrayOfString &iy_aux_vars, const Index &ilos, const Index &nf)
Definition: rte.cc:1473
void mueller_stokes2modif(Matrix &Cm, const Index &stokes_dim)
mueller_stokes2modif
Definition: rte.cc:1893
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:1305
constexpr Numeric SPEED_OF_LIGHT
Definition: rte.cc:31
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
Definition: rte.cc:828
void pos2true_latlon(Numeric &lat, Numeric &lon, const Index &atmosphere_dim, const ConstVectorView &lat_grid, const ConstVectorView &lat_true, const ConstVectorView &lon_true, const ConstVectorView &pos)
Determines the true alt and lon for an "ARTS position".
Definition: rte.cc:1914
Numeric dotprod_with_los(const ConstVectorView &los, const Numeric &u, const Numeric &v, const Numeric &w, const Index &atmosphere_dim)
Calculates the dot product between a field and a LOS.
Definition: rte.cc:651
void apply_iy_unit2(Tensor3View J, const ConstMatrixView &iy, const String &iy_unit, const ConstVectorView &f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
Largely as apply_iy_unit but operates on jacobian data.
Definition: rte.cc:162
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:1755
void apply_iy_unit(MatrixView iy, const String &iy_unit, const ConstVectorView &f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
Performs conversion from radiance to other units, as well as applies refractive index to fulfill the ...
Definition: rte.cc:89
void get_iy(Workspace &ws, Matrix &iy, const Index &cloudbox_on, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda)
Basic call of iy_main_agenda.
Definition: rte.cc:672
void adjust_los(VectorView los, const Index &atmosphere_dim)
Ensures that the zenith and azimuth angles of a line-of-sight vector are inside defined ranges.
Definition: rte.cc:66
void defocusing_general_sub(Workspace &ws, Vector &pos, Vector &rte_los, Index &background, const Vector &rte_pos, const Numeric &lo0, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Verbosity &verbosity)
Just to avoid duplicatuion of code in defocusing_general.
Definition: rte.cc:294
void defocusing_sat2sat(Workspace &ws, Numeric &dlf, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Ppath &ppath, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Numeric &dza, const Verbosity &verbosity)
Calculates defocusing for limb measurements between two satellites.
Definition: rte.cc:535
void mueller_modif2stokes(Matrix &Cs, const Index &stokes_dim)
mueller_modif2stokes
Definition: rte.cc:1845
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const String &iy_unit, const Tensor3 &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:713
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, const ConstVectorView &f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, const ConstMatrixView &ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
Definition: rte.cc:1051
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, Matrix &geo_pos_matrix, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
Performs calculations for one measurement block, on iy-level.
Definition: rte.cc:1591
void mueller_rotation(Matrix &L, const Index &stokes_dim, const Numeric &rotangle)
mueller_rotation
Definition: rte.cc:1866
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:2086
void bending_angle1d(Numeric &alpha, const Ppath &ppath)
Calculates the bending angle for a 1D atmosphere.
Definition: rte.cc:249
void rtmethods_jacobian_finalisation(Workspace &ws, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &atmosphere_dim, const Ppath &ppath, const Vector &ppvar_p, const Vector &ppvar_t, const Matrix &ppvar_vmr, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &jac_species_i)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
Definition: rte.cc:1948
void mirror_los(Vector &los_mirrored, const ConstVectorView &los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:1792
void yCalc_mblock_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfVector &iyb_aux_array, Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Matrix &y_geo, Matrix &jacobian, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity, const Index &mblock_index, const Index &n1y, const Index &j_analytical_do)
Performs calculations for one measurement block, on y-level.
Definition: rte.cc:2128
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &ppath_f_grid, const Vector &ppath_magnetic_field, const Vector &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const Vector &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition: rte.cc:1116
constexpr Numeric TEMP_0_C
Definition: rte.cc:32
void muellersparse_rotation(Sparse &H, const Index &stokes_dim, const Numeric &rotangle)
muellersparse_rotation
Definition: rte.cc:1816
Vector get_stepwise_f_partials(const ConstVectorView &line_of_sight, const ConstVectorView &f_grid, const Jacobian::Atm wind_type, const Index &atmosphere_dim)
Computes the ratio that a partial derivative with regards to frequency relates to the wind of come co...
Definition: rte.cc:1154
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
Returns the "range" of y corresponding to a measurement block.
Definition: rte.cc:1095
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_line_of_sight, const Index &lte, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
Definition: rte.cc:38
void get_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:1101
void defocusing_general(Workspace &ws, Numeric &dlf, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Ppath &ppath, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Numeric &dza, const Verbosity &verbosity)
Defocusing for arbitrary geometry (zenith angle part only)
Definition: rte.cc:416
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Definition: rte.cc:947
Declaration of functions in rte.cc.
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
Converts atmospheric a grid position to weights for interpolation of a field defined ONLY inside the ...
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Header file for special_interp.cc.
EnergyLevelMap InterpToGridPos(Index atmosphere_dim, const ArrayOfGridPos &p, const ArrayOfGridPos &lat, const ArrayOfGridPos &lon) const
Structure to store a grid position.
Definition: interpolation.h:56
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
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath_struct.h:55
String background
Radiative background.
Definition: ppath_struct.h:25
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath_struct.h:33
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath_struct.h:53
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath_struct.h:45
Vector start_pos
Start position.
Definition: ppath_struct.h:27
Vector lstep
The length between ppath points.
Definition: ppath_struct.h:39
Numeric start_lstep
Length between sensor and atmospheric boundary.
Definition: ppath_struct.h:31
Numeric constant
The propagation path constant (only used for 1D)
Definition: ppath_struct.h:23
Vector r
Radius of each ppath point.
Definition: ppath_struct.h:37
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath_struct.h:51
Index dim
Atmospheric dimensionality.
Definition: ppath_struct.h:19
Vector end_pos
End position.
Definition: ppath_struct.h:41
Vector start_los
Start line-of-sight.
Definition: ppath_struct.h:29
Vector nreal
The real part of the refractive index at each path position.
Definition: ppath_struct.h:47
Vector end_los
End line-of-sight.
Definition: ppath_struct.h:43
#define u
#define v
#define w
#define a
#define b