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