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