ARTS 2.5.11 (git: 725533f0)
m_jacobian.cc
Go to the documentation of this file.
1
9/*===========================================================================
10 === External declarations
11 ===========================================================================*/
12
13#include <cmath>
14#include <string>
15#include "absorption.h"
16#include "arts.h"
17#include "auto_md.h"
18#include "check_input.h"
19#include "cloudbox.h"
20#include "arts_conversions.h"
21#include "interp.h"
22#include "jacobian.h"
23#include "m_xml.h"
24#include "math_funcs.h"
25#include "matpack_math.h"
26#include "messages.h"
27#include "physics_funcs.h"
28#include "rte.h"
29
30/*===========================================================================
31 === The methods, with general methods first followed by the Add/Calc method
32 === pairs for each retrieval quantity.
33 ===========================================================================*/
34
35//----------------------------------------------------------------------------
36// Basic methods:
37//----------------------------------------------------------------------------
38
39/* Workspace method: Doxygen documentation will be auto-generated */
40void jacobianCalcDoNothing(Matrix& jacobian _U_,
41 const Index& mblock_index _U_,
42 const Vector& iyb _U_,
43 const Vector& yb _U_,
44 const Verbosity&) {
45 /* Nothing to do here for the analytical case, this function just exists
46 to satisfy the required inputs and outputs of the jacobian_agenda */
47}
48
49/* Workspace method: Doxygen documentation will be auto-generated */
51 Index& jacobian_do,
52 Agenda& jacobian_agenda,
53 const ArrayOfRetrievalQuantity& jacobian_quantities,
54 const Verbosity& verbosity) {
55 // Make sure that the array is not empty
56 if (jacobian_quantities.empty())
57 throw runtime_error(
58 "No retrieval quantities has been added to *jacobian_quantities*.");
59
60 jacobian_agenda.check(ws_in, verbosity);
61 jacobian_do = 1;
62}
63
64/* Workspace method: Doxygen documentation will be auto-generated */
66 ArrayOfRetrievalQuantity& jacobian_quantities,
67 Agenda& jacobian_agenda,
68 const Verbosity&) {
69 jacobian_quantities.resize(0);
70 jacobian_agenda = Agenda{ws};
71 jacobian_agenda.set_name("jacobian_agenda");
72}
73
74/* Workspace method: Doxygen documentation will be auto-generated */
76 Index& jacobian_do,
77 Agenda& jacobian_agenda,
78 ArrayOfRetrievalQuantity& jacobian_quantities,
79 const Verbosity& verbosity) {
80 jacobian_do = 0;
81 jacobianInit(ws, jacobian_quantities, jacobian_agenda, verbosity);
82}
83
84//----------------------------------------------------------------------------
85// Absorption species:
86//----------------------------------------------------------------------------
87
88/* Workspace method: Doxygen documentation will be auto-generated */
91 Agenda& jacobian_agenda,
92 const Index& atmosphere_dim,
93 const Vector& p_grid,
94 const Vector& lat_grid,
95 const Vector& lon_grid,
96 const Vector& rq_p_grid,
97 const Vector& rq_lat_grid,
98 const Vector& rq_lon_grid,
99 const String& species,
100 const String& mode,
101 const Index& for_species_tag,
102 const Verbosity& verbosity) {
105
107 if (not for_species_tag) {
108 ArrayOfSpeciesTag test(species);
109 if (test.nelem() not_eq 1)
110 throw std::runtime_error(
111 "Trying to add a species as a species tag of multiple species.\n"
112 "This is not supported. Please give just a single species instead.\n"
113 "Otherwise consider if you intended for_species_tag to be evaluated true.\n");
114 qi = QuantumIdentifier(test[0].Isotopologue());
115 }
116
117 // Check that this species is not already included in the jacobian.
118 for (Index it = 0; it < jq.nelem(); it++) {
119 ARTS_USER_ERROR_IF (jq[it] == Jacobian::Special::ArrayOfSpeciesTagVMR && jq[it].Subtag() == species,
120 "The gas species:\n", species, "\nis already included in *jacobian_quantities*.")
121 ARTS_USER_ERROR_IF (jq[it] == Jacobian::Line::VMR and jq[it].QuantumIdentity() == qi,
122 "The atmospheric species of:\n", species, "\nis already included in *jacobian_quantities*\n"
123 "as: ", jq[it])
124 }
125
126 // Check retrieval grids, here we just check the length of the grids
127 // vs. the atmosphere dimension
128 ArrayOfVector grids(atmosphere_dim);
129 {
130 ostringstream os;
131 if (!check_retrieval_grids(grids,
132 os,
133 p_grid,
134 lat_grid,
135 lon_grid,
136 rq_p_grid,
137 rq_lat_grid,
138 rq_lon_grid,
139 "retrieval pressure grid",
140 "retrieval latitude grid",
141 "retrievallongitude_grid",
142 atmosphere_dim))
143 throw runtime_error(os.str());
144 }
145
146 // Check that mode is correct
147 if (mode != "vmr" && mode != "nd" && mode != "rel" && mode != "rh" &&
148 mode != "q") {
149 throw runtime_error(
150 "The retrieval mode can only be \"vmr\", \"nd\", "
151 "\"rel\", \"rh\" or \"q\".");
152 }
153 if ((mode == "rh" || mode == "q") && species.substr(0, 3) != "H2O") {
154 throw runtime_error(
155 "Retrieval modes \"rh\" and \"q\" can only be applied "
156 "on species starting with H2O.");
157 }
158
159 // Create the new retrieval quantity
161 rq.Subtag(species);
162 rq.Mode(mode);
163 rq.Grids(grids);
164 if (for_species_tag == 0) {
165 rq.Target(Jacobian::Target(Jacobian::Line::VMR, qi, qi.Species()));
166 } else {
167 ArrayOfSpeciesTag test(species);
168 rq.Target(Jacobian::Target(Jacobian::Special::ArrayOfSpeciesTagVMR, test));
169 }
170 rq.Target().perturbation = 0.001;
171
172 // Add it to the *jacobian_quantities*
173 jq.push_back(rq);
174
175 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
176}
177
178//----------------------------------------------------------------------------
179// Frequency shift
180//----------------------------------------------------------------------------
181
182/* Workspace method: Doxygen documentation will be auto-generated */
184 ArrayOfRetrievalQuantity& jacobian_quantities,
185 Agenda& jacobian_agenda,
186 const Vector& f_grid,
187 const Numeric& df,
188 const Verbosity&) {
189 // Check that this jacobian type is not already included.
190 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
191 if (jacobian_quantities[it] == Jacobian::Sensor::FrequencyShift) {
192 ostringstream os;
193 os << "Fit of frequency shift is already included in\n"
194 << "*jacobian_quantities*.";
195 throw runtime_error(os.str());
196 }
197 }
198
199 // Checks of frequencies
200 if (df <= 0) throw runtime_error("The argument *df* must be > 0.");
201 if (df > 1e6)
202 throw runtime_error("The argument *df* is not allowed to exceed 1 MHz.");
203 const Index nf = f_grid.nelem();
204 if (nf < 2)
205 throw runtime_error(
206 "Frequency shifts and *f_grid* of length 1 can "
207 "not be combined.");
208 const Numeric maxdf = f_grid[nf - 1] - f_grid[nf - 2];
209 if (df > maxdf) {
210 ostringstream os;
211 os << "The value of *df* is too big with respect to spacing of "
212 << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
213 << "between the two last elements of *f_grid*.\n"
214 << "This spacing is : " << maxdf / 1e3 << " kHz\n"
215 << "The value of df is: " << df / 1e3 << " kHz";
216 throw runtime_error(os.str());
217 }
218
219 // Create the new retrieval quantity
221 rq.Target() = Jacobian::Target(Jacobian::Sensor::FrequencyShift);
222 rq.Mode("");
223 rq.Target().perturbation = df;
224
225 // Dummy vector of length 1
226 Vector grid(1, 0);
227 ArrayOfVector grids(1, grid);
228 rq.Grids(grids);
229
230 // Add it to the *jacobian_quantities*
231 jacobian_quantities.push_back(rq);
232
233 // Add corresponding calculation method to the jacobian agenda
234 jacobian_agenda.append("jacobianCalcFreqShift", "");
235}
236
237/* Workspace method: Doxygen documentation will be auto-generated */
238void jacobianCalcFreqShift(Matrix& jacobian,
239 const Index& mblock_index,
240 const Vector& iyb,
241 const Vector& yb,
242 const Index& stokes_dim,
243 const Vector& f_grid,
244 const Matrix& mblock_dlos,
245 const Sparse& sensor_response,
246 const ArrayOfRetrievalQuantity& jacobian_quantities,
247 const Verbosity&) {
248 // Set some useful (and needed) variables.
250 ArrayOfIndex ji;
251
252 // Find the retrieval quantity related to this method.
253 // This works since the combined MainTag and Subtag is individual.
254 bool found = false;
255 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
256 if (jacobian_quantities[n] == Jacobian::Sensor::FrequencyShift) {
257 bool any_affine;
258 ArrayOfArrayOfIndex jacobian_indices;
260 jacobian_indices, any_affine, jacobian_quantities, true);
261 //
262 found = true;
263 rq = jacobian_quantities[n];
264 ji = jacobian_indices[n];
265 }
266 }
267 if (!found) {
268 throw runtime_error(
269 "There is no such frequency retrieval quantity defined.\n");
270 }
271
272 // Check that sensor_response is consistent with yb and iyb
273 //
274 if (sensor_response.nrows() != yb.nelem())
275 throw runtime_error("Mismatch in size between *sensor_response* and *yb*.");
276 if (sensor_response.ncols() != iyb.nelem())
277 throw runtime_error(
278 "Mismatch in size between *sensor_response* and *iyb*.");
279
280 // Get disturbed (part of) y
281 //
282 const Index n1y = sensor_response.nrows();
283 Vector dy(n1y);
284 {
285 const Index nf2 = f_grid.nelem();
286 const Index nlos2 = mblock_dlos.nrows();
287 const Index niyb = nf2 * nlos2 * stokes_dim;
288
289 // Interpolation weights
290 //
291 constexpr Index porder = 3;
292 //
293 Vector fg_new = f_grid, iyb2(niyb);
294 fg_new += rq.Target().perturbation;
295 //
296 const auto lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<porder>>(fg_new, f_grid);
297 const auto itw = interpweights(lag);
298
299 // Do interpolation
300 for (Index ilos = 0; ilos < nlos2; ilos++) {
301 const Index row0 = ilos * nf2 * stokes_dim;
302
303 for (Index is = 0; is < stokes_dim; is++) {
304 iyb2[Range(row0 + is, nf2, stokes_dim)] = reinterp(iyb[Range(row0 + is, nf2, stokes_dim)], itw, lag);
305 }
306 }
307
308 // Determine difference
309 //
310 mult(dy, sensor_response, iyb2);
311 //
312 for (Index i = 0; i < n1y; i++) {
313 dy[i] = (dy[i] - yb[i]) / rq.Target().perturbation;
314 }
315 }
316
317 //--- Set jacobian ---
318 ARTS_ASSERT(rq.Grids()[0].nelem() == 1);
319 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
320 jacobian(rowind, ji[0]) = dy;
321}
322
323//----------------------------------------------------------------------------
324// Frequency stretch
325//----------------------------------------------------------------------------
326
327/* Workspace method: Doxygen documentation will be auto-generated */
329 ArrayOfRetrievalQuantity& jacobian_quantities,
330 Agenda& jacobian_agenda,
331 const Vector& f_grid,
332 const Numeric& df,
333 const Verbosity&) {
334 // Check that this jacobian type is not already included.
335 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
336 if (jacobian_quantities[it] == Jacobian::Sensor::FrequencyStretch) {
337 ostringstream os;
338 os << "Fit of frequency stretch is already included in\n"
339 << "*jacobian_quantities*.";
340 throw runtime_error(os.str());
341 }
342 }
343
344 // Checks of df
345 if (df <= 0) throw runtime_error("The argument *df* must be > 0.");
346 if (df > 1e6)
347 throw runtime_error("The argument *df* is not allowed to exceed 1 MHz.");
348 const Index nf = f_grid.nelem();
349 const Numeric maxdf = f_grid[nf - 1] - f_grid[nf - 2];
350 if (df > maxdf) {
351 ostringstream os;
352 os << "The value of *df* is too big with respect to spacing of "
353 << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
354 << "between the two last elements of *f_grid*.\n"
355 << "This spacing is : " << maxdf / 1e3 << " kHz\n"
356 << "The value of df is: " << df / 1e3 << " kHz";
357 throw runtime_error(os.str());
358 }
359
360 // Create the new retrieval quantity
362 rq.Target() = Jacobian::Target(Jacobian::Sensor::FrequencyStretch);
363 rq.Mode("");
364 rq.Target().perturbation=df;
365
366 // Dummy vector of length 1
367 Vector grid(1, 0);
368 ArrayOfVector grids(1, grid);
369 rq.Grids(grids);
370
371 // Add it to the *jacobian_quantities*
372 jacobian_quantities.push_back(rq);
373
374 // Add corresponding calculation method to the jacobian agenda
375 jacobian_agenda.append("jacobianCalcFreqStretch", "");
376}
377
378/* Workspace method: Doxygen documentation will be auto-generated */
380 Matrix& jacobian,
381 const Index& mblock_index,
382 const Vector& iyb,
383 const Vector& yb,
384 const Index& stokes_dim,
385 const Vector& f_grid,
386 const Matrix& mblock_dlos,
387 const Sparse& sensor_response,
388 const ArrayOfIndex& sensor_response_pol_grid,
389 const Vector& sensor_response_f_grid,
390 const Matrix& sensor_response_dlos_grid,
391 const ArrayOfRetrievalQuantity& jacobian_quantities,
392 const Verbosity&) {
393 // The code here is close to identical to the one for Shift. The main
394 // difference is that dy is weighted with poly_order 1 basis function.
395
396 // Set some useful (and needed) variables.
398 ArrayOfIndex ji;
399
400 // Find the retrieval quantity related to this method.
401 // This works since the combined MainTag and Subtag is individual.
402 bool found = false;
403 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
404 if (jacobian_quantities[n] == Jacobian::Sensor::FrequencyStretch) {
405 bool any_affine;
406 ArrayOfArrayOfIndex jacobian_indices;
408 jacobian_indices, any_affine, jacobian_quantities, true);
409 //
410 found = true;
411 rq = jacobian_quantities[n];
412 ji = jacobian_indices[n];
413 }
414 }
415 if (!found) {
416 throw runtime_error(
417 "There is no such frequency retrieval quantity defined.\n");
418 }
419
420 // Check that sensor_response is consistent with yb and iyb
421 //
422 if (sensor_response.nrows() != yb.nelem())
423 throw runtime_error("Mismatch in size between *sensor_response* and *yb*.");
424 if (sensor_response.ncols() != iyb.nelem())
425 throw runtime_error(
426 "Mismatch in size between *sensor_response* and *iyb*.");
427
428 // Get disturbed (part of) y
429 //
430 const Index n1y = sensor_response.nrows();
431 Vector dy(n1y);
432 {
433 const Index nf2 = f_grid.nelem();
434 const Index nlos2 = mblock_dlos.nrows();
435 const Index niyb = nf2 * nlos2 * stokes_dim;
436
437 // Interpolation weights
438 //
439 constexpr Index porder = 3;
440 //
441 Vector fg_new = f_grid, iyb2(niyb);
442 //
443 fg_new += rq.Target().perturbation;
444 //
445 const auto lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<porder>>(fg_new, f_grid);
446 const auto itw = interpweights(lag);
447
448 // Do interpolation
449 for (Index ilos = 0; ilos < nlos2; ilos++) {
450 const Index row0 = ilos * nf2 * stokes_dim;
451
452 for (Index is = 0; is < stokes_dim; is++) {
453 iyb2[Range(row0 + is, nf2, stokes_dim)] = reinterp(iyb[Range(row0 + is, nf2, stokes_dim)], itw, lag);
454 }
455 }
456
457 // Determine difference
458 //
459 mult(dy, sensor_response, iyb2);
460 //
461 for (Index i = 0; i < n1y; i++) {
462 dy[i] = (dy[i] - yb[i]) / rq.Target().perturbation;
463 }
464
465 // dy above corresponds now to shift. Convert to stretch:
466 //
467 Vector w;
468 polynomial_basis_func(w, sensor_response_f_grid, 1);
469 //
470 const Index nf = sensor_response_f_grid.nelem();
471 const Index npol = sensor_response_pol_grid.nelem();
472 const Index nlos = sensor_response_dlos_grid.nrows();
473 //
474 for (Index l = 0; l < nlos; l++) {
475 for (Index f = 0; f < nf; f++) {
476 const Index row1 = (l * nf + f) * npol;
477 for (Index p = 0; p < npol; p++) {
478 dy[row1 + p] *= w[f];
479 }
480 }
481 }
482 }
483
484 //--- Set jacobians ---
485 ARTS_ASSERT(rq.Grids()[0].nelem() == 1);
486 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
487 jacobian(rowind, ji[0]) = dy;
488}
489
490//----------------------------------------------------------------------------
491// Pointing:
492//----------------------------------------------------------------------------
493
494/* Workspace method: Doxygen documentation will be auto-generated */
496 ArrayOfRetrievalQuantity& jacobian_quantities,
497 Agenda& jacobian_agenda,
498 const Matrix& sensor_pos,
499 const ArrayOfTime& sensor_time,
500 const Index& poly_order,
501 const String& calcmode,
502 const Numeric& dza,
503 const Verbosity&) {
504 // Check that poly_order is -1 or positive
505 if (poly_order < -1)
506 throw runtime_error(
507 "The polynomial order has to be positive or -1 for gitter.");
508
509 // Check that this jacobian type is not already included.
510 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
511 if (jacobian_quantities[it].Target().isPointing()) {
512 ostringstream os;
513 os << "Fit of zenith angle pointing off-set is already included in\n"
514 << "*jacobian_quantities*.";
515 throw runtime_error(os.str());
516 }
517 }
518
519 // Checks of dza
520 if (dza <= 0) throw runtime_error("The argument *dza* must be > 0.");
521 if (dza > 0.1)
522 throw runtime_error("The argument *dza* is not allowed to exceed 0.1 deg.");
523
524 // Check that sensor_time is consistent with sensor_pos
525 if (sensor_time.nelem() != sensor_pos.nrows()) {
526 ostringstream os;
527 os << "The WSV *sensor_time* must be defined for every "
528 << "measurement block.\n";
529 throw runtime_error(os.str());
530 }
531
532 // Do not allow that *poly_order* is not too large compared to *sensor_time*
533 if (poly_order > sensor_time.nelem() - 1) {
534 throw runtime_error(
535 "The polynomial order can not be >= length of *sensor_time*.");
536 }
537
538 // Create the new retrieval quantity
540 if (calcmode == "recalc") {
541 rq.Target() = Jacobian::Target(Jacobian::Sensor::PointingZenithRecalc);
542 jacobian_agenda.append("jacobianCalcPointingZaRecalc", "");
543 } else if (calcmode == "interp") {
544 rq.Target() = Jacobian::Target(Jacobian::Sensor::PointingZenithInterp);
545 jacobian_agenda.append("jacobianCalcPointingZaInterp", "");
546 } else
547 throw runtime_error(
548 R"(Possible choices for *calcmode* are "recalc" and "interp".)");
549 rq.Target().perturbation = dza;
550
551 // To store the value or the polynomial order, create a vector with length
552 // poly_order+1, in case of gitter set the size of the grid vector to be the
553 // number of measurement blocks, all elements set to -1.
554 Vector grid=uniform_grid(0, poly_order + 1, 1);
555 if (poly_order == -1) {
556 grid.resize(sensor_pos.nrows());
557 grid = -1.0;
558 }
559 ArrayOfVector grids(1, grid);
560 rq.Grids(grids);
561
562 // Add it to the *jacobian_quantities*
563 jacobian_quantities.push_back(rq);
564}
565
566/* Workspace method: Doxygen documentation will be auto-generated */
568 Matrix& jacobian,
569 const Index& mblock_index,
570 const Vector& iyb,
571 const Vector& yb _U_,
572 const Index& stokes_dim,
573 const Vector& f_grid,
574 const Matrix& DEBUG_ONLY(sensor_los),
575 const Matrix& mblock_dlos,
576 const Sparse& sensor_response,
577 const ArrayOfTime& sensor_time,
578 const ArrayOfRetrievalQuantity& jacobian_quantities,
579 const Verbosity&) {
580 if (mblock_dlos.nrows() < 2)
581 throw runtime_error(
582 "The method demands that *mblock_dlos* has "
583 "more than one row.");
584
585 if (!(is_increasing(mblock_dlos(joker, 0)) ||
586 is_decreasing(mblock_dlos(joker, 0))))
587 throw runtime_error(
588 "The method demands that the zenith angles in "
589 "*mblock_dlos* are sorted (increasing or decreasing).");
590
591 // Set some useful variables.
593 ArrayOfIndex ji;
594
595 // Find the retrieval quantity related to this method.
596 // This works since the combined MainTag and Subtag is individual.
597 bool found = false;
598 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
599 if (jacobian_quantities[n] == Jacobian::Sensor::PointingZenithInterp) {
600 bool any_affine;
601 ArrayOfArrayOfIndex jacobian_indices;
603 jacobian_indices, any_affine, jacobian_quantities, true);
604 //
605 found = true;
606 rq = jacobian_quantities[n];
607 ji = jacobian_indices[n];
608 }
609 }
610 if (!found) {
611 throw runtime_error(
612 "There is no such pointing retrieval quantity defined.\n");
613 }
614
615 // Get "dy", by inter/extra-polation of existing iyb
616 //
617 const Index n1y = sensor_response.nrows();
618 Vector dy(n1y);
619 {
620 // Sizes
621 const Index nf = f_grid.nelem();
622 const Index nza = mblock_dlos.nrows();
623
624 // Shifted zenith angles
625 Vector za1{mblock_dlos(joker, 0)};
626 za1 -= rq.Target().perturbation;
627 Vector za2{mblock_dlos(joker, 0)};
628 za2 += rq.Target().perturbation;
629
630 // Find interpolation weights
631 ArrayOfGridPos gp1(nza), gp2(nza);
632 gridpos(
633 gp1, mblock_dlos(joker, 0), za1, 1e6); // Note huge extrapolation!
634 gridpos(
635 gp2, mblock_dlos(joker, 0), za2, 1e6); // Note huge extrapolation!
636 Matrix itw1(nza, 2), itw2(nza, 2);
637 interpweights(itw1, gp1);
638 interpweights(itw2, gp2);
639
640 // Make interpolation (for all azimuth angles, frequencies and Stokes)
641 //
642 Vector iyb1(iyb.nelem()), iyb2(iyb.nelem());
643 //
644 for (Index iza = 0; iza < nza; iza++) {
645 for (Index iv = 0; iv < nf; iv++) {
646 for (Index is = 0; is < stokes_dim; is++) {
647 const Range r(iv * stokes_dim + is, nza, nf * stokes_dim);
648 interp(iyb1[r], itw1, iyb[r], gp1);
649 interp(iyb2[r], itw2, iyb[r], gp2);
650 }
651 }
652 }
653
654 // Apply sensor and take difference
655 //
656 Vector y1(n1y), y2(n1y);
657 mult(y1, sensor_response, iyb1);
658 mult(y2, sensor_response, iyb2);
659 //
660 for (Index i = 0; i < n1y; i++) {
661 dy[i] = (y2[i] - y1[i]) / (2 * rq.Target().perturbation);
662 }
663 }
664
665 //--- Create jacobians ---
666
667 const Index lg = rq.Grids()[0].nelem();
668 const Index it = ji[0];
669 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
670 const Index row0 = rowind.offset;
671
672 // Handle pointing "jitter" seperately
673 if (rq.Grids()[0][0] == -1) // Not all values are set here,
674 { // but should already have been
675 ARTS_ASSERT(lg == sensor_los.nrows()); // set to 0
676 ARTS_ASSERT(rq.Grids()[0][mblock_index] == -1);
677 jacobian(rowind, it + mblock_index) = dy;
678 }
679
680 // Polynomial representation
681 else {
682 Vector w;
683 for (Index c = 0; c < lg; c++) {
684 ARTS_ASSERT(Numeric(c) == rq.Grids()[0][c]);
685 //
686 polynomial_basis_func(w, time_vector(sensor_time), c);
687 //
688 for (Index i = 0; i < n1y; i++) {
689 jacobian(row0 + i, it + c) = w[mblock_index] * dy[i];
690 }
691 }
692 }
693}
694
695/* Workspace method: Doxygen documentation will be auto-generated */
697 Workspace& ws,
698 Matrix& jacobian,
699 const Index& mblock_index,
700 const Vector& iyb _U_,
701 const Vector& yb,
702 const Index& atmosphere_dim,
703 const EnergyLevelMap& nlte_field,
704 const Index& cloudbox_on,
705 const Index& stokes_dim,
706 const Vector& f_grid,
707 const Matrix& sensor_pos,
708 const Matrix& sensor_los,
709 const Matrix& transmitter_pos,
710 const Matrix& mblock_dlos,
711 const Sparse& sensor_response,
712 const ArrayOfTime& sensor_time,
713 const String& iy_unit,
714 const Agenda& iy_main_agenda,
715 const ArrayOfRetrievalQuantity& jacobian_quantities,
716 const Verbosity& verbosity) {
717 // Set some useful variables.
719 ArrayOfIndex ji;
720
721 // Find the retrieval quantity related to this method.
722 // This works since the combined MainTag and Subtag is individual.
723 bool found = false;
724 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
725 if (jacobian_quantities[n] == Jacobian::Sensor::PointingZenithRecalc) {
726 bool any_affine;
727 ArrayOfArrayOfIndex jacobian_indices;
729 jacobian_indices, any_affine, jacobian_quantities, true);
730 //
731 found = true;
732 rq = jacobian_quantities[n];
733 ji = jacobian_indices[n];
734 }
735 }
736 if (!found) {
737 throw runtime_error(
738 "There is no such pointing retrieval quantity defined.\n");
739 }
740
741 // Get "dy", by calling iyb_calc with shifted sensor_los.
742 //
743 const Index n1y = sensor_response.nrows();
744 Vector dy(n1y);
745 {
746 Vector iyb2;
747 Matrix los = sensor_los;
748 Matrix geo_pos;
749 ArrayOfVector iyb_aux;
750 ArrayOfMatrix diyb_dx;
751
752 los(joker, 0) += rq.Target().perturbation;
753
754 iyb_calc(ws,
755 iyb2,
756 iyb_aux,
757 diyb_dx,
758 geo_pos,
759 mblock_index,
760 atmosphere_dim,
761 nlte_field,
762 cloudbox_on,
763 stokes_dim,
764 f_grid,
765 sensor_pos,
766 los,
767 transmitter_pos,
768 mblock_dlos,
769 iy_unit,
770 iy_main_agenda,
771 0,
775 verbosity);
776
777 // Apply sensor and take difference
778 //
779 mult(dy, sensor_response, iyb2);
780 //
781 for (Index i = 0; i < n1y; i++) {
782 dy[i] = (dy[i] - yb[i]) / rq.Target().perturbation;
783 }
784 }
785
786 //--- Create jacobians ---
787
788 const Index lg = rq.Grids()[0].nelem();
789 const Index it = ji[0];
790 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
791 const Index row0 = rowind.offset;
792
793 // Handle "jitter" seperately
794 if (rq.Grids()[0][0] == -1) // Not all values are set here,
795 { // but should already have been
796 ARTS_ASSERT(lg == sensor_los.nrows()); // set to 0
797 ARTS_ASSERT(rq.Grids()[0][mblock_index] == -1);
798 jacobian(rowind, it + mblock_index) = dy;
799 }
800
801 // Polynomial representation
802 else {
803 Vector w;
804 for (Index c = 0; c < lg; c++) {
805 ARTS_ASSERT(Numeric(c) == rq.Grids()[0][c]);
806 //
807 polynomial_basis_func(w, time_vector(sensor_time), c);
808 //
809 for (Index i = 0; i < n1y; i++) {
810 jacobian(row0 + i, it + c) = w[mblock_index] * dy[i];
811 }
812 }
813 }
814}
815
816//----------------------------------------------------------------------------
817// Polynomial baseline fits:
818//----------------------------------------------------------------------------
819
820/* Workspace method: Doxygen documentation will be auto-generated */
823 Agenda& jacobian_agenda,
824 const ArrayOfIndex& sensor_response_pol_grid,
825 const Matrix& sensor_response_dlos_grid,
826 const Matrix& sensor_pos,
827 const Index& poly_order,
828 const Index& no_pol_variation,
829 const Index& no_los_variation,
830 const Index& no_mblock_variation,
831 const Verbosity&) {
832 // Check that poly_order is >= 0
833 if (poly_order < 0)
834 throw runtime_error("The polynomial order has to be >= 0.");
835
836 // Check that polyfit is not already included in the jacobian.
837 for (Index it = 0; it < jq.nelem(); it++) {
838 if (jq[it] == Jacobian::Sensor::Polyfit) {
839 ostringstream os;
840 os << "Polynomial baseline fit is already included in\n"
841 << "*jacobian_quantities*.";
842 throw runtime_error(os.str());
843 }
844 }
845
846 // "Grids"
847 //
848 // Grid dimensions correspond here to
849 // 1: polynomial order
850 // 2: polarisation
851 // 3: viewing direction
852 // 4: measurement block
853 //
854 ArrayOfVector grids(4);
855 //
856 if (no_pol_variation)
857 grids[1] = Vector(1, 1);
858 else
859 grids[1] = uniform_grid(0, sensor_response_pol_grid.nelem(), 1);
860 if (no_los_variation)
861 grids[2] = Vector(1, 1);
862 else
863 grids[2] = uniform_grid(0, sensor_response_dlos_grid.nrows(), 1);
864 if (no_mblock_variation)
865 grids[3] = Vector(1, 1);
866 else
867 grids[3] = uniform_grid(0, sensor_pos.nrows(), 1);
868
869 // Create the new retrieval quantity
871 rq.Target() = Jacobian::Target(Jacobian::Sensor::Polyfit);
872 rq.Mode("");
873 rq.Target().perturbation = 0;
874
875 // Each polynomial coeff. is treated as a retrieval quantity
876 //
877 for (Index i = 0; i <= poly_order; i++) {
878 ostringstream sstr;
879 sstr << "Coefficient " << i;
880 rq.Subtag(sstr.str());
881
882 // Grid is a scalar, use polynomial coeff.
883 grids[0] = Vector(1, (Numeric)i);
884 rq.Grids(grids);
885
886 // Add it to the *jacobian_quantities*
887 jq.push_back(rq);
888
889 // Add pointing method to the jacobian agenda
890 jacobian_agenda.append("jacobianCalcPolyfit", i);
891 }
892}
893
894/* Workspace method: Doxygen documentation will be auto-generated */
895void jacobianCalcPolyfit(Matrix& jacobian,
896 const Index& mblock_index,
897 const Vector& iyb _U_,
898 const Vector& yb _U_,
899 const Sparse& sensor_response,
900 const ArrayOfIndex& sensor_response_pol_grid,
901 const Vector& sensor_response_f_grid,
902 const Matrix& sensor_response_dlos_grid,
903 const ArrayOfRetrievalQuantity& jacobian_quantities,
904 const Index& poly_coeff,
905 const Verbosity&) {
906 // Find the retrieval quantity related to this method
908 ArrayOfIndex ji;
909 bool found = false;
910 Index iq;
911 ostringstream sstr;
912 sstr << "Coefficient " << poly_coeff;
913 for (iq = 0; iq < jacobian_quantities.nelem() && !found; iq++) {
914 if (jacobian_quantities[iq] == Jacobian::Sensor::Polyfit &&
915 jacobian_quantities[iq].Subtag() == sstr.str()) {
916 found = true;
917 break;
918 }
919 }
920 if (!found) {
921 throw runtime_error(
922 "There is no Polyfit jacobian defined, in general "
923 "or for the selected polynomial coefficient.\n");
924 }
925
926 // Size and check of sensor_response
927 //
928 const Index nf = sensor_response_f_grid.nelem();
929 const Index npol = sensor_response_pol_grid.nelem();
930 const Index nlos = sensor_response_dlos_grid.nrows();
931
932 // Make a vector with values to distribute over *jacobian*
933 //
934 Vector w;
935 //
936 polynomial_basis_func(w, sensor_response_f_grid, poly_coeff);
937
938 // Fill J
939 //
940 ArrayOfArrayOfIndex jacobian_indices;
941 {
942 bool any_affine;
943 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
944 }
945 //
946 ArrayOfVector jg = jacobian_quantities[iq].Grids();
947 const Index n1 = jg[1].nelem();
948 const Index n2 = jg[2].nelem();
949 const Index n3 = jg[3].nelem();
950 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
951 const Index row4 = rowind.offset;
952 Index col4 = jacobian_indices[iq][0];
953
954 if (n3 > 1) {
955 col4 += mblock_index * n2 * n1;
956 }
957
958 for (Index l = 0; l < nlos; l++) {
959 const Index row3 = row4 + l * nf * npol;
960 const Index col3 = col4 + l * n1;
961
962 for (Index f = 0; f < nf; f++) {
963 const Index row2 = row3 + f * npol;
964
965 for (Index p = 0; p < npol; p++) {
966 Index col1 = col3;
967 if (n1 > 1) {
968 col1 += p;
969 }
970
971 jacobian(row2 + p, col1) = w[f];
972 }
973 }
974 }
975}
976
977//----------------------------------------------------------------------------
978// Scattering species:
979//----------------------------------------------------------------------------
980
981/* Workspace method: Doxygen documentation will be auto-generated */
984 Agenda& jacobian_agenda,
985 const Index& atmosphere_dim,
986 const Vector& p_grid,
987 const Vector& lat_grid,
988 const Vector& lon_grid,
989 const Vector& rq_p_grid,
990 const Vector& rq_lat_grid,
991 const Vector& rq_lon_grid,
992 const String& species,
993 const String& quantity,
994 const Verbosity& verbosity) {
997
998 // Check that this species+quantity combination is not already included in
999 // the jacobian.
1000 for (Index it = 0; it < jq.nelem(); it++) {
1001 if (jq[it] == Jacobian::Special::ScatteringString && jq[it].Subtag() == species &&
1002 jq[it].SubSubtag() == quantity) {
1003 ostringstream os;
1004 os << "The combintaion of\n scattering species: " << species
1005 << "\n retrieval quantity: " << quantity
1006 << "\nis already included in *jacobian_quantities*.";
1007 throw runtime_error(os.str());
1008 }
1009 }
1010
1011 // Check retrieval grids, here we just check the length of the grids
1012 // vs. the atmosphere dimension
1013 ArrayOfVector grids(atmosphere_dim);
1014 {
1015 ostringstream os;
1016 if (!check_retrieval_grids(grids,
1017 os,
1018 p_grid,
1019 lat_grid,
1020 lon_grid,
1021 rq_p_grid,
1022 rq_lat_grid,
1023 rq_lon_grid,
1024 "retrieval pressure grid",
1025 "retrieval latitude grid",
1026 "retrievallongitude_grid",
1027 atmosphere_dim))
1028 throw runtime_error(os.str());
1029 }
1030
1031 // Create the new retrieval quantity
1033 rq.Target() = Jacobian::Target(Jacobian::Special::ScatteringString, species);
1034 rq.Subtag(species);
1035 rq.SubSubtag(quantity);
1036 rq.Grids(grids);
1037
1038 // Add it to the *jacobian_quantities*
1039 jq.push_back(rq);
1040
1041 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1042}
1043
1044//----------------------------------------------------------------------------
1045// Sinusoidal baseline fits:
1046//----------------------------------------------------------------------------
1047
1048/* Workspace method: Doxygen documentation will be auto-generated */
1051 Agenda& jacobian_agenda,
1052 const ArrayOfIndex& sensor_response_pol_grid,
1053 const Matrix& sensor_response_dlos_grid,
1054 const Matrix& sensor_pos,
1055 const Vector& period_lengths,
1056 const Index& no_pol_variation,
1057 const Index& no_los_variation,
1058 const Index& no_mblock_variation,
1059 const Verbosity&) {
1060 const Index np = period_lengths.nelem();
1061
1062 // Check that poly_order is >= 0
1063 if (np == 0) throw runtime_error("No sinusoidal periods has benn given.");
1064
1065 // Check that polyfit is not already included in the jacobian.
1066 for (Index it = 0; it < jq.nelem(); it++) {
1067 if (jq[it] == Jacobian::Sensor::Sinefit) {
1068 ostringstream os;
1069 os << "Polynomial baseline fit is already included in\n"
1070 << "*jacobian_quantities*.";
1071 throw runtime_error(os.str());
1072 }
1073 }
1074
1075 // "Grids"
1076 //
1077 // Grid dimensions correspond here to
1078 // 1: polynomial order
1079 // 2: polarisation
1080 // 3: viewing direction
1081 // 4: measurement block
1082 //
1083 ArrayOfVector grids(4);
1084 //
1085 if (no_pol_variation)
1086 grids[1] = Vector(1, 1);
1087 else
1088 grids[1] = uniform_grid(0, sensor_response_pol_grid.nelem(), 1);
1089 if (no_los_variation)
1090 grids[2] = Vector(1, 1);
1091 else
1092 grids[2] = uniform_grid(0, sensor_response_dlos_grid.nrows(), 1);
1093 if (no_mblock_variation)
1094 grids[3] = Vector(1, 1);
1095 else
1096 grids[3] = uniform_grid(0, sensor_pos.nrows(), 1);
1097
1098 // Create the new retrieval quantity
1100 rq.Target() = Jacobian::Target(Jacobian::Sensor::Sinefit);
1101 rq.Mode("");
1102 rq.Target().perturbation = 0;
1103
1104 // Each sinefit coeff. pair is treated as a retrieval quantity
1105 //
1106 for (Index i = 0; i < np; i++) {
1107 ostringstream sstr;
1108 sstr << "Period " << i;
1109 rq.Subtag(sstr.str());
1110
1111 // "Grid" has length 2, set to period length
1112 grids[0] = Vector(2, period_lengths[i]);
1113 rq.Grids(grids);
1114
1115 // Add it to the *jacobian_quantities*
1116 jq.push_back(rq);
1117
1118 // Add pointing method to the jacobian agenda
1119 jacobian_agenda.append("jacobianCalcSinefit", i);
1120 }
1121}
1122
1123/* Workspace method: Doxygen documentation will be auto-generated */
1124void jacobianCalcSinefit(Matrix& jacobian,
1125 const Index& mblock_index,
1126 const Vector& iyb _U_,
1127 const Vector& yb _U_,
1128 const Sparse& sensor_response,
1129 const ArrayOfIndex& sensor_response_pol_grid,
1130 const Vector& sensor_response_f_grid,
1131 const Matrix& sensor_response_dlos_grid,
1132 const ArrayOfRetrievalQuantity& jacobian_quantities,
1133 const Index& period_index,
1134 const Verbosity&) {
1135 // Find the retrieval quantity related to this method
1137 ArrayOfIndex ji;
1138 bool found = false;
1139 Index iq;
1140 ostringstream sstr;
1141 sstr << "Period " << period_index;
1142 for (iq = 0; iq < jacobian_quantities.nelem() && !found; iq++) {
1143 if (jacobian_quantities[iq] == Jacobian::Sensor::Sinefit &&
1144 jacobian_quantities[iq].Subtag() == sstr.str()) {
1145 found = true;
1146 break;
1147 }
1148 }
1149 if (!found) {
1150 throw runtime_error(
1151 "There is no Sinefit jacobian defined, in general "
1152 "or for the selected period length.\n");
1153 }
1154
1155 // Size and check of sensor_response
1156 //
1157 const Index nf = sensor_response_f_grid.nelem();
1158 const Index npol = sensor_response_pol_grid.nelem();
1159 const Index nlos = sensor_response_dlos_grid.nrows();
1160
1161 // Make vectors with values to distribute over *jacobian*
1162 //
1163 // (period length stored in grid 0)
1164 ArrayOfVector jg = jacobian_quantities[iq].Grids();
1165 //
1166 Vector s(nf), c(nf);
1167 //
1168 for (Index f = 0; f < nf; f++) {
1169 Numeric a = (sensor_response_f_grid[f] - sensor_response_f_grid[0]) *
1170 Constant::two_pi / jg[0][0];
1171 s[f] = sin(a);
1172 c[f] = cos(a);
1173 }
1174
1175 // Fill J
1176 //
1177 ArrayOfArrayOfIndex jacobian_indices;
1178 {
1179 bool any_affine;
1180 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
1181 }
1182 //
1183 const Index n1 = jg[1].nelem();
1184 const Index n2 = jg[2].nelem();
1185 const Index n3 = jg[3].nelem();
1186 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
1187 const Index row4 = rowind.offset;
1188 Index col4 = jacobian_indices[iq][0];
1189
1190 if (n3 > 1) {
1191 col4 += mblock_index * n2 * n1 * 2;
1192 }
1193
1194 for (Index l = 0; l < nlos; l++) {
1195 const Index row3 = row4 + l * nf * npol;
1196 const Index col3 = col4 + l * n1 * 2;
1197
1198 for (Index f = 0; f < nf; f++) {
1199 const Index row2 = row3 + f * npol;
1200
1201 for (Index p = 0; p < npol; p++) {
1202 Index col1 = col3;
1203 if (n1 > 1) {
1204 col1 += p * 2;
1205 }
1206
1207 jacobian(row2 + p, col1) = s[f];
1208 jacobian(row2 + p, col1 + 1) = c[f];
1209 }
1210 }
1211 }
1212}
1213
1214//----------------------------------------------------------------------------
1215// Surface quantities
1216//----------------------------------------------------------------------------
1217
1218/* Workspace method: Doxygen documentation will be auto-generated */
1221 Agenda& jacobian_agenda,
1222 const Index& atmosphere_dim,
1223 const Vector& lat_grid,
1224 const Vector& lon_grid,
1225 const Vector& rq_lat_grid,
1226 const Vector& rq_lon_grid,
1227 const String& quantity,
1228 const Verbosity& verbosity) {
1231
1232 // Check that this species is not already included in the jacobian.
1233 for (Index it = 0; it < jq.nelem(); it++) {
1234 if (jq[it] == Jacobian::Special::SurfaceString && jq[it].Subtag() == quantity) {
1235 ostringstream os;
1236 os << quantity << " is already included as a surface variable "
1237 << "in *jacobian_quantities*.";
1238 throw runtime_error(os.str());
1239 }
1240 }
1241
1242 // Check retrieval grids, here we just check the length of the grids
1243 // vs. the atmosphere dimension
1244 ArrayOfVector grids(max(atmosphere_dim - 1, Index(1)));
1245 {
1246 ostringstream os;
1247 if (!check_retrieval_grids(grids,
1248 os,
1249 lat_grid,
1250 lon_grid,
1251 rq_lat_grid,
1252 rq_lon_grid,
1253 "retrieval latitude grid",
1254 "retrievallongitude_grid",
1255 atmosphere_dim))
1256 throw runtime_error(os.str());
1257 }
1258
1259 // Create the new retrieval quantity
1261 rq.Target() = Jacobian::Target(Jacobian::Special::SurfaceString, quantity);
1262 rq.Subtag(quantity);
1263 rq.Grids(grids);
1264
1265 // Add it to the *jacobian_quantities*
1266 jq.push_back(rq);
1267
1268 // Add dummy
1269 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1270}
1271
1272//----------------------------------------------------------------------------
1273// Temperatures (atmospheric):
1274//----------------------------------------------------------------------------
1275
1276/* Workspace method: Doxygen documentation will be auto-generated */
1279 Agenda& jacobian_agenda,
1280 const Index& atmosphere_dim,
1281 const Vector& p_grid,
1282 const Vector& lat_grid,
1283 const Vector& lon_grid,
1284 const Vector& rq_p_grid,
1285 const Vector& rq_lat_grid,
1286 const Vector& rq_lon_grid,
1287 const String& hse,
1288 const Verbosity& verbosity) {
1290
1291 // Check that temperature is not already included in the jacobian.
1292 // We only check the main tag.
1293 for (Index it = 0; it < jq.nelem(); it++) {
1294 if (jq[it] == Jacobian::Atm::Temperature) {
1295 ostringstream os;
1296 os << "Temperature is already included in *jacobian_quantities*.";
1297 throw runtime_error(os.str());
1298 }
1299 }
1300
1301 // Check retrieval grids, here we just check the length of the grids
1302 // vs. the atmosphere dimension
1303 ArrayOfVector grids(atmosphere_dim);
1304 {
1305 ostringstream os;
1306 if (!check_retrieval_grids(grids,
1307 os,
1308 p_grid,
1309 lat_grid,
1310 lon_grid,
1311 rq_p_grid,
1312 rq_lat_grid,
1313 rq_lon_grid,
1314 "retrieval pressure grid",
1315 "retrieval latitude grid",
1316 "retrievallongitude_grid",
1317 atmosphere_dim))
1318 throw runtime_error(os.str());
1319 }
1320
1321 // Set subtag
1322 String subtag;
1323 if (hse == "on") {
1324 subtag = "HSE on";
1325 } else if (hse == "off") {
1326 subtag = "HSE off";
1327 } else {
1328 ostringstream os;
1329 os << "The keyword for hydrostatic equilibrium can only be set to\n"
1330 << "\"on\" or \"off\"\n";
1331 throw runtime_error(os.str());
1332 }
1333
1334 // Create the new retrieval quantity
1336 rq.Subtag(subtag);
1337 rq.Grids(grids);
1338 rq.Target(Jacobian::Target(Jacobian::Atm::Temperature));
1339 rq.Target().perturbation = 0.1;
1340
1341 // Add it to the *jacobian_quantities*
1342 jq.push_back(rq);
1343
1344 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1345}
1346
1347//----------------------------------------------------------------------------
1348// Winds:
1349//----------------------------------------------------------------------------
1350
1351/* Workspace method: Doxygen documentation will be auto-generated */
1354 Agenda& jacobian_agenda,
1355 const Index& atmosphere_dim,
1356 const Vector& p_grid,
1357 const Vector& lat_grid,
1358 const Vector& lon_grid,
1359 const Vector& rq_p_grid,
1360 const Vector& rq_lat_grid,
1361 const Vector& rq_lon_grid,
1362 const String& component,
1363 const Numeric& dfrequency,
1364 const Verbosity& verbosity) {
1367
1368 // Create the new retrieval quantity
1369 const auto opt = Options::toWindMagJacobianOrThrow(component);
1371 switch (opt) {
1372 case Options::WindMagJacobian::u:
1373 rq.Target(Jacobian::Target(Jacobian::Atm::WindU));
1374 break;
1375 case Options::WindMagJacobian::v:
1376 rq.Target(Jacobian::Target(Jacobian::Atm::WindV));
1377 break;
1378 case Options::WindMagJacobian::w:
1379 rq.Target(Jacobian::Target(Jacobian::Atm::WindW));
1380 break;
1381 case Options::WindMagJacobian::strength:
1382 rq.Target(Jacobian::Target(Jacobian::Atm::WindMagnitude));
1383 break;
1384 case Options::WindMagJacobian::FINAL:
1385 ARTS_ASSERT(false, "This error should be caught earlier")
1386 }
1387
1388 // Check that this species is not already included in the jacobian.
1389 for (Index it = 0; it < jq.nelem(); it++) {
1390 ARTS_USER_ERROR_IF (jq[it].Target().sameTargetType(rq.Target()),
1391 "The wind component:\n", component, "\n"
1392 "is already included in *jacobian_quantities*.")
1393 }
1394
1395 // Check retrieval grids, here we just check the length of the grids
1396 // vs. the atmosphere dimension
1397 ArrayOfVector grids(atmosphere_dim);
1398 {
1399 ostringstream os;
1400 if (!check_retrieval_grids(grids,
1401 os,
1402 p_grid,
1403 lat_grid,
1404 lon_grid,
1405 rq_p_grid,
1406 rq_lat_grid,
1407 rq_lon_grid,
1408 "retrieval pressure grid",
1409 "retrieval latitude grid",
1410 "retrievallongitude_grid",
1411 atmosphere_dim))
1412 throw runtime_error(os.str());
1413 }
1414
1415 rq.Subtag(component); // nb. This should be possible to remove...
1416 rq.Grids(grids);
1417 rq.Target().perturbation = dfrequency;
1418
1419 // Add it to the *jacobian_quantities*
1420 jq.push_back(rq);
1421
1422 out3 << " Calculations done by propagation matrix expression.\n";
1423 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1424}
1425
1426//----------------------------------------------------------------------------
1427// Magnetic field:
1428//----------------------------------------------------------------------------
1429
1430/* Workspace method: Doxygen documentation will be auto-generated */
1433 Agenda& jacobian_agenda,
1434 const Index& atmosphere_dim,
1435 const Vector& p_grid,
1436 const Vector& lat_grid,
1437 const Vector& lon_grid,
1438 const Vector& rq_p_grid,
1439 const Vector& rq_lat_grid,
1440 const Vector& rq_lon_grid,
1441 const String& component,
1442 const Numeric& dB,
1443 const Verbosity& verbosity) {
1446
1447 // Create the new retrieval quantity
1448 const auto opt = Options::toWindMagJacobianOrThrow(component);
1450 switch (opt) {
1451 case Options::WindMagJacobian::u:
1452 rq.Target(Jacobian::Target(Jacobian::Atm::MagneticU));
1453 break;
1454 case Options::WindMagJacobian::v:
1455 rq.Target(Jacobian::Target(Jacobian::Atm::MagneticV));
1456 break;
1457 case Options::WindMagJacobian::w:
1458 rq.Target(Jacobian::Target(Jacobian::Atm::MagneticW));
1459 break;
1460 case Options::WindMagJacobian::strength:
1461 rq.Target(Jacobian::Target(Jacobian::Atm::MagneticMagnitude));
1462 break;
1463 case Options::WindMagJacobian::FINAL:
1464 ARTS_ASSERT(false, "This error should be caught earlier")
1465 }
1466
1467 // Check that this species is not already included in the jacobian.
1468 for (Index it = 0; it < jq.nelem(); it++) {
1469 ARTS_USER_ERROR_IF (jq[it].Target().sameTargetType(rq.Target()),
1470 "The magnetic component:\n", component, "\n"
1471 "is already included in *jacobian_quantities*.")
1472 }
1473
1474 // Check retrieval grids, here we just check the length of the grids
1475 // vs. the atmosphere dimension
1476 ArrayOfVector grids(atmosphere_dim);
1477 {
1478 ostringstream os;
1479 if (!check_retrieval_grids(grids,
1480 os,
1481 p_grid,
1482 lat_grid,
1483 lon_grid,
1484 rq_p_grid,
1485 rq_lat_grid,
1486 rq_lon_grid,
1487 "retrieval pressure grid",
1488 "retrieval latitude grid",
1489 "retrievallongitude_grid",
1490 atmosphere_dim))
1491 throw runtime_error(os.str());
1492 }
1493
1494 rq.Grids(grids);
1495 rq.Target().perturbation = dB;
1496
1497 // Add it to the *jacobian_quantities*
1498 jq.push_back(rq);
1499
1500 // Add gas species method to the jacobian agenda
1501 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1502}
1503
1504//----------------------------------------------------------------------------
1505// Catalog parameters:
1506//----------------------------------------------------------------------------
1507
1508/* Workspace method: Doxygen documentation will be auto-generated */
1511 Agenda& jacobian_agenda,
1512 const QuantumIdentifier& line_identity,
1513 const String& species,
1514 const String& variable,
1515 const String& coefficient,
1516 const Verbosity& verbosity) {
1518
1519 const auto jpt = select_derivativeLineShape(variable, coefficient);
1520
1521 out3 << "Attempting to create RT tag for " << line_identity << " " << variable
1522 << " " << coefficient << " for ";
1523
1524 // Create the quantity
1526 rq.Mode(species);
1527 rq.Grids(ArrayOfVector(0, Vector(0)));
1528
1529 // Map the species
1530 if (species == LineShape::self_broadening) {
1531 rq.Target(Jacobian::Target(jpt, line_identity, line_identity.Species()));
1532 } else if (species == LineShape::bath_broadening) {
1533 rq.Target(Jacobian::Target(jpt, line_identity, Species::Species::Bath));
1534 } else {
1535 rq.Target(Jacobian::Target(jpt, line_identity, SpeciesTag(species).Spec()));
1536 }
1537 out3 << species << ' ' << rq.Target().species_id << "\n";
1538
1539 // Test this is not a copy
1540 for (auto& q : jq)
1541 if (q.HasSameInternalsAs(rq))
1542 throw std::runtime_error("Error with copies of the quantities");
1543
1544 // Append and do housekeeping
1545 jq.push_back(rq);
1546 out3 << "Creation was successful!\n";
1547 jacobian_agenda.append("jacobianCalcDoNothing",
1548 TokVal()); // old code activation
1549}
1550
1551/* Workspace method: Doxygen documentation will be auto-generated */
1553 Workspace& ws,
1555 Agenda& jacobian_agenda,
1556 const ArrayOfQuantumIdentifier& line_identities,
1557 const ArrayOfString& species,
1558 const ArrayOfString& variables,
1559 const ArrayOfString& coefficients,
1560 const Verbosity& verbosity) {
1561 if (not(line_identities.nelem() or species.nelem() or variables.nelem() or
1562 coefficients.nelem()))
1563 throw std::runtime_error("Must have at least 1-long lists for all GINs");
1564
1565 ArrayOfString vars;
1566 if (variables[0] == "ALL")
1567 vars = ArrayOfString(LineShape::enumstrs::VariableNames);
1568 else
1569 vars = variables;
1570
1571 ArrayOfString coeffs;
1572 if (coefficients[0] == "ALL")
1573 coeffs = ArrayOfString(Options::enumstrs::LineShapeCoeffNames);
1574 else
1575 coeffs = coefficients;
1576
1577 for (auto& l : line_identities)
1578 for (auto& s : species)
1579 for (auto& v : vars)
1580 for (auto& c : coeffs)
1582 ws, jq, jacobian_agenda, l, s, v, c, verbosity);
1583}
1584
1585/* Workspace method: Doxygen documentation will be auto-generated */
1588 Agenda& jacobian_agenda,
1589 const QuantumIdentifier& catalog_identity,
1590 const String& catalog_parameter,
1591 const Verbosity& verbosity) {
1593
1594 // Create the new retrieval quantity
1595 const auto opt = Options::toBasicCatParamJacobianOrThrow(catalog_parameter);
1597 switch (opt) {
1598 case Options::BasicCatParamJacobian::LineCenter:
1599 rq.Target(Jacobian::Target(Jacobian::Line::Center, catalog_identity, catalog_identity.Species()));
1600 break;
1601 case Options::BasicCatParamJacobian::LineStrength:
1602 rq.Target(Jacobian::Target(Jacobian::Line::Strength, catalog_identity, catalog_identity.Species()));
1603 break;
1604 case Options::BasicCatParamJacobian::FINAL:
1605 ARTS_ASSERT(false, "This error should be caught earlier")
1606 }
1607
1608 // Check that this is not already included in the jacobian.
1609 for (Index it = 0; it < jq.nelem(); it++) {
1610 ARTS_USER_ERROR_IF (rq.Target().sameTargetType(jq[it].Target()),
1611 "The catalog identifier:\n",
1612 catalog_identity, " for ID: ", catalog_identity, "\n"
1613 "is already included in jacobian_quantities*.")
1614 }
1615
1616 rq.Grids(ArrayOfVector(0, Vector(0)));
1617
1618 // Add it to the *jacobian_quantities*
1619 jq.push_back(rq);
1620
1621 out3 << " Calculations done by propagation matrix expressions.\n";
1622
1623 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1624}
1625
1626/* Workspace method: Doxygen documentation will be auto-generated */
1628 Workspace& ws,
1630 Agenda& jacobian_agenda,
1631 const ArrayOfQuantumIdentifier& catalog_identities,
1632 const ArrayOfString& catalog_parameters,
1633 const Verbosity& verbosity) {
1635
1636 out2 << " Adding " << catalog_identities.nelem() * catalog_parameters.nelem()
1637 << " expression(s) to the Jacobian calculations.\n";
1638
1639 for (auto& qi : catalog_identities)
1640 for (auto& param : catalog_parameters)
1642 ws, jq, jacobian_agenda, qi, param, verbosity);
1643}
1644
1645//----------------------------------------------------------------------------
1646// NLTE temperature:
1647//----------------------------------------------------------------------------
1648
1649/* Workspace method: Doxygen documentation will be auto-generated */
1652 Agenda& jacobian_agenda,
1653 const Index& atmosphere_dim,
1654 const Vector& p_grid,
1655 const Vector& lat_grid,
1656 const Vector& lon_grid,
1657 const Vector& rq_p_grid,
1658 const Vector& rq_lat_grid,
1659 const Vector& rq_lon_grid,
1660 const QuantumIdentifier& energy_level_identity,
1661 const Numeric& dx,
1662 const Verbosity& verbosity) {
1664
1665 // Check that this species is not already included in the jacobian.
1666 for (Index it = 0; it < jq.nelem(); it++) {
1667 if (jq[it] == Jacobian::Line::NLTE and
1668 jq[it].QuantumIdentity() == energy_level_identity) {
1669 ostringstream os;
1670 os << "The NLTE identifier:\n"
1671 << energy_level_identity << "\nis already included in "
1672 << "*jacobian_quantities*.";
1673 throw std::runtime_error(os.str());
1674 }
1675 }
1676
1677 // Check retrieval grids, here we just check the length of the grids
1678 // vs. the atmosphere dimension
1679 ArrayOfVector grids(atmosphere_dim);
1680 {
1681 ostringstream os;
1682 if (not check_retrieval_grids(grids,
1683 os,
1684 p_grid,
1685 lat_grid,
1686 lon_grid,
1687 rq_p_grid,
1688 rq_lat_grid,
1689 rq_lon_grid,
1690 "retrieval pressure grid",
1691 "retrieval latitude grid",
1692 "retrievallongitude_grid",
1693 atmosphere_dim))
1694 throw runtime_error(os.str());
1695 }
1696
1697 // Create the new retrieval quantity
1699 rq.Target(Jacobian::Target(Jacobian::Line::NLTE, energy_level_identity, energy_level_identity.Species()));
1700 rq.Target().perturbation = dx;
1701 rq.Grids(grids);
1702
1703 // Add it to the *jacobian_quantities*
1704 jq.push_back(rq);
1705
1706 out3 << " Calculations done by propagation matrix expressions.\n";
1707
1708 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1709}
1710
1713 Agenda& jacobian_agenda,
1714 const Index& atmosphere_dim,
1715 const Vector& p_grid,
1716 const Vector& lat_grid,
1717 const Vector& lon_grid,
1718 const Vector& rq_p_grid,
1719 const Vector& rq_lat_grid,
1720 const Vector& rq_lon_grid,
1721 const ArrayOfQuantumIdentifier& energy_level_identities,
1722 const Numeric& dx,
1723 const Verbosity& verbosity) {
1724 for (const auto& qi : energy_level_identities)
1725 jacobianAddNLTE(ws,
1726 jq,
1727 jacobian_agenda,
1728 atmosphere_dim,
1729 p_grid,
1730 lat_grid,
1731 lon_grid,
1732 rq_p_grid,
1733 rq_lat_grid,
1734 rq_lon_grid,
1735 qi,
1736 dx,
1737 verbosity);
1738}
1739
1740/* Workspace method: Doxygen documentation will be auto-generated */
1743 Agenda& jacobian_agenda,
1744 const Index& atmosphere_dim,
1745 const Vector& p_grid,
1746 const Vector& lat_grid,
1747 const Vector& lon_grid,
1748 const Vector& rq_p_grid,
1749 const Vector& rq_lat_grid,
1750 const Vector& rq_lon_grid,
1751 const String& species,
1752 const Verbosity& verbosity) {
1755
1756 // Check retrieval grids, here we just check the length of the grids
1757 // vs. the atmosphere dimension
1758 ArrayOfVector grids(atmosphere_dim);
1759 {
1760 ostringstream os;
1761 if (!check_retrieval_grids(grids,
1762 os,
1763 p_grid,
1764 lat_grid,
1765 lon_grid,
1766 rq_p_grid,
1767 rq_lat_grid,
1768 rq_lon_grid,
1769 "retrieval pressure grid",
1770 "retrieval latitude grid",
1771 "retrievallongitude_grid",
1772 atmosphere_dim))
1773 throw runtime_error(os.str());
1774 }
1775
1776 // Create the new retrieval quantity
1778 rq.Grids(grids);
1779
1780 // Make sure modes are valid and complain if they are repeated
1781 if (species == "electrons") {
1782 for (Index it = 0; it < jq.nelem(); it++) {
1783 if (jq[it] == Jacobian::Atm::Electrons) {
1784 ostringstream os;
1785 os << "Electrons are already included in *jacobian_quantities*.";
1786 throw std::runtime_error(os.str());
1787 }
1788 }
1789 rq.Target(Jacobian::Target(Jacobian::Atm::Electrons));
1790
1791 } else if (species == "particulates") {
1792 for (Index it = 0; it < jq.nelem(); it++) {
1793 if (jq[it] == Jacobian::Atm::Particulates) {
1794 ostringstream os;
1795 os << "Particulates are already included in *jacobian_quantities*.";
1796 throw std::runtime_error(os.str());
1797 }
1798 }
1799 rq.Target(Jacobian::Target(Jacobian::Atm::Particulates));
1800
1801 } else {
1802 ostringstream os;
1803 os << "Unknown special species jacobian: \"" << species
1804 << "\"\nPlease see *jacobianAddSpecialSpecies* for viable options.";
1805 throw std::runtime_error(os.str());
1806 }
1807
1808 // Add it to the *jacobian_quantities*
1809 jq.push_back(rq);
1810
1811 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1812}
1813
1814//----------------------------------------------------------------------------
1815// Adjustments and transformations
1816//----------------------------------------------------------------------------
1817
1818/* Workspace method: Doxygen documentation will be auto-generated */
1820 Matrix& jacobian,
1821 const ArrayOfRetrievalQuantity& jacobian_quantities,
1822 const Vector& x,
1823 const Verbosity&) {
1824 // For flexibility inside inversion_iteration_agenda, we should accept empty
1825 // Jacobian
1826 if (jacobian.empty()) {
1827 return;
1828 }
1829
1830 // Adjustments
1831 //
1832 // Unfortunately the adjustment requires both range indices and the
1833 // untransformed x, which makes things a bit messy
1834 bool vars_init = false;
1836 Vector x0;
1837 //
1838 for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
1839 if (jacobian_quantities[q].Target().isSpeciesVMR() &&
1840 jacobian_quantities[q].Mode() == "rel") {
1841 if (!vars_init) {
1842 bool any_affine;
1843 jac_ranges_indices(jis0, any_affine, jacobian_quantities, true);
1844 x0 = x;
1845 transform_x_back(x0, jacobian_quantities);
1846 vars_init = true;
1847 }
1848 for (Index i = jis0[q][0]; i <= jis0[q][1]; i++) {
1849 if (x[i] != 1) {
1850 jacobian(joker, i) /= x[i];
1851 }
1852 }
1853 }
1854 }
1855
1856 // Transformations
1857 transform_jacobian(jacobian, x, jacobian_quantities);
1858}
1859
1860/* Workspace method: Doxygen documentation will be auto-generated */
1862 const Matrix& transformation_matrix,
1863 const Vector& offset_vector,
1864 const Verbosity& /*v*/
1865) {
1866 if (jqs.empty()) {
1867 runtime_error(
1868 "Jacobian quantities is empty, so there is nothing to add the "
1869 "transformation to.");
1870 }
1871
1872 Index nelem = jqs.back().Grids().nelem();
1873
1874 if (!(nelem == transformation_matrix.nrows())) {
1875 runtime_error(
1876 "Dimension of transformation matrix incompatible with retrieval grids.");
1877 }
1878 if (!(nelem == offset_vector.nelem())) {
1879 runtime_error(
1880 "Dimension of offset vector incompatible with retrieval grids.");
1881 }
1882
1883 jqs.back().SetTransformationMatrix(Matrix{transpose(transformation_matrix)});
1884 jqs.back().SetOffsetVector(offset_vector);
1885}
1886
1887/* Workspace method: Doxygen documentation will be auto-generated */
1889 const String& transformation_func,
1890 const Numeric& z_min,
1891 const Numeric& z_max,
1892 const Verbosity& /*v*/
1893) {
1894 if (jqs.empty())
1895 throw runtime_error(
1896 "Jacobian quantities is empty, so there is nothing to add the "
1897 "transformation to.");
1898
1899 if (transformation_func == "none") {
1900 jqs.back().SetTransformationFunc("");
1901 return;
1902 }
1903
1904 Vector pars;
1905
1906 if (transformation_func == "atanh") {
1907 if (z_max <= z_min)
1908 throw runtime_error(
1909 "For option atanh, the GIN *z_max* must be set and be > z_min.");
1910 pars.resize(2);
1911 pars[0] = z_min;
1912 pars[1] = z_max;
1913 } else if (transformation_func == "log" || transformation_func == "log10") {
1914 pars.resize(1);
1915 pars[0] = z_min;
1916 } else {
1917 ostringstream os;
1918 os << "Valid options for *transformation_func* are:\n"
1919 << "\"none\", \"log\", \"log10\" and \"atanh\"\n"
1920 << "But found: \"" << transformation_func << "\"";
1921 throw runtime_error(os.str());
1922 }
1923
1924 jqs.back().SetTransformationFunc(transformation_func);
1925 jqs.back().SetTFuncParameters(pars);
1926}
1927
1928//----------------------------------------------------------------------------
1929// Methods for doing perturbations
1930//----------------------------------------------------------------------------
1931
1932/* Workspace method: Doxygen documentation will be auto-generated */
1933void AtmFieldPerturb(Tensor3& perturbed_field,
1934 const Index& atmosphere_dim,
1935 const Vector& p_grid,
1936 const Vector& lat_grid,
1937 const Vector& lon_grid,
1938 const Tensor3& original_field,
1939 const Vector& p_ret_grid,
1940 const Vector& lat_ret_grid,
1941 const Vector& lon_ret_grid,
1942 const Index& pert_index,
1943 const Numeric& pert_size,
1944 const String& pert_mode,
1945 const Verbosity&) {
1946 // Input checks (more below)
1947 chk_atm_field("original_field",
1948 original_field,
1949 atmosphere_dim,
1950 p_grid,
1951 lat_grid,
1952 lon_grid,
1953 false );
1954
1955 // Pack retrieval grids into an ArrayOfVector
1956 ArrayOfVector ret_grids(atmosphere_dim);
1957 ret_grids[0] = p_ret_grid;
1958 if (atmosphere_dim>1){
1959 ret_grids[1] = lat_ret_grid;
1960 if (atmosphere_dim>2){
1961 ret_grids[2] = lon_ret_grid;
1962 }
1963 }
1964
1965 // Find mapping from retrieval grids to atmospheric grids
1966 ArrayOfGridPos gp_p, gp_lat, gp_lon;
1967 Index n_p, n_lat, n_lon;
1969 gp_lat,
1970 gp_lon,
1971 n_p,
1972 n_lat,
1973 n_lon,
1974 ret_grids,
1975 atmosphere_dim,
1976 p_grid,
1977 lat_grid,
1978 lon_grid);
1979
1980 // Now we can chec *pert_index*
1981 if (pert_index<0){
1982 throw runtime_error("Bad *pert_index*. It is negative.");
1983 }
1984 const Index n_tot = n_p * n_lat * n_lon;
1985 if (pert_index >= n_tot){
1986 throw runtime_error("Bad *pert_index*. It is too high with respect "
1987 "to length of retrieval grids.");
1988 }
1989
1990 // Create x-vector that matches perturbation
1991 Vector x(n_tot);
1992 if (pert_mode == "absolute" ){
1993 x = 0;
1994 x[pert_index] = pert_size;
1995 }
1996 else if (pert_mode == "relative" ){
1997 x = 1;
1998 x[pert_index] += pert_size;
1999 }
2000 else{
2001 throw runtime_error("Bad *pert_mode*. Allowed choices are: "
2002 """absolute"" and ""relative"".");
2003 }
2004
2005 // Map x to a perturbation defined at atmospheric grids
2006 Tensor3 x3d(n_p, n_lat, n_lon), pert(n_p, n_lat, n_lon);
2007 reshape(x3d, x);
2008 regrid_atmfield_by_gp_oem(pert, atmosphere_dim, x3d, gp_p, gp_lat, gp_lon);
2009
2010 // Init perturbed_field, if not equal to original_field
2011 if (&perturbed_field != &original_field) {
2012 perturbed_field = original_field;
2013 }
2014
2015 // Apply perturbation
2016 if (pert_mode == "absolute" ){
2017 perturbed_field += pert;
2018 }
2019 else{
2020 perturbed_field *= pert;
2021 }
2022}
2023
2024/* Workspace method: Doxygen documentation will be auto-generated */
2025void AtmFieldPerturbAtmGrids(Tensor3& perturbed_field,
2026 const Index& atmosphere_dim,
2027 const Vector& p_grid,
2028 const Vector& lat_grid,
2029 const Vector& lon_grid,
2030 const Tensor3& original_field,
2031 const Index& pert_index,
2032 const Numeric& pert_size,
2033 const String& pert_mode,
2034 const Verbosity&) {
2035 // Some sizes
2036 const Index n_p = p_grid.nelem();
2037 const Index n_lat = atmosphere_dim<2 ? 1 : lat_grid.nelem();
2038 const Index n_lon = atmosphere_dim<3 ? 1 : lon_grid.nelem();
2039
2040 // Check input
2041 chk_atm_field("original_field",
2042 original_field,
2043 atmosphere_dim,
2044 p_grid,
2045 lat_grid,
2046 lon_grid,
2047 false );
2048 if (pert_index<0){
2049 throw runtime_error("Bad *pert_index*. It is negative.");
2050 }
2051 if (pert_index >= n_p * n_lat * n_lon){
2052 throw runtime_error("Bad *pert_index*. It is too high with respect "
2053 "to length of atmospheric grids.");
2054 }
2055
2056 // Determine indexes with respect to atmospheric grids
2057 Index tot_index = pert_index;
2058 const Index lon_index = atmosphere_dim<3 ? 0 : tot_index / (n_lat * n_p);
2059 tot_index -= lon_index * n_lat * n_p;
2060 const Index lat_index = atmosphere_dim<2 ? 0 : tot_index / n_p;
2061 tot_index -= lat_index * n_p;
2062 const Index p_index = tot_index;
2063
2064 // Init perturbed_field, if not equal to original_field
2065 if (&perturbed_field != &original_field) {
2066 perturbed_field = original_field;
2067 }
2068
2069 // Perturb
2070 if (pert_mode == "absolute" ){
2071 perturbed_field(p_index,
2072 atmosphere_dim>1 ? lat_index : 0,
2073 atmosphere_dim>2 ? lon_index : 0) += pert_size;
2074 }
2075 else if (pert_mode == "relative"){
2076 perturbed_field(p_index,
2077 atmosphere_dim>1 ? lat_index : 0,
2078 atmosphere_dim>2 ? lon_index : 0) *= 1 + pert_size;
2079 }
2080 else{
2081 throw runtime_error("Bad *pert_mode*. Allowed choices are: "
2082 """absolute"" and ""relative"".");
2083 }
2084}
2085
2086/* Workspace method: Doxygen documentation will be auto-generated */
2088 const Index& atmosphere_dim,
2089 const Vector& p_grid,
2090 const Vector& lat_grid,
2091 const Vector& lon_grid,
2092 const Verbosity&) {
2093 const Index n_p = p_grid.nelem();
2094 const Index n_lat = atmosphere_dim<2 ? 1 : lat_grid.nelem();
2095 const Index n_lon = atmosphere_dim<3 ? 1 : lon_grid.nelem();
2096
2097 n = n_p * n_lat * n_lon;
2098}
2099
2100/* Workspace method: Doxygen documentation will be auto-generated */
2101void jacobianFromTwoY(Matrix& jacobian,
2102 const Vector& y_pert,
2103 const Vector& y,
2104 const Numeric& pert_size,
2105 const Verbosity&) {
2106 const Index n = y.nelem();
2107 if( y_pert.nelem() != n ){
2108 throw runtime_error("Inconsistency in length of *y_pert* and *y*.");
2109 }
2110 jacobian = ExhaustiveConstMatrixView{y_pert};
2111 jacobian -= ExhaustiveConstMatrixView{y};
2112 jacobian /= pert_size;
2113}
2114
2115/* Workspace method: Doxygen documentation will be auto-generated */
2116void jacobianFromYbatch(Matrix& jacobian,
2117 const ArrayOfVector& ybatch,
2118 const Vector& y,
2119 const Numeric& pert_size,
2120 const Verbosity&) {
2121 const Index n = y.nelem();
2122 const Index l = ybatch.nelem();
2123 if (l>0){
2124 if( ybatch[0].nelem() != n )
2125 throw runtime_error("Inconsistency in length of y and ybatch[0].");
2126 }
2127 jacobian.resize(n,l);
2128 for (Index i=0; i<l; i++) {
2129 jacobian(joker,i) = ybatch[i];
2130 jacobian(joker,i) -= y;
2131 }
2132 jacobian /= pert_size;
2133}
2134
2135/* Workspace method: Doxygen documentation will be auto-generated */
2136void particle_bulkprop_fieldPerturb(Tensor4& particle_bulkprop_field,
2137 const Index& atmosphere_dim,
2138 const Vector& p_grid,
2139 const Vector& lat_grid,
2140 const Vector& lon_grid,
2141 const ArrayOfString& particle_bulkprop_names,
2142 const String& particle_type,
2143 const Vector& p_ret_grid,
2144 const Vector& lat_ret_grid,
2145 const Vector& lon_ret_grid,
2146 const Index& pert_index,
2147 const Numeric& pert_size,
2148 const String& pert_mode,
2149 const Verbosity& verbosity) {
2150 // Locate particle_type among particle_bulkprop_names
2151 Index iq = find_first(particle_bulkprop_names, particle_type);
2152 if (iq < 0) {
2153 ostringstream os;
2154 os << "Could not find " << particle_type << " in *particle_bulkprop_names*.\n";
2155 throw std::runtime_error(os.str());
2156 }
2157
2158 Tensor3 original_field, perturbed_field;
2159 original_field = particle_bulkprop_field(iq,joker,joker,joker);
2160 AtmFieldPerturb(perturbed_field,
2161 atmosphere_dim,
2162 p_grid,
2163 lat_grid,
2164 lon_grid,
2165 original_field,
2166 p_ret_grid,
2167 lat_ret_grid,
2168 lon_ret_grid,
2169 pert_index,
2170 pert_size,
2171 pert_mode,
2172 verbosity);
2173 particle_bulkprop_field(iq,joker,joker,joker) = perturbed_field;
2174}
2175
2176/* Workspace method: Doxygen documentation will be auto-generated */
2177void particle_bulkprop_fieldPerturbAtmGrids(Tensor4& particle_bulkprop_field,
2178 const Index& atmosphere_dim,
2179 const Vector& p_grid,
2180 const Vector& lat_grid,
2181 const Vector& lon_grid,
2182 const ArrayOfString& particle_bulkprop_names,
2183 const String& particle_type,
2184 const Index& pert_index,
2185 const Numeric& pert_size,
2186 const String& pert_mode,
2187 const Verbosity& verbosity) {
2188 // Locate particle_type among particle_bulkprop_names
2189 Index iq = find_first(particle_bulkprop_names, particle_type);
2190 if (iq < 0) {
2191 ostringstream os;
2192 os << "Could not find " << particle_type << " in *particle_bulkprop_names*.\n";
2193 throw std::runtime_error(os.str());
2194 }
2195
2196 Tensor3 original_field, perturbed_field;
2197 original_field = particle_bulkprop_field(iq,joker,joker,joker);
2198 AtmFieldPerturbAtmGrids(perturbed_field,
2199 atmosphere_dim,
2200 p_grid,
2201 lat_grid,
2202 lon_grid,
2203 original_field,
2204 pert_index,
2205 pert_size,
2206 pert_mode,
2207 verbosity);
2208 particle_bulkprop_field(iq,joker,joker,joker) = perturbed_field;
2209}
2210
2211/* Workspace method: Doxygen documentation will be auto-generated */
2212void vmr_fieldPerturb(Tensor4& vmr_field,
2213 const Index& atmosphere_dim,
2214 const Vector& p_grid,
2215 const Vector& lat_grid,
2216 const Vector& lon_grid,
2217 const ArrayOfArrayOfSpeciesTag& abs_species,
2218 const String& species,
2219 const Vector& p_ret_grid,
2220 const Vector& lat_ret_grid,
2221 const Vector& lon_ret_grid,
2222 const Index& pert_index,
2223 const Numeric& pert_size,
2224 const String& pert_mode,
2225 const Verbosity& verbosity) {
2226 // Locate vmr_species among abs_species
2227 Index iq = -1;
2228 for (Index i = 0; i < abs_species.nelem(); i++) {
2229 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
2230 iq = i;
2231 break;
2232 }
2233 }
2234 if (iq < 0) {
2235 ostringstream os;
2236 os << "Could not find " << species << " in *abs_species*.\n";
2237 throw std::runtime_error(os.str());
2238 }
2239
2240 Tensor3 original_field, perturbed_field;
2241 original_field = vmr_field(iq,joker,joker,joker);
2242 AtmFieldPerturb(perturbed_field,
2243 atmosphere_dim,
2244 p_grid,
2245 lat_grid,
2246 lon_grid,
2247 original_field,
2248 p_ret_grid,
2249 lat_ret_grid,
2250 lon_ret_grid,
2251 pert_index,
2252 pert_size,
2253 pert_mode,
2254 verbosity);
2255 vmr_field(iq,joker,joker,joker) = perturbed_field;
2256}
2257
2258/* Workspace method: Doxygen documentation will be auto-generated */
2259void vmr_fieldPerturbAtmGrids(Tensor4& vmr_field,
2260 const Index& atmosphere_dim,
2261 const Vector& p_grid,
2262 const Vector& lat_grid,
2263 const Vector& lon_grid,
2264 const ArrayOfArrayOfSpeciesTag& abs_species,
2265 const String& species,
2266 const Index& pert_index,
2267 const Numeric& pert_size,
2268 const String& pert_mode,
2269 const Verbosity& verbosity) {
2270 // Locate vmr_species among abs_species
2271 Index iq = -1;
2272 for (Index i = 0; i < abs_species.nelem(); i++) {
2273 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
2274 iq = i;
2275 break;
2276 }
2277 }
2278 if (iq < 0) {
2279 ostringstream os;
2280 os << "Could not find " << species << " in *abs_species*.\n";
2281 throw std::runtime_error(os.str());
2282 }
2283
2284 Tensor3 original_field, perturbed_field;
2285 original_field = vmr_field(iq,joker,joker,joker);
2286 AtmFieldPerturbAtmGrids(perturbed_field,
2287 atmosphere_dim,
2288 p_grid,
2289 lat_grid,
2290 lon_grid,
2291 original_field,
2292 pert_index,
2293 pert_size,
2294 pert_mode,
2295 verbosity);
2296 vmr_field(iq,joker,joker,joker) = perturbed_field;
2297}
2298
2299
Declarations required for the calculation of absorption coefficients.
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition array.h:173
base max(const Array< base > &x)
Max function.
Definition array.h:128
Array< ArrayOfIndex > ArrayOfArrayOfIndex
Definition array.h:121
The global header file for ARTS.
Common ARTS conversions.
Vector time_vector(const ArrayOfTime &times)
Converts from each Time to seconds and returns as Vector.
Definition artstime.cc:149
void chk_atm_field(const String &x_name, ConstTensor3View x, const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, const bool &chk_lat90)
chk_atm_field (simple fields)
The Agenda class.
void append(const String &methodname, const TokVal &keywordvalue)
Appends methods to an agenda.
void check(Workspace &ws_in, const Verbosity &verbosity)
Checks consistency of an agenda.
void set_name(const String &nname)
Set agenda name.
This can be used to make arrays out of anything.
Definition array.h:31
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition jacobian.h:310
const String & Mode() const
Returns the mode.
Definition jacobian.h:379
Jacobian::Target & Target()
Get the Jacobian Target.
Definition jacobian.h:416
const String & SubSubtag() const
Returns the sub-sub-tag.
Definition jacobian.h:364
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition jacobian.h:393
const String & Subtag() const
Returns the sub-tag.
Definition jacobian.h:350
Workspace class.
Internal cloudbox functions.
#define _U_
Definition config.h:177
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define DEBUG_ONLY(...)
Definition debug.h:55
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
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 polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition jacobian.cc:831
void transform_jacobian(Matrix &jacobian, const Vector x, const ArrayOfRetrievalQuantity &jqs)
Applies both functional and affine transformations.
Definition jacobian.cc:74
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition jacobian.cc:216
bool check_retrieval_grids(ArrayOfVector &grids, ostringstream &os, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &p_retr, const Vector &lat_retr, const Vector &lon_retr, const String &p_retr_name, const String &lat_retr_name, const String &lon_retr_name, const Index &dim)
Check that the retrieval grids are defined for each atmosphere dim.
Definition jacobian.cc:624
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
Definition jacobian.cc:29
Routines for setting up the jacobian.
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition jacobian.h:529
Jacobian::Line select_derivativeLineShape(const String &var, const String &coeff)
Return the derivative type based on string input.
void jacobianAdjustAndTransform(Matrix &jacobian, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Verbosity &)
WORKSPACE METHOD: jacobianAdjustAndTransform.
void jacobianAddWind(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dfrequency, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddWind.
void jacobianAddSpecialSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddSpecialSpecies.
void jacobianCalcFreqStretch(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &mblock_dlos, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: jacobianCalcFreqStretch.
void jacobianClose(Workspace &ws_in, Index &jacobian_do, Agenda &jacobian_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianClose.
Definition m_jacobian.cc:50
void vmr_fieldPerturb(Tensor4 &vmr_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Vector &p_ret_grid, const Vector &lat_ret_grid, const Vector &lon_ret_grid, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: vmr_fieldPerturb.
void jacobianCalcPolyfit(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &poly_coeff, const Verbosity &)
WORKSPACE METHOD: jacobianCalcPolyfit.
void jacobianAddMagField(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dB, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddMagField.
void jacobianCalcPointingZaRecalc(Workspace &ws, Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const ArrayOfTime &sensor_time, const String &iy_unit, const Agenda &iy_main_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianCalcPointingZaRecalc.
void jacobianFromTwoY(Matrix &jacobian, const Vector &y_pert, const Vector &y, const Numeric &pert_size, const Verbosity &)
WORKSPACE METHOD: jacobianFromTwoY.
void AtmFieldPerturb(Tensor3 &perturbed_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &original_field, const Vector &p_ret_grid, const Vector &lat_ret_grid, const Vector &lon_ret_grid, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &)
WORKSPACE METHOD: AtmFieldPerturb.
void jacobianFromYbatch(Matrix &jacobian, const ArrayOfVector &ybatch, const Vector &y, const Numeric &pert_size, const Verbosity &)
WORKSPACE METHOD: jacobianFromYbatch.
void jacobianAddSurfaceQuantity(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddSurfaceQuantity.
void jacobianAddShapeCatalogParameter(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const QuantumIdentifier &line_identity, const String &species, const String &variable, const String &coefficient, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddShapeCatalogParameter.
void jacobianInit(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Verbosity &)
WORKSPACE METHOD: jacobianInit.
Definition m_jacobian.cc:65
void jacobianAddFreqStretch(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqStretch.
void jacobianAddPointingZa(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Matrix &sensor_pos, const ArrayOfTime &sensor_time, const Index &poly_order, const String &calcmode, const Numeric &dza, const Verbosity &)
WORKSPACE METHOD: jacobianAddPointingZa.
void particle_bulkprop_fieldPerturb(Tensor4 &particle_bulkprop_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfString &particle_bulkprop_names, const String &particle_type, const Vector &p_ret_grid, const Vector &lat_ret_grid, const Vector &lon_ret_grid, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: particle_bulkprop_fieldPerturb.
void jacobianAddNLTEs(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const ArrayOfQuantumIdentifier &energy_level_identities, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddNLTEs.
void vmr_fieldPerturbAtmGrids(Tensor4 &vmr_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: vmr_fieldPerturbAtmGrids.
void jacobianSetAffineTransformation(ArrayOfRetrievalQuantity &jqs, const Matrix &transformation_matrix, const Vector &offset_vector, const Verbosity &)
WORKSPACE METHOD: jacobianSetAffineTransformation.
void jacobianAddNLTE(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const QuantumIdentifier &energy_level_identity, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddNLTE.
void IndexNumberOfAtmosphericPoints(Index &n, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Verbosity &)
WORKSPACE METHOD: IndexNumberOfAtmosphericPoints.
void jacobianCalcPointingZaInterp(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &DEBUG_ONLY(sensor_los), const Matrix &mblock_dlos, const Sparse &sensor_response, const ArrayOfTime &sensor_time, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
void jacobianCalcDoNothing(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Verbosity &)
WORKSPACE METHOD: jacobianCalcDoNothing.
Definition m_jacobian.cc:40
void jacobianAddSinefit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Vector &period_lengths, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &)
WORKSPACE METHOD: jacobianAddSinefit.
void jacobianAddTemperature(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &hse, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddTemperature.
void particle_bulkprop_fieldPerturbAtmGrids(Tensor4 &particle_bulkprop_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfString &particle_bulkprop_names, const String &particle_type, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: particle_bulkprop_fieldPerturbAtmGrids.
void jacobianCalcFreqShift(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &mblock_dlos, const Sparse &sensor_response, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: jacobianCalcFreqShift.
void jacobianAddBasicCatalogParameter(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const QuantumIdentifier &catalog_identity, const String &catalog_parameter, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddBasicCatalogParameter.
void jacobianAddAbsSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Index &for_species_tag, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddAbsSpecies.
Definition m_jacobian.cc:89
void jacobianAddShapeCatalogParameters(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfQuantumIdentifier &line_identities, const ArrayOfString &species, const ArrayOfString &variables, const ArrayOfString &coefficients, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddShapeCatalogParameters.
void jacobianAddScatSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddScatSpecies.
void jacobianAddFreqShift(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqShift.
void AtmFieldPerturbAtmGrids(Tensor3 &perturbed_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &original_field, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &)
WORKSPACE METHOD: AtmFieldPerturbAtmGrids.
void jacobianSetFuncTransformation(ArrayOfRetrievalQuantity &jqs, const String &transformation_func, const Numeric &z_min, const Numeric &z_max, const Verbosity &)
WORKSPACE METHOD: jacobianSetFuncTransformation.
void jacobianAddBasicCatalogParameters(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfQuantumIdentifier &catalog_identities, const ArrayOfString &catalog_parameters, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddBasicCatalogParameters.
void jacobianOff(Workspace &ws, Index &jacobian_do, Agenda &jacobian_agenda, ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianOff.
Definition m_jacobian.cc:75
void jacobianAddPolyfit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Index &poly_order, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &)
WORKSPACE METHOD: jacobianAddPolyfit.
void jacobianCalcSinefit(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &period_index, const Verbosity &)
WORKSPACE METHOD: jacobianCalcSinefit.
Workspace methods and template functions for supergeneric XML IO.
void reshape(Tensor3View X, ConstVectorView x)
reshape
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition messages.h:189
#define CREATE_OUT2
Definition messages.h:188
Array< String > ArrayOfString
An array of Strings.
Definition mystring.h:202
constexpr Numeric two_pi
Two times pi.
constexpr std::string_view bath_broadening
Name for bath broadening in printing and reading user input.
constexpr std::string_view self_broadening
Name for self broadening in printing and reading user input.
This file contains declerations of functions of physical character.
Quantum::Number::GlobalState QuantumIdentifier
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, Matrix &geo_pos_matrix, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
Performs calculations for one measurement block, on iy-level.
Definition rte.cc:1591
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
Returns the "range" of y corresponding to a measurement block.
Definition rte.cc:1095
Declaration of functions in rte.cc.
void regrid_atmfield_by_gp_oem(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of atmospheric field OEM-type.
void get_gp_rq_to_atmgrids(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, Index &n_p, Index &n_lat, Index &n_lon, const ArrayOfVector &ret_grids, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
Species::Tag SpeciesTag
Holds all information required for individual partial derivatives.
Definition jacobian.h:93
Numeric perturbation
Perturbations for methods where theoretical computations are impossible or plain slow.
Definition jacobian.h:110
bool sameTargetType(const Target &other) const noexcept
Definition jacobian.h:177
Species::Species species_id
Species ID for line parameters.
Definition jacobian.h:122
A logical struct for global quantum numbers with species identifiers.
Species::Species Species() const noexcept
#define v
#define w
#define a
#define c