ARTS 2.5.0 (git: 9ee3ac6c)
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().AtmType(), 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 MatrixView J,
1170 Tensor3View dJ_dx,
1171 const PropagationMatrix& K,
1172 const StokesVector& a,
1173 const StokesVector& S,
1174 const ArrayOfPropagationMatrix& dK_dx,
1175 const ArrayOfStokesVector& da_dx,
1176 const ArrayOfStokesVector& dS_dx,
1177 const ConstVectorView& B,
1178 const ConstVectorView& dB_dT,
1179 const ArrayOfRetrievalQuantity& jacobian_quantities,
1180 const bool& jacobian_do) {
1181 const Index nq = jacobian_quantities.nelem();
1182 const Index nf = K.NumberOfFrequencies();
1183 const Index ns = K.StokesDimensions();
1184
1185 dJ_dx = 0;
1186 for (Index i1 = 0; i1 < nf; i1++) {
1187 Matrix invK(ns, ns);
1188 Vector j(ns);
1189
1190 if (K.IsRotational(i1)) {
1191 dJ_dx(joker, i1, joker) = 0;
1192 J(i1, joker) = 0;
1193 continue;
1194 }
1195
1196 K.MatrixInverseAtPosition(invK, i1);
1197
1198 // Set a B to j
1199 j = a.VectorAtPosition(i1);
1200 ;
1201 j *= B[i1];
1202
1203 // Add S to j
1204 if (not S.IsEmpty()) j += S.VectorAtPosition(i1);
1205 ;
1206
1207 // Compute J = K^-1 (a B + S)
1208 mult(J(i1, joker), invK, j);
1209
1210 // Compute dJ = K^-1((da B + a dB + dS) - dK K^-1(a B + S))
1211 if (jacobian_do)
1212 //FOR_ANALYTICAL_JACOBIANS_DO
1213 //(
1214 for (Index iq = 0; iq < nq; iq++) {
1215 if (not(jacobian_quantities[iq] == Jacobian::Type::Sensor) and not(jacobian_quantities[iq] == Jacobian::Special::SurfaceString)) {
1216 Matrix dk(ns, ns), tmp_matrix(ns, ns);
1217 Vector dj(ns, 0), tmp(ns);
1218
1219 // Control parameters for special jacobians that are computed elsewhere
1220 //const bool has_dk = (not dK_dx[iq].IsEmpty()); // currently always
1221 //const bool has_ds = (not dS_dx[iq].IsEmpty()); // evaluate as true
1222 const bool has_dt =
1223 (jacobian_quantities[iq] == Jacobian::Atm::Temperature);
1224
1225 // Sets the -K^-1 dK/dx K^-1 (a B + S) term
1226 //if(has_dk)
1227 //{
1228 dK_dx[iq].MatrixAtPosition(dk, i1);
1229 mult(tmp, dk, J(i1, joker));
1230
1231 dj = da_dx[iq].VectorAtPosition(i1);
1232 dj *= B[i1];
1233
1234 dj -= tmp;
1235
1236 // Adds a dB to dj
1237 if (has_dt) {
1238 tmp = a.VectorAtPosition(i1);
1239 tmp *= dB_dT[i1];
1240 dj += tmp;
1241 }
1242
1243 // Adds dS to dj
1244 //if(has_ds)
1245 dj += dS_dx[iq].VectorAtPosition(i1);
1246
1247 mult(dJ_dx(iq, i1, joker), invK, dj);
1248 //}
1249 }
1250 // don't need that anymore now that we zero dJ_dx at the very beginning
1251 //else
1252 //{
1253 // dJ_dx(iq, i1, joker) = 0;
1254 //}
1255 }
1256 //)
1257 }
1258
1259 if (nq) dJ_dx *= 0.5;
1260}
1261
1263 const ConstVectorView& f_grid,
1264 const ConstVectorView& ppath_wind,
1265 const ConstVectorView& ppath_line_of_sight,
1266 const Numeric& rte_alonglos_v,
1267 const Index& atmosphere_dim) {
1268 Numeric v_doppler = rte_alonglos_v;
1269
1270 if (ppath_wind[0] not_eq 0 or ppath_wind[1] not_eq 0 or
1271 ppath_wind[2] not_eq 0)
1272 v_doppler += dotprod_with_los(ppath_line_of_sight,
1273 ppath_wind[0],
1274 ppath_wind[1],
1275 ppath_wind[2],
1276 atmosphere_dim);
1277 ppath_f_grid = f_grid;
1278
1279 if (v_doppler not_eq 0) ppath_f_grid *= 1 - v_doppler / SPEED_OF_LIGHT;
1280}
1281
1283 const ConstVectorView& f_grid,
1284 const Jacobian::Atm wind_type,
1285 const Index& atmosphere_dim) {
1286 // Doppler relevant velocity
1287 Numeric dv_doppler_dx = 0.0;
1288
1289 Vector deriv(f_grid);
1290
1291 switch (wind_type) {
1293 dv_doppler_dx = 1.0;
1294 break;
1296 dv_doppler_dx =
1297 (dotprod_with_los(line_of_sight, 1, 0, 0, atmosphere_dim));
1298 break;
1300 dv_doppler_dx =
1301 (dotprod_with_los(line_of_sight, 0, 1, 0, atmosphere_dim));
1302 break;
1304 dv_doppler_dx =
1305 (dotprod_with_los(line_of_sight, 0, 0, 1, atmosphere_dim));
1306 break;
1307 default:
1308 ARTS_ASSERT(false, "Not allowed to call this function without a wind parameter as wind_type");
1309 break;
1310 }
1311
1312 deriv *= - dv_doppler_dx / Constant::c;
1313 return deriv;
1314}
1315
1317 StokesVector& ap,
1319 ArrayOfStokesVector& dap_dx,
1321 const ArrayOfRetrievalQuantity& jacobian_quantities,
1322 const ConstMatrixView& ppath_1p_pnd, // the ppath_pnd at this ppath point
1323 const ArrayOfMatrix&
1324 ppath_dpnd_dx, // the full ppath_dpnd_dx, ie all ppath points
1325 const Index ppath_1p_id,
1326 const ArrayOfArrayOfSingleScatteringData& scat_data,
1327 const ConstVectorView& ppath_line_of_sight,
1328 const ConstVectorView& ppath_temperature,
1329 const Index& atmosphere_dim,
1330 const bool& jacobian_do) {
1331 const Index nf = Kp.NumberOfFrequencies(), stokes_dim = Kp.StokesDimensions();
1332
1333 //StokesVector da_aux(nf, stokes_dim);
1334 //PropagationMatrix dK_aux(nf, stokes_dim);
1335
1337
1338 // Direction of outgoing scattered radiation (which is reversed to
1339 // LOS). Only used for extracting scattering properties.
1340 Vector dir;
1341 mirror_los(dir, ppath_line_of_sight, atmosphere_dim);
1342 Matrix dir_array(1, 2, 0.);
1343 dir_array(0, joker) = dir;
1344
1345 ArrayOfArrayOfTensor5 ext_mat_Nse;
1346 ArrayOfArrayOfTensor4 abs_vec_Nse;
1347 ArrayOfArrayOfIndex ptypes_Nse;
1348 Matrix t_ok;
1349 ArrayOfTensor5 ext_mat_ssbulk;
1350 ArrayOfTensor4 abs_vec_ssbulk;
1351 ArrayOfIndex ptype_ssbulk;
1352 Tensor5 ext_mat_bulk;
1353 Tensor4 abs_vec_bulk;
1354 Index ptype_bulk;
1355
1356 // get per-scat-elem data here. and fold with pnd.
1357 // keep per-scat-elem data to fold with dpnd_dx further down in
1358 // analyt-jac-loop.
1359 opt_prop_NScatElems(ext_mat_Nse,
1360 abs_vec_Nse,
1361 ptypes_Nse,
1362 t_ok,
1363 scat_data,
1364 stokes_dim,
1365 ppath_temperature,
1366 dir_array,
1367 -1);
1368
1369 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
1370 abs_vec_ssbulk,
1371 ptype_ssbulk,
1372 ext_mat_Nse,
1373 abs_vec_Nse,
1374 ptypes_Nse,
1375 ppath_1p_pnd,
1376 t_ok);
1377 opt_prop_Bulk(ext_mat_bulk,
1378 abs_vec_bulk,
1379 ptype_bulk,
1380 ext_mat_ssbulk,
1381 abs_vec_ssbulk,
1382 ptype_ssbulk);
1383
1384 const Index nf_ssd = abs_vec_bulk.nbooks(); // number of freqs in extracted
1385 // optprops. if 1, we need to
1386 // duplicate the ext/abs output.
1387
1388 for (Index iv = 0; iv < nf; iv++) {
1389 if (nf_ssd > 1) {
1390 ap.SetAtPosition(abs_vec_bulk(iv, 0, 0, joker), iv);
1391 Kp.SetAtPosition(ext_mat_bulk(iv, 0, 0, joker, joker), iv);
1392 } else {
1393 ap.SetAtPosition(abs_vec_bulk(0, 0, 0, joker), iv);
1394 Kp.SetAtPosition(ext_mat_bulk(0, 0, 0, joker, joker), iv);
1395 }
1396 }
1397
1398 if (jacobian_do)
1400 if (ppath_dpnd_dx[iq].empty()) {
1401 dap_dx[iq].SetZero();
1402 dKp_dx[iq].SetZero();
1403 } else {
1404 // check, whether we have any non-zero ppath_dpnd_dx in this
1405 // pnd-affecting x? might speed up things a little bit.
1406 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
1407 abs_vec_ssbulk,
1408 ptype_ssbulk,
1409 ext_mat_Nse,
1410 abs_vec_Nse,
1411 ptypes_Nse,
1412 ppath_dpnd_dx[iq](joker, Range(ppath_1p_id, 1)),
1413 t_ok);
1414 opt_prop_Bulk(ext_mat_bulk,
1415 abs_vec_bulk,
1416 ptype_bulk,
1417 ext_mat_ssbulk,
1418 abs_vec_ssbulk,
1419 ptype_ssbulk);
1420 for (Index iv = 0; iv < nf; iv++) {
1421 if (nf_ssd > 1) {
1422 dap_dx[iq].SetAtPosition(abs_vec_bulk(iv, 0, 0, joker), iv);
1423 dKp_dx[iq].SetAtPosition(ext_mat_bulk(iv, 0, 0, joker, joker),
1424 iv);
1425 } else {
1426 dap_dx[iq].SetAtPosition(abs_vec_bulk(0, 0, 0, joker), iv);
1427 dKp_dx[iq].SetAtPosition(ext_mat_bulk(0, 0, 0, joker, joker), iv);
1428 }
1429 }
1430 })
1431}
1432
1434 StokesVector& Sp,
1435 ArrayOfStokesVector& dSp_dx,
1436 const ArrayOfRetrievalQuantity& jacobian_quantities,
1437 const ConstVectorView& ppath_1p_pnd, // the ppath_pnd at this ppath point
1438 const ArrayOfMatrix&
1439 ppath_dpnd_dx, // the full ppath_dpnd_dx, ie all ppath points
1440 const Index ppath_1p_id,
1441 const ArrayOfArrayOfSingleScatteringData& scat_data,
1442 const ConstTensor7View& cloudbox_field,
1443 const ConstVectorView& za_grid,
1444 const ConstVectorView& aa_grid,
1445 const ConstMatrixView& ppath_line_of_sight,
1446 const GridPos& ppath_pressure,
1447 const Vector& temperature,
1448 const Index& atmosphere_dim,
1449 const bool& jacobian_do,
1450 const Index& t_interp_order) {
1451 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
1452 "This function handles so far only 1D atmospheres.");
1453
1454 const Index nf = Sp.NumberOfFrequencies();
1455 const Index stokes_dim = Sp.StokesDimensions();
1456 const Index ne = ppath_1p_pnd.nelem();
1457 ARTS_ASSERT(TotalNumberOfElements(scat_data) == ne);
1458 const Index nza = za_grid.nelem();
1459 const Index naa = aa_grid.nelem();
1460 const Index nq = jacobian_do ? jacobian_quantities.nelem() : 0;
1461
1462 // interpolate incident field to this ppath point (no need to do this
1463 // separately per scatelem)
1464 GridPos gp_p;
1465 gridpos_copy(gp_p, ppath_pressure);
1466 Vector itw_p(2);
1467 interpweights(itw_p, gp_p);
1468 Tensor3 inc_field(nf, nza, stokes_dim, 0.);
1469 for (Index iv = 0; iv < nf; iv++) {
1470 for (Index iza = 0; iza < nza; iza++) {
1471 for (Index i = 0; i < stokes_dim; i++) {
1472 inc_field(iv, iza, i) =
1473 interp(itw_p, cloudbox_field(iv, joker, 0, 0, iza, 0, i), gp_p);
1474 }
1475 }
1476 }
1477
1478 // create matrix of incident directions (flat representation of the
1479 // za_grid * aa_grid matrix)
1480 Matrix idir(nza * naa, 2);
1481 Index ia = 0;
1482 for (Index iza = 0; iza < nza; iza++) {
1483 for (Index iaa = 0; iaa < naa; iaa++) {
1484 idir(ia, 0) = za_grid[iza];
1485 idir(ia, 1) = aa_grid[iaa];
1486 ia++;
1487 }
1488 }
1489
1490 // setting prop (aka scattered) direction
1491 Matrix pdir(1, 2);
1492 //if( ppath_line_of_sight.ncols()==2 )
1493 // pdir(0,joker) = ppath_line_of_sight;
1494 //else // 1D case only (currently the only handled case). azimuth not defined.
1495 {
1496 pdir(0, 0) = ppath_line_of_sight(0, 0);
1497 pdir(0, 1) = 0.;
1498 }
1499
1500 // some further variables needed for pha_mat extraction
1501 Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
1502 Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
1503 Tensor6 pha_mat_1se(nf_ssd, 1, 1, nza * naa, stokes_dim, stokes_dim);
1504 Vector t_ok(1);
1505 Index ptype;
1506 Tensor3 scat_source_1se(ne, nf, stokes_dim, 0.);
1507
1508 Index ise_flat = 0;
1509 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++) {
1510 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
1511 // determine whether we have some valid pnd for this
1512 // scatelem (in pnd or dpnd)
1513 Index val_pnd = 0;
1514 if (ppath_1p_pnd[ise_flat] != 0) {
1515 val_pnd = 1;
1516 } else if (jacobian_do) {
1517 for (Index iq = 0; (!val_pnd) && (iq < nq); iq++) {
1518 if ((not(jacobian_quantities[iq] == Jacobian::Type::Sensor) and not(jacobian_quantities[iq] == Jacobian::Special::SurfaceString)) &&
1519 !ppath_dpnd_dx[iq].empty() &&
1520 ppath_dpnd_dx[iq](ise_flat, ppath_1p_id) != 0) {
1521 val_pnd = 1;
1522 }
1523 }
1524 }
1525
1526 if (val_pnd) {
1527 pha_mat_1ScatElem(pha_mat_1se,
1528 ptype,
1529 t_ok,
1530 scat_data[i_ss][i_se],
1531 temperature,
1532 pdir,
1533 idir,
1534 0,
1535 t_interp_order);
1536 ARTS_USER_ERROR_IF (t_ok[0] == 0,
1537 "Interpolation error for (flat-array) scattering "
1538 "element #", ise_flat, "\n"
1539 "at location/temperature point #", ppath_1p_id, "\n")
1540
1541 Index this_iv = 0;
1542 for (Index iv = 0; iv < nf; iv++) {
1543 if (!duplicate_freqs) {
1544 this_iv = iv;
1545 }
1546 Tensor3 product_fields(nza, naa, stokes_dim, 0.);
1547
1548 ia = 0;
1549 for (Index iza = 0; iza < nza; iza++) {
1550 for (Index iaa = 0; iaa < naa; iaa++) {
1551 for (Index i = 0; i < stokes_dim; i++) {
1552 for (Index j = 0; j < stokes_dim; j++) {
1553 product_fields(iza, iaa, i) +=
1554 pha_mat_1se(this_iv, 0, 0, ia, i, j) *
1555 inc_field(iv, iza, j);
1556 }
1557 }
1558 ia++;
1559 }
1560 }
1561
1562 for (Index i = 0; i < stokes_dim; i++) {
1563 scat_source_1se(ise_flat, iv, i) = AngIntegrate_trapezoid(
1564 product_fields(joker, joker, i), za_grid, aa_grid);
1565 }
1566 } // for iv
1567 } // if val_pnd
1568
1569 ise_flat++;
1570
1571 } // for i_se
1572 } // for i_ss
1573
1574 for (Index iv = 0; iv < nf; iv++) {
1575 Vector scat_source(stokes_dim, 0.);
1576 for (ise_flat = 0; ise_flat < ne; ise_flat++) {
1577 for (Index i = 0; i < stokes_dim; i++) {
1578 scat_source[i] +=
1579 scat_source_1se(ise_flat, iv, i) * ppath_1p_pnd[ise_flat];
1580 }
1581 }
1582
1583 Sp.SetAtPosition(scat_source, iv);
1584
1585 if (jacobian_do) {
1587 if (ppath_dpnd_dx[iq].empty()) { dSp_dx[iq].SetZero(); } else {
1588 scat_source = 0.;
1589 for (ise_flat = 0; ise_flat < ne; ise_flat++) {
1590 for (Index i = 0; i < stokes_dim; i++) {
1591 scat_source[i] += scat_source_1se(ise_flat, iv, i) *
1592 ppath_dpnd_dx[iq](ise_flat, ppath_1p_id);
1593 dSp_dx[iq].SetAtPosition(scat_source, iv);
1594 }
1595 }
1596 })
1597 }
1598 } // for iv
1599}
1600
1603 Tensor3View T,
1604 Tensor4View dT_close_dx,
1605 Tensor4View dT_far_dx,
1606 const ConstTensor3View& cumulative_transmission_close,
1607 const PropagationMatrix& K_close,
1608 const PropagationMatrix& K_far,
1609 const ArrayOfPropagationMatrix& dK_close_dx,
1610 const ArrayOfPropagationMatrix& dK_far_dx,
1611 const Numeric& ppath_distance,
1612 const bool& first_level,
1613 const Numeric& dppath_distance_dT_HSE_close,
1614 const Numeric& dppath_distance_dT_HSE_far,
1615 const Index& temperature_derivative_position_if_hse_is_active) {
1616 // Frequency counter
1617 const Index nf = K_close.NumberOfFrequencies();
1618 const Index stokes_dim = T.ncols();
1619
1620 if (first_level) {
1621 if (stokes_dim > 1)
1622 for (Index iv = 0; iv < nf; iv++)
1624 else
1626 } else {
1627 // Compute the transmission of the layer between close and far
1628 if (not dK_close_dx.nelem())
1629 compute_transmission_matrix(T, ppath_distance, K_close, K_far);
1630 else
1632 T,
1633 dT_close_dx,
1634 dT_far_dx,
1635 ppath_distance,
1636 K_close,
1637 K_far,
1638 dK_close_dx,
1639 dK_far_dx,
1640 dppath_distance_dT_HSE_close,
1641 dppath_distance_dT_HSE_far,
1642 temperature_derivative_position_if_hse_is_active);
1643
1644 // Cumulate transmission
1645 if (stokes_dim > 1)
1646 for (Index iv = 0; iv < nf; iv++)
1648 cumulative_transmission_close(iv, joker, joker),
1649 T(iv, joker, joker));
1650 else {
1651 cumulative_transmission = cumulative_transmission_close;
1653 }
1654 }
1655}
1656
1657void iyb_calc_body(bool& failed,
1658 String& fail_msg,
1659 ArrayOfArrayOfMatrix& iy_aux_array,
1660 Workspace& ws,
1661 Ppath& ppath,
1662 Vector& iyb,
1663 ArrayOfMatrix& diyb_dx,
1664 const Index& mblock_index,
1665 const Index& atmosphere_dim,
1666 const EnergyLevelMap& nlte_field,
1667 const Index& cloudbox_on,
1668 const Index& stokes_dim,
1669 const ConstMatrixView& sensor_pos,
1670 const ConstMatrixView& sensor_los,
1671 const ConstMatrixView& transmitter_pos,
1672 const ConstMatrixView& mblock_dlos_grid,
1673 const String& iy_unit,
1674 const Agenda& iy_main_agenda,
1675 const Index& j_analytical_do,
1676 const ArrayOfRetrievalQuantity& jacobian_quantities,
1677 const ArrayOfArrayOfIndex& jacobian_indices,
1678 const ConstVectorView& f_grid,
1679 const ArrayOfString& iy_aux_vars,
1680 const Index& ilos,
1681 const Index& nf) {
1682 // The try block here is necessary to correctly handle
1683 // exceptions inside the parallel region.
1684 try {
1685 //--- LOS of interest
1686 //
1687 Vector los(sensor_los.ncols());
1688 //
1689 los = sensor_los(mblock_index, joker);
1690 if (mblock_dlos_grid.ncols() == 1) {
1691 los[0] += mblock_dlos_grid(ilos, 0);
1692 adjust_los(los, atmosphere_dim);
1693 } else {
1694 add_za_aa(los[0],
1695 los[1],
1696 los[0],
1697 los[1],
1698 mblock_dlos_grid(ilos, 0),
1699 mblock_dlos_grid(ilos, 1));
1700 }
1701
1702 //--- rtp_pos 1 and 2
1703 //
1704 Vector rtp_pos, rtp_pos2(0);
1705 //
1706 rtp_pos = sensor_pos(mblock_index, joker);
1707 if (transmitter_pos.nrows()) {
1708 rtp_pos2 = transmitter_pos(mblock_index, joker);
1709 }
1710
1711 // Calculate iy and associated variables
1712 //
1713 Matrix iy;
1714 ArrayOfTensor3 diy_dx;
1715 Tensor3 iy_transmittance(0, 0, 0);
1716 const Index iy_agenda_call1 = 1;
1717 const Index iy_id =
1718 (Index)1e6 * (mblock_index + 1) + (Index)1e3 * (ilos + 1);
1719 //
1721 iy,
1722 iy_aux_array[ilos],
1723 ppath,
1724 diy_dx,
1725 iy_agenda_call1,
1726 iy_transmittance,
1727 iy_aux_vars,
1728 iy_id,
1729 iy_unit,
1730 cloudbox_on,
1731 j_analytical_do,
1732 f_grid,
1733 nlte_field,
1734 rtp_pos,
1735 los,
1736 rtp_pos2,
1737 iy_main_agenda);
1738
1739 // Start row in iyb etc. for present LOS
1740 //
1741 const Index row0 = ilos * nf * stokes_dim;
1742
1743 // Jacobian part
1744 //
1745 if (j_analytical_do) {
1747 for (Index ip = 0;
1748 ip < jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1;
1749 ip++) {
1750 for (Index is = 0; is < stokes_dim; is++) {
1751 diyb_dx[iq](Range(row0 + is, nf, stokes_dim), ip) =
1752 diy_dx[iq](ip, joker, is);
1753 }
1754 })
1755 }
1756
1757 // iy : copy to iyb
1758 for (Index is = 0; is < stokes_dim; is++) {
1759 iyb[Range(row0 + is, nf, stokes_dim)] = iy(joker, is);
1760 }
1761
1762 } // End try
1763
1764 catch (const std::exception& e) {
1765#pragma omp critical(iyb_calc_fail)
1766 {
1767 fail_msg = e.what();
1768 failed = true;
1769 }
1770 }
1771}
1772
1774 Vector& iyb,
1775 ArrayOfVector& iyb_aux,
1776 ArrayOfMatrix& diyb_dx,
1777 Matrix& geo_pos_matrix,
1778 const Index& mblock_index,
1779 const Index& atmosphere_dim,
1780 const EnergyLevelMap& nlte_field,
1781 const Index& cloudbox_on,
1782 const Index& stokes_dim,
1783 const ConstVectorView& f_grid,
1784 const ConstMatrixView& sensor_pos,
1785 const ConstMatrixView& sensor_los,
1786 const ConstMatrixView& transmitter_pos,
1787 const ConstMatrixView& mblock_dlos_grid,
1788 const String& iy_unit,
1789 const Agenda& iy_main_agenda,
1790 const Agenda& geo_pos_agenda,
1791 const Index& j_analytical_do,
1792 const ArrayOfRetrievalQuantity& jacobian_quantities,
1793 const ArrayOfArrayOfIndex& jacobian_indices,
1794 const ArrayOfString& iy_aux_vars,
1795 const Verbosity& verbosity) {
1797
1798 // Sizes
1799 const Index nf = f_grid.nelem();
1800 const Index nlos = mblock_dlos_grid.nrows();
1801 const Index niyb = nf * nlos * stokes_dim;
1802 // Set up size of containers for data of 1 measurement block.
1803 // (can not be made below due to parallalisation)
1804 iyb.resize(niyb);
1805 //
1806 if (j_analytical_do) {
1807 diyb_dx.resize(jacobian_indices.nelem());
1808 FOR_ANALYTICAL_JACOBIANS_DO2(diyb_dx[iq].resize(
1809 niyb, jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1);)
1810 } else {
1811 diyb_dx.resize(0);
1812 }
1813 // Assume that geo_pos_agenda returns empty geo_pos.
1814 geo_pos_matrix.resize(nlos, 5);
1815 geo_pos_matrix = NAN;
1816
1817 // For iy_aux we don't know the number of quantities, and we have to store
1818 // all outout
1819 ArrayOfArrayOfMatrix iy_aux_array(nlos);
1820
1821 // We have to make a local copy of the Workspace and the agendas because
1822 // only non-reference types can be declared firstprivate in OpenMP
1823 Workspace l_ws(ws);
1824 Agenda l_iy_main_agenda(iy_main_agenda);
1825 Agenda l_geo_pos_agenda(geo_pos_agenda);
1826
1827 String fail_msg;
1828 bool failed = false;
1829 if (nlos >= arts_omp_get_max_threads() || nlos * 10 >= nf) {
1830 out3 << " Parallelizing los loop (" << nlos << " iterations, " << nf
1831 << " frequencies)\n";
1832
1833 // Start of actual calculations
1834#pragma omp parallel for if (!arts_omp_in_parallel()) \
1835 firstprivate(l_ws, l_iy_main_agenda, l_geo_pos_agenda)
1836 for (Index ilos = 0; ilos < nlos; ilos++) {
1837 // Skip remaining iterations if an error occurred
1838 if (failed) continue;
1839
1840 Ppath ppath;
1841 iyb_calc_body(failed,
1842 fail_msg,
1843 iy_aux_array,
1844 l_ws,
1845 ppath,
1846 iyb,
1847 diyb_dx,
1848 mblock_index,
1849 atmosphere_dim,
1850 nlte_field,
1851 cloudbox_on,
1852 stokes_dim,
1853 sensor_pos,
1854 sensor_los,
1855 transmitter_pos,
1856 mblock_dlos_grid,
1857 iy_unit,
1858 l_iy_main_agenda,
1859 j_analytical_do,
1860 jacobian_quantities,
1861 jacobian_indices,
1862 f_grid,
1863 iy_aux_vars,
1864 ilos,
1865 nf);
1866
1867 // Skip remaining iterations if an error occurred
1868 if (failed) continue;
1869
1870 // Note that this code is found in two places inside the function
1871 Vector geo_pos;
1872 try {
1873 geo_pos_agendaExecute(l_ws, geo_pos, ppath, l_geo_pos_agenda);
1874 if (geo_pos.nelem()) {
1875 ARTS_USER_ERROR_IF (geo_pos.nelem() != 5,
1876 "Wrong size of *geo_pos* obtained from *geo_pos_agenda*.\n"
1877 "The length of *geo_pos* must be zero or five.");
1878
1879 geo_pos_matrix(ilos, joker) = geo_pos;
1880 }
1881 } catch (const std::exception& e) {
1882#pragma omp critical(iyb_calc_fail)
1883 {
1884 fail_msg = e.what();
1885 failed = true;
1886 }
1887 }
1888 }
1889 } else {
1890 out3 << " Not parallelizing los loop (" << nlos << " iterations, " << nf
1891 << " frequencies)\n";
1892
1893 for (Index ilos = 0; ilos < nlos; ilos++) {
1894 // Skip remaining iterations if an error occurred
1895 if (failed) continue;
1896
1897 Ppath ppath;
1898 iyb_calc_body(failed,
1899 fail_msg,
1900 iy_aux_array,
1901 l_ws,
1902 ppath,
1903 iyb,
1904 diyb_dx,
1905 mblock_index,
1906 atmosphere_dim,
1907 nlte_field,
1908 cloudbox_on,
1909 stokes_dim,
1910 sensor_pos,
1911 sensor_los,
1912 transmitter_pos,
1913 mblock_dlos_grid,
1914 iy_unit,
1915 l_iy_main_agenda,
1916 j_analytical_do,
1917 jacobian_quantities,
1918 jacobian_indices,
1919 f_grid,
1920 iy_aux_vars,
1921 ilos,
1922 nf);
1923
1924 // Skip remaining iterations if an error occurred
1925 if (failed) continue;
1926
1927 // Note that this code is found in two places inside the function
1928 Vector geo_pos;
1929 try {
1930 geo_pos_agendaExecute(l_ws, geo_pos, ppath, l_geo_pos_agenda);
1931 if (geo_pos.nelem()) {
1932 ARTS_USER_ERROR_IF (geo_pos.nelem() != 5,
1933 "Wrong size of *geo_pos* obtained from *geo_pos_agenda*.\n"
1934 "The length of *geo_pos* must be zero or five.");
1935
1936 geo_pos_matrix(ilos, joker) = geo_pos;
1937 }
1938 } catch (const std::exception& e) {
1939#pragma omp critical(iyb_calc_fail)
1940 {
1941 fail_msg = e.what();
1942 failed = true;
1943 }
1944 }
1945 }
1946 }
1947
1948 ARTS_USER_ERROR_IF (failed,
1949 "Run-time error in function: iyb_calc\n", fail_msg);
1950
1951 // Compile iyb_aux
1952 //
1953 const Index nq = iy_aux_array[0].nelem();
1954 iyb_aux.resize(nq);
1955 //
1956 for (Index q = 0; q < nq; q++) {
1957 iyb_aux[q].resize(niyb);
1958 //
1959 for (Index ilos = 0; ilos < nlos; ilos++) {
1960 const Index row0 = ilos * nf * stokes_dim;
1961 for (Index iv = 0; iv < nf; iv++) {
1962 const Index row1 = row0 + iv * stokes_dim;
1963 const Index i1 = min(iv, iy_aux_array[ilos][q].nrows() - 1);
1964 for (Index is = 0; is < stokes_dim; is++) {
1965 Index i2 = min(is, iy_aux_array[ilos][q].ncols() - 1);
1966 iyb_aux[q][row1 + is] = iy_aux_array[ilos][q](i1, i2);
1967 }
1968 }
1969 }
1970 }
1971}
1972
1973void iy_transmittance_mult(Tensor3& iy_trans_total,
1974 const ConstTensor3View& iy_trans_old,
1975 const ConstTensor3View& iy_trans_new) {
1976 const Index nf = iy_trans_old.npages();
1977 const Index ns = iy_trans_old.ncols();
1978
1979 ARTS_ASSERT(ns == iy_trans_old.nrows());
1980 ARTS_ASSERT(nf == iy_trans_new.npages());
1981 ARTS_ASSERT(ns == iy_trans_new.nrows());
1982 ARTS_ASSERT(ns == iy_trans_new.ncols());
1983
1984 iy_trans_total.resize(nf, ns, ns);
1985
1986 for (Index iv = 0; iv < nf; iv++) {
1987 mult(iy_trans_total(iv, joker, joker),
1988 iy_trans_old(iv, joker, joker),
1989 iy_trans_new(iv, joker, joker));
1990 }
1991}
1992
1994 const ConstTensor3View& iy_trans,
1995 const ConstMatrixView& iy_old) {
1996 const Index nf = iy_trans.npages();
1997 const Index ns = iy_trans.ncols();
1998
1999 ARTS_ASSERT(ns == iy_trans.nrows());
2000 ARTS_ASSERT(nf == iy_old.nrows());
2001 ARTS_ASSERT(ns == iy_old.ncols());
2002
2003 iy_new.resize(nf, ns);
2004
2005 for (Index iv = 0; iv < nf; iv++) {
2006 mult(iy_new(iv, joker), iy_trans(iv, joker, joker), iy_old(iv, joker));
2007 }
2008}
2009
2010void mirror_los(Vector& los_mirrored,
2011 const ConstVectorView& los,
2012 const Index& atmosphere_dim) {
2013 los_mirrored.resize(2);
2014 //
2015 if (atmosphere_dim == 1) {
2016 los_mirrored[0] = 180 - los[0];
2017 los_mirrored[1] = 180;
2018 } else if (atmosphere_dim == 2) {
2019 los_mirrored[0] = 180 - fabs(los[0]);
2020 if (los[0] >= 0) {
2021 los_mirrored[1] = 180;
2022 } else {
2023 los_mirrored[1] = 0;
2024 }
2025 } else if (atmosphere_dim == 3) {
2026 los_mirrored[0] = 180 - los[0];
2027 los_mirrored[1] = los[1] + 180;
2028 if (los_mirrored[1] > 180) {
2029 los_mirrored[1] -= 360;
2030 }
2031 }
2032}
2033
2035 Numeric& lon,
2036 const Index& atmosphere_dim,
2037 const ConstVectorView& lat_grid,
2038 const ConstVectorView& lat_true,
2039 const ConstVectorView& lon_true,
2040 const ConstVectorView& pos) {
2041 ARTS_ASSERT(pos.nelem() == atmosphere_dim);
2042
2043 if (atmosphere_dim == 1) {
2044 ARTS_ASSERT(lat_true.nelem() == 1);
2045 ARTS_ASSERT(lon_true.nelem() == 1);
2046 //
2047 lat = lat_true[0];
2048 lon = lon_true[0];
2049 }
2050
2051 else if (atmosphere_dim == 2) {
2052 ARTS_ASSERT(lat_true.nelem() == lat_grid.nelem());
2053 ARTS_ASSERT(lon_true.nelem() == lat_grid.nelem());
2054 GridPos gp;
2055 Vector itw(2);
2056 gridpos(gp, lat_grid, pos[1]);
2057 interpweights(itw, gp);
2058 lat = interp(itw, lat_true, gp);
2059 lon = interp(itw, lon_true, gp);
2060 }
2061
2062 else {
2063 lat = pos[1];
2064 lon = pos[2];
2065 }
2066}
2067
2069 Workspace& ws,
2070 ArrayOfTensor3& diy_dx,
2071 ArrayOfTensor3& diy_dpath,
2072 const Index& ns,
2073 const Index& nf,
2074 const Index& np,
2075 const Index& atmosphere_dim,
2076 const Ppath& ppath,
2077 const Vector& ppvar_p,
2078 const Vector& ppvar_t,
2079 const Matrix& ppvar_vmr,
2080 const Index& iy_agenda_call1,
2081 const Tensor3& iy_transmittance,
2082 const Agenda& water_p_eq_agenda,
2083 const ArrayOfRetrievalQuantity& jacobian_quantities,
2084 const ArrayOfIndex& jac_species_i) {
2085 // Weight with iy_transmittance
2086 if (!iy_agenda_call1) {
2087 Matrix X, Y;
2088 //
2090 Y.resize(ns, diy_dpath[iq].npages());
2091 for (Index iv = 0; iv < nf; iv++) {
2092 X = transpose(diy_dpath[iq](joker, iv, joker));
2093 mult(Y, iy_transmittance(iv, joker, joker), X);
2094 diy_dpath[iq](joker, iv, joker) = transpose(Y);
2095 })
2096 }
2097
2098 // Handle abs species retrieval units, both internally and impact on T-jacobian
2099 //
2100 Tensor3 water_p_eq(0, 0, 0);
2101 //
2102 // Conversion for abs species itself
2103 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
2104 // Let x be VMR, and z the selected retrieval unit.
2105 // We have then that diy/dz = diy/dx * dx/dz
2106 //
2107 if (not(jacobian_quantities[iq] == Jacobian::Type::Sensor) and
2108 not(jacobian_quantities[iq] == Jacobian::Special::SurfaceString) and
2109 jac_species_i[iq] >= 0) {
2110 if (jacobian_quantities[iq].Mode() == "vmr") {
2111 }
2112
2113 else if (jacobian_quantities[iq].Mode() == "rel") {
2114 // Here x = vmr*z
2115 for (Index ip = 0; ip < np; ip++) {
2116 diy_dpath[iq](ip, joker, joker) *= ppvar_vmr(jac_species_i[iq], ip);
2117 }
2118 }
2119
2120 else if (jacobian_quantities[iq].Mode() == "nd") {
2121 // Here x = z/nd_tot
2122 for (Index ip = 0; ip < np; ip++) {
2123 diy_dpath[iq](ip, joker, joker) /=
2124 number_density(ppvar_p[ip], ppvar_t[ip]);
2125 }
2126 }
2127
2128 else if (jacobian_quantities[iq].Mode() == "rh") {
2129 // Here x = (p_sat/p) * z
2130 Tensor3 t_data(ppvar_t.nelem(), 1, 1);
2131 t_data(joker, 0, 0) = ppvar_t;
2132 water_p_eq_agendaExecute(ws, water_p_eq, t_data, water_p_eq_agenda);
2133 for (Index ip = 0; ip < np; ip++) {
2134 diy_dpath[iq](ip, joker, joker) *= water_p_eq(ip, 0, 0) / ppvar_p[ip];
2135 }
2136 }
2137
2138 else if (jacobian_quantities[iq].Mode() == "q") {
2139 // Here we use the approximation of x = z/0.622
2140 diy_dpath[iq](joker, joker, joker) /= 0.622;
2141 }
2142
2143 else {
2144 ARTS_ASSERT(0);
2145 }
2146 }
2147 }
2148
2149 // Correction of temperature Jacobian
2150 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
2151 // Let a be unit for abs species, and iy = f(T,a(T))
2152 // We have then that diy/dT = df/dT + df/da*da/dT
2153 // diy_dpath holds already df/dT. Remains is to add
2154 // df/da*da/dT for which abs species having da/dT != 0
2155 // This is only true for "nd" and "rh"
2156 //
2157 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
2158 // Loop abs species, again
2159 for (Index ia = 0; ia < jacobian_quantities.nelem(); ia++) {
2160 if (jac_species_i[ia] >= 0) {
2161 if (jacobian_quantities[ia].Mode() == "nd") {
2162 for (Index ip = 0; ip < np; ip++) {
2163 Matrix ddterm = diy_dpath[ia](ip, joker, joker);
2164 ddterm *= ppvar_vmr(jac_species_i[ia], ip) *
2165 (number_density(ppvar_p[ip], ppvar_t[ip] + 1) -
2166 number_density(ppvar_p[ip], ppvar_t[ip]));
2167 diy_dpath[iq](ip, joker, joker) += ddterm;
2168 }
2169 } else if (jacobian_quantities[ia].Mode() == "rh") {
2170 Tensor3 t_data(ppvar_t.nelem(), 1, 1);
2171 t_data(joker, 0, 0) = ppvar_t;
2172 // Calculate water sat. pressure if not already done
2173 if (water_p_eq.npages() == 0) {
2175 ws, water_p_eq, t_data, water_p_eq_agenda);
2176 }
2177 // Sat.pressure for +1K
2178 Tensor3 water_p_eq1K;
2179 t_data(joker, 0, 0) += 1;
2181 ws, water_p_eq1K, t_data, water_p_eq_agenda);
2182
2183 for (Index ip = 0; ip < np; ip++) {
2184 const Numeric p_eq = water_p_eq(ip, 0, 0);
2185 const Numeric p_eq1K = water_p_eq1K(ip, 0, 0);
2186 Matrix ddterm = diy_dpath[ia](ip, joker, joker);
2187 ddterm *= ppvar_vmr(jac_species_i[ia], ip) *
2188 (ppvar_p[ip] / pow(p_eq, 2.0)) * (p_eq1K - p_eq);
2189 diy_dpath[iq](ip, joker, joker) += ddterm;
2190 }
2191 }
2192 }
2193 }
2194 }
2195 }
2196
2197 // Map to retrieval grids
2199 jacobian_quantities[iq],
2200 diy_dpath[iq],
2201 atmosphere_dim,
2202 ppath,
2203 ppvar_p);)
2204}
2205
2207 Matrix& iy,
2208 ArrayOfTensor3& diy_dx,
2209 Tensor3& ppvar_iy,
2210 const Index& ns,
2211 const Index& np,
2212 const Vector& f_grid,
2213 const Ppath& ppath,
2214 const ArrayOfRetrievalQuantity& jacobian_quantities,
2215 const Index& j_analytical_do,
2216 const String& iy_unit) {
2217 // Determine refractive index to use for the n2 radiance law
2218 Numeric n = 1.0; // First guess is that sensor is in space
2219 //
2220 if (ppath.end_lstep == 0) // If true, sensor inside the atmosphere
2221 {
2222 n = ppath.nreal[np - 1];
2223 }
2224
2225 // Polarisation index variable
2226 ArrayOfIndex i_pol(ns);
2227 for (Index is = 0; is < ns; is++) {
2228 i_pol[is] = is + 1;
2229 }
2230
2231 // Jacobian part (must be converted to Tb before iy for PlanckBT)
2232 //
2233 if (j_analytical_do) {
2235 apply_iy_unit2(diy_dx[iq], iy, iy_unit, f_grid, n, i_pol);)
2236 }
2237
2238 // iy
2239 apply_iy_unit(iy, iy_unit, f_grid, n, i_pol);
2240
2241 // ppvar_iy
2242 for (Index ip = 0; ip < ppath.np; ip++) {
2244 ppvar_iy(joker, joker, ip), iy_unit, f_grid, ppath.nreal[ip], i_pol);
2245 }
2246}
2247
2248void yCalc_mblock_loop_body(bool& failed,
2249 String& fail_msg,
2250 ArrayOfArrayOfVector& iyb_aux_array,
2251 Workspace& ws,
2252 Vector& y,
2253 Vector& y_f,
2254 ArrayOfIndex& y_pol,
2255 Matrix& y_pos,
2256 Matrix& y_los,
2257 Matrix& y_geo,
2258 Matrix& jacobian,
2259 const Index& atmosphere_dim,
2260 const EnergyLevelMap& nlte_field,
2261 const Index& cloudbox_on,
2262 const Index& stokes_dim,
2263 const Vector& f_grid,
2264 const Matrix& sensor_pos,
2265 const Matrix& sensor_los,
2266 const Matrix& transmitter_pos,
2267 const Matrix& mblock_dlos_grid,
2268 const Sparse& sensor_response,
2269 const Vector& sensor_response_f,
2270 const ArrayOfIndex& sensor_response_pol,
2271 const Matrix& sensor_response_dlos,
2272 const String& iy_unit,
2273 const Agenda& iy_main_agenda,
2274 const Agenda& geo_pos_agenda,
2275 const Agenda& jacobian_agenda,
2276 const Index& jacobian_do,
2277 const ArrayOfRetrievalQuantity& jacobian_quantities,
2278 const ArrayOfArrayOfIndex& jacobian_indices,
2279 const ArrayOfString& iy_aux_vars,
2280 const Verbosity& verbosity,
2281 const Index& mblock_index,
2282 const Index& n1y,
2283 const Index& j_analytical_do) {
2284 try {
2285 // Calculate monochromatic pencil beam data for 1 measurement block
2286 //
2287 Vector iyb, iyb_error, yb(n1y);
2288 ArrayOfMatrix diyb_dx;
2289 Matrix geo_pos_matrix;
2290 //
2291 iyb_calc(ws,
2292 iyb,
2293 iyb_aux_array[mblock_index],
2294 diyb_dx,
2295 geo_pos_matrix,
2296 mblock_index,
2297 atmosphere_dim,
2298 nlte_field,
2299 cloudbox_on,
2300 stokes_dim,
2301 f_grid,
2302 sensor_pos,
2303 sensor_los,
2304 transmitter_pos,
2305 mblock_dlos_grid,
2306 iy_unit,
2307 iy_main_agenda,
2308 geo_pos_agenda,
2309 j_analytical_do,
2310 jacobian_quantities,
2311 jacobian_indices,
2312 iy_aux_vars,
2313 verbosity);
2314
2315 // Apply sensor response matrix on iyb, and put into y
2316 //
2317 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
2318 const Index row0 = rowind.get_start();
2319 //
2320 mult(yb, sensor_response, iyb);
2321 //
2322 y[rowind] = yb; // *yb* also used below, as input to jacobian_agenda
2323
2324 // Fill information variables. And search for NaNs in *y*.
2325 //
2326 for (Index i = 0; i < n1y; i++) {
2327 const Index ii = row0 + i;
2329 "One or several NaNs found in *y*.");
2330 y_f[ii] = sensor_response_f[i];
2331 y_pol[ii] = sensor_response_pol[i];
2332 y_pos(ii, joker) = sensor_pos(mblock_index, joker);
2333 y_los(ii, joker) = sensor_los(mblock_index, joker);
2334 y_los(ii, 0) += sensor_response_dlos(i, 0);
2335 if (sensor_response_dlos.ncols() > 1) {
2336 y_los(ii, 1) += sensor_response_dlos(i, 1);
2337 }
2338 }
2339
2340 // Apply sensor response matrix on diyb_dx, and put into jacobian
2341 // (that is, analytical jacobian part)
2342 //
2343 if (j_analytical_do) {
2345 mult(jacobian(rowind,
2346 Range(jacobian_indices[iq][0],
2347 jacobian_indices[iq][1] -
2348 jacobian_indices[iq][0] + 1)),
2349 sensor_response,
2350 diyb_dx[iq]);)
2351 }
2352
2353 // Calculate remaining parts of *jacobian*
2354 //
2355 if (jacobian_do) {
2357 ws, jacobian, mblock_index, iyb, yb, jacobian_agenda);
2358 }
2359
2360 // Handle geo-positioning
2361 if (!std::isnan(geo_pos_matrix(0, 0))) // No data are flagged as NaN
2362 {
2363 // We set geo_pos based on the max value in sensor_response
2364 const Index nfs = f_grid.nelem() * stokes_dim;
2365 for (Index i = 0; i < n1y; i++) {
2366 Index jmax = -1;
2367 Numeric rmax = -99e99;
2368 for (Index j = 0; j < sensor_response.ncols(); j++) {
2369 if (sensor_response(i, j) > rmax) {
2370 rmax = sensor_response(i, j);
2371 jmax = j;
2372 }
2373 }
2374 const auto jhit = Index(floor(jmax / nfs));
2375 y_geo(row0 + i, joker) = geo_pos_matrix(jhit, joker);
2376 }
2377 }
2378 }
2379
2380 catch (const std::exception& e) {
2381#pragma omp critical(yCalc_fail)
2382 {
2383 fail_msg = e.what();
2384 failed = true;
2385 }
2386 }
2387}
2388
2390 const Vector& f_grid,
2391 const Numeric& ze_tref,
2392 const Numeric& k2) {
2393 const Index nf = f_grid.nelem();
2394
2395 ARTS_ASSERT(fac.nelem() == nf);
2396
2397 // Refractive index for water (if needed)
2398 Matrix complex_n(0, 0);
2399 if (k2 <= 0) {
2400 complex_n_water_liebe93(complex_n, f_grid, ze_tref);
2401 }
2402
2403 // Common conversion factor
2404 const Numeric a = 4e18 / (PI * PI * PI * PI);
2405
2406 for (Index iv = 0; iv < nf; iv++) {
2407 // Reference dielectric factor.
2408 Numeric K2;
2409 if (k2 >= 0) {
2410 K2 = k2;
2411 } else {
2412 Complex n(complex_n(iv, 0), complex_n(iv, 1));
2413 Complex n2 = n * n;
2414 Complex K = (n2 - Numeric(1.0)) / (n2 + Numeric(2.0));
2415 Numeric absK = abs(K);
2416 K2 = absK * absK;
2417 }
2418
2419 // Wavelength
2420 Numeric la = SPEED_OF_LIGHT / f_grid[iv];
2421
2422 fac[iv] = a * la * la * la * la / K2;
2423 }
2424}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:39
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Index nrows
Index ncols
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:23295
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:23658
void geo_pos_agendaExecute(Workspace &ws, Vector &geo_pos, const Ppath &ppath, const Agenda &input_agenda)
Definition: auto_md.cc:23214
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:22824
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:23905
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:23520
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:24619
void iy_surface_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, ArrayOfTensor4 &dsurface_rmatrix_dx, 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:23700
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:44
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
A constant view of a Matrix.
Definition: matpackI.h:1014
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Tensor3.
Definition: matpackIII.h:132
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
A constant view of a Tensor4.
Definition: matpackIV.h:133
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
A constant view of a Tensor7.
Definition: matpackVII.h:147
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
EnergyLevelMap InterpToGridPos(Index atmosphere_dim, const ArrayOfGridPos &p, const ArrayOfGridPos &lat, const ArrayOfGridPos &lon) const
The MatrixView class.
Definition: matpackI.h:1125
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
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.
bool IsEmpty() const
Asks if the class is empty.
void MatrixInverseAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Return the matrix inverse at the position.
bool IsRotational(const Index iv=0, const Index iz=0, const Index ia=0) const
False if diagonal element is non-zero in internal Matrix representation.
The range class.
Definition: matpackI.h:165
constexpr Index get_start() const ARTS_NOEXCEPT
Returns the start index of the range.
Definition: matpackI.h:332
The Sparse class.
Definition: matpackII.h:67
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
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)
VectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0)
Get a vectorview to the position.
The Tensor3View class.
Definition: matpackIII.h:239
The Tensor3 class.
Definition: matpackIII.h:339
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:658
The Tensor4View class.
Definition: matpackIV.h:284
The Tensor4 class.
Definition: matpackIV.h:421
The Tensor5 class.
Definition: matpackV.h:506
The Tensor6 class.
Definition: matpackVI.h:1088
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
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
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:560
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition: jacobian.h:571
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:563
#define abs(x)
#define q
#define ns
#define min(a, b)
#define dx
#define max(a, b)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2235
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:290
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
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:57
WindMagnitude
Definition: jacobian.h:58
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:39
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.
void compute_transmission_matrix_and_derivative(Tensor3View T, Tensor4View dT_dx_upper_level, Tensor4View dT_dx_lower_level, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const ArrayOfPropagationMatrix &dupper_level_dx, const ArrayOfPropagationMatrix &dlower_level_dx, const Numeric &dr_dTu, const Numeric &dr_dTl, const Index it, const Index iz, const Index ia)
void compute_transmission_matrix(Tensor3View T, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
#define u
#define v
#define w
#define a
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:71
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:1316
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:2389
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_magnetic_field, const ConstVectorView &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const ConstVectorView &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition: rte.cc:1131
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const ConstTensor3View &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const ConstVectorView &rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const ConstVectorView &f_grid, const String &iy_unit, const ConstTensor3View &surface_props_data, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Verbosity &verbosity)
Determines iy of the "background" of a propgation path.
Definition: rte.cc:728
void get_stepwise_scattersky_source(StokesVector &Sp, ArrayOfStokesVector &dSp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstTensor7View &cloudbox_field, const ConstVectorView &za_grid, const ConstVectorView &aa_grid, const ConstMatrixView &ppath_line_of_sight, const GridPos &ppath_pressure, const Vector &temperature, const Index &atmosphere_dim, const bool &jacobian_do, const Index &t_interp_order)
Calculates the stepwise scattering source terms.
Definition: rte.cc:1433
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
Definition: rte.cc:843
const Numeric TEMP_0_C
void get_stepwise_frequency_grid(VectorView ppath_f_grid, const ConstVectorView &f_grid, const ConstVectorView &ppath_wind, const ConstVectorView &ppath_line_of_sight, const Numeric &rte_alonglos_v, const Index &atmosphere_dim)
Inverse of get_stepwise_f_partials.
Definition: rte.cc:1262
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:2034
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:2248
void iy_transmittance_mult(Tensor3 &iy_trans_total, const ConstTensor3View &iy_trans_old, const ConstTensor3View &iy_trans_new)
Multiplicates iy_transmittance with transmissions.
Definition: rte.cc:1973
void 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 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:1657
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 rtmethods_unit_conversion(Matrix &iy, ArrayOfTensor3 &diy_dx, Tensor3 &ppvar_iy, const Index &ns, const Index &np, const Vector &f_grid, const Ppath &ppath, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &j_analytical_do, const String &iy_unit)
This function handles the unit conversion to be done at the end of some radiative transfer WSMs.
Definition: rte.cc:2206
const Numeric SPEED_OF_LIGHT
void get_stepwise_transmission_matrix(Tensor3View cumulative_transmission, Tensor3View T, Tensor4View dT_close_dx, Tensor4View dT_far_dx, const ConstTensor3View &cumulative_transmission_close, const PropagationMatrix &K_close, const PropagationMatrix &K_far, const ArrayOfPropagationMatrix &dK_close_dx, const ArrayOfPropagationMatrix &dK_far_dx, const Numeric &ppath_distance, const bool &first_level, const Numeric &dppath_distance_dT_HSE_close, const Numeric &dppath_distance_dT_HSE_far, const Index &temperature_derivative_position_if_hse_is_active)
Computes layer transmission matrix and cumulative transmission.
Definition: rte.cc:1601
void bending_angle1d(Numeric &alpha, const Ppath &ppath)
Calculates the bending angle for a 1D atmosphere.
Definition: rte.cc:266
void get_stepwise_effective_source(MatrixView J, Tensor3View dJ_dx, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK_dx, const ArrayOfStokesVector &da_dx, const ArrayOfStokesVector &dS_dx, const ConstVectorView &B, const ConstVectorView &dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Gets the effective source at propagation path point.
Definition: rte.cc:1168
void rtmethods_jacobian_finalisation(Workspace &ws, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &atmosphere_dim, const Ppath &ppath, const Vector &ppvar_p, const Vector &ppvar_t, const Matrix &ppvar_vmr, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &jac_species_i)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
Definition: rte.cc:2068
void 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:2010
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:1773
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:1282
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.
Structure to store a grid position.
Definition: interpolation.h:73
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Matrix los
Line-of-sight at each ppath point.
Definition: ppath.h:66
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath.h:86
String background
Radiative background.
Definition: ppath.h:56
Index np
Number of points describing the ppath.
Definition: ppath.h:52
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath.h:64
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath.h:84
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath.h:76
Vector start_pos
Start position.
Definition: ppath.h:58
Vector lstep
The length between ppath points.
Definition: ppath.h:70
Numeric start_lstep
Length between sensor and atmospheric boundary.
Definition: ppath.h:62
Numeric constant
The propagation path constant (only used for 1D)
Definition: ppath.h:54
Vector r
Radius of each ppath point.
Definition: ppath.h:68
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath.h:82
Index dim
Atmospheric dimensionality.
Definition: ppath.h:50
Vector end_pos
End position.
Definition: ppath.h:72
Vector start_los
Start line-of-sight.
Definition: ppath.h:60
Vector nreal
The real part of the refractive index at each path position.
Definition: ppath.h:78
Vector end_los
End line-of-sight.
Definition: ppath.h:74
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.