ARTS 2.5.0 (git: 9ee3ac6c)
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 "constants.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, verbosity);
78 jacobian_do = 1;
79}
80
81/* Workspace method: Doxygen documentation will be auto-generated */
82void jacobianInit(ArrayOfRetrievalQuantity& jacobian_quantities,
83 Agenda& jacobian_agenda,
84 const Verbosity&) {
85 jacobian_quantities.resize(0);
86 jacobian_agenda = Agenda();
87 jacobian_agenda.set_name("jacobian_agenda");
88}
89
90/* Workspace method: Doxygen documentation will be auto-generated */
91void jacobianOff(Index& jacobian_do,
92 Agenda& jacobian_agenda,
93 ArrayOfRetrievalQuantity& jacobian_quantities,
94 const Verbosity& verbosity) {
95 jacobian_do = 0;
96 jacobianInit(jacobian_quantities, jacobian_agenda, verbosity);
97}
98
99//----------------------------------------------------------------------------
100// Absorption species:
101//----------------------------------------------------------------------------
102
103/* Workspace method: Doxygen documentation will be auto-generated */
106 Agenda& jacobian_agenda,
107 const Index& atmosphere_dim,
108 const Vector& p_grid,
109 const Vector& lat_grid,
110 const Vector& lon_grid,
111 const Vector& rq_p_grid,
112 const Vector& rq_lat_grid,
113 const Vector& rq_lon_grid,
114 const String& species,
115 const String& mode,
116 const Index& for_species_tag,
117 const Verbosity& verbosity) {
120
122 if (not for_species_tag) {
123 ArrayOfSpeciesTag test(species);
124 if (test.nelem() not_eq 1)
125 throw std::runtime_error(
126 "Trying to add a species as a species tag of multiple species.\n"
127 "This is not supported. Please give just a single species instead.\n"
128 "Otherwise consider if you intended for_species_tag to be evaluated true.\n");
129 qi = QuantumIdentifier(test[0].Isotopologue(), Quantum::IdentifierType::All);
130 }
131
132 // Check that this species is not already included in the jacobian.
133 for (Index it = 0; it < jq.nelem(); it++) {
134 ARTS_USER_ERROR_IF (jq[it] == Jacobian::Special::ArrayOfSpeciesTagVMR && jq[it].Subtag() == species,
135 "The gas species:\n", species, "\nis already included in *jacobian_quantities*.")
136 ARTS_USER_ERROR_IF (jq[it] == Jacobian::Line::VMR and jq[it].QuantumIdentity() == qi,
137 "The atmospheric species of:\n", species, "\nis already included in *jacobian_quantities*\n"
138 "as: ", jq[it])
139 }
140
141 // Check retrieval grids, here we just check the length of the grids
142 // vs. the atmosphere dimension
143 ArrayOfVector grids(atmosphere_dim);
144 {
145 ostringstream os;
146 if (!check_retrieval_grids(grids,
147 os,
148 p_grid,
149 lat_grid,
150 lon_grid,
151 rq_p_grid,
152 rq_lat_grid,
153 rq_lon_grid,
154 "retrieval pressure grid",
155 "retrieval latitude grid",
156 "retrievallongitude_grid",
157 atmosphere_dim))
158 throw runtime_error(os.str());
159 }
160
161 // Check that mode is correct
162 if (mode != "vmr" && mode != "nd" && mode != "rel" && mode != "rh" &&
163 mode != "q") {
164 throw runtime_error(
165 "The retrieval mode can only be \"vmr\", \"nd\", "
166 "\"rel\", \"rh\" or \"q\".");
167 }
168 if ((mode == "rh" || mode == "q") && species.substr(0, 3) != "H2O") {
169 throw runtime_error(
170 "Retrieval modes \"rh\" and \"q\" can only be applied "
171 "on species starting with H2O.");
172 }
173
174 // Create the new retrieval quantity
176 rq.Subtag(species);
177 rq.Mode(mode);
178 rq.Grids(grids);
179 if (for_species_tag == 0) {
180 rq.Target(Jacobian::Target(Jacobian::Line::VMR, qi, qi.Species()));
181 } else {
182 ArrayOfSpeciesTag test(species);
183 rq.Target(Jacobian::Target(Jacobian::Special::ArrayOfSpeciesTagVMR, test));
184 }
185 rq.Target().Perturbation(0.001);
186
187 // Add it to the *jacobian_quantities*
188 jq.push_back(rq);
189
190 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
191}
192
193//----------------------------------------------------------------------------
194// Frequency shift
195//----------------------------------------------------------------------------
196
197/* Workspace method: Doxygen documentation will be auto-generated */
199 ArrayOfRetrievalQuantity& jacobian_quantities,
200 Agenda& jacobian_agenda,
201 const Vector& f_grid,
202 const Numeric& df,
203 const Verbosity&) {
204 // Check that this jacobian type is not already included.
205 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
206 if (jacobian_quantities[it] == Jacobian::Sensor::FrequencyShift) {
207 ostringstream os;
208 os << "Fit of frequency shift is already included in\n"
209 << "*jacobian_quantities*.";
210 throw runtime_error(os.str());
211 }
212 }
213
214 // Checks of frequencies
215 if (df <= 0) throw runtime_error("The argument *df* must be > 0.");
216 if (df > 1e6)
217 throw runtime_error("The argument *df* is not allowed to exceed 1 MHz.");
218 const Index nf = f_grid.nelem();
219 if (nf < 2)
220 throw runtime_error(
221 "Frequency shifts and *f_grid* of length 1 can "
222 "not be combined.");
223 const Numeric maxdf = f_grid[nf - 1] - f_grid[nf - 2];
224 if (df > maxdf) {
225 ostringstream os;
226 os << "The value of *df* is too big with respect to spacing of "
227 << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
228 << "between the two last elements of *f_grid*.\n"
229 << "This spacing is : " << maxdf / 1e3 << " kHz\n"
230 << "The value of df is: " << df / 1e3 << " kHz";
231 throw runtime_error(os.str());
232 }
233
234 // Create the new retrieval quantity
236 rq.Target() = Jacobian::Target(Jacobian::Sensor::FrequencyShift);
237 rq.Mode("");
238 rq.Target().Perturbation(df);
239
240 // Dummy vector of length 1
241 Vector grid(1, 0);
242 ArrayOfVector grids(1, grid);
243 rq.Grids(grids);
244
245 // Add it to the *jacobian_quantities*
246 jacobian_quantities.push_back(rq);
247
248 // Add corresponding calculation method to the jacobian agenda
249 jacobian_agenda.append("jacobianCalcFreqShift", "");
250}
251
252/* Workspace method: Doxygen documentation will be auto-generated */
254 const Index& mblock_index,
255 const Vector& iyb,
256 const Vector& yb,
257 const Index& stokes_dim,
258 const Vector& f_grid,
259 const Matrix& mblock_dlos_grid,
260 const Sparse& sensor_response,
261 const ArrayOfRetrievalQuantity& jacobian_quantities,
262 const Verbosity&) {
263 // Set some useful (and needed) variables.
265 ArrayOfIndex ji;
266
267 // Find the retrieval quantity related to this method.
268 // This works since the combined MainTag and Subtag is individual.
269 bool found = false;
270 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
271 if (jacobian_quantities[n] == Jacobian::Sensor::FrequencyShift) {
272 bool any_affine;
273 ArrayOfArrayOfIndex jacobian_indices;
275 jacobian_indices, any_affine, jacobian_quantities, true);
276 //
277 found = true;
278 rq = jacobian_quantities[n];
279 ji = jacobian_indices[n];
280 }
281 }
282 if (!found) {
283 throw runtime_error(
284 "There is no such frequency retrieval quantity defined.\n");
285 }
286
287 // Check that sensor_response is consistent with yb and iyb
288 //
289 if (sensor_response.nrows() != yb.nelem())
290 throw runtime_error("Mismatch in size between *sensor_response* and *yb*.");
291 if (sensor_response.ncols() != iyb.nelem())
292 throw runtime_error(
293 "Mismatch in size between *sensor_response* and *iyb*.");
294
295 // Get disturbed (part of) y
296 //
297 const Index n1y = sensor_response.nrows();
298 Vector dy(n1y);
299 {
300 const Index nf2 = f_grid.nelem();
301 const Index nlos2 = mblock_dlos_grid.nrows();
302 const Index niyb = nf2 * nlos2 * stokes_dim;
303
304 // Interpolation weights
305 //
306 constexpr Index porder = 3;
307 //
308 Vector fg_new = f_grid, iyb2(niyb);
309 fg_new += rq.Target().Perturbation();
310 //
311 const auto lag = Interpolation::FixedLagrangeVector<porder>(fg_new, f_grid);
312 const auto itw = interpweights(lag);
313
314 // Do interpolation
315 for (Index ilos = 0; ilos < nlos2; ilos++) {
316 const Index row0 = ilos * nf2 * stokes_dim;
317
318 for (Index is = 0; is < stokes_dim; is++) {
319 iyb2[Range(row0 + is, nf2, stokes_dim)] = reinterp(iyb[Range(row0 + is, nf2, stokes_dim)], itw, lag);
320 }
321 }
322
323 // Determine difference
324 //
325 mult(dy, sensor_response, iyb2);
326 //
327 for (Index i = 0; i < n1y; i++) {
328 dy[i] = (dy[i] - yb[i]) / rq.Target().Perturbation();
329 }
330 }
331
332 //--- Set jacobian ---
333 ARTS_ASSERT(rq.Grids()[0].nelem() == 1);
334 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
335 jacobian(rowind, ji[0]) = dy;
336}
337
338//----------------------------------------------------------------------------
339// Frequency stretch
340//----------------------------------------------------------------------------
341
342/* Workspace method: Doxygen documentation will be auto-generated */
344 ArrayOfRetrievalQuantity& jacobian_quantities,
345 Agenda& jacobian_agenda,
346 const Vector& f_grid,
347 const Numeric& df,
348 const Verbosity&) {
349 // Check that this jacobian type is not already included.
350 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
351 if (jacobian_quantities[it] == Jacobian::Sensor::FrequencyStretch) {
352 ostringstream os;
353 os << "Fit of frequency stretch is already included in\n"
354 << "*jacobian_quantities*.";
355 throw runtime_error(os.str());
356 }
357 }
358
359 // Checks of df
360 if (df <= 0) throw runtime_error("The argument *df* must be > 0.");
361 if (df > 1e6)
362 throw runtime_error("The argument *df* is not allowed to exceed 1 MHz.");
363 const Index nf = f_grid.nelem();
364 const Numeric maxdf = f_grid[nf - 1] - f_grid[nf - 2];
365 if (df > maxdf) {
366 ostringstream os;
367 os << "The value of *df* is too big with respect to spacing of "
368 << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
369 << "between the two last elements of *f_grid*.\n"
370 << "This spacing is : " << maxdf / 1e3 << " kHz\n"
371 << "The value of df is: " << df / 1e3 << " kHz";
372 throw runtime_error(os.str());
373 }
374
375 // Create the new retrieval quantity
377 rq.Target() = Jacobian::Target(Jacobian::Sensor::FrequencyStretch);
378 rq.Mode("");
379 rq.Target().Perturbation(df);
380
381 // Dummy vector of length 1
382 Vector grid(1, 0);
383 ArrayOfVector grids(1, grid);
384 rq.Grids(grids);
385
386 // Add it to the *jacobian_quantities*
387 jacobian_quantities.push_back(rq);
388
389 // Add corresponding calculation method to the jacobian agenda
390 jacobian_agenda.append("jacobianCalcFreqStretch", "");
391}
392
393/* Workspace method: Doxygen documentation will be auto-generated */
395 Matrix& jacobian,
396 const Index& mblock_index,
397 const Vector& iyb,
398 const Vector& yb,
399 const Index& stokes_dim,
400 const Vector& f_grid,
401 const Matrix& mblock_dlos_grid,
402 const Sparse& sensor_response,
403 const ArrayOfIndex& sensor_response_pol_grid,
404 const Vector& sensor_response_f_grid,
405 const Matrix& sensor_response_dlos_grid,
406 const ArrayOfRetrievalQuantity& jacobian_quantities,
407 const Verbosity&) {
408 // The code here is close to identical to the one for Shift. The main
409 // difference is that dy is weighted with poly_order 1 basis function.
410
411 // Set some useful (and needed) variables.
413 ArrayOfIndex ji;
414
415 // Find the retrieval quantity related to this method.
416 // This works since the combined MainTag and Subtag is individual.
417 bool found = false;
418 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
419 if (jacobian_quantities[n] == Jacobian::Sensor::FrequencyStretch) {
420 bool any_affine;
421 ArrayOfArrayOfIndex jacobian_indices;
423 jacobian_indices, any_affine, jacobian_quantities, true);
424 //
425 found = true;
426 rq = jacobian_quantities[n];
427 ji = jacobian_indices[n];
428 }
429 }
430 if (!found) {
431 throw runtime_error(
432 "There is no such frequency retrieval quantity defined.\n");
433 }
434
435 // Check that sensor_response is consistent with yb and iyb
436 //
437 if (sensor_response.nrows() != yb.nelem())
438 throw runtime_error("Mismatch in size between *sensor_response* and *yb*.");
439 if (sensor_response.ncols() != iyb.nelem())
440 throw runtime_error(
441 "Mismatch in size between *sensor_response* and *iyb*.");
442
443 // Get disturbed (part of) y
444 //
445 const Index n1y = sensor_response.nrows();
446 Vector dy(n1y);
447 {
448 const Index nf2 = f_grid.nelem();
449 const Index nlos2 = mblock_dlos_grid.nrows();
450 const Index niyb = nf2 * nlos2 * stokes_dim;
451
452 // Interpolation weights
453 //
454 constexpr Index porder = 3;
455 //
456 Vector fg_new = f_grid, iyb2(niyb);
457 //
458 fg_new += rq.Target().Perturbation();
459 //
460 const auto lag = Interpolation::FixedLagrangeVector<porder>(fg_new, f_grid);
461 const auto itw = interpweights(lag);
462
463 // Do interpolation
464 for (Index ilos = 0; ilos < nlos2; ilos++) {
465 const Index row0 = ilos * nf2 * stokes_dim;
466
467 for (Index is = 0; is < stokes_dim; is++) {
468 iyb2[Range(row0 + is, nf2, stokes_dim)] = reinterp(iyb[Range(row0 + is, nf2, stokes_dim)], itw, lag);
469 }
470 }
471
472 // Determine difference
473 //
474 mult(dy, sensor_response, iyb2);
475 //
476 for (Index i = 0; i < n1y; i++) {
477 dy[i] = (dy[i] - yb[i]) / rq.Target().Perturbation();
478 }
479
480 // dy above corresponds now to shift. Convert to stretch:
481 //
482 Vector w;
483 polynomial_basis_func(w, sensor_response_f_grid, 1);
484 //
485 const Index nf = sensor_response_f_grid.nelem();
486 const Index npol = sensor_response_pol_grid.nelem();
487 const Index nlos = sensor_response_dlos_grid.nrows();
488 //
489 for (Index l = 0; l < nlos; l++) {
490 for (Index f = 0; f < nf; f++) {
491 const Index row1 = (l * nf + f) * npol;
492 for (Index p = 0; p < npol; p++) {
493 dy[row1 + p] *= w[f];
494 }
495 }
496 }
497 }
498
499 //--- Set jacobians ---
500 ARTS_ASSERT(rq.Grids()[0].nelem() == 1);
501 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
502 jacobian(rowind, ji[0]) = dy;
503}
504
505//----------------------------------------------------------------------------
506// Pointing:
507//----------------------------------------------------------------------------
508
509/* Workspace method: Doxygen documentation will be auto-generated */
511 ArrayOfRetrievalQuantity& jacobian_quantities,
512 Agenda& jacobian_agenda,
513 const Matrix& sensor_pos,
514 const ArrayOfTime& sensor_time,
515 const Index& poly_order,
516 const String& calcmode,
517 const Numeric& dza,
518 const Verbosity&) {
519 // Check that poly_order is -1 or positive
520 if (poly_order < -1)
521 throw runtime_error(
522 "The polynomial order has to be positive or -1 for gitter.");
523
524 // Check that this jacobian type is not already included.
525 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
526 if (jacobian_quantities[it].Target().isPointing()) {
527 ostringstream os;
528 os << "Fit of zenith angle pointing off-set is already included in\n"
529 << "*jacobian_quantities*.";
530 throw runtime_error(os.str());
531 }
532 }
533
534 // Checks of dza
535 if (dza <= 0) throw runtime_error("The argument *dza* must be > 0.");
536 if (dza > 0.1)
537 throw runtime_error("The argument *dza* is not allowed to exceed 0.1 deg.");
538
539 // Check that sensor_time is consistent with sensor_pos
540 if (sensor_time.nelem() != sensor_pos.nrows()) {
541 ostringstream os;
542 os << "The WSV *sensor_time* must be defined for every "
543 << "measurement block.\n";
544 throw runtime_error(os.str());
545 }
546
547 // Do not allow that *poly_order* is not too large compared to *sensor_time*
548 if (poly_order > sensor_time.nelem() - 1) {
549 throw runtime_error(
550 "The polynomial order can not be >= length of *sensor_time*.");
551 }
552
553 // Create the new retrieval quantity
555 if (calcmode == "recalc") {
556 rq.Target() = Jacobian::Target(Jacobian::Sensor::PointingZenithRecalc);
557 jacobian_agenda.append("jacobianCalcPointingZaRecalc", "");
558 } else if (calcmode == "interp") {
559 rq.Target() = Jacobian::Target(Jacobian::Sensor::PointingZenithInterp);
560 jacobian_agenda.append("jacobianCalcPointingZaInterp", "");
561 } else
562 throw runtime_error(
563 "Possible choices for *calcmode* are \"recalc\" and \"interp\".");
564 rq.Target().Perturbation(dza);
565
566 // To store the value or the polynomial order, create a vector with length
567 // poly_order+1, in case of gitter set the size of the grid vector to be the
568 // number of measurement blocks, all elements set to -1.
569 Vector grid(0, poly_order + 1, 1);
570 if (poly_order == -1) {
571 grid.resize(sensor_pos.nrows());
572 grid = -1.0;
573 }
574 ArrayOfVector grids(1, grid);
575 rq.Grids(grids);
576
577 // Add it to the *jacobian_quantities*
578 jacobian_quantities.push_back(rq);
579}
580
581/* Workspace method: Doxygen documentation will be auto-generated */
583 Matrix& jacobian,
584 const Index& mblock_index,
585 const Vector& iyb,
586 const Vector& yb _U_,
587 const Index& stokes_dim,
588 const Vector& f_grid,
589 const Matrix& DEBUG_ONLY(sensor_los),
590 const Matrix& mblock_dlos_grid,
591 const Sparse& sensor_response,
592 const ArrayOfTime& sensor_time,
593 const ArrayOfRetrievalQuantity& jacobian_quantities,
594 const Verbosity&) {
595 if (mblock_dlos_grid.nrows() < 2)
596 throw runtime_error(
597 "The method demands that *mblock_dlos_grid* has "
598 "more than one row.");
599
600 if (!(is_increasing(mblock_dlos_grid(joker, 0)) ||
601 is_decreasing(mblock_dlos_grid(joker, 0))))
602 throw runtime_error(
603 "The method demands that the zenith angles in "
604 "*mblock_dlos_grid* are sorted (increasing or decreasing).");
605
606 // Set some useful variables.
608 ArrayOfIndex ji;
609
610 // Find the retrieval quantity related to this method.
611 // This works since the combined MainTag and Subtag is individual.
612 bool found = false;
613 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
614 if (jacobian_quantities[n] == Jacobian::Sensor::PointingZenithInterp) {
615 bool any_affine;
616 ArrayOfArrayOfIndex jacobian_indices;
618 jacobian_indices, any_affine, jacobian_quantities, true);
619 //
620 found = true;
621 rq = jacobian_quantities[n];
622 ji = jacobian_indices[n];
623 }
624 }
625 if (!found) {
626 throw runtime_error(
627 "There is no such pointing retrieval quantity defined.\n");
628 }
629
630 // Get "dy", by inter/extra-polation of existing iyb
631 //
632 const Index n1y = sensor_response.nrows();
633 Vector dy(n1y);
634 {
635 // Sizes
636 const Index nf = f_grid.nelem();
637 const Index nza = mblock_dlos_grid.nrows();
638
639 // Shifted zenith angles
640 Vector za1 = mblock_dlos_grid(joker, 0);
641 za1 -= rq.Target().Perturbation();
642 Vector za2 = mblock_dlos_grid(joker, 0);
643 za2 += rq.Target().Perturbation();
644
645 // Find interpolation weights
646 ArrayOfGridPos gp1(nza), gp2(nza);
647 gridpos(
648 gp1, mblock_dlos_grid(joker, 0), za1, 1e6); // Note huge extrapolation!
649 gridpos(
650 gp2, mblock_dlos_grid(joker, 0), za2, 1e6); // Note huge extrapolation!
651 Matrix itw1(nza, 2), itw2(nza, 2);
652 interpweights(itw1, gp1);
653 interpweights(itw2, gp2);
654
655 // Make interpolation (for all azimuth angles, frequencies and Stokes)
656 //
657 Vector iyb1(iyb.nelem()), iyb2(iyb.nelem());
658 //
659 for (Index iza = 0; iza < nza; iza++) {
660 for (Index iv = 0; iv < nf; iv++) {
661 for (Index is = 0; is < stokes_dim; is++) {
662 const Range r(iv * stokes_dim + is, nza, nf * stokes_dim);
663 interp(iyb1[r], itw1, iyb[r], gp1);
664 interp(iyb2[r], itw2, iyb[r], gp2);
665 }
666 }
667 }
668
669 // Apply sensor and take difference
670 //
671 Vector y1(n1y), y2(n1y);
672 mult(y1, sensor_response, iyb1);
673 mult(y2, sensor_response, iyb2);
674 //
675 for (Index i = 0; i < n1y; i++) {
676 dy[i] = (y2[i] - y1[i]) / (2 * rq.Target().Perturbation());
677 }
678 }
679
680 //--- Create jacobians ---
681
682 const Index lg = rq.Grids()[0].nelem();
683 const Index it = ji[0];
684 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
685 const Index row0 = rowind.get_start();
686
687 // Handle pointing "jitter" seperately
688 if (rq.Grids()[0][0] == -1) // Not all values are set here,
689 { // but should already have been
690 ARTS_ASSERT(lg == sensor_los.nrows()); // set to 0
691 ARTS_ASSERT(rq.Grids()[0][mblock_index] == -1);
692 jacobian(rowind, it + mblock_index) = dy;
693 }
694
695 // Polynomial representation
696 else {
697 Vector w;
698 for (Index c = 0; c < lg; c++) {
699 ARTS_ASSERT(Numeric(c) == rq.Grids()[0][c]);
700 //
701 polynomial_basis_func(w, time_vector(sensor_time), c);
702 //
703 for (Index i = 0; i < n1y; i++) {
704 jacobian(row0 + i, it + c) = w[mblock_index] * dy[i];
705 }
706 }
707 }
708}
709
710/* Workspace method: Doxygen documentation will be auto-generated */
712 Workspace& ws,
713 Matrix& jacobian,
714 const Index& mblock_index,
715 const Vector& iyb _U_,
716 const Vector& yb,
717 const Index& atmosphere_dim,
718 const EnergyLevelMap& nlte_field,
719 const Index& cloudbox_on,
720 const Index& stokes_dim,
721 const Vector& f_grid,
722 const Matrix& sensor_pos,
723 const Matrix& sensor_los,
724 const Matrix& transmitter_pos,
725 const Matrix& mblock_dlos_grid,
726 const Sparse& sensor_response,
727 const ArrayOfTime& sensor_time,
728 const String& iy_unit,
729 const Agenda& iy_main_agenda,
730 const Agenda& geo_pos_agenda,
731 const ArrayOfRetrievalQuantity& jacobian_quantities,
732 const Verbosity& verbosity) {
733 // Set some useful variables.
735 ArrayOfIndex ji;
736
737 // Find the retrieval quantity related to this method.
738 // This works since the combined MainTag and Subtag is individual.
739 bool found = false;
740 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
741 if (jacobian_quantities[n] == Jacobian::Sensor::PointingZenithRecalc) {
742 bool any_affine;
743 ArrayOfArrayOfIndex jacobian_indices;
745 jacobian_indices, any_affine, jacobian_quantities, true);
746 //
747 found = true;
748 rq = jacobian_quantities[n];
749 ji = jacobian_indices[n];
750 }
751 }
752 if (!found) {
753 throw runtime_error(
754 "There is no such pointing retrieval quantity defined.\n");
755 }
756
757 // Get "dy", by calling iyb_calc with shifted sensor_los.
758 //
759 const Index n1y = sensor_response.nrows();
760 Vector dy(n1y);
761 {
762 Vector iyb2;
763 Matrix los = sensor_los;
764 Matrix geo_pos;
765 ArrayOfVector iyb_aux;
766 ArrayOfMatrix diyb_dx;
767
768 los(joker, 0) += rq.Target().Perturbation();
769
770 iyb_calc(ws,
771 iyb2,
772 iyb_aux,
773 diyb_dx,
774 geo_pos,
775 mblock_index,
776 atmosphere_dim,
777 nlte_field,
778 cloudbox_on,
779 stokes_dim,
780 f_grid,
781 sensor_pos,
782 los,
783 transmitter_pos,
784 mblock_dlos_grid,
785 iy_unit,
786 iy_main_agenda,
787 geo_pos_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) {
1390 rq.Target(Jacobian::Target(Jacobian::Atm::WindU));
1391 break;
1393 rq.Target(Jacobian::Target(Jacobian::Atm::WindV));
1394 break;
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) {
1469 rq.Target(Jacobian::Target(Jacobian::Atm::MagneticU));
1470 break;
1472 rq.Target(Jacobian::Target(Jacobian::Atm::MagneticV));
1473 break;
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 if (line_identity.type not_eq Quantum::IdentifierType::Transition)
1537 throw std::runtime_error("Identity has to identify a line");
1538
1539 const auto jpt = select_derivativeLineShape(variable, coefficient);
1540
1541 out3 << "Attempting to create RT tag for " << line_identity << " " << variable
1542 << " " << coefficient << " for ";
1543
1544 // Create the quantity
1546 rq.Mode(species);
1547 rq.Grids(ArrayOfVector(0, Vector(0)));
1548
1549 // Map the species
1550 if (species == LineShape::self_broadening) {
1551 rq.Target(Jacobian::Target(jpt, line_identity, line_identity.Species()));
1552 } else if (species == LineShape::bath_broadening) {
1553 rq.Target(Jacobian::Target(jpt, line_identity, Species::Species::Bath));
1554 } else {
1555 rq.Target(Jacobian::Target(jpt, line_identity, SpeciesTag(species).Spec()));
1556 }
1557 out3 << species << ' ' << rq.Target().LineSpecies() << "\n";
1558
1559 // Test this is not a copy
1560 for (auto& q : jq)
1561 if (q.HasSameInternalsAs(rq))
1562 throw std::runtime_error("Error with copies of the quantities");
1563
1564 // Append and do housekeeping
1565 jq.push_back(rq);
1566 out3 << "Creation was successful!\n";
1567 jacobian_agenda.append("jacobianCalcDoNothing",
1568 TokVal()); // old code activation
1569}
1570
1571/* Workspace method: Doxygen documentation will be auto-generated */
1573 Workspace& ws,
1575 Agenda& jacobian_agenda,
1576 const ArrayOfQuantumIdentifier& line_identities,
1577 const ArrayOfString& species,
1578 const ArrayOfString& variables,
1579 const ArrayOfString& coefficients,
1580 const Verbosity& verbosity) {
1581 if (not(line_identities.nelem() or species.nelem() or variables.nelem() or
1582 coefficients.nelem()))
1583 throw std::runtime_error("Must have at least 1-long lists for all GINs");
1584
1585 ArrayOfString vars;
1586 if (variables[0] == "ALL")
1587 vars = ArrayOfString(LineShape::enumstrs::VariableNames);
1588 else
1589 vars = variables;
1590
1591 ArrayOfString coeffs;
1592 if (coefficients[0] == "ALL")
1593 coeffs = ArrayOfString(Options::enumstrs::LineShapeCoeffNames);
1594 else
1595 coeffs = coefficients;
1596
1597 for (auto& l : line_identities)
1598 for (auto& s : species)
1599 for (auto& v : vars)
1600 for (auto& c : coeffs)
1602 ws, jq, jacobian_agenda, l, s, v, c, verbosity);
1603}
1604
1605/* Workspace method: Doxygen documentation will be auto-generated */
1608 Agenda& jacobian_agenda,
1609 const QuantumIdentifier& catalog_identity,
1610 const String& catalog_parameter,
1611 const Verbosity& verbosity) {
1613
1614 // Create the new retrieval quantity
1615 const auto opt = Options::toBasicCatParamJacobianOrThrow(catalog_parameter);
1617 switch (opt) {
1618 case Options::BasicCatParamJacobian::LineCenter:
1619 rq.Target(Jacobian::Target(Jacobian::Line::Center, catalog_identity, catalog_identity.Species()));
1620 break;
1622 rq.Target(Jacobian::Target(Jacobian::Line::Strength, catalog_identity, catalog_identity.Species()));
1623 break;
1624 case Options::BasicCatParamJacobian::FINAL:
1625 ARTS_ASSERT(false, "This error should be caught earlier")
1626 }
1627
1628 // Check that this is not already included in the jacobian.
1629 for (Index it = 0; it < jq.nelem(); it++) {
1630 ARTS_USER_ERROR_IF (rq.Target().sameTargetType(jq[it].Target()),
1631 "The catalog identifier:\n",
1632 catalog_identity, " for ID: ", catalog_identity, "\n"
1633 "is already included in jacobian_quantities*.")
1634 }
1635
1636 rq.Grids(ArrayOfVector(0, Vector(0)));
1637
1638 // Add it to the *jacobian_quantities*
1639 jq.push_back(rq);
1640
1641 out3 << " Calculations done by propagation matrix expressions.\n";
1642
1643 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1644}
1645
1646/* Workspace method: Doxygen documentation will be auto-generated */
1648 Workspace& ws,
1650 Agenda& jacobian_agenda,
1651 const ArrayOfQuantumIdentifier& catalog_identities,
1652 const ArrayOfString& catalog_parameters,
1653 const Verbosity& verbosity) {
1655
1656 out2 << " Adding " << catalog_identities.nelem() * catalog_parameters.nelem()
1657 << " expression(s) to the Jacobian calculations.\n";
1658
1659 for (auto& qi : catalog_identities)
1660 for (auto& param : catalog_parameters)
1662 ws, jq, jacobian_agenda, qi, param, verbosity);
1663}
1664
1665//----------------------------------------------------------------------------
1666// NLTE temperature:
1667//----------------------------------------------------------------------------
1668
1669/* Workspace method: Doxygen documentation will be auto-generated */
1672 Agenda& jacobian_agenda,
1673 const Index& atmosphere_dim,
1674 const Vector& p_grid,
1675 const Vector& lat_grid,
1676 const Vector& lon_grid,
1677 const Vector& rq_p_grid,
1678 const Vector& rq_lat_grid,
1679 const Vector& rq_lon_grid,
1680 const QuantumIdentifier& energy_level_identity,
1681 const Numeric& dx,
1682 const Verbosity& verbosity) {
1684
1685 // Check that this species is not already included in the jacobian.
1686 for (Index it = 0; it < jq.nelem(); it++) {
1687 if (jq[it] == Jacobian::Line::NLTE and
1688 jq[it].QuantumIdentity() == energy_level_identity) {
1689 ostringstream os;
1690 os << "The NLTE identifier:\n"
1691 << energy_level_identity << "\nis already included in "
1692 << "*jacobian_quantities*.";
1693 throw std::runtime_error(os.str());
1694 }
1695 }
1696
1697 // Check retrieval grids, here we just check the length of the grids
1698 // vs. the atmosphere dimension
1699 ArrayOfVector grids(atmosphere_dim);
1700 {
1701 ostringstream os;
1702 if (not check_retrieval_grids(grids,
1703 os,
1704 p_grid,
1705 lat_grid,
1706 lon_grid,
1707 rq_p_grid,
1708 rq_lat_grid,
1709 rq_lon_grid,
1710 "retrieval pressure grid",
1711 "retrieval latitude grid",
1712 "retrievallongitude_grid",
1713 atmosphere_dim))
1714 throw runtime_error(os.str());
1715 }
1716
1717 // Create the new retrieval quantity
1719 rq.Target(Jacobian::Target(Jacobian::Line::NLTE, energy_level_identity, energy_level_identity.Species()));
1720 rq.Target().Perturbation(dx);
1721 rq.Grids(grids);
1722
1723 // Add it to the *jacobian_quantities*
1724 jq.push_back(rq);
1725
1726 out3 << " Calculations done by propagation matrix expressions.\n";
1727
1728 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1729}
1730
1733 Agenda& jacobian_agenda,
1734 const Index& atmosphere_dim,
1735 const Vector& p_grid,
1736 const Vector& lat_grid,
1737 const Vector& lon_grid,
1738 const Vector& rq_p_grid,
1739 const Vector& rq_lat_grid,
1740 const Vector& rq_lon_grid,
1741 const ArrayOfQuantumIdentifier& energy_level_identities,
1742 const Numeric& dx,
1743 const Verbosity& verbosity) {
1744 for (const auto& qi : energy_level_identities)
1745 jacobianAddNLTE(ws,
1746 jq,
1747 jacobian_agenda,
1748 atmosphere_dim,
1749 p_grid,
1750 lat_grid,
1751 lon_grid,
1752 rq_p_grid,
1753 rq_lat_grid,
1754 rq_lon_grid,
1755 qi,
1756 dx,
1757 verbosity);
1758}
1759
1760/* Workspace method: Doxygen documentation will be auto-generated */
1763 Agenda& jacobian_agenda,
1764 const Index& atmosphere_dim,
1765 const Vector& p_grid,
1766 const Vector& lat_grid,
1767 const Vector& lon_grid,
1768 const Vector& rq_p_grid,
1769 const Vector& rq_lat_grid,
1770 const Vector& rq_lon_grid,
1771 const String& species,
1772 const Verbosity& verbosity) {
1775
1776 // Check retrieval grids, here we just check the length of the grids
1777 // vs. the atmosphere dimension
1778 ArrayOfVector grids(atmosphere_dim);
1779 {
1780 ostringstream os;
1781 if (!check_retrieval_grids(grids,
1782 os,
1783 p_grid,
1784 lat_grid,
1785 lon_grid,
1786 rq_p_grid,
1787 rq_lat_grid,
1788 rq_lon_grid,
1789 "retrieval pressure grid",
1790 "retrieval latitude grid",
1791 "retrievallongitude_grid",
1792 atmosphere_dim))
1793 throw runtime_error(os.str());
1794 }
1795
1796 // Create the new retrieval quantity
1798 rq.Grids(grids);
1799
1800 // Make sure modes are valid and complain if they are repeated
1801 if (species == "electrons") {
1802 for (Index it = 0; it < jq.nelem(); it++) {
1803 if (jq[it] == Jacobian::Atm::Electrons) {
1804 ostringstream os;
1805 os << "Electrons are already included in *jacobian_quantities*.";
1806 throw std::runtime_error(os.str());
1807 }
1808 }
1809 rq.Target(Jacobian::Target(Jacobian::Atm::Electrons));
1810
1811 } else if (species == "particulates") {
1812 for (Index it = 0; it < jq.nelem(); it++) {
1813 if (jq[it] == Jacobian::Atm::Particulates) {
1814 ostringstream os;
1815 os << "Particulates are already included in *jacobian_quantities*.";
1816 throw std::runtime_error(os.str());
1817 }
1818 }
1819 rq.Target(Jacobian::Target(Jacobian::Atm::Particulates));
1820
1821 } else {
1822 ostringstream os;
1823 os << "Unknown special species jacobian: \"" << species
1824 << "\"\nPlease see *jacobianAddSpecialSpecies* for viable options.";
1825 throw std::runtime_error(os.str());
1826 }
1827
1828 // Add it to the *jacobian_quantities*
1829 jq.push_back(rq);
1830
1831 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1832}
1833
1834//----------------------------------------------------------------------------
1835// Adjustments and transformations
1836//----------------------------------------------------------------------------
1837
1838/* Workspace method: Doxygen documentation will be auto-generated */
1840 Matrix& jacobian,
1841 const ArrayOfRetrievalQuantity& jacobian_quantities,
1842 const Vector& x,
1843 const Verbosity&) {
1844 // For flexibility inside inversion_iteration_agenda, we should accept empty
1845 // Jacobian
1846 if (jacobian.empty()) {
1847 return;
1848 }
1849
1850 // Adjustments
1851 //
1852 // Unfortunately the adjustment requires both range indices and the
1853 // untransformed x, which makes things a bit messy
1854 bool vars_init = false;
1856 Vector x0;
1857 //
1858 for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
1859 if (jacobian_quantities[q].Target().isSpeciesVMR() &&
1860 jacobian_quantities[q].Mode() == "rel") {
1861 if (!vars_init) {
1862 bool any_affine;
1863 jac_ranges_indices(jis0, any_affine, jacobian_quantities, true);
1864 x0 = x;
1865 transform_x_back(x0, jacobian_quantities);
1866 vars_init = true;
1867 }
1868 for (Index i = jis0[q][0]; i <= jis0[q][1]; i++) {
1869 if (x[i] != 1) {
1870 jacobian(joker, i) /= x[i];
1871 }
1872 }
1873 }
1874 }
1875
1876 // Transformations
1877 transform_jacobian(jacobian, x, jacobian_quantities);
1878}
1879
1880/* Workspace method: Doxygen documentation will be auto-generated */
1882 const Matrix& transformation_matrix,
1883 const Vector& offset_vector,
1884 const Verbosity& /*v*/
1885) {
1886 if (jqs.empty()) {
1887 runtime_error(
1888 "Jacobian quantities is empty, so there is nothing to add the "
1889 "transformation to.");
1890 }
1891
1892 Index nelem = jqs.back().Grids().nelem();
1893
1894 if (!(nelem == transformation_matrix.nrows())) {
1895 runtime_error(
1896 "Dimension of transformation matrix incompatible with retrieval grids.");
1897 }
1898 if (!(nelem == offset_vector.nelem())) {
1899 runtime_error(
1900 "Dimension of offset vector incompatible with retrieval grids.");
1901 }
1902
1903 jqs.back().SetTransformationMatrix(transpose(transformation_matrix));
1904 jqs.back().SetOffsetVector(offset_vector);
1905}
1906
1907/* Workspace method: Doxygen documentation will be auto-generated */
1909 const String& transformation_func,
1910 const Numeric& z_min,
1911 const Numeric& z_max,
1912 const Verbosity& /*v*/
1913) {
1914 if (jqs.empty())
1915 throw runtime_error(
1916 "Jacobian quantities is empty, so there is nothing to add the "
1917 "transformation to.");
1918
1919 if (transformation_func == "none") {
1920 jqs.back().SetTransformationFunc("");
1921 return;
1922 }
1923
1924 Vector pars;
1925
1926 if (transformation_func == "atanh") {
1927 if (z_max <= z_min)
1928 throw runtime_error(
1929 "For option atanh, the GIN *z_max* must be set and be > z_min.");
1930 pars.resize(2);
1931 pars[0] = z_min;
1932 pars[1] = z_max;
1933 } else if (transformation_func == "log" || transformation_func == "log10") {
1934 pars.resize(1);
1935 pars[0] = z_min;
1936 } else {
1937 ostringstream os;
1938 os << "Valid options for *transformation_func* are:\n"
1939 << "\"none\", \"log\", \"log10\" and \"atanh\"\n"
1940 << "But found: \"" << transformation_func << "\"";
1941 throw runtime_error(os.str());
1942 }
1943
1944 jqs.back().SetTransformationFunc(transformation_func);
1945 jqs.back().SetTFuncParameters(pars);
1946}
1947
1948//----------------------------------------------------------------------------
1949// Methods for doing perturbations
1950//----------------------------------------------------------------------------
1951
1952/* Workspace method: Doxygen documentation will be auto-generated */
1953void AtmFieldPerturb(Tensor3& perturbed_field,
1954 const Index& atmosphere_dim,
1955 const Vector& p_grid,
1956 const Vector& lat_grid,
1957 const Vector& lon_grid,
1958 const Tensor3& original_field,
1959 const Vector& p_ret_grid,
1960 const Vector& lat_ret_grid,
1961 const Vector& lon_ret_grid,
1962 const Index& pert_index,
1963 const Numeric& pert_size,
1964 const String& pert_mode,
1965 const Verbosity&) {
1966 // Input checks (more below)
1967 chk_atm_field("original_field",
1968 original_field,
1969 atmosphere_dim,
1970 p_grid,
1971 lat_grid,
1972 lon_grid,
1973 false );
1974
1975 // Pack retrieval grids into an ArrayOfVector
1976 ArrayOfVector ret_grids(atmosphere_dim);
1977 ret_grids[0] = p_ret_grid;
1978 if (atmosphere_dim>1){
1979 ret_grids[1] = lat_ret_grid;
1980 if (atmosphere_dim>2){
1981 ret_grids[2] = lon_ret_grid;
1982 }
1983 }
1984
1985 // Find mapping from retrieval grids to atmospheric grids
1986 ArrayOfGridPos gp_p, gp_lat, gp_lon;
1987 Index n_p, n_lat, n_lon;
1989 gp_lat,
1990 gp_lon,
1991 n_p,
1992 n_lat,
1993 n_lon,
1994 ret_grids,
1995 atmosphere_dim,
1996 p_grid,
1997 lat_grid,
1998 lon_grid);
1999
2000 // Now we can chec *pert_index*
2001 if (pert_index<0){
2002 throw runtime_error("Bad *pert_index*. It is negative.");
2003 }
2004 const Index n_tot = n_p * n_lat * n_lon;
2005 if (pert_index >= n_tot){
2006 throw runtime_error("Bad *pert_index*. It is too high with respect "
2007 "to length of retrieval grids.");
2008 }
2009
2010 // Create x-vector that matches perturbation
2011 Vector x(n_tot);
2012 if (pert_mode == "absolute" ){
2013 x = 0;
2014 x[pert_index] = pert_size;
2015 }
2016 else if (pert_mode == "relative" ){
2017 x = 1;
2018 x[pert_index] += pert_size;
2019 }
2020 else{
2021 throw runtime_error("Bad *pert_mode*. Allowed choices are: "
2022 """absolute"" and ""relative"".");
2023 }
2024
2025 // Map x to a perturbation defined at atmospheric grids
2026 Tensor3 x3d(n_p, n_lat, n_lon), pert(n_p, n_lat, n_lon);
2027 reshape(x3d, x);
2028 regrid_atmfield_by_gp_oem(pert, atmosphere_dim, x3d, gp_p, gp_lat, gp_lon);
2029
2030 // Init perturbed_field, if not equal to original_field
2031 if (&perturbed_field != &original_field) {
2032 perturbed_field = original_field;
2033 }
2034
2035 // Apply perturbation
2036 if (pert_mode == "absolute" ){
2037 perturbed_field += pert;
2038 }
2039 else{
2040 perturbed_field *= pert;
2041 }
2042}
2043
2044/* Workspace method: Doxygen documentation will be auto-generated */
2045void AtmFieldPerturbAtmGrids(Tensor3& perturbed_field,
2046 const Index& atmosphere_dim,
2047 const Vector& p_grid,
2048 const Vector& lat_grid,
2049 const Vector& lon_grid,
2050 const Tensor3& original_field,
2051 const Index& pert_index,
2052 const Numeric& pert_size,
2053 const String& pert_mode,
2054 const Verbosity&) {
2055 // Some sizes
2056 const Index n_p = p_grid.nelem();
2057 const Index n_lat = atmosphere_dim<2 ? 1 : lat_grid.nelem();
2058 const Index n_lon = atmosphere_dim<3 ? 1 : lon_grid.nelem();
2059
2060 // Check input
2061 chk_atm_field("original_field",
2062 original_field,
2063 atmosphere_dim,
2064 p_grid,
2065 lat_grid,
2066 lon_grid,
2067 false );
2068 if (pert_index<0){
2069 throw runtime_error("Bad *pert_index*. It is negative.");
2070 }
2071 if (pert_index >= n_p * n_lat * n_lon){
2072 throw runtime_error("Bad *pert_index*. It is too high with respect "
2073 "to length of atmospheric grids.");
2074 }
2075
2076 // Determine indexes with respect to atmospheric grids
2077 Index tot_index = pert_index;
2078 const Index lon_index = atmosphere_dim<3 ? 0 : tot_index / (n_lat * n_p);
2079 tot_index -= lon_index * n_lat * n_p;
2080 const Index lat_index = atmosphere_dim<2 ? 0 : tot_index / n_p;
2081 tot_index -= lat_index * n_p;
2082 const Index p_index = tot_index;
2083
2084 // Init perturbed_field, if not equal to original_field
2085 if (&perturbed_field != &original_field) {
2086 perturbed_field = original_field;
2087 }
2088
2089 // Perturb
2090 if (pert_mode == "absolute" ){
2091 perturbed_field(p_index,
2092 atmosphere_dim>1 ? lat_index : 0,
2093 atmosphere_dim>2 ? lon_index : 0) += pert_size;
2094 }
2095 else if (pert_mode == "relative"){
2096 perturbed_field(p_index,
2097 atmosphere_dim>1 ? lat_index : 0,
2098 atmosphere_dim>2 ? lon_index : 0) *= 1 + pert_size;
2099 }
2100 else{
2101 throw runtime_error("Bad *pert_mode*. Allowed choices are: "
2102 """absolute"" and ""relative"".");
2103 }
2104}
2105
2106/* Workspace method: Doxygen documentation will be auto-generated */
2108 const Index& atmosphere_dim,
2109 const Vector& p_grid,
2110 const Vector& lat_grid,
2111 const Vector& lon_grid,
2112 const Verbosity&) {
2113 const Index n_p = p_grid.nelem();
2114 const Index n_lat = atmosphere_dim<2 ? 1 : lat_grid.nelem();
2115 const Index n_lon = atmosphere_dim<3 ? 1 : lon_grid.nelem();
2116
2117 n = n_p * n_lat * n_lon;
2118}
2119
2120/* Workspace method: Doxygen documentation will be auto-generated */
2122 const Vector& y_pert,
2123 const Vector& y,
2124 const Numeric& pert_size,
2125 const Verbosity&) {
2126 const Index n = y.nelem();
2127 if( y_pert.nelem() != n ){
2128 throw runtime_error("Inconsistency in length of *y_pert* and *y*.");
2129 }
2130 jacobian = y_pert;
2131 jacobian -= y;
2132 jacobian /= pert_size;
2133}
2134
2135/* Workspace method: Doxygen documentation will be auto-generated */
2137 const ArrayOfVector& ybatch,
2138 const Vector& y,
2139 const Numeric& pert_size,
2140 const Verbosity&) {
2141 const Index n = y.nelem();
2142 const Index l = ybatch.nelem();
2143 if (l>0){
2144 if( ybatch[0].nelem() != n )
2145 throw runtime_error("Inconsistency in length of y and ybatch[0].");
2146 }
2147 jacobian.resize(n,l);
2148 for (Index i=0; i<l; i++) {
2149 jacobian(joker,i) = ybatch[i];
2150 jacobian(joker,i) -= y;
2151 }
2152 jacobian /= pert_size;
2153}
2154
2155/* Workspace method: Doxygen documentation will be auto-generated */
2156void particle_bulkprop_fieldPerturb(Tensor4& particle_bulkprop_field,
2157 const Index& atmosphere_dim,
2158 const Vector& p_grid,
2159 const Vector& lat_grid,
2160 const Vector& lon_grid,
2161 const ArrayOfString& particle_bulkprop_names,
2162 const String& particle_type,
2163 const Vector& p_ret_grid,
2164 const Vector& lat_ret_grid,
2165 const Vector& lon_ret_grid,
2166 const Index& pert_index,
2167 const Numeric& pert_size,
2168 const String& pert_mode,
2169 const Verbosity& verbosity) {
2170 // Locate particle_type among particle_bulkprop_names
2171 Index iq = find_first(particle_bulkprop_names, particle_type);
2172 if (iq < 0) {
2173 ostringstream os;
2174 os << "Could not find " << particle_type << " in *particle_bulkprop_names*.\n";
2175 throw std::runtime_error(os.str());
2176 }
2177
2178 Tensor3 original_field, perturbed_field;
2179 original_field = particle_bulkprop_field(iq,joker,joker,joker);
2180 AtmFieldPerturb(perturbed_field,
2181 atmosphere_dim,
2182 p_grid,
2183 lat_grid,
2184 lon_grid,
2185 original_field,
2186 p_ret_grid,
2187 lat_ret_grid,
2188 lon_ret_grid,
2189 pert_index,
2190 pert_size,
2191 pert_mode,
2192 verbosity);
2193 particle_bulkprop_field(iq,joker,joker,joker) = perturbed_field;
2194}
2195
2196/* Workspace method: Doxygen documentation will be auto-generated */
2198 const Index& atmosphere_dim,
2199 const Vector& p_grid,
2200 const Vector& lat_grid,
2201 const Vector& lon_grid,
2202 const ArrayOfString& particle_bulkprop_names,
2203 const String& particle_type,
2204 const Index& pert_index,
2205 const Numeric& pert_size,
2206 const String& pert_mode,
2207 const Verbosity& verbosity) {
2208 // Locate particle_type among particle_bulkprop_names
2209 Index iq = find_first(particle_bulkprop_names, particle_type);
2210 if (iq < 0) {
2211 ostringstream os;
2212 os << "Could not find " << particle_type << " in *particle_bulkprop_names*.\n";
2213 throw std::runtime_error(os.str());
2214 }
2215
2216 Tensor3 original_field, perturbed_field;
2217 original_field = particle_bulkprop_field(iq,joker,joker,joker);
2218 AtmFieldPerturbAtmGrids(perturbed_field,
2219 atmosphere_dim,
2220 p_grid,
2221 lat_grid,
2222 lon_grid,
2223 original_field,
2224 pert_index,
2225 pert_size,
2226 pert_mode,
2227 verbosity);
2228 particle_bulkprop_field(iq,joker,joker,joker) = perturbed_field;
2229}
2230
2231/* Workspace method: Doxygen documentation will be auto-generated */
2233 const Index& atmosphere_dim,
2234 const Vector& p_grid,
2235 const Vector& lat_grid,
2236 const Vector& lon_grid,
2237 const ArrayOfArrayOfSpeciesTag& abs_species,
2238 const String& species,
2239 const Vector& p_ret_grid,
2240 const Vector& lat_ret_grid,
2241 const Vector& lon_ret_grid,
2242 const Index& pert_index,
2243 const Numeric& pert_size,
2244 const String& pert_mode,
2245 const Verbosity& verbosity) {
2246 // Locate vmr_species among abs_species
2247 Index iq = -1;
2248 for (Index i = 0; i < abs_species.nelem(); i++) {
2249 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
2250 iq = i;
2251 break;
2252 }
2253 }
2254 if (iq < 0) {
2255 ostringstream os;
2256 os << "Could not find " << species << " in *abs_species*.\n";
2257 throw std::runtime_error(os.str());
2258 }
2259
2260 Tensor3 original_field, perturbed_field;
2261 original_field = vmr_field(iq,joker,joker,joker);
2262 AtmFieldPerturb(perturbed_field,
2263 atmosphere_dim,
2264 p_grid,
2265 lat_grid,
2266 lon_grid,
2267 original_field,
2268 p_ret_grid,
2269 lat_ret_grid,
2270 lon_ret_grid,
2271 pert_index,
2272 pert_size,
2273 pert_mode,
2274 verbosity);
2275 vmr_field(iq,joker,joker,joker) = perturbed_field;
2276}
2277
2278/* Workspace method: Doxygen documentation will be auto-generated */
2280 const Index& atmosphere_dim,
2281 const Vector& p_grid,
2282 const Vector& lat_grid,
2283 const Vector& lon_grid,
2284 const ArrayOfArrayOfSpeciesTag& abs_species,
2285 const String& species,
2286 const Index& pert_index,
2287 const Numeric& pert_size,
2288 const String& pert_mode,
2289 const Verbosity& verbosity) {
2290 // Locate vmr_species among abs_species
2291 Index iq = -1;
2292 for (Index i = 0; i < abs_species.nelem(); i++) {
2293 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
2294 iq = i;
2295 break;
2296 }
2297 }
2298 if (iq < 0) {
2299 ostringstream os;
2300 os << "Could not find " << species << " in *abs_species*.\n";
2301 throw std::runtime_error(os.str());
2302 }
2303
2304 Tensor3 original_field, perturbed_field;
2305 original_field = vmr_field(iq,joker,joker,joker);
2306 AtmFieldPerturbAtmGrids(perturbed_field,
2307 atmosphere_dim,
2308 p_grid,
2309 lat_grid,
2310 lon_grid,
2311 original_field,
2312 pert_index,
2313 pert_size,
2314 pert_mode,
2315 verbosity);
2316 vmr_field(iq,joker,joker,joker) = perturbed_field;
2317}
2318
2319
Declarations required for the calculation of absorption coefficients.
Array< ArrayOfIndex > ArrayOfArrayOfIndex
Definition: array.h:44
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:292
Array< Vector > ArrayOfVector
An array of vectors.
Definition: array.h:57
The global header file for ARTS.
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:44
void append(const String &methodname, const TokVal &keywordvalue)
Appends methods to an agenda.
Definition: agenda_class.cc:58
void check(Workspace &ws, const Verbosity &verbosity)
Checks consistency of an agenda.
Definition: agenda_class.cc:80
void set_name(const String &nname)
Set agenda name.
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:430
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
The range class.
Definition: matpackI.h:165
constexpr Index get_start() const ARTS_NOEXCEPT
Returns the start index of the range.
Definition: matpackI.h:332
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:345
const String & Mode() const
Returns the mode.
Definition: jacobian.h:414
Jacobian::Target & Target()
Get the Jacobian Target.
Definition: jacobian.h:451
const String & SubSubtag() const
Returns the sub-sub-tag.
Definition: jacobian.h:399
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:428
const String & Subtag() const
Returns the sub-tag.
Definition: jacobian.h:385
The Sparse class.
Definition: matpackII.h:67
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
The Tensor3 class.
Definition: matpackIII.h:339
The Tensor4 class.
Definition: matpackIV.h:421
This stores arbitrary token values and remembers the type.
Definition: token.h:42
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
Internal cloudbox functions.
#define _U_
Definition: config.h:180
Constants of physical expressions as constexpr.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define DEBUG_ONLY(...)
Definition: debug.h:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#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:847
void transform_jacobian(Matrix &jacobian, const Vector x, const ArrayOfRetrievalQuantity &jqs)
Applies both functional and affine transformations.
Definition: jacobian.cc:90
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition: jacobian.cc:232
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:640
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:45
Routines for setting up the jacobian.
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:560
#define q
#define dx
#define max(a, b)
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:215
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
void jacobianAdjustAndTransform(Matrix &jacobian, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Verbosity &)
WORKSPACE METHOD: jacobianAdjustAndTransform.
Definition: m_jacobian.cc:1839
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:253
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 Agenda &geo_pos_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianCalcPointingZaRecalc.
Definition: m_jacobian.cc:711
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:1761
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:2232
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:2121
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:1953
void jacobianFromYbatch(Matrix &jacobian, const ArrayOfVector &ybatch, const Vector &y, const Numeric &pert_size, const Verbosity &)
WORKSPACE METHOD: jacobianFromYbatch.
Definition: m_jacobian.cc:2136
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 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:343
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:510
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:2156
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:1731
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:2279
void jacobianSetAffineTransformation(ArrayOfRetrievalQuantity &jqs, const Matrix &transformation_matrix, const Vector &offset_vector, const Verbosity &)
WORKSPACE METHOD: jacobianSetAffineTransformation.
Definition: m_jacobian.cc:1881
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:1670
void jacobianClose(Workspace &ws, Index &jacobian_do, Agenda &jacobian_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianClose.
Definition: m_jacobian.cc:67
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:2107
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:2197
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:1606
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:104
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:1572
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 jacobianOff(Index &jacobian_do, Agenda &jacobian_agenda, ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianOff.
Definition: m_jacobian.cc:91
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:582
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:198
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:2045
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:1908
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:1647
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 jacobianInit(ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Verbosity &)
WORKSPACE METHOD: jacobianInit.
Definition: m_jacobian.cc:82
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:394
Workspace methods and template functions for supergeneric XML IO.
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:735
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:207
#define CREATE_OUT2
Definition: messages.h:206
Array< String > ArrayOfString
An array of Strings.
Definition: mystring.h:290
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
Particulates Polyfit
Definition: jacobian.h:86
Particulates FrequencyStretch
Definition: jacobian.h:85
Temperature
Definition: jacobian.h:57
Particulates FrequencyShift
Definition: jacobian.h:84
Particulates Sinefit
Definition: jacobian.h:87
Particulates PointingZenithInterp
Definition: jacobian.h:88
MagneticMagnitude
Definition: jacobian.h:59
WindMagnitude
Definition: jacobian.h:58
X3 LineStrength
Definition: constants.h:563
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
This file contains declerations of functions of physical character.
#define u
#define v
#define w
#define a
#define c
Quantum::Identifier QuantumIdentifier
Definition: quantum.h:471
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, Matrix &geo_pos_matrix, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const ConstVectorView &f_grid, const ConstMatrixView &sensor_pos, const ConstMatrixView &sensor_los, const ConstMatrixView &transmitter_pos, const ConstMatrixView &mblock_dlos_grid, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
Performs calculations for one measurement block, on iy-level.
Definition: rte.cc:1773
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
Returns the "range" of y corresponding to a measurement block.
Definition: rte.cc:1110
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:99