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