ARTS 2.5.4 (git: 7d04b88e)
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_in, verbosity);
78 jacobian_do = 1;
79}
80
81/* Workspace method: Doxygen documentation will be auto-generated */
83 ArrayOfRetrievalQuantity& jacobian_quantities,
84 Agenda& jacobian_agenda,
85 const Verbosity&) {
86 jacobian_quantities.resize(0);
87 jacobian_agenda = Agenda{ws};
88 jacobian_agenda.set_name("jacobian_agenda");
89}
90
91/* Workspace method: Doxygen documentation will be auto-generated */
93 Index& jacobian_do,
94 Agenda& jacobian_agenda,
95 ArrayOfRetrievalQuantity& jacobian_quantities,
96 const Verbosity& verbosity) {
97 jacobian_do = 0;
98 jacobianInit(ws, jacobian_quantities, jacobian_agenda, verbosity);
99}
100
101//----------------------------------------------------------------------------
102// Absorption species:
103//----------------------------------------------------------------------------
104
105/* Workspace method: Doxygen documentation will be auto-generated */
108 Agenda& jacobian_agenda,
109 const Index& atmosphere_dim,
110 const Vector& p_grid,
111 const Vector& lat_grid,
112 const Vector& lon_grid,
113 const Vector& rq_p_grid,
114 const Vector& rq_lat_grid,
115 const Vector& rq_lon_grid,
116 const String& species,
117 const String& mode,
118 const Index& for_species_tag,
119 const Verbosity& verbosity) {
122
124 if (not for_species_tag) {
125 ArrayOfSpeciesTag test(species);
126 if (test.nelem() not_eq 1)
127 throw std::runtime_error(
128 "Trying to add a species as a species tag of multiple species.\n"
129 "This is not supported. Please give just a single species instead.\n"
130 "Otherwise consider if you intended for_species_tag to be evaluated true.\n");
131 qi = QuantumIdentifier(test[0].Isotopologue());
132 }
133
134 // Check that this species is not already included in the jacobian.
135 for (Index it = 0; it < jq.nelem(); it++) {
136 ARTS_USER_ERROR_IF (jq[it] == Jacobian::Special::ArrayOfSpeciesTagVMR && jq[it].Subtag() == species,
137 "The gas species:\n", species, "\nis already included in *jacobian_quantities*.")
138 ARTS_USER_ERROR_IF (jq[it] == Jacobian::Line::VMR and jq[it].QuantumIdentity() == qi,
139 "The atmospheric species of:\n", species, "\nis already included in *jacobian_quantities*\n"
140 "as: ", jq[it])
141 }
142
143 // Check retrieval grids, here we just check the length of the grids
144 // vs. the atmosphere dimension
145 ArrayOfVector grids(atmosphere_dim);
146 {
147 ostringstream os;
148 if (!check_retrieval_grids(grids,
149 os,
150 p_grid,
151 lat_grid,
152 lon_grid,
153 rq_p_grid,
154 rq_lat_grid,
155 rq_lon_grid,
156 "retrieval pressure grid",
157 "retrieval latitude grid",
158 "retrievallongitude_grid",
159 atmosphere_dim))
160 throw runtime_error(os.str());
161 }
162
163 // Check that mode is correct
164 if (mode != "vmr" && mode != "nd" && mode != "rel" && mode != "rh" &&
165 mode != "q") {
166 throw runtime_error(
167 "The retrieval mode can only be \"vmr\", \"nd\", "
168 "\"rel\", \"rh\" or \"q\".");
169 }
170 if ((mode == "rh" || mode == "q") && species.substr(0, 3) != "H2O") {
171 throw runtime_error(
172 "Retrieval modes \"rh\" and \"q\" can only be applied "
173 "on species starting with H2O.");
174 }
175
176 // Create the new retrieval quantity
178 rq.Subtag(species);
179 rq.Mode(mode);
180 rq.Grids(grids);
181 if (for_species_tag == 0) {
182 rq.Target(Jacobian::Target(Jacobian::Line::VMR, qi, qi.Species()));
183 } else {
184 ArrayOfSpeciesTag test(species);
186 }
187 rq.Target().perturbation = 0.001;
188
189 // Add it to the *jacobian_quantities*
190 jq.push_back(rq);
191
192 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
193}
194
195//----------------------------------------------------------------------------
196// Frequency shift
197//----------------------------------------------------------------------------
198
199/* Workspace method: Doxygen documentation will be auto-generated */
201 ArrayOfRetrievalQuantity& jacobian_quantities,
202 Agenda& jacobian_agenda,
203 const Vector& f_grid,
204 const Numeric& df,
205 const Verbosity&) {
206 // Check that this jacobian type is not already included.
207 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
208 if (jacobian_quantities[it] == Jacobian::Sensor::FrequencyShift) {
209 ostringstream os;
210 os << "Fit of frequency shift is already included in\n"
211 << "*jacobian_quantities*.";
212 throw runtime_error(os.str());
213 }
214 }
215
216 // Checks of frequencies
217 if (df <= 0) throw runtime_error("The argument *df* must be > 0.");
218 if (df > 1e6)
219 throw runtime_error("The argument *df* is not allowed to exceed 1 MHz.");
220 const Index nf = f_grid.nelem();
221 if (nf < 2)
222 throw runtime_error(
223 "Frequency shifts and *f_grid* of length 1 can "
224 "not be combined.");
225 const Numeric maxdf = f_grid[nf - 1] - f_grid[nf - 2];
226 if (df > maxdf) {
227 ostringstream os;
228 os << "The value of *df* is too big with respect to spacing of "
229 << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
230 << "between the two last elements of *f_grid*.\n"
231 << "This spacing is : " << maxdf / 1e3 << " kHz\n"
232 << "The value of df is: " << df / 1e3 << " kHz";
233 throw runtime_error(os.str());
234 }
235
236 // Create the new retrieval quantity
238 rq.Target() = Jacobian::Target(Jacobian::Sensor::FrequencyShift);
239 rq.Mode("");
240 rq.Target().perturbation = df;
241
242 // Dummy vector of length 1
243 Vector grid(1, 0);
244 ArrayOfVector grids(1, grid);
245 rq.Grids(grids);
246
247 // Add it to the *jacobian_quantities*
248 jacobian_quantities.push_back(rq);
249
250 // Add corresponding calculation method to the jacobian agenda
251 jacobian_agenda.append("jacobianCalcFreqShift", "");
252}
253
254/* Workspace method: Doxygen documentation will be auto-generated */
256 const Index& mblock_index,
257 const Vector& iyb,
258 const Vector& yb,
259 const Index& stokes_dim,
260 const Vector& f_grid,
261 const Matrix& mblock_dlos_grid,
262 const Sparse& sensor_response,
263 const ArrayOfRetrievalQuantity& jacobian_quantities,
264 const Verbosity&) {
265 // Set some useful (and needed) variables.
267 ArrayOfIndex ji;
268
269 // Find the retrieval quantity related to this method.
270 // This works since the combined MainTag and Subtag is individual.
271 bool found = false;
272 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
273 if (jacobian_quantities[n] == Jacobian::Sensor::FrequencyShift) {
274 bool any_affine;
275 ArrayOfArrayOfIndex jacobian_indices;
277 jacobian_indices, any_affine, jacobian_quantities, true);
278 //
279 found = true;
280 rq = jacobian_quantities[n];
281 ji = jacobian_indices[n];
282 }
283 }
284 if (!found) {
285 throw runtime_error(
286 "There is no such frequency retrieval quantity defined.\n");
287 }
288
289 // Check that sensor_response is consistent with yb and iyb
290 //
291 if (sensor_response.nrows() != yb.nelem())
292 throw runtime_error("Mismatch in size between *sensor_response* and *yb*.");
293 if (sensor_response.ncols() != iyb.nelem())
294 throw runtime_error(
295 "Mismatch in size between *sensor_response* and *iyb*.");
296
297 // Get disturbed (part of) y
298 //
299 const Index n1y = sensor_response.nrows();
300 Vector dy(n1y);
301 {
302 const Index nf2 = f_grid.nelem();
303 const Index nlos2 = mblock_dlos_grid.nrows();
304 const Index niyb = nf2 * nlos2 * stokes_dim;
305
306 // Interpolation weights
307 //
308 constexpr Index porder = 3;
309 //
310 Vector fg_new = f_grid, iyb2(niyb);
311 fg_new += rq.Target().perturbation;
312 //
313 const auto lag = Interpolation::FixedLagrangeVector<porder>(fg_new, f_grid);
314 const auto itw = interpweights(lag);
315
316 // Do interpolation
317 for (Index ilos = 0; ilos < nlos2; ilos++) {
318 const Index row0 = ilos * nf2 * stokes_dim;
319
320 for (Index is = 0; is < stokes_dim; is++) {
321 iyb2[Range(row0 + is, nf2, stokes_dim)] = reinterp(iyb[Range(row0 + is, nf2, stokes_dim)], itw, lag);
322 }
323 }
324
325 // Determine difference
326 //
327 mult(dy, sensor_response, iyb2);
328 //
329 for (Index i = 0; i < n1y; i++) {
330 dy[i] = (dy[i] - yb[i]) / rq.Target().perturbation;
331 }
332 }
333
334 //--- Set jacobian ---
335 ARTS_ASSERT(rq.Grids()[0].nelem() == 1);
336 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
337 jacobian(rowind, ji[0]) = dy;
338}
339
340//----------------------------------------------------------------------------
341// Frequency stretch
342//----------------------------------------------------------------------------
343
344/* Workspace method: Doxygen documentation will be auto-generated */
346 ArrayOfRetrievalQuantity& jacobian_quantities,
347 Agenda& jacobian_agenda,
348 const Vector& f_grid,
349 const Numeric& df,
350 const Verbosity&) {
351 // Check that this jacobian type is not already included.
352 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
353 if (jacobian_quantities[it] == Jacobian::Sensor::FrequencyStretch) {
354 ostringstream os;
355 os << "Fit of frequency stretch is already included in\n"
356 << "*jacobian_quantities*.";
357 throw runtime_error(os.str());
358 }
359 }
360
361 // Checks of df
362 if (df <= 0) throw runtime_error("The argument *df* must be > 0.");
363 if (df > 1e6)
364 throw runtime_error("The argument *df* is not allowed to exceed 1 MHz.");
365 const Index nf = f_grid.nelem();
366 const Numeric maxdf = f_grid[nf - 1] - f_grid[nf - 2];
367 if (df > maxdf) {
368 ostringstream os;
369 os << "The value of *df* is too big with respect to spacing of "
370 << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
371 << "between the two last elements of *f_grid*.\n"
372 << "This spacing is : " << maxdf / 1e3 << " kHz\n"
373 << "The value of df is: " << df / 1e3 << " kHz";
374 throw runtime_error(os.str());
375 }
376
377 // Create the new retrieval quantity
379 rq.Target() = Jacobian::Target(Jacobian::Sensor::FrequencyStretch);
380 rq.Mode("");
381 rq.Target().perturbation=df;
382
383 // Dummy vector of length 1
384 Vector grid(1, 0);
385 ArrayOfVector grids(1, grid);
386 rq.Grids(grids);
387
388 // Add it to the *jacobian_quantities*
389 jacobian_quantities.push_back(rq);
390
391 // Add corresponding calculation method to the jacobian agenda
392 jacobian_agenda.append("jacobianCalcFreqStretch", "");
393}
394
395/* Workspace method: Doxygen documentation will be auto-generated */
397 Matrix& jacobian,
398 const Index& mblock_index,
399 const Vector& iyb,
400 const Vector& yb,
401 const Index& stokes_dim,
402 const Vector& f_grid,
403 const Matrix& mblock_dlos_grid,
404 const Sparse& sensor_response,
405 const ArrayOfIndex& sensor_response_pol_grid,
406 const Vector& sensor_response_f_grid,
407 const Matrix& sensor_response_dlos_grid,
408 const ArrayOfRetrievalQuantity& jacobian_quantities,
409 const Verbosity&) {
410 // The code here is close to identical to the one for Shift. The main
411 // difference is that dy is weighted with poly_order 1 basis function.
412
413 // Set some useful (and needed) variables.
415 ArrayOfIndex ji;
416
417 // Find the retrieval quantity related to this method.
418 // This works since the combined MainTag and Subtag is individual.
419 bool found = false;
420 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
421 if (jacobian_quantities[n] == Jacobian::Sensor::FrequencyStretch) {
422 bool any_affine;
423 ArrayOfArrayOfIndex jacobian_indices;
425 jacobian_indices, any_affine, jacobian_quantities, true);
426 //
427 found = true;
428 rq = jacobian_quantities[n];
429 ji = jacobian_indices[n];
430 }
431 }
432 if (!found) {
433 throw runtime_error(
434 "There is no such frequency retrieval quantity defined.\n");
435 }
436
437 // Check that sensor_response is consistent with yb and iyb
438 //
439 if (sensor_response.nrows() != yb.nelem())
440 throw runtime_error("Mismatch in size between *sensor_response* and *yb*.");
441 if (sensor_response.ncols() != iyb.nelem())
442 throw runtime_error(
443 "Mismatch in size between *sensor_response* and *iyb*.");
444
445 // Get disturbed (part of) y
446 //
447 const Index n1y = sensor_response.nrows();
448 Vector dy(n1y);
449 {
450 const Index nf2 = f_grid.nelem();
451 const Index nlos2 = mblock_dlos_grid.nrows();
452 const Index niyb = nf2 * nlos2 * stokes_dim;
453
454 // Interpolation weights
455 //
456 constexpr Index porder = 3;
457 //
458 Vector fg_new = f_grid, iyb2(niyb);
459 //
460 fg_new += rq.Target().perturbation;
461 //
462 const auto lag = Interpolation::FixedLagrangeVector<porder>(fg_new, f_grid);
463 const auto itw = interpweights(lag);
464
465 // Do interpolation
466 for (Index ilos = 0; ilos < nlos2; ilos++) {
467 const Index row0 = ilos * nf2 * stokes_dim;
468
469 for (Index is = 0; is < stokes_dim; is++) {
470 iyb2[Range(row0 + is, nf2, stokes_dim)] = reinterp(iyb[Range(row0 + is, nf2, stokes_dim)], itw, lag);
471 }
472 }
473
474 // Determine difference
475 //
476 mult(dy, sensor_response, iyb2);
477 //
478 for (Index i = 0; i < n1y; i++) {
479 dy[i] = (dy[i] - yb[i]) / rq.Target().perturbation;
480 }
481
482 // dy above corresponds now to shift. Convert to stretch:
483 //
484 Vector w;
485 polynomial_basis_func(w, sensor_response_f_grid, 1);
486 //
487 const Index nf = sensor_response_f_grid.nelem();
488 const Index npol = sensor_response_pol_grid.nelem();
489 const Index nlos = sensor_response_dlos_grid.nrows();
490 //
491 for (Index l = 0; l < nlos; l++) {
492 for (Index f = 0; f < nf; f++) {
493 const Index row1 = (l * nf + f) * npol;
494 for (Index p = 0; p < npol; p++) {
495 dy[row1 + p] *= w[f];
496 }
497 }
498 }
499 }
500
501 //--- Set jacobians ---
502 ARTS_ASSERT(rq.Grids()[0].nelem() == 1);
503 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
504 jacobian(rowind, ji[0]) = dy;
505}
506
507//----------------------------------------------------------------------------
508// Pointing:
509//----------------------------------------------------------------------------
510
511/* Workspace method: Doxygen documentation will be auto-generated */
513 ArrayOfRetrievalQuantity& jacobian_quantities,
514 Agenda& jacobian_agenda,
515 const Matrix& sensor_pos,
516 const ArrayOfTime& sensor_time,
517 const Index& poly_order,
518 const String& calcmode,
519 const Numeric& dza,
520 const Verbosity&) {
521 // Check that poly_order is -1 or positive
522 if (poly_order < -1)
523 throw runtime_error(
524 "The polynomial order has to be positive or -1 for gitter.");
525
526 // Check that this jacobian type is not already included.
527 for (Index it = 0; it < jacobian_quantities.nelem(); it++) {
528 if (jacobian_quantities[it].Target().isPointing()) {
529 ostringstream os;
530 os << "Fit of zenith angle pointing off-set is already included in\n"
531 << "*jacobian_quantities*.";
532 throw runtime_error(os.str());
533 }
534 }
535
536 // Checks of dza
537 if (dza <= 0) throw runtime_error("The argument *dza* must be > 0.");
538 if (dza > 0.1)
539 throw runtime_error("The argument *dza* is not allowed to exceed 0.1 deg.");
540
541 // Check that sensor_time is consistent with sensor_pos
542 if (sensor_time.nelem() != sensor_pos.nrows()) {
543 ostringstream os;
544 os << "The WSV *sensor_time* must be defined for every "
545 << "measurement block.\n";
546 throw runtime_error(os.str());
547 }
548
549 // Do not allow that *poly_order* is not too large compared to *sensor_time*
550 if (poly_order > sensor_time.nelem() - 1) {
551 throw runtime_error(
552 "The polynomial order can not be >= length of *sensor_time*.");
553 }
554
555 // Create the new retrieval quantity
557 if (calcmode == "recalc") {
558 rq.Target() = Jacobian::Target(Jacobian::Sensor::PointingZenithRecalc);
559 jacobian_agenda.append("jacobianCalcPointingZaRecalc", "");
560 } else if (calcmode == "interp") {
561 rq.Target() = Jacobian::Target(Jacobian::Sensor::PointingZenithInterp);
562 jacobian_agenda.append("jacobianCalcPointingZaInterp", "");
563 } else
564 throw runtime_error(
565 R"(Possible choices for *calcmode* are "recalc" and "interp".)");
566 rq.Target().perturbation = dza;
567
568 // To store the value or the polynomial order, create a vector with length
569 // poly_order+1, in case of gitter set the size of the grid vector to be the
570 // number of measurement blocks, all elements set to -1.
571 Vector grid(0, poly_order + 1, 1);
572 if (poly_order == -1) {
573 grid.resize(sensor_pos.nrows());
574 grid = -1.0;
575 }
576 ArrayOfVector grids(1, grid);
577 rq.Grids(grids);
578
579 // Add it to the *jacobian_quantities*
580 jacobian_quantities.push_back(rq);
581}
582
583/* Workspace method: Doxygen documentation will be auto-generated */
585 Matrix& jacobian,
586 const Index& mblock_index,
587 const Vector& iyb,
588 const Vector& yb _U_,
589 const Index& stokes_dim,
590 const Vector& f_grid,
591 const Matrix& DEBUG_ONLY(sensor_los),
592 const Matrix& mblock_dlos_grid,
593 const Sparse& sensor_response,
594 const ArrayOfTime& sensor_time,
595 const ArrayOfRetrievalQuantity& jacobian_quantities,
596 const Verbosity&) {
597 if (mblock_dlos_grid.nrows() < 2)
598 throw runtime_error(
599 "The method demands that *mblock_dlos_grid* has "
600 "more than one row.");
601
602 if (!(is_increasing(mblock_dlos_grid(joker, 0)) ||
603 is_decreasing(mblock_dlos_grid(joker, 0))))
604 throw runtime_error(
605 "The method demands that the zenith angles in "
606 "*mblock_dlos_grid* are sorted (increasing or decreasing).");
607
608 // Set some useful variables.
610 ArrayOfIndex ji;
611
612 // Find the retrieval quantity related to this method.
613 // This works since the combined MainTag and Subtag is individual.
614 bool found = false;
615 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
616 if (jacobian_quantities[n] == Jacobian::Sensor::PointingZenithInterp) {
617 bool any_affine;
618 ArrayOfArrayOfIndex jacobian_indices;
620 jacobian_indices, any_affine, jacobian_quantities, true);
621 //
622 found = true;
623 rq = jacobian_quantities[n];
624 ji = jacobian_indices[n];
625 }
626 }
627 if (!found) {
628 throw runtime_error(
629 "There is no such pointing retrieval quantity defined.\n");
630 }
631
632 // Get "dy", by inter/extra-polation of existing iyb
633 //
634 const Index n1y = sensor_response.nrows();
635 Vector dy(n1y);
636 {
637 // Sizes
638 const Index nf = f_grid.nelem();
639 const Index nza = mblock_dlos_grid.nrows();
640
641 // Shifted zenith angles
642 Vector za1 = mblock_dlos_grid(joker, 0);
643 za1 -= rq.Target().perturbation;
644 Vector za2 = mblock_dlos_grid(joker, 0);
645 za2 += rq.Target().perturbation;
646
647 // Find interpolation weights
648 ArrayOfGridPos gp1(nza), gp2(nza);
649 gridpos(
650 gp1, mblock_dlos_grid(joker, 0), za1, 1e6); // Note huge extrapolation!
651 gridpos(
652 gp2, mblock_dlos_grid(joker, 0), za2, 1e6); // Note huge extrapolation!
653 Matrix itw1(nza, 2), itw2(nza, 2);
654 interpweights(itw1, gp1);
655 interpweights(itw2, gp2);
656
657 // Make interpolation (for all azimuth angles, frequencies and Stokes)
658 //
659 Vector iyb1(iyb.nelem()), iyb2(iyb.nelem());
660 //
661 for (Index iza = 0; iza < nza; iza++) {
662 for (Index iv = 0; iv < nf; iv++) {
663 for (Index is = 0; is < stokes_dim; is++) {
664 const Range r(iv * stokes_dim + is, nza, nf * stokes_dim);
665 interp(iyb1[r], itw1, iyb[r], gp1);
666 interp(iyb2[r], itw2, iyb[r], gp2);
667 }
668 }
669 }
670
671 // Apply sensor and take difference
672 //
673 Vector y1(n1y), y2(n1y);
674 mult(y1, sensor_response, iyb1);
675 mult(y2, sensor_response, iyb2);
676 //
677 for (Index i = 0; i < n1y; i++) {
678 dy[i] = (y2[i] - y1[i]) / (2 * rq.Target().perturbation);
679 }
680 }
681
682 //--- Create jacobians ---
683
684 const Index lg = rq.Grids()[0].nelem();
685 const Index it = ji[0];
686 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
687 const Index row0 = rowind.get_start();
688
689 // Handle pointing "jitter" seperately
690 if (rq.Grids()[0][0] == -1) // Not all values are set here,
691 { // but should already have been
692 ARTS_ASSERT(lg == sensor_los.nrows()); // set to 0
693 ARTS_ASSERT(rq.Grids()[0][mblock_index] == -1);
694 jacobian(rowind, it + mblock_index) = dy;
695 }
696
697 // Polynomial representation
698 else {
699 Vector w;
700 for (Index c = 0; c < lg; c++) {
701 ARTS_ASSERT(Numeric(c) == rq.Grids()[0][c]);
702 //
703 polynomial_basis_func(w, time_vector(sensor_time), c);
704 //
705 for (Index i = 0; i < n1y; i++) {
706 jacobian(row0 + i, it + c) = w[mblock_index] * dy[i];
707 }
708 }
709 }
710}
711
712/* Workspace method: Doxygen documentation will be auto-generated */
714 Workspace& ws,
715 Matrix& jacobian,
716 const Index& mblock_index,
717 const Vector& iyb _U_,
718 const Vector& yb,
719 const Index& atmosphere_dim,
720 const EnergyLevelMap& nlte_field,
721 const Index& cloudbox_on,
722 const Index& stokes_dim,
723 const Vector& f_grid,
724 const Matrix& sensor_pos,
725 const Matrix& sensor_los,
726 const Matrix& transmitter_pos,
727 const Matrix& mblock_dlos_grid,
728 const Sparse& sensor_response,
729 const ArrayOfTime& sensor_time,
730 const String& iy_unit,
731 const Agenda& iy_main_agenda,
732 const ArrayOfRetrievalQuantity& jacobian_quantities,
733 const Verbosity& verbosity) {
734 // Set some useful variables.
736 ArrayOfIndex ji;
737
738 // Find the retrieval quantity related to this method.
739 // This works since the combined MainTag and Subtag is individual.
740 bool found = false;
741 for (Index n = 0; n < jacobian_quantities.nelem() && !found; n++) {
742 if (jacobian_quantities[n] == Jacobian::Sensor::PointingZenithRecalc) {
743 bool any_affine;
744 ArrayOfArrayOfIndex jacobian_indices;
746 jacobian_indices, any_affine, jacobian_quantities, true);
747 //
748 found = true;
749 rq = jacobian_quantities[n];
750 ji = jacobian_indices[n];
751 }
752 }
753 if (!found) {
754 throw runtime_error(
755 "There is no such pointing retrieval quantity defined.\n");
756 }
757
758 // Get "dy", by calling iyb_calc with shifted sensor_los.
759 //
760 const Index n1y = sensor_response.nrows();
761 Vector dy(n1y);
762 {
763 Vector iyb2;
764 Matrix los = sensor_los;
765 Matrix geo_pos;
766 ArrayOfVector iyb_aux;
767 ArrayOfMatrix diyb_dx;
768
769 los(joker, 0) += rq.Target().perturbation;
770
771 iyb_calc(ws,
772 iyb2,
773 iyb_aux,
774 diyb_dx,
775 geo_pos,
776 mblock_index,
777 atmosphere_dim,
778 nlte_field,
779 cloudbox_on,
780 stokes_dim,
781 f_grid,
782 sensor_pos,
783 los,
784 transmitter_pos,
785 mblock_dlos_grid,
786 iy_unit,
787 iy_main_agenda,
788 0,
792 verbosity);
793
794 // Apply sensor and take difference
795 //
796 mult(dy, sensor_response, iyb2);
797 //
798 for (Index i = 0; i < n1y; i++) {
799 dy[i] = (dy[i] - yb[i]) / rq.Target().perturbation;
800 }
801 }
802
803 //--- Create jacobians ---
804
805 const Index lg = rq.Grids()[0].nelem();
806 const Index it = ji[0];
807 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
808 const Index row0 = rowind.get_start();
809
810 // Handle "jitter" seperately
811 if (rq.Grids()[0][0] == -1) // Not all values are set here,
812 { // but should already have been
813 ARTS_ASSERT(lg == sensor_los.nrows()); // set to 0
814 ARTS_ASSERT(rq.Grids()[0][mblock_index] == -1);
815 jacobian(rowind, it + mblock_index) = dy;
816 }
817
818 // Polynomial representation
819 else {
820 Vector w;
821 for (Index c = 0; c < lg; c++) {
822 ARTS_ASSERT(Numeric(c) == rq.Grids()[0][c]);
823 //
824 polynomial_basis_func(w, time_vector(sensor_time), c);
825 //
826 for (Index i = 0; i < n1y; i++) {
827 jacobian(row0 + i, it + c) = w[mblock_index] * dy[i];
828 }
829 }
830 }
831}
832
833//----------------------------------------------------------------------------
834// Polynomial baseline fits:
835//----------------------------------------------------------------------------
836
837/* Workspace method: Doxygen documentation will be auto-generated */
840 Agenda& jacobian_agenda,
841 const ArrayOfIndex& sensor_response_pol_grid,
842 const Matrix& sensor_response_dlos_grid,
843 const Matrix& sensor_pos,
844 const Index& poly_order,
845 const Index& no_pol_variation,
846 const Index& no_los_variation,
847 const Index& no_mblock_variation,
848 const Verbosity&) {
849 // Check that poly_order is >= 0
850 if (poly_order < 0)
851 throw runtime_error("The polynomial order has to be >= 0.");
852
853 // Check that polyfit is not already included in the jacobian.
854 for (Index it = 0; it < jq.nelem(); it++) {
855 if (jq[it] == Jacobian::Sensor::Polyfit) {
856 ostringstream os;
857 os << "Polynomial baseline fit is already included in\n"
858 << "*jacobian_quantities*.";
859 throw runtime_error(os.str());
860 }
861 }
862
863 // "Grids"
864 //
865 // Grid dimensions correspond here to
866 // 1: polynomial order
867 // 2: polarisation
868 // 3: viewing direction
869 // 4: measurement block
870 //
871 ArrayOfVector grids(4);
872 //
873 if (no_pol_variation)
874 grids[1] = Vector(1, 1);
875 else
876 grids[1] = Vector(0, sensor_response_pol_grid.nelem(), 1);
877 if (no_los_variation)
878 grids[2] = Vector(1, 1);
879 else
880 grids[2] = Vector(0, sensor_response_dlos_grid.nrows(), 1);
881 if (no_mblock_variation)
882 grids[3] = Vector(1, 1);
883 else
884 grids[3] = Vector(0, sensor_pos.nrows(), 1);
885
886 // Create the new retrieval quantity
888 rq.Target() = Jacobian::Target(Jacobian::Sensor::Polyfit);
889 rq.Mode("");
890 rq.Target().perturbation = 0;
891
892 // Each polynomial coeff. is treated as a retrieval quantity
893 //
894 for (Index i = 0; i <= poly_order; i++) {
895 ostringstream sstr;
896 sstr << "Coefficient " << i;
897 rq.Subtag(sstr.str());
898
899 // Grid is a scalar, use polynomial coeff.
900 grids[0] = Vector(1, (Numeric)i);
901 rq.Grids(grids);
902
903 // Add it to the *jacobian_quantities*
904 jq.push_back(rq);
905
906 // Add pointing method to the jacobian agenda
907 jacobian_agenda.append("jacobianCalcPolyfit", i);
908 }
909}
910
911/* Workspace method: Doxygen documentation will be auto-generated */
913 const Index& mblock_index,
914 const Vector& iyb _U_,
915 const Vector& yb _U_,
916 const Sparse& sensor_response,
917 const ArrayOfIndex& sensor_response_pol_grid,
918 const Vector& sensor_response_f_grid,
919 const Matrix& sensor_response_dlos_grid,
920 const ArrayOfRetrievalQuantity& jacobian_quantities,
921 const Index& poly_coeff,
922 const Verbosity&) {
923 // Find the retrieval quantity related to this method
925 ArrayOfIndex ji;
926 bool found = false;
927 Index iq;
928 ostringstream sstr;
929 sstr << "Coefficient " << poly_coeff;
930 for (iq = 0; iq < jacobian_quantities.nelem() && !found; iq++) {
931 if (jacobian_quantities[iq] == Jacobian::Sensor::Polyfit &&
932 jacobian_quantities[iq].Subtag() == sstr.str()) {
933 found = true;
934 break;
935 }
936 }
937 if (!found) {
938 throw runtime_error(
939 "There is no Polyfit jacobian defined, in general "
940 "or for the selected polynomial coefficient.\n");
941 }
942
943 // Size and check of sensor_response
944 //
945 const Index nf = sensor_response_f_grid.nelem();
946 const Index npol = sensor_response_pol_grid.nelem();
947 const Index nlos = sensor_response_dlos_grid.nrows();
948
949 // Make a vector with values to distribute over *jacobian*
950 //
951 Vector w;
952 //
953 polynomial_basis_func(w, sensor_response_f_grid, poly_coeff);
954
955 // Fill J
956 //
957 ArrayOfArrayOfIndex jacobian_indices;
958 {
959 bool any_affine;
960 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
961 }
962 //
963 ArrayOfVector jg = jacobian_quantities[iq].Grids();
964 const Index n1 = jg[1].nelem();
965 const Index n2 = jg[2].nelem();
966 const Index n3 = jg[3].nelem();
967 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
968 const Index row4 = rowind.get_start();
969 Index col4 = jacobian_indices[iq][0];
970
971 if (n3 > 1) {
972 col4 += mblock_index * n2 * n1;
973 }
974
975 for (Index l = 0; l < nlos; l++) {
976 const Index row3 = row4 + l * nf * npol;
977 const Index col3 = col4 + l * n1;
978
979 for (Index f = 0; f < nf; f++) {
980 const Index row2 = row3 + f * npol;
981
982 for (Index p = 0; p < npol; p++) {
983 Index col1 = col3;
984 if (n1 > 1) {
985 col1 += p;
986 }
987
988 jacobian(row2 + p, col1) = w[f];
989 }
990 }
991 }
992}
993
994//----------------------------------------------------------------------------
995// Scattering species:
996//----------------------------------------------------------------------------
997
998/* Workspace method: Doxygen documentation will be auto-generated */
1001 Agenda& jacobian_agenda,
1002 const Index& atmosphere_dim,
1003 const Vector& p_grid,
1004 const Vector& lat_grid,
1005 const Vector& lon_grid,
1006 const Vector& rq_p_grid,
1007 const Vector& rq_lat_grid,
1008 const Vector& rq_lon_grid,
1009 const String& species,
1010 const String& quantity,
1011 const Verbosity& verbosity) {
1014
1015 // Check that this species+quantity combination is not already included in
1016 // the jacobian.
1017 for (Index it = 0; it < jq.nelem(); it++) {
1018 if (jq[it] == Jacobian::Special::ScatteringString && jq[it].Subtag() == species &&
1019 jq[it].SubSubtag() == quantity) {
1020 ostringstream os;
1021 os << "The combintaion of\n scattering species: " << species
1022 << "\n retrieval quantity: " << quantity
1023 << "\nis already included in *jacobian_quantities*.";
1024 throw runtime_error(os.str());
1025 }
1026 }
1027
1028 // Check retrieval grids, here we just check the length of the grids
1029 // vs. the atmosphere dimension
1030 ArrayOfVector grids(atmosphere_dim);
1031 {
1032 ostringstream os;
1033 if (!check_retrieval_grids(grids,
1034 os,
1035 p_grid,
1036 lat_grid,
1037 lon_grid,
1038 rq_p_grid,
1039 rq_lat_grid,
1040 rq_lon_grid,
1041 "retrieval pressure grid",
1042 "retrieval latitude grid",
1043 "retrievallongitude_grid",
1044 atmosphere_dim))
1045 throw runtime_error(os.str());
1046 }
1047
1048 // Create the new retrieval quantity
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);
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) {
1391 break;
1394 break;
1397 break;
1398 case Options::WindMagJacobian::strength:
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) {
1470 break;
1473 break;
1476 break;
1477 case Options::WindMagJacobian::strength:
1479 break;
1480 case Options::WindMagJacobian::FINAL:
1481 ARTS_ASSERT(false, "This error should be caught earlier")
1482 }
1483
1484 // Check that this species is not already included in the jacobian.
1485 for (Index it = 0; it < jq.nelem(); it++) {
1486 ARTS_USER_ERROR_IF (jq[it].Target().sameTargetType(rq.Target()),
1487 "The magnetic component:\n", component, "\n"
1488 "is already included in *jacobian_quantities*.")
1489 }
1490
1491 // Check retrieval grids, here we just check the length of the grids
1492 // vs. the atmosphere dimension
1493 ArrayOfVector grids(atmosphere_dim);
1494 {
1495 ostringstream os;
1496 if (!check_retrieval_grids(grids,
1497 os,
1498 p_grid,
1499 lat_grid,
1500 lon_grid,
1501 rq_p_grid,
1502 rq_lat_grid,
1503 rq_lon_grid,
1504 "retrieval pressure grid",
1505 "retrieval latitude grid",
1506 "retrievallongitude_grid",
1507 atmosphere_dim))
1508 throw runtime_error(os.str());
1509 }
1510
1511 rq.Grids(grids);
1512 rq.Target().perturbation = dB;
1513
1514 // Add it to the *jacobian_quantities*
1515 jq.push_back(rq);
1516
1517 // Add gas species method to the jacobian agenda
1518 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1519}
1520
1521//----------------------------------------------------------------------------
1522// Catalog parameters:
1523//----------------------------------------------------------------------------
1524
1525/* Workspace method: Doxygen documentation will be auto-generated */
1528 Agenda& jacobian_agenda,
1529 const QuantumIdentifier& line_identity,
1530 const String& species,
1531 const String& variable,
1532 const String& coefficient,
1533 const Verbosity& verbosity) {
1535
1536 const auto jpt = select_derivativeLineShape(variable, coefficient);
1537
1538 out3 << "Attempting to create RT tag for " << line_identity << " " << variable
1539 << " " << coefficient << " for ";
1540
1541 // Create the quantity
1543 rq.Mode(species);
1544 rq.Grids(ArrayOfVector(0, Vector(0)));
1545
1546 // Map the species
1547 if (species == LineShape::self_broadening) {
1548 rq.Target(Jacobian::Target(jpt, line_identity, line_identity.Species()));
1549 } else if (species == LineShape::bath_broadening) {
1550 rq.Target(Jacobian::Target(jpt, line_identity, Species::Species::Bath));
1551 } else {
1552 rq.Target(Jacobian::Target(jpt, line_identity, SpeciesTag(species).Spec()));
1553 }
1554 out3 << species << ' ' << rq.Target().species_id << "\n";
1555
1556 // Test this is not a copy
1557 for (auto& q : jq)
1558 if (q.HasSameInternalsAs(rq))
1559 throw std::runtime_error("Error with copies of the quantities");
1560
1561 // Append and do housekeeping
1562 jq.push_back(rq);
1563 out3 << "Creation was successful!\n";
1564 jacobian_agenda.append("jacobianCalcDoNothing",
1565 TokVal()); // old code activation
1566}
1567
1568/* Workspace method: Doxygen documentation will be auto-generated */
1570 Workspace& ws,
1572 Agenda& jacobian_agenda,
1573 const ArrayOfQuantumIdentifier& line_identities,
1574 const ArrayOfString& species,
1575 const ArrayOfString& variables,
1576 const ArrayOfString& coefficients,
1577 const Verbosity& verbosity) {
1578 if (not(line_identities.nelem() or species.nelem() or variables.nelem() or
1579 coefficients.nelem()))
1580 throw std::runtime_error("Must have at least 1-long lists for all GINs");
1581
1582 ArrayOfString vars;
1583 if (variables[0] == "ALL")
1584 vars = ArrayOfString(LineShape::enumstrs::VariableNames);
1585 else
1586 vars = variables;
1587
1588 ArrayOfString coeffs;
1589 if (coefficients[0] == "ALL")
1590 coeffs = ArrayOfString(Options::enumstrs::LineShapeCoeffNames);
1591 else
1592 coeffs = coefficients;
1593
1594 for (auto& l : line_identities)
1595 for (auto& s : species)
1596 for (auto& v : vars)
1597 for (auto& c : coeffs)
1599 ws, jq, jacobian_agenda, l, s, v, c, verbosity);
1600}
1601
1602/* Workspace method: Doxygen documentation will be auto-generated */
1605 Agenda& jacobian_agenda,
1606 const QuantumIdentifier& catalog_identity,
1607 const String& catalog_parameter,
1608 const Verbosity& verbosity) {
1610
1611 // Create the new retrieval quantity
1612 const auto opt = Options::toBasicCatParamJacobianOrThrow(catalog_parameter);
1614 switch (opt) {
1615 case Options::BasicCatParamJacobian::LineCenter:
1616 rq.Target(Jacobian::Target(Jacobian::Line::Center, catalog_identity, catalog_identity.Species()));
1617 break;
1618 case Options::BasicCatParamJacobian::LineStrength:
1619 rq.Target(Jacobian::Target(Jacobian::Line::Strength, catalog_identity, catalog_identity.Species()));
1620 break;
1621 case Options::BasicCatParamJacobian::FINAL:
1622 ARTS_ASSERT(false, "This error should be caught earlier")
1623 }
1624
1625 // Check that this is not already included in the jacobian.
1626 for (Index it = 0; it < jq.nelem(); it++) {
1627 ARTS_USER_ERROR_IF (rq.Target().sameTargetType(jq[it].Target()),
1628 "The catalog identifier:\n",
1629 catalog_identity, " for ID: ", catalog_identity, "\n"
1630 "is already included in jacobian_quantities*.")
1631 }
1632
1633 rq.Grids(ArrayOfVector(0, Vector(0)));
1634
1635 // Add it to the *jacobian_quantities*
1636 jq.push_back(rq);
1637
1638 out3 << " Calculations done by propagation matrix expressions.\n";
1639
1640 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1641}
1642
1643/* Workspace method: Doxygen documentation will be auto-generated */
1645 Workspace& ws,
1647 Agenda& jacobian_agenda,
1648 const ArrayOfQuantumIdentifier& catalog_identities,
1649 const ArrayOfString& catalog_parameters,
1650 const Verbosity& verbosity) {
1652
1653 out2 << " Adding " << catalog_identities.nelem() * catalog_parameters.nelem()
1654 << " expression(s) to the Jacobian calculations.\n";
1655
1656 for (auto& qi : catalog_identities)
1657 for (auto& param : catalog_parameters)
1659 ws, jq, jacobian_agenda, qi, param, verbosity);
1660}
1661
1662//----------------------------------------------------------------------------
1663// NLTE temperature:
1664//----------------------------------------------------------------------------
1665
1666/* Workspace method: Doxygen documentation will be auto-generated */
1669 Agenda& jacobian_agenda,
1670 const Index& atmosphere_dim,
1671 const Vector& p_grid,
1672 const Vector& lat_grid,
1673 const Vector& lon_grid,
1674 const Vector& rq_p_grid,
1675 const Vector& rq_lat_grid,
1676 const Vector& rq_lon_grid,
1677 const QuantumIdentifier& energy_level_identity,
1678 const Numeric& dx,
1679 const Verbosity& verbosity) {
1681
1682 // Check that this species is not already included in the jacobian.
1683 for (Index it = 0; it < jq.nelem(); it++) {
1684 if (jq[it] == Jacobian::Line::NLTE and
1685 jq[it].QuantumIdentity() == energy_level_identity) {
1686 ostringstream os;
1687 os << "The NLTE identifier:\n"
1688 << energy_level_identity << "\nis already included in "
1689 << "*jacobian_quantities*.";
1690 throw std::runtime_error(os.str());
1691 }
1692 }
1693
1694 // Check retrieval grids, here we just check the length of the grids
1695 // vs. the atmosphere dimension
1696 ArrayOfVector grids(atmosphere_dim);
1697 {
1698 ostringstream os;
1699 if (not check_retrieval_grids(grids,
1700 os,
1701 p_grid,
1702 lat_grid,
1703 lon_grid,
1704 rq_p_grid,
1705 rq_lat_grid,
1706 rq_lon_grid,
1707 "retrieval pressure grid",
1708 "retrieval latitude grid",
1709 "retrievallongitude_grid",
1710 atmosphere_dim))
1711 throw runtime_error(os.str());
1712 }
1713
1714 // Create the new retrieval quantity
1716 rq.Target(Jacobian::Target(Jacobian::Line::NLTE, energy_level_identity, energy_level_identity.Species()));
1717 rq.Target().perturbation = dx;
1718 rq.Grids(grids);
1719
1720 // Add it to the *jacobian_quantities*
1721 jq.push_back(rq);
1722
1723 out3 << " Calculations done by propagation matrix expressions.\n";
1724
1725 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1726}
1727
1730 Agenda& jacobian_agenda,
1731 const Index& atmosphere_dim,
1732 const Vector& p_grid,
1733 const Vector& lat_grid,
1734 const Vector& lon_grid,
1735 const Vector& rq_p_grid,
1736 const Vector& rq_lat_grid,
1737 const Vector& rq_lon_grid,
1738 const ArrayOfQuantumIdentifier& energy_level_identities,
1739 const Numeric& dx,
1740 const Verbosity& verbosity) {
1741 for (const auto& qi : energy_level_identities)
1742 jacobianAddNLTE(ws,
1743 jq,
1744 jacobian_agenda,
1745 atmosphere_dim,
1746 p_grid,
1747 lat_grid,
1748 lon_grid,
1749 rq_p_grid,
1750 rq_lat_grid,
1751 rq_lon_grid,
1752 qi,
1753 dx,
1754 verbosity);
1755}
1756
1757/* Workspace method: Doxygen documentation will be auto-generated */
1760 Agenda& jacobian_agenda,
1761 const Index& atmosphere_dim,
1762 const Vector& p_grid,
1763 const Vector& lat_grid,
1764 const Vector& lon_grid,
1765 const Vector& rq_p_grid,
1766 const Vector& rq_lat_grid,
1767 const Vector& rq_lon_grid,
1768 const String& species,
1769 const Verbosity& verbosity) {
1772
1773 // Check retrieval grids, here we just check the length of the grids
1774 // vs. the atmosphere dimension
1775 ArrayOfVector grids(atmosphere_dim);
1776 {
1777 ostringstream os;
1778 if (!check_retrieval_grids(grids,
1779 os,
1780 p_grid,
1781 lat_grid,
1782 lon_grid,
1783 rq_p_grid,
1784 rq_lat_grid,
1785 rq_lon_grid,
1786 "retrieval pressure grid",
1787 "retrieval latitude grid",
1788 "retrievallongitude_grid",
1789 atmosphere_dim))
1790 throw runtime_error(os.str());
1791 }
1792
1793 // Create the new retrieval quantity
1795 rq.Grids(grids);
1796
1797 // Make sure modes are valid and complain if they are repeated
1798 if (species == "electrons") {
1799 for (Index it = 0; it < jq.nelem(); it++) {
1800 if (jq[it] == Jacobian::Atm::Electrons) {
1801 ostringstream os;
1802 os << "Electrons are already included in *jacobian_quantities*.";
1803 throw std::runtime_error(os.str());
1804 }
1805 }
1807
1808 } else if (species == "particulates") {
1809 for (Index it = 0; it < jq.nelem(); it++) {
1810 if (jq[it] == Jacobian::Atm::Particulates) {
1811 ostringstream os;
1812 os << "Particulates are already included in *jacobian_quantities*.";
1813 throw std::runtime_error(os.str());
1814 }
1815 }
1816 rq.Target(Jacobian::Target(Jacobian::Atm::Particulates));
1817
1818 } else {
1819 ostringstream os;
1820 os << "Unknown special species jacobian: \"" << species
1821 << "\"\nPlease see *jacobianAddSpecialSpecies* for viable options.";
1822 throw std::runtime_error(os.str());
1823 }
1824
1825 // Add it to the *jacobian_quantities*
1826 jq.push_back(rq);
1827
1828 jacobian_agenda.append("jacobianCalcDoNothing", TokVal());
1829}
1830
1831//----------------------------------------------------------------------------
1832// Adjustments and transformations
1833//----------------------------------------------------------------------------
1834
1835/* Workspace method: Doxygen documentation will be auto-generated */
1837 Matrix& jacobian,
1838 const ArrayOfRetrievalQuantity& jacobian_quantities,
1839 const Vector& x,
1840 const Verbosity&) {
1841 // For flexibility inside inversion_iteration_agenda, we should accept empty
1842 // Jacobian
1843 if (jacobian.empty()) {
1844 return;
1845 }
1846
1847 // Adjustments
1848 //
1849 // Unfortunately the adjustment requires both range indices and the
1850 // untransformed x, which makes things a bit messy
1851 bool vars_init = false;
1853 Vector x0;
1854 //
1855 for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
1856 if (jacobian_quantities[q].Target().isSpeciesVMR() &&
1857 jacobian_quantities[q].Mode() == "rel") {
1858 if (!vars_init) {
1859 bool any_affine;
1860 jac_ranges_indices(jis0, any_affine, jacobian_quantities, true);
1861 x0 = x;
1862 transform_x_back(x0, jacobian_quantities);
1863 vars_init = true;
1864 }
1865 for (Index i = jis0[q][0]; i <= jis0[q][1]; i++) {
1866 if (x[i] != 1) {
1867 jacobian(joker, i) /= x[i];
1868 }
1869 }
1870 }
1871 }
1872
1873 // Transformations
1874 transform_jacobian(jacobian, x, jacobian_quantities);
1875}
1876
1877/* Workspace method: Doxygen documentation will be auto-generated */
1879 const Matrix& transformation_matrix,
1880 const Vector& offset_vector,
1881 const Verbosity& /*v*/
1882) {
1883 if (jqs.empty()) {
1884 runtime_error(
1885 "Jacobian quantities is empty, so there is nothing to add the "
1886 "transformation to.");
1887 }
1888
1889 Index nelem = jqs.back().Grids().nelem();
1890
1891 if (!(nelem == transformation_matrix.nrows())) {
1892 runtime_error(
1893 "Dimension of transformation matrix incompatible with retrieval grids.");
1894 }
1895 if (!(nelem == offset_vector.nelem())) {
1896 runtime_error(
1897 "Dimension of offset vector incompatible with retrieval grids.");
1898 }
1899
1900 jqs.back().SetTransformationMatrix(transpose(transformation_matrix));
1901 jqs.back().SetOffsetVector(offset_vector);
1902}
1903
1904/* Workspace method: Doxygen documentation will be auto-generated */
1906 const String& transformation_func,
1907 const Numeric& z_min,
1908 const Numeric& z_max,
1909 const Verbosity& /*v*/
1910) {
1911 if (jqs.empty())
1912 throw runtime_error(
1913 "Jacobian quantities is empty, so there is nothing to add the "
1914 "transformation to.");
1915
1916 if (transformation_func == "none") {
1917 jqs.back().SetTransformationFunc("");
1918 return;
1919 }
1920
1921 Vector pars;
1922
1923 if (transformation_func == "atanh") {
1924 if (z_max <= z_min)
1925 throw runtime_error(
1926 "For option atanh, the GIN *z_max* must be set and be > z_min.");
1927 pars.resize(2);
1928 pars[0] = z_min;
1929 pars[1] = z_max;
1930 } else if (transformation_func == "log" || transformation_func == "log10") {
1931 pars.resize(1);
1932 pars[0] = z_min;
1933 } else {
1934 ostringstream os;
1935 os << "Valid options for *transformation_func* are:\n"
1936 << "\"none\", \"log\", \"log10\" and \"atanh\"\n"
1937 << "But found: \"" << transformation_func << "\"";
1938 throw runtime_error(os.str());
1939 }
1940
1941 jqs.back().SetTransformationFunc(transformation_func);
1942 jqs.back().SetTFuncParameters(pars);
1943}
1944
1945//----------------------------------------------------------------------------
1946// Methods for doing perturbations
1947//----------------------------------------------------------------------------
1948
1949/* Workspace method: Doxygen documentation will be auto-generated */
1950void AtmFieldPerturb(Tensor3& perturbed_field,
1951 const Index& atmosphere_dim,
1952 const Vector& p_grid,
1953 const Vector& lat_grid,
1954 const Vector& lon_grid,
1955 const Tensor3& original_field,
1956 const Vector& p_ret_grid,
1957 const Vector& lat_ret_grid,
1958 const Vector& lon_ret_grid,
1959 const Index& pert_index,
1960 const Numeric& pert_size,
1961 const String& pert_mode,
1962 const Verbosity&) {
1963 // Input checks (more below)
1964 chk_atm_field("original_field",
1965 original_field,
1966 atmosphere_dim,
1967 p_grid,
1968 lat_grid,
1969 lon_grid,
1970 false );
1971
1972 // Pack retrieval grids into an ArrayOfVector
1973 ArrayOfVector ret_grids(atmosphere_dim);
1974 ret_grids[0] = p_ret_grid;
1975 if (atmosphere_dim>1){
1976 ret_grids[1] = lat_ret_grid;
1977 if (atmosphere_dim>2){
1978 ret_grids[2] = lon_ret_grid;
1979 }
1980 }
1981
1982 // Find mapping from retrieval grids to atmospheric grids
1983 ArrayOfGridPos gp_p, gp_lat, gp_lon;
1984 Index n_p, n_lat, n_lon;
1986 gp_lat,
1987 gp_lon,
1988 n_p,
1989 n_lat,
1990 n_lon,
1991 ret_grids,
1992 atmosphere_dim,
1993 p_grid,
1994 lat_grid,
1995 lon_grid);
1996
1997 // Now we can chec *pert_index*
1998 if (pert_index<0){
1999 throw runtime_error("Bad *pert_index*. It is negative.");
2000 }
2001 const Index n_tot = n_p * n_lat * n_lon;
2002 if (pert_index >= n_tot){
2003 throw runtime_error("Bad *pert_index*. It is too high with respect "
2004 "to length of retrieval grids.");
2005 }
2006
2007 // Create x-vector that matches perturbation
2008 Vector x(n_tot);
2009 if (pert_mode == "absolute" ){
2010 x = 0;
2011 x[pert_index] = pert_size;
2012 }
2013 else if (pert_mode == "relative" ){
2014 x = 1;
2015 x[pert_index] += pert_size;
2016 }
2017 else{
2018 throw runtime_error("Bad *pert_mode*. Allowed choices are: "
2019 """absolute"" and ""relative"".");
2020 }
2021
2022 // Map x to a perturbation defined at atmospheric grids
2023 Tensor3 x3d(n_p, n_lat, n_lon), pert(n_p, n_lat, n_lon);
2024 reshape(x3d, x);
2025 regrid_atmfield_by_gp_oem(pert, atmosphere_dim, x3d, gp_p, gp_lat, gp_lon);
2026
2027 // Init perturbed_field, if not equal to original_field
2028 if (&perturbed_field != &original_field) {
2029 perturbed_field = original_field;
2030 }
2031
2032 // Apply perturbation
2033 if (pert_mode == "absolute" ){
2034 perturbed_field += pert;
2035 }
2036 else{
2037 perturbed_field *= pert;
2038 }
2039}
2040
2041/* Workspace method: Doxygen documentation will be auto-generated */
2042void AtmFieldPerturbAtmGrids(Tensor3& perturbed_field,
2043 const Index& atmosphere_dim,
2044 const Vector& p_grid,
2045 const Vector& lat_grid,
2046 const Vector& lon_grid,
2047 const Tensor3& original_field,
2048 const Index& pert_index,
2049 const Numeric& pert_size,
2050 const String& pert_mode,
2051 const Verbosity&) {
2052 // Some sizes
2053 const Index n_p = p_grid.nelem();
2054 const Index n_lat = atmosphere_dim<2 ? 1 : lat_grid.nelem();
2055 const Index n_lon = atmosphere_dim<3 ? 1 : lon_grid.nelem();
2056
2057 // Check input
2058 chk_atm_field("original_field",
2059 original_field,
2060 atmosphere_dim,
2061 p_grid,
2062 lat_grid,
2063 lon_grid,
2064 false );
2065 if (pert_index<0){
2066 throw runtime_error("Bad *pert_index*. It is negative.");
2067 }
2068 if (pert_index >= n_p * n_lat * n_lon){
2069 throw runtime_error("Bad *pert_index*. It is too high with respect "
2070 "to length of atmospheric grids.");
2071 }
2072
2073 // Determine indexes with respect to atmospheric grids
2074 Index tot_index = pert_index;
2075 const Index lon_index = atmosphere_dim<3 ? 0 : tot_index / (n_lat * n_p);
2076 tot_index -= lon_index * n_lat * n_p;
2077 const Index lat_index = atmosphere_dim<2 ? 0 : tot_index / n_p;
2078 tot_index -= lat_index * n_p;
2079 const Index p_index = tot_index;
2080
2081 // Init perturbed_field, if not equal to original_field
2082 if (&perturbed_field != &original_field) {
2083 perturbed_field = original_field;
2084 }
2085
2086 // Perturb
2087 if (pert_mode == "absolute" ){
2088 perturbed_field(p_index,
2089 atmosphere_dim>1 ? lat_index : 0,
2090 atmosphere_dim>2 ? lon_index : 0) += pert_size;
2091 }
2092 else if (pert_mode == "relative"){
2093 perturbed_field(p_index,
2094 atmosphere_dim>1 ? lat_index : 0,
2095 atmosphere_dim>2 ? lon_index : 0) *= 1 + pert_size;
2096 }
2097 else{
2098 throw runtime_error("Bad *pert_mode*. Allowed choices are: "
2099 """absolute"" and ""relative"".");
2100 }
2101}
2102
2103/* Workspace method: Doxygen documentation will be auto-generated */
2105 const Index& atmosphere_dim,
2106 const Vector& p_grid,
2107 const Vector& lat_grid,
2108 const Vector& lon_grid,
2109 const Verbosity&) {
2110 const Index n_p = p_grid.nelem();
2111 const Index n_lat = atmosphere_dim<2 ? 1 : lat_grid.nelem();
2112 const Index n_lon = atmosphere_dim<3 ? 1 : lon_grid.nelem();
2113
2114 n = n_p * n_lat * n_lon;
2115}
2116
2117/* Workspace method: Doxygen documentation will be auto-generated */
2119 const Vector& y_pert,
2120 const Vector& y,
2121 const Numeric& pert_size,
2122 const Verbosity&) {
2123 const Index n = y.nelem();
2124 if( y_pert.nelem() != n ){
2125 throw runtime_error("Inconsistency in length of *y_pert* and *y*.");
2126 }
2127 jacobian = y_pert;
2128 jacobian -= y;
2129 jacobian /= pert_size;
2130}
2131
2132/* Workspace method: Doxygen documentation will be auto-generated */
2134 const ArrayOfVector& ybatch,
2135 const Vector& y,
2136 const Numeric& pert_size,
2137 const Verbosity&) {
2138 const Index n = y.nelem();
2139 const Index l = ybatch.nelem();
2140 if (l>0){
2141 if( ybatch[0].nelem() != n )
2142 throw runtime_error("Inconsistency in length of y and ybatch[0].");
2143 }
2144 jacobian.resize(n,l);
2145 for (Index i=0; i<l; i++) {
2146 jacobian(joker,i) = ybatch[i];
2147 jacobian(joker,i) -= y;
2148 }
2149 jacobian /= pert_size;
2150}
2151
2152/* Workspace method: Doxygen documentation will be auto-generated */
2153void particle_bulkprop_fieldPerturb(Tensor4& particle_bulkprop_field,
2154 const Index& atmosphere_dim,
2155 const Vector& p_grid,
2156 const Vector& lat_grid,
2157 const Vector& lon_grid,
2158 const ArrayOfString& particle_bulkprop_names,
2159 const String& particle_type,
2160 const Vector& p_ret_grid,
2161 const Vector& lat_ret_grid,
2162 const Vector& lon_ret_grid,
2163 const Index& pert_index,
2164 const Numeric& pert_size,
2165 const String& pert_mode,
2166 const Verbosity& verbosity) {
2167 // Locate particle_type among particle_bulkprop_names
2168 Index iq = find_first(particle_bulkprop_names, particle_type);
2169 if (iq < 0) {
2170 ostringstream os;
2171 os << "Could not find " << particle_type << " in *particle_bulkprop_names*.\n";
2172 throw std::runtime_error(os.str());
2173 }
2174
2175 Tensor3 original_field, perturbed_field;
2176 original_field = particle_bulkprop_field(iq,joker,joker,joker);
2177 AtmFieldPerturb(perturbed_field,
2178 atmosphere_dim,
2179 p_grid,
2180 lat_grid,
2181 lon_grid,
2182 original_field,
2183 p_ret_grid,
2184 lat_ret_grid,
2185 lon_ret_grid,
2186 pert_index,
2187 pert_size,
2188 pert_mode,
2189 verbosity);
2190 particle_bulkprop_field(iq,joker,joker,joker) = perturbed_field;
2191}
2192
2193/* Workspace method: Doxygen documentation will be auto-generated */
2195 const Index& atmosphere_dim,
2196 const Vector& p_grid,
2197 const Vector& lat_grid,
2198 const Vector& lon_grid,
2199 const ArrayOfString& particle_bulkprop_names,
2200 const String& particle_type,
2201 const Index& pert_index,
2202 const Numeric& pert_size,
2203 const String& pert_mode,
2204 const Verbosity& verbosity) {
2205 // Locate particle_type among particle_bulkprop_names
2206 Index iq = find_first(particle_bulkprop_names, particle_type);
2207 if (iq < 0) {
2208 ostringstream os;
2209 os << "Could not find " << particle_type << " in *particle_bulkprop_names*.\n";
2210 throw std::runtime_error(os.str());
2211 }
2212
2213 Tensor3 original_field, perturbed_field;
2214 original_field = particle_bulkprop_field(iq,joker,joker,joker);
2215 AtmFieldPerturbAtmGrids(perturbed_field,
2216 atmosphere_dim,
2217 p_grid,
2218 lat_grid,
2219 lon_grid,
2220 original_field,
2221 pert_index,
2222 pert_size,
2223 pert_mode,
2224 verbosity);
2225 particle_bulkprop_field(iq,joker,joker,joker) = perturbed_field;
2226}
2227
2228/* Workspace method: Doxygen documentation will be auto-generated */
2230 const Index& atmosphere_dim,
2231 const Vector& p_grid,
2232 const Vector& lat_grid,
2233 const Vector& lon_grid,
2234 const ArrayOfArrayOfSpeciesTag& abs_species,
2235 const String& species,
2236 const Vector& p_ret_grid,
2237 const Vector& lat_ret_grid,
2238 const Vector& lon_ret_grid,
2239 const Index& pert_index,
2240 const Numeric& pert_size,
2241 const String& pert_mode,
2242 const Verbosity& verbosity) {
2243 // Locate vmr_species among abs_species
2244 Index iq = -1;
2245 for (Index i = 0; i < abs_species.nelem(); i++) {
2246 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
2247 iq = i;
2248 break;
2249 }
2250 }
2251 if (iq < 0) {
2252 ostringstream os;
2253 os << "Could not find " << species << " in *abs_species*.\n";
2254 throw std::runtime_error(os.str());
2255 }
2256
2257 Tensor3 original_field, perturbed_field;
2258 original_field = vmr_field(iq,joker,joker,joker);
2259 AtmFieldPerturb(perturbed_field,
2260 atmosphere_dim,
2261 p_grid,
2262 lat_grid,
2263 lon_grid,
2264 original_field,
2265 p_ret_grid,
2266 lat_ret_grid,
2267 lon_ret_grid,
2268 pert_index,
2269 pert_size,
2270 pert_mode,
2271 verbosity);
2272 vmr_field(iq,joker,joker,joker) = perturbed_field;
2273}
2274
2275/* Workspace method: Doxygen documentation will be auto-generated */
2277 const Index& atmosphere_dim,
2278 const Vector& p_grid,
2279 const Vector& lat_grid,
2280 const Vector& lon_grid,
2281 const ArrayOfArrayOfSpeciesTag& abs_species,
2282 const String& species,
2283 const Index& pert_index,
2284 const Numeric& pert_size,
2285 const String& pert_mode,
2286 const Verbosity& verbosity) {
2287 // Locate vmr_species among abs_species
2288 Index iq = -1;
2289 for (Index i = 0; i < abs_species.nelem(); i++) {
2290 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
2291 iq = i;
2292 break;
2293 }
2294 }
2295 if (iq < 0) {
2296 ostringstream os;
2297 os << "Could not find " << species << " in *abs_species*.\n";
2298 throw std::runtime_error(os.str());
2299 }
2300
2301 Tensor3 original_field, perturbed_field;
2302 original_field = vmr_field(iq,joker,joker,joker);
2303 AtmFieldPerturbAtmGrids(perturbed_field,
2304 atmosphere_dim,
2305 p_grid,
2306 lat_grid,
2307 lon_grid,
2308 original_field,
2309 pert_index,
2310 pert_size,
2311 pert_mode,
2312 verbosity);
2313 vmr_field(iq,joker,joker,joker) = perturbed_field;
2314}
2315
2316
Declarations required for the calculation of absorption coefficients.
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:190
Array< ArrayOfIndex > ArrayOfArrayOfIndex
Definition: array.h:138
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:69
void append(const String &methodname, const TokVal &keywordvalue)
Appends methods to an agenda.
Definition: agenda_class.cc:87
void check(Workspace &ws_in, const Verbosity &verbosity)
Checks consistency of an agenda.
void set_name(const String &nname)
Set agenda name.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
bool empty() const noexcept
Definition: matpackI.h:1080
Index nrows() const noexcept
Definition: matpackI.h:1075
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:554
The Matrix class.
Definition: matpackI.h:1283
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1063
The range class.
Definition: matpackI.h:177
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:358
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:324
const String & Mode() const
Returns the mode.
Definition: jacobian.h:393
Jacobian::Target & Target()
Get the Jacobian Target.
Definition: jacobian.h:430
const String & SubSubtag() const
Returns the sub-sub-tag.
Definition: jacobian.h:378
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:407
const String & Subtag() const
Returns the sub-tag.
Definition: jacobian.h:364
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4 class.
Definition: matpackIV.h:429
Definition: tokval.h:33
The Vector class.
Definition: matpackI.h:922
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
Workspace class.
Definition: workspace_ng.h:53
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:538
#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:218
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:255
void jacobianAdjustAndTransform(Matrix &jacobian, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Verbosity &)
WORKSPACE METHOD: jacobianAdjustAndTransform.
Definition: m_jacobian.cc:1836
void jacobianCalcFreqShift(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: jacobianCalcFreqShift.
Definition: m_jacobian.cc:255
void jacobianAddWind(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dfrequency, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddWind.
Definition: m_jacobian.cc:1369
void jacobianCalcPointingZaRecalc(Workspace &ws, Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const ArrayOfTime &sensor_time, const String &iy_unit, const Agenda &iy_main_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianCalcPointingZaRecalc.
Definition: m_jacobian.cc:713
void jacobianAddSpecialSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddSpecialSpecies.
Definition: m_jacobian.cc:1758
void jacobianClose(Workspace &ws_in, Index &jacobian_do, Agenda &jacobian_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianClose.
Definition: m_jacobian.cc:67
void vmr_fieldPerturb(Tensor4 &vmr_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Vector &p_ret_grid, const Vector &lat_ret_grid, const Vector &lon_ret_grid, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: vmr_fieldPerturb.
Definition: m_jacobian.cc:2229
void jacobianCalcPolyfit(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &poly_coeff, const Verbosity &)
WORKSPACE METHOD: jacobianCalcPolyfit.
Definition: m_jacobian.cc:912
void jacobianAddMagField(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dB, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddMagField.
Definition: m_jacobian.cc:1448
void jacobianFromTwoY(Matrix &jacobian, const Vector &y_pert, const Vector &y, const Numeric &pert_size, const Verbosity &)
WORKSPACE METHOD: jacobianFromTwoY.
Definition: m_jacobian.cc:2118
void AtmFieldPerturb(Tensor3 &perturbed_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &original_field, const Vector &p_ret_grid, const Vector &lat_ret_grid, const Vector &lon_ret_grid, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &)
WORKSPACE METHOD: AtmFieldPerturb.
Definition: m_jacobian.cc:1950
void jacobianFromYbatch(Matrix &jacobian, const ArrayOfVector &ybatch, const Vector &y, const Numeric &pert_size, const Verbosity &)
WORKSPACE METHOD: jacobianFromYbatch.
Definition: m_jacobian.cc:2133
void jacobianAddSurfaceQuantity(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddSurfaceQuantity.
Definition: m_jacobian.cc:1236
void jacobianAddShapeCatalogParameter(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const QuantumIdentifier &line_identity, const String &species, const String &variable, const String &coefficient, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddShapeCatalogParameter.
Definition: m_jacobian.cc:1526
void jacobianInit(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Verbosity &)
WORKSPACE METHOD: jacobianInit.
Definition: m_jacobian.cc:82
void jacobianAddFreqStretch(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqStretch.
Definition: m_jacobian.cc:345
void jacobianAddPointingZa(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Matrix &sensor_pos, const ArrayOfTime &sensor_time, const Index &poly_order, const String &calcmode, const Numeric &dza, const Verbosity &)
WORKSPACE METHOD: jacobianAddPointingZa.
Definition: m_jacobian.cc:512
void particle_bulkprop_fieldPerturb(Tensor4 &particle_bulkprop_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfString &particle_bulkprop_names, const String &particle_type, const Vector &p_ret_grid, const Vector &lat_ret_grid, const Vector &lon_ret_grid, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: particle_bulkprop_fieldPerturb.
Definition: m_jacobian.cc:2153
void jacobianAddNLTEs(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const ArrayOfQuantumIdentifier &energy_level_identities, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddNLTEs.
Definition: m_jacobian.cc:1728
void vmr_fieldPerturbAtmGrids(Tensor4 &vmr_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: vmr_fieldPerturbAtmGrids.
Definition: m_jacobian.cc:2276
void jacobianSetAffineTransformation(ArrayOfRetrievalQuantity &jqs, const Matrix &transformation_matrix, const Vector &offset_vector, const Verbosity &)
WORKSPACE METHOD: jacobianSetAffineTransformation.
Definition: m_jacobian.cc:1878
void jacobianAddNLTE(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const QuantumIdentifier &energy_level_identity, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddNLTE.
Definition: m_jacobian.cc:1667
void IndexNumberOfAtmosphericPoints(Index &n, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Verbosity &)
WORKSPACE METHOD: IndexNumberOfAtmosphericPoints.
Definition: m_jacobian.cc:2104
void jacobianCalcDoNothing(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Verbosity &)
WORKSPACE METHOD: jacobianCalcDoNothing.
Definition: m_jacobian.cc:57
void jacobianAddSinefit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Vector &period_lengths, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &)
WORKSPACE METHOD: jacobianAddSinefit.
Definition: m_jacobian.cc:1066
void jacobianAddTemperature(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &hse, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddTemperature.
Definition: m_jacobian.cc:1294
void particle_bulkprop_fieldPerturbAtmGrids(Tensor4 &particle_bulkprop_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfString &particle_bulkprop_names, const String &particle_type, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &verbosity)
WORKSPACE METHOD: particle_bulkprop_fieldPerturbAtmGrids.
Definition: m_jacobian.cc:2194
void jacobianAddBasicCatalogParameter(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const QuantumIdentifier &catalog_identity, const String &catalog_parameter, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddBasicCatalogParameter.
Definition: m_jacobian.cc:1603
void jacobianAddAbsSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Index &for_species_tag, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddAbsSpecies.
Definition: m_jacobian.cc:106
void jacobianAddShapeCatalogParameters(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfQuantumIdentifier &line_identities, const ArrayOfString &species, const ArrayOfString &variables, const ArrayOfString &coefficients, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddShapeCatalogParameters.
Definition: m_jacobian.cc:1569
void jacobianAddScatSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddScatSpecies.
Definition: m_jacobian.cc:999
void jacobianCalcPointingZaInterp(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &DEBUG_ONLY(sensor_los), const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const ArrayOfTime &sensor_time, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
Definition: m_jacobian.cc:584
void jacobianAddFreqShift(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqShift.
Definition: m_jacobian.cc:200
void AtmFieldPerturbAtmGrids(Tensor3 &perturbed_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &original_field, const Index &pert_index, const Numeric &pert_size, const String &pert_mode, const Verbosity &)
WORKSPACE METHOD: AtmFieldPerturbAtmGrids.
Definition: m_jacobian.cc:2042
void jacobianSetFuncTransformation(ArrayOfRetrievalQuantity &jqs, const String &transformation_func, const Numeric &z_min, const Numeric &z_max, const Verbosity &)
WORKSPACE METHOD: jacobianSetFuncTransformation.
Definition: m_jacobian.cc:1905
void jacobianAddBasicCatalogParameters(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfQuantumIdentifier &catalog_identities, const ArrayOfString &catalog_parameters, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddBasicCatalogParameters.
Definition: m_jacobian.cc:1644
void jacobianOff(Workspace &ws, Index &jacobian_do, Agenda &jacobian_agenda, ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianOff.
Definition: m_jacobian.cc:92
void jacobianAddPolyfit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Index &poly_order, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &)
WORKSPACE METHOD: jacobianAddPolyfit.
Definition: m_jacobian.cc:838
void jacobianCalcSinefit(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &period_index, const Verbosity &)
WORKSPACE METHOD: jacobianCalcSinefit.
Definition: m_jacobian.cc:1141
void jacobianCalcFreqStretch(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: jacobianCalcFreqStretch.
Definition: m_jacobian.cc:396
Workspace methods and template functions for supergeneric XML IO.
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:758
Array< Vector > ArrayOfVector
An array of vectors.
Definition: matpackI.h:1439
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:206
#define CREATE_OUT2
Definition: messages.h:205
Array< String > ArrayOfString
An array of Strings.
Definition: mystring.h:218
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 ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ArrayOfSpeciesTagVMR
Definition: jacobian.h:101
Temperature
Definition: jacobian.h:58
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ScatteringString
Definition: jacobian.h:102
MagneticMagnitude
Definition: jacobian.h:60
WindMagnitude
Definition: jacobian.h:59
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::Number::GlobalState QuantumIdentifier
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, Matrix &geo_pos_matrix, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const ConstVectorView &f_grid, const ConstMatrixView &sensor_pos, const ConstMatrixView &sensor_los, const ConstMatrixView &transmitter_pos, const ConstMatrixView &mblock_dlos_grid, const String &iy_unit, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
Performs calculations for one measurement block, on iy-level.
Definition: rte.cc:1608
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:1112
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:101
Holds all information required for individual partial derivatives.
Definition: jacobian.h:107
Numeric perturbation
Perturbations for methods where theoretical computations are impossible or plain slow.
Definition: jacobian.h:124
bool sameTargetType(const Target &other) const noexcept
Definition: jacobian.h:191
Species::Species species_id
Species ID for line parameters.
Definition: jacobian.h:136
A logical struct for global quantum numbers with species identifiers.
Species::Species Species() const noexcept
The Sparse class.
Definition: matpackII.h:72
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69