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