ARTS 2.5.4 (git: 7d04b88e)
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
26#include "jacobian.h"
27#include "arts.h"
28#include "auto_md.h"
29#include "check_input.h"
30#include "global_data.h"
31#include "lin_alg.h"
32#include "physics_funcs.h"
33#include "rte.h"
34#include "special_interp.h"
35
36extern const Numeric NAT_LOG_TEN;
37extern const Numeric PI;
38
39ostream& operator<<(ostream& os, const RetrievalQuantity& ot) {
40 return os << "\n Target = " << ot.Target()
41 << "\n Sub tag = " << ot.Subtag()
42 << "\n Mode = " << ot.Mode();
43}
44
46 bool& any_affine,
47 const ArrayOfRetrievalQuantity& jqs,
48 const bool& before_affine) {
49 jis.resize(jqs.nelem());
50
51 any_affine = false;
52
53 // Indices before affine transformation
54 if (before_affine) {
55 for (Index i = 0; i < jqs.nelem(); ++i) {
56 jis[i] = ArrayOfIndex(2);
57 if (i > 0) {
58 jis[i][0] = jis[i - 1][1] + 1;
59 } else {
60 jis[i][0] = 0;
61 }
62 const RetrievalQuantity& jq = jqs[i];
63 jis[i][1] = jis[i][0] + jq.nelem() - 1;
64 if (jq.HasAffine()) {
65 any_affine = true;
66 }
67 }
68 }
69
70 // After affine transformation
71 else {
72 for (Index i = 0; i < jqs.nelem(); ++i) {
73 jis[i] = ArrayOfIndex(2);
74 if (i > 0) {
75 jis[i][0] = jis[i - 1][1] + 1;
76 } else {
77 jis[i][0] = 0;
78 }
79 const RetrievalQuantity& jq = jqs[i];
80 if (jq.HasAffine()) {
81 jis[i][1] = jis[i][0] + jq.TransformationMatrix().ncols() - 1;
82 any_affine = true;
83 } else {
84 jis[i][1] = jis[i][0] + jq.nelem() - 1;
85 }
86 }
87 }
88}
89
91 const Vector x,
92 const ArrayOfRetrievalQuantity& jqs) {
93 // Range indices without affine trans
95 bool any_affine;
96 //
97 jac_ranges_indices(jis, any_affine, jqs, true);
98
99 Vector x_t(x);
100 transform_x_back(x_t, jqs, false);
101
102 // Apply functional transformations
103 for (Index i = 0; i < jqs.nelem(); ++i) {
104 const RetrievalQuantity& jq = jqs[i];
105 const String tfun = jq.TransformationFunc();
106 // Remember to add new functions also to transform_jacobian and transform_x_back
107 if (tfun == "") {
108 // Nothing to do
109 } else if (tfun == "log") {
110 for (Index c = jis[i][0]; c <= jis[i][1]; ++c) {
111 jacobian(joker, c) *= exp(x_t[c]);
112 }
113 } else if (tfun == "log10") {
114 for (Index c = jis[i][0]; c <= jis[i][1]; ++c) {
115 jacobian(joker, c) *= NAT_LOG_TEN * pow(10.0, x_t[c]);
116 }
117 } else if (tfun == "atanh") {
118 const Vector& pars = jq.TFuncParameters();
119 for (Index c = jis[i][0]; c <= jis[i][1]; ++c) {
120 jacobian(joker, c) *=
121 2 * (pars[1] - pars[0]) / pow(exp(-x_t[c]) + exp(x_t[c]), 2.0);
122 }
123 } else {
124 ARTS_ASSERT(0);
125 }
126 }
127
128 // Apply affine transformations
129 if (any_affine) {
131 jac_ranges_indices(jis_t, any_affine, jqs);
132
133 Matrix jacobian_t(jacobian.nrows(), jis_t.back()[1] + 1);
134
135 for (Index i = 0; i < jqs.nelem(); ++i) {
136 const RetrievalQuantity& jq = jqs[i];
137 Index col_start = jis[i][0];
138 Index col_extent = jis[i][1] - jis[i][0] + 1;
139 Range col_range(col_start, col_extent);
140 Index col_start_t = jis_t[i][0];
141 Index col_extent_t = jis_t[i][1] - jis_t[i][0] + 1;
142 Range col_range_t(col_start_t, col_extent_t);
143 if (jq.HasAffine()) {
144 mult(jacobian_t(joker, col_range_t),
145 jacobian(joker, col_range),
147 } else {
148 jacobian_t(joker, col_range_t) = jacobian(joker, col_range);
149 }
150 }
151 swap(jacobian_t, jacobian);
152 }
153}
154
156 // Range indices without affine trans
158 bool any_affine;
159 //
160 jac_ranges_indices(jis, any_affine, jqs, true);
161
162 // Apply functional transformations
163 for (Index i = 0; i < jqs.nelem(); ++i) {
164 const RetrievalQuantity& jq = jqs[i];
165 const String tfun = jq.TransformationFunc();
166 // Remember to add new functions also to transform_jacobian and transform_x_back
167 if (tfun == "") {
168 // Nothing to do
169 } else if (tfun == "log") {
170 const Vector& pars = jq.TFuncParameters();
171 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
172 ARTS_USER_ERROR_IF (x[r] <= pars[0],
173 "log-transformation selected for retrieval quantity with\n"
174 "index ", i, " (0-based), but at least one value <= z_min\n"
175 "found for this quantity. This is not allowed.")
176 x[r] = log(x[r] - pars[0]);
177 }
178 } else if (tfun == "log10") {
179 const Vector& pars = jq.TFuncParameters();
180 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
181 ARTS_USER_ERROR_IF (x[r] <= 0,
182 "log10-transformation selected for retrieval quantity with\n"
183 "index ", i, " (0-based), but at least one value <= z_min\n"
184 "found for this quantity. This is not allowed.")
185 x[r] = log10(x[r] - pars[0]);
186 }
187 } else if (tfun == "atanh") {
188 const Vector& pars = jq.TFuncParameters();
189 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
190 ARTS_USER_ERROR_IF (x[r] <= pars[0],
191 "atanh-transformation selected for retrieval quantity with\n"
192 "index ", i, " (0-based), but at least one value <= z_min\n"
193 "found for this quantity. This is not allowed.")
194 ARTS_USER_ERROR_IF (x[r] >= pars[1],
195 "atanh-transformation selected for retrieval quantity with\n"
196 "index ", i, " (0-based), but at least one value is\n"
197 ">= z_max. This is not allowed.")
198 x[r] = atanh(2 * (x[r] - pars[0]) / (pars[1] - pars[0]) - 1);
199 }
200 } else {
201 ARTS_ASSERT(0);
202 }
203 }
204
205 // Apply affine transformations
206 if (any_affine) {
208 jac_ranges_indices(jis_t, any_affine, jqs);
209
210 Vector x_t(jis_t.back()[1] + 1);
211
212 for (Index i = 0; i < jqs.nelem(); ++i) {
213 const RetrievalQuantity& jq = jqs[i];
214 Index col_start = jis[i][0];
215 Index col_extent = jis[i][1] - jis[i][0] + 1;
216 Range col_range(col_start, col_extent);
217 Index col_start_t = jis_t[i][0];
218 Index col_extent_t = jis_t[i][1] - jis_t[i][0] + 1;
219 Range col_range_t(col_start_t, col_extent_t);
220 if (jq.HasAffine()) {
221 Vector t(x[col_range]);
222 t -= jq.OffsetVector();
223 mult(x_t[col_range_t], transpose(jq.TransformationMatrix()), t);
224 } else {
225 x_t[col_range_t] = x[col_range];
226 }
227 }
228 swap(x, x_t);
229 }
230}
231
233 const ArrayOfRetrievalQuantity& jqs,
234 bool revert_functional_transforms) {
235 // Range indices without affine trans
237 bool any_affine;
238 //
239 jac_ranges_indices(jis, any_affine, jqs, true);
240
241 // Revert affine transformations
242 // Apply affine transformations
243 if (any_affine) {
245 jac_ranges_indices(jis_t, any_affine, jqs);
246
247 Vector x(jis.back()[1] + 1);
248
249 for (Index i = 0; i < jqs.nelem(); ++i) {
250 const RetrievalQuantity& jq = jqs[i];
251 Index col_start = jis[i][0];
252 Index col_extent = jis[i][1] - jis[i][0] + 1;
253 Range col_range(col_start, col_extent);
254 Index col_start_t = jis_t[i][0];
255 Index col_extent_t = jis_t[i][1] - jis_t[i][0] + 1;
256 Range col_range_t(col_start_t, col_extent_t);
257 if (jq.HasAffine()) {
258 mult(x[col_range], jq.TransformationMatrix(), x_t[col_range_t]);
259 x[col_range] += jq.OffsetVector();
260 } else {
261 x[col_range] = x_t[col_range_t];
262 }
263 }
264 swap(x_t, x);
265 }
266
267 if (revert_functional_transforms) {
268 // Revert functional transformations
269 for (Index i = 0; i < jqs.nelem(); ++i) {
270 const RetrievalQuantity& jq = jqs[i];
271 const String tfun = jq.TransformationFunc();
272 // Remember to add new functions also to transform_jacobian and transform_x_back
273 if (tfun == "") {
274 // Nothing to do
275 } else if (tfun == "log") {
276 const Vector& pars = jq.TFuncParameters();
277 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
278 x_t[r] = pars[0] + exp(x_t[r]);
279 }
280 } else if (tfun == "log10") {
281 const Vector& pars = jq.TFuncParameters();
282 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
283 x_t[r] = pars[0] + pow(10.0, x_t[r]);
284 }
285 } else if (tfun == "atanh") {
286 const Vector& pars = jq.TFuncParameters();
287 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
288 x_t[r] = pars[0] + ((pars[1] - pars[0]) / 2) * (1 + tanh(x_t[r]));
289 }
290 } else {
291 ARTS_ASSERT(0);
292 }
293 }
294 }
295}
296
297/*===========================================================================
298 === Help sub-functions to handle analytical jacobians (in alphabetical order)
299 ===========================================================================*/
300
301// Small help function, to make the code below cleaner
303 ConstMatrixView diy_dq,
304 const Numeric& w) {
305 for (Index irow = 0; irow < diy_dx.nrows(); irow++) {
306 for (Index icol = 0; icol < diy_dx.ncols(); icol++) {
307 diy_dx(irow, icol) += w * diy_dq(irow, icol);
308 }
309 }
310}
311
313 const RetrievalQuantity& jacobian_quantity,
314 ConstTensor3View diy_dpath,
315 const Index& atmosphere_dim,
316 const Ppath& ppath,
317 ConstVectorView ppath_p) {
318 ARTS_ASSERT(jacobian_quantity.Grids().nelem() == atmosphere_dim or jacobian_quantity.Grids().empty());
319
320 // We want here an extrapolation to infinity ->
321 // extremly high extrapolation factor
322 const Numeric extpolfac = 1.0e99;
323
324 if (ppath.np > 1) // Otherwise nothing to do here
325 {
326 // Pressure
327 Index nr1 = jacobian_quantity.Grids().empty() ? 0 : jacobian_quantity.Grids()[0].nelem();
328 ArrayOfGridPos gp_p(ppath.np);
329 if (nr1 > 1) {
330 p2gridpos(gp_p, jacobian_quantity.Grids()[0], ppath_p, extpolfac);
332 } else {
333 gp4length1grid(gp_p);
334 }
335
336 // Latitude
337 Index nr2 = 1;
338 ArrayOfGridPos gp_lat;
339 if (atmosphere_dim > 1) {
340 gp_lat.resize(ppath.np);
341 nr2 = jacobian_quantity.Grids().empty() ? 0 : jacobian_quantity.Grids()[1].nelem();
342 if (nr2 > 1) {
343 gridpos(gp_lat,
344 jacobian_quantity.Grids()[1],
345 ppath.pos(joker, 1),
346 extpolfac);
348 } else {
349 gp4length1grid(gp_lat);
350 }
351 }
352
353 // Longitude
354 ArrayOfGridPos gp_lon;
355 if (atmosphere_dim > 2) {
356 Index nr3 = jacobian_quantity.Grids().empty() ? 0 : jacobian_quantity.Grids()[2].nelem();
357 gp_lon.resize(ppath.np);
358 if (nr3 > 1) {
359 gridpos(gp_lon,
360 jacobian_quantity.Grids()[2],
361 ppath.pos(joker, 2),
362 extpolfac);
364 } else {
365 gp4length1grid(gp_lon);
366 }
367 }
368
369 //- 1D
370 if (atmosphere_dim == 1) {
371 for (Index ip = 0; ip < ppath.np; ip++) {
372 if (gp_p[ip].fd[1] > 0) {
373 from_dpath_to_dx(diy_dx(gp_p[ip].idx, joker, joker),
374 diy_dpath(ip, joker, joker),
375 gp_p[ip].fd[1]);
376 }
377 if (gp_p[ip].fd[0] > 0) {
378 from_dpath_to_dx(diy_dx(gp_p[ip].idx + 1, joker, joker),
379 diy_dpath(ip, joker, joker),
380 gp_p[ip].fd[0]);
381 }
382 }
383 }
384
385 //- 2D
386 else if (atmosphere_dim == 2) {
387 for (Index ip = 0; ip < ppath.np; ip++) {
388 Index ix = nr1 * gp_lat[ip].idx + gp_p[ip].idx;
389 // Low lat, low p
390 if (gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[1] > 0)
391 from_dpath_to_dx(diy_dx(ix, joker, joker),
392 diy_dpath(ip, joker, joker),
393 gp_lat[ip].fd[1] * gp_p[ip].fd[1]);
394 // Low lat, high p
395 if (gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[0] > 0)
396 from_dpath_to_dx(diy_dx(ix + 1, joker, joker),
397 diy_dpath(ip, joker, joker),
398 gp_lat[ip].fd[1] * gp_p[ip].fd[0]);
399 // High lat, low p
400 if (gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[1] > 0)
401 from_dpath_to_dx(diy_dx(ix + nr1, joker, joker),
402 diy_dpath(ip, joker, joker),
403 gp_lat[ip].fd[0] * gp_p[ip].fd[1]);
404 // High lat, high p
405 if (gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[0] > 0)
406 from_dpath_to_dx(diy_dx(ix + nr1 + 1, joker, joker),
407 diy_dpath(ip, joker, joker),
408 gp_lat[ip].fd[0] * gp_p[ip].fd[0]);
409 }
410 }
411
412 //- 3D
413 else if (atmosphere_dim == 3) {
414 for (Index ip = 0; ip < ppath.np; ip++) {
415 Index ix =
416 nr2 * nr1 * gp_lon[ip].idx + nr1 * gp_lat[ip].idx + gp_p[ip].idx;
417 // Low lon, low lat, low p
418 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[1] > 0)
420 diy_dx(ix, joker, joker),
421 diy_dpath(ip, joker, joker),
422 gp_lon[ip].fd[1] * gp_lat[ip].fd[1] * gp_p[ip].fd[1]);
423 // Low lon, low lat, high p
424 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[0] > 0)
426 diy_dx(ix + 1, joker, joker),
427 diy_dpath(ip, joker, joker),
428 gp_lon[ip].fd[1] * gp_lat[ip].fd[1] * gp_p[ip].fd[0]);
429 // Low lon, high lat, low p
430 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[1] > 0)
432 diy_dx(ix + nr1, joker, joker),
433 diy_dpath(ip, joker, joker),
434 gp_lon[ip].fd[1] * gp_lat[ip].fd[0] * gp_p[ip].fd[1]);
435 // Low lon, high lat, high p
436 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[0] > 0)
438 diy_dx(ix + nr1 + 1, joker, joker),
439 diy_dpath(ip, joker, joker),
440 gp_lon[ip].fd[1] * gp_lat[ip].fd[0] * gp_p[ip].fd[0]);
441
442 // Increase *ix* (to be valid for high lon level)
443 ix += nr2 * nr1;
444
445 // High lon, low lat, low p
446 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[1] > 0)
448 diy_dx(ix, joker, joker),
449 diy_dpath(ip, joker, joker),
450 gp_lon[ip].fd[0] * gp_lat[ip].fd[1] * gp_p[ip].fd[1]);
451 // High lon, low lat, high p
452 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[0] > 0)
454 diy_dx(ix + 1, joker, joker),
455 diy_dpath(ip, joker, joker),
456 gp_lon[ip].fd[0] * gp_lat[ip].fd[1] * gp_p[ip].fd[0]);
457 // High lon, high lat, low p
458 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[1] > 0)
460 diy_dx(ix + nr1, joker, joker),
461 diy_dpath(ip, joker, joker),
462 gp_lon[ip].fd[0] * gp_lat[ip].fd[0] * gp_p[ip].fd[1]);
463 // High lon, high lat, high p
464 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[0] > 0)
466 diy_dx(ix + nr1 + 1, joker, joker),
467 diy_dpath(ip, joker, joker),
468 gp_lon[ip].fd[0] * gp_lat[ip].fd[0] * gp_p[ip].fd[0]);
469 }
470 }
471 }
472}
473
475 const RetrievalQuantity& jacobian_quantity,
476 ConstMatrixView diy_dpos,
477 const Index& atmosphere_dim,
478 ConstVectorView rtp_pos) {
479 ARTS_ASSERT(jacobian_quantity.Grids().nelem() ==
480 max(atmosphere_dim - 1, Index(1)));
481 ARTS_ASSERT(rtp_pos.nelem() == atmosphere_dim);
482
483 // We want here an extrapolation to infinity ->
484 // extremly high extrapolation factor
485 const Numeric extpolfac = 1.0e99;
486
487 // Handle 1D separately
488 if (atmosphere_dim == 1) {
489 diy_dx(0, joker, joker) = diy_dpos;
490 return;
491 }
492
493 // Latitude
494 Index nr1 = 1;
495 ArrayOfGridPos gp_lat;
496 {
497 gp_lat.resize(1);
498 nr1 = jacobian_quantity.Grids()[0].nelem();
499 if (nr1 > 1) {
500 gridpos(gp_lat,
501 jacobian_quantity.Grids()[0],
502 Vector(1, rtp_pos[1]),
503 extpolfac);
505 } else {
506 gp4length1grid(gp_lat);
507 }
508 }
509
510 // Longitude
511 ArrayOfGridPos gp_lon;
512 if (atmosphere_dim > 2) {
513 gp_lon.resize(1);
514 if (jacobian_quantity.Grids()[1].nelem() > 1) {
515 gridpos(gp_lon,
516 jacobian_quantity.Grids()[1],
517 Vector(1, rtp_pos[2]),
518 extpolfac);
520 } else {
521 gp4length1grid(gp_lon);
522 }
523 }
524
525 //- 2D
526 if (atmosphere_dim == 2) {
527 if (gp_lat[0].fd[1] > 0) {
528 from_dpath_to_dx(diy_dx(gp_lat[0].idx, joker, joker),
529 diy_dpos(joker, joker),
530 gp_lat[0].fd[1]);
531 }
532 if (gp_lat[0].fd[0] > 0) {
533 from_dpath_to_dx(diy_dx(gp_lat[0].idx + 1, joker, joker),
534 diy_dpos(joker, joker),
535 gp_lat[0].fd[0]);
536 }
537 }
538 //- 3D
539 else {
540 Index ix = nr1 * gp_lon[0].idx + gp_lat[0].idx;
541 // Low lon, low lat
542 if (gp_lon[0].fd[1] > 0 && gp_lat[0].fd[1] > 0)
543 from_dpath_to_dx(diy_dx(ix, joker, joker),
544 diy_dpos(joker, joker),
545 gp_lon[0].fd[1] * gp_lat[0].fd[1]);
546 // Low lon, high lat
547 if (gp_lon[0].fd[1] > 0 && gp_lat[0].fd[0] > 0)
548 from_dpath_to_dx(diy_dx(ix + 1, joker, joker),
549 diy_dpos(joker, joker),
550 gp_lon[0].fd[1] * gp_lat[0].fd[0]);
551 // High lon, low lat
552 if (gp_lon[0].fd[0] > 0 && gp_lat[0].fd[1] > 0)
553 from_dpath_to_dx(diy_dx(ix + nr1, joker, joker),
554 diy_dpos(joker, joker),
555 gp_lon[0].fd[0] * gp_lat[0].fd[1]);
556 // High lon, high lat
557 if (gp_lon[0].fd[0] > 0 && gp_lat[0].fd[0] > 0)
558 from_dpath_to_dx(diy_dx(ix + nr1 + 1, joker, joker),
559 diy_dpos(joker, joker),
560 gp_lon[0].fd[0] * gp_lat[0].fd[0]);
561 }
562}
563
565 const ArrayOfArrayOfSpeciesTag& abs_species) {
566 ArrayOfIndex aoi(jacobian_quantities.nelem(), -1);
567
569 if (jacobian_quantities[iq] == Jacobian::Line::VMR) {
570 auto p = std::find_if(abs_species.cbegin(), abs_species.cend(),
571 [qid=jacobian_quantities[iq].QuantumIdentity()](auto& specs){
572 return std::any_of(specs.cbegin(), specs.cend(),
573 [qid](auto& spec){return qid.Isotopologue() == spec.Isotopologue();});
574 });
575 if (p not_eq abs_species.cend()) {
576 aoi[iq] = Index(abs_species.cend() - p);
577 } else {
579 "Could not find ",
580 jacobian_quantities[iq].Subtag(),
581 " in species of abs_species.\n")
582 }
583 } else if (jacobian_quantities[iq] == Jacobian::Special::ArrayOfSpeciesTagVMR) {
584 ArrayOfSpeciesTag atag(jacobian_quantities[iq].Subtag());
585 aoi[iq] = chk_contains("abs_species", abs_species, atag);
586 } else if (jacobian_quantities[iq] == Jacobian::Atm::Particulates) {
587 aoi[iq] = -9999;
588 } else if (jacobian_quantities[iq] == Jacobian::Atm::Electrons) {
589 aoi[iq] = -9999;
590 }
591 )
592
593 return aoi;
594}
595
596ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity& jacobian_quantities, Index np, Index nf, Index ns, bool active) {
597 ArrayOfTensor3 diy_dpath(jacobian_quantities.nelem());
598
599 const Index nn = active ? np * nf : nf;
600 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq] = Tensor3(np, nn, ns, 0.0);)
601
602 return diy_dpath;
603}
604
605ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity& jacobian_quantities, Index np, Index nf, Index ns, bool active) {
606 ArrayOfTensor3 diy_dx(jacobian_quantities.nelem());
607
608 bool any_affine;
609 ArrayOfArrayOfIndex jacobian_indices;
610 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
611
612 const Index nn = active ? np * nf : nf;
613 FOR_ANALYTICAL_JACOBIANS_DO2(diy_dx[iq] = Tensor3(jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1, nn, ns, 0.0);)
614
615 return diy_dx;
616}
617
619 const ArrayOfString& scat_species,
620 const bool cloudbox_on) {
621 ArrayOfIndex aoi(jacobian_quantities.nelem(), -1);
622
624 if (cloudbox_on and jacobian_quantities[iq] == Jacobian::Special::ScatteringString) {
625 aoi[iq] = find_first(scat_species, jacobian_quantities[iq].Subtag());
626 ARTS_USER_ERROR_IF (aoi[iq] < 0,
627 "Jacobian quantity with index ", iq, " refers to\n"
628 " ", jacobian_quantities[iq].Subtag(),
629 "\nbut this species could not be found in *scat_species*.")
630 }
631 )
632
633 return aoi;
634}
635
636/*===========================================================================
637 === Other functions, in alphabetical order
638 ===========================================================================*/
639
641 ostringstream& os,
642 const Vector& p_grid,
643 const Vector& lat_grid,
644 const Vector& lon_grid,
645 const Vector& p_retr,
646 const Vector& lat_retr,
647 const Vector& lon_retr,
648 const String& p_retr_name,
649 const String& lat_retr_name,
650 const String& lon_retr_name,
651 const Index& dim) {
652 ARTS_ASSERT(grids.nelem() == dim);
653
654 if (p_retr.nelem() == 0) {
655 os << "The grid vector *" << p_retr_name << "* is empty,"
656 << " at least one pressure level\n"
657 << "should be specified.";
658 return false;
659 } else if (!is_decreasing(p_retr)) {
660 os << "The pressure grid vector *" << p_retr_name << "* is not a\n"
661 << "strictly decreasing vector, which is required.";
662 return false;
663 } else if (p_grid.nelem() == 1 and p_grid.nelem() == p_retr.nelem()) {
664 if (p_grid[0] not_eq p_retr[0]) {
665 os << "Mismatching 1-long grids for " << p_retr_name;
666 return false;
667 }
668
669 // Necessary repeat but grids are OK
670 grids[0] = p_retr;
671 } else if (log(p_retr[0]) > 1.5 * log(p_grid[0]) - 0.5 * log(p_grid[1]) ||
672 log(p_retr[p_retr.nelem() - 1]) <
673 1.5 * log(p_grid[p_grid.nelem() - 1]) -
674 0.5 * log(p_grid[p_grid.nelem() - 2])) {
675 os << "The grid vector *" << p_retr_name << "* is not covered by the\n"
676 << "corresponding atmospheric grid.";
677 return false;
678 } else {
679 // Pressure grid ok, add it to grids
680 grids[0] = p_retr;
681 }
682
683 if (dim >= 2) {
684 // If 2D and 3D atmosphere, check latitude grid
685 if (lat_retr.nelem() == 0) {
686 os << "The grid vector *" << lat_retr_name << "* is empty,"
687 << " at least one latitude\n"
688 << "should be specified for a 2D/3D atmosphere.";
689 return false;
690 } else if (!is_increasing(lat_retr)) {
691 os << "The latitude grid vector *" << lat_retr_name << "* is not a\n"
692 << "strictly increasing vector, which is required.";
693 return false;
694 } else if (lat_grid.nelem() == 1 and lat_grid.nelem() == lat_retr.nelem()) {
695 if (lat_grid[0] not_eq lat_retr[0]) {
696 os << "Mismatching 1-long grids for " << lat_retr_name;
697 return false;
698 }
699
700 // Necessary repeat but grids are OK
701 grids[1] = lat_retr;
702 } else if (lat_retr[0] < 1.5 * lat_grid[0] - 0.5 * lat_grid[1] ||
703 lat_retr[lat_retr.nelem() - 1] >
704 1.5 * lat_grid[lat_grid.nelem() - 1] -
705 0.5 * lat_grid[lat_grid.nelem() - 2]) {
706 os << "The grid vector *" << lat_retr_name << "* is not covered by the\n"
707 << "corresponding atmospheric grid.";
708 return false;
709 } else {
710 // Latitude grid ok, add it to grids
711 grids[1] = lat_retr;
712 }
713 if (dim == 3) {
714 // For 3D atmosphere check longitude grid
715 if (lon_retr.nelem() == 0) {
716 os << "The grid vector *" << lon_retr_name << "* is empty,"
717 << " at least one longitude\n"
718 << "should be specified for a 3D atmosphere.";
719 return false;
720 } else if (!is_increasing(lon_retr)) {
721 os << "The longitude grid vector *" << lon_retr_name << "* is not a\n"
722 << "strictly increasing vector, which is required.";
723 return false;
724 } else if (lon_grid.nelem() == 1 and
725 lon_grid.nelem() == lon_retr.nelem()) {
726 if (lon_grid[0] not_eq lon_retr[0]) {
727 os << "Mismatching 1-long grids for " << lon_retr_name;
728 return false;
729 }
730
731 // Necessary repeat but grids are OK
732 grids[2] = lon_retr;
733 } else if (lon_retr[0] < 1.5 * lon_grid[0] - 0.5 * lon_grid[1] ||
734 lon_retr[lon_retr.nelem() - 1] >
735 1.5 * lon_grid[lon_grid.nelem() - 1] -
736 0.5 * lon_grid[lon_grid.nelem() - 2]) {
737 os << "The grid vector *" << lon_retr_name
738 << "* is not covered by the\n"
739 << "corresponding atmospheric grid.";
740 return false;
741 } else {
742 // Longitude grid ok, add it to grids
743 grids[2] = lon_retr;
744 }
745 }
746 }
747 return true;
748}
749
751 ostringstream& os,
752 const Vector& lat_grid,
753 const Vector& lon_grid,
754 const Vector& lat_retr,
755 const Vector& lon_retr,
756 const String& lat_retr_name,
757 const String& lon_retr_name,
758 const Index& dim) {
759 ARTS_ASSERT(grids.nelem() == max(dim - 1, Index(1)));
760
761 if (dim == 1) {
762 // Here we only need to create a length 1 dummy grid
763 grids[0].resize(1);
764 grids[0][0] = 0;
765 }
766
767 if (dim >= 2) {
768 // If 2D and 3D atmosphere, check latitude grid
769 if (lat_retr.nelem() == 0) {
770 os << "The grid vector *" << lat_retr_name << "* is empty,"
771 << " at least one latitude\n"
772 << "should be specified for a 2D/3D atmosphere.";
773 return false;
774 } else if (!is_increasing(lat_retr)) {
775 os << "The latitude grid vector *" << lat_retr_name << "* is not a\n"
776 << "strictly increasing vector, which is required.";
777 return false;
778 } else if (lat_grid.nelem() == 1 and lat_grid.nelem() == lat_retr.nelem()) {
779 if (lat_grid[0] not_eq lat_retr[0]) {
780 os << "Mismatching 1-long grids for " << lat_retr_name;
781 return false;
782 }
783
784 // Necessary repeat but grids are OK
785 grids[0] = lat_retr;
786 } else if (lat_retr[0] < 1.5 * lat_grid[0] - 0.5 * lat_grid[1] ||
787 lat_retr[lat_retr.nelem() - 1] >
788 1.5 * lat_grid[lat_grid.nelem() - 1] -
789 0.5 * lat_grid[lat_grid.nelem() - 2]) {
790 os << "The grid vector *" << lat_retr_name << "* is not covered by the\n"
791 << "corresponding atmospheric grid.";
792 return false;
793 } else {
794 // Latitude grid ok, add it to grids
795 grids[0] = lat_retr;
796 }
797
798 if (dim == 3) {
799 // For 3D atmosphere check longitude grid
800 if (lon_retr.nelem() == 0) {
801 os << "The grid vector *" << lon_retr_name << "* is empty,"
802 << " at least one longitude\n"
803 << "should be specified for a 3D atmosphere.";
804 return false;
805 } else if (!is_increasing(lon_retr)) {
806 os << "The longitude grid vector *" << lon_retr_name << "* is not a\n"
807 << "strictly increasing vector, which is required.";
808 return false;
809 } else if (lon_grid.nelem() == 1 and
810 lon_grid.nelem() == lon_retr.nelem()) {
811 if (lon_grid[0] not_eq lon_retr[0]) {
812 os << "Mismatching 1-long grids for " << lon_retr_name;
813 return false;
814 }
815
816 // Necessary repeat but grids are OK
817 grids[1] = lon_retr;
818 } else if (lon_retr[0] < 1.5 * lon_grid[0] - 0.5 * lon_grid[1] ||
819 lon_retr[lon_retr.nelem() - 1] >
820 1.5 * lon_grid[lon_grid.nelem() - 1] -
821 0.5 * lon_grid[lon_grid.nelem() - 2]) {
822 os << "The grid vector *" << lon_retr_name
823 << "* is not covered by the\n"
824 << "corresponding atmospheric grid.";
825 return false;
826 } else {
827 // Longitude grid ok, add it to grids
828 grids[1] = lon_retr;
829 }
830 }
831 }
832 return true;
833}
834
836 for (Index i = 0; i < gp.nelem(); i++) {
837 if (gp[i].fd[0] < 0) {
838 gp[i].fd[0] = 0;
839 gp[i].fd[1] = 1;
840 } else if (gp[i].fd[0] > 1) {
841 gp[i].fd[0] = 1;
842 gp[i].fd[1] = 0;
843 }
844 }
845}
846
848 const Vector& x,
849 const Index& poly_coeff) {
850 const Index l = x.nelem();
851
852 ARTS_ASSERT(l > poly_coeff);
853
854 if (b.nelem() != l) b.resize(l);
855
856 if (poly_coeff == 0) {
857 b = 1.0;
858 } else {
859 const Numeric xmin = min(x);
860 const Numeric dx = 0.5 * (max(x) - xmin);
861 //
862 for (Index i = 0; i < l; i++) {
863 b[i] = (x[i] - xmin) / dx - 1.0;
864 b[i] = pow(b[i], int(poly_coeff));
865 }
866 //
867 b -= mean(b);
868 }
869}
870
871void calcBaselineFit(Vector& y_baseline,
872 const Vector& x,
873 const Index& mblock_index,
874 const Sparse& sensor_response,
875 const ArrayOfIndex& sensor_response_pol_grid,
876 const Vector& sensor_response_f_grid,
877 const Matrix& sensor_response_dlos_grid,
878 const RetrievalQuantity& rq,
879 const Index rq_index,
880 const ArrayOfArrayOfIndex& jacobian_indices) {
881 bool is_sine_fit = false;
882 if (rq == Jacobian::Sensor::Polyfit) {
883 is_sine_fit = false;
884 } else if (rq == Jacobian::Sensor::Sinefit) {
885 is_sine_fit = true;
886 } else {
888 "Retrieval quantity is neither a polynomial or a sine "
889 " baseline fit.");
890 }
891
892 // Size and check of sensor_response
893 //
894 const Index nf = sensor_response_f_grid.nelem();
895 const Index npol = sensor_response_pol_grid.nelem();
896 const Index nlos = sensor_response_dlos_grid.nrows();
897
898 // Evaluate basis functions for fits.
899 Vector w, s, c;
900 if (is_sine_fit) {
901 s.resize(nf);
902 c.resize(nf);
903 Numeric period = rq.Grids()[0][0];
904 for (Index f = 0; f < nf; f++) {
905 Numeric a = (sensor_response_f_grid[f] - sensor_response_f_grid[0]) * 2 *
906 PI / period;
907 s[f] = sin(a);
908 c[f] = cos(a);
909 }
910 } else {
911 Numeric poly_coeff = rq.Grids()[0][0];
913 w, sensor_response_f_grid, static_cast<Index>(poly_coeff));
914 }
915
916 // Compute baseline
917 ArrayOfVector jg = rq.Grids();
918 const Index n1 = jg[1].nelem();
919 const Index n2 = jg[2].nelem();
920 const Index n3 = jg[3].nelem();
921 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
922 const Index row4 = rowind.get_start();
923 Index col4 = jacobian_indices[rq_index][0];
924
925 if (n3 > 1) {
926 col4 += mblock_index * n2 * n1;
927 }
928
929 for (Index l = 0; l < nlos; l++) {
930 const Index row3 = row4 + l * nf * npol;
931 const Index col3 = col4 + l * n1 * (is_sine_fit ? 2 : 1);
932
933 for (Index f = 0; f < nf; f++) {
934 const Index row2 = row3 + f * npol;
935
936 for (Index p = 0; p < npol; p++) {
937 Index col1 = col3;
938 if (n1 > 1) {
939 col1 += p;
940 }
941 if (is_sine_fit) {
942 y_baseline[row2 + p] += x[col1] * s[f] + x[col1 + 1] * c[f];
943 } else {
944 y_baseline[row2 + p] += w[f] * x[col1];
945 }
946 }
947 }
948 }
949}
950
952 const String& unit,
953 const Numeric& vmr,
954 const Numeric& p,
955 const Numeric& t) {
956 if (unit == "rel" || unit == "logrel") {
957 x = 1;
958 } else if (unit == "vmr") {
959 if (vmr == 0) {
960 x = 0;
961 return;
962 }
963 x = 1 / vmr;
964 } else if (unit == "nd") {
965 if (vmr == 0) {
966 x = 0;
967 return;
968 }
969 x = 1 / (vmr * number_density(p, t));
970 } else {
972 "Allowed options for gas species jacobians are "
973 "\"rel\", \"vmr\" and \"nd\".\nYou have selected: ",
974 unit, '\n')
975 }
976}
977
979 const String& unit,
980 const Numeric& vmr,
981 const Numeric& p,
982 const Numeric& t) {
983 if (unit == "rel" || unit == "logrel") {
984 x = vmr;
985 } else if (unit == "vmr") {
986 x = 1;
987 } else if (unit == "nd") {
988 x = 1 / number_density(p, t);
989 } else {
991 "Allowed options for gas species jacobians are "
992 "\"rel\", \"vmr\" and \"nd\".\nYou have selected: ",
993 unit, '\n')
994 }
995}
996
997//======================================================================
998// Propmat partials descriptions
999//======================================================================
1000
1001bool is_wind_parameter(const RetrievalQuantity& t) noexcept {
1002 return t.Target().isWind();
1003}
1004
1006 return t.Target().isWind() or t.Target().isFrequency();
1007}
1008
1011}
1012
1013bool is_nlte_parameter(const RetrievalQuantity& t) noexcept {
1014 return t == Jacobian::Line::NLTE;
1015}
1016
1017#define ISLINESHAPETYPE(X) \
1018 bool is_pressure_broadening_##X(const RetrievalQuantity& t) noexcept { \
1019 return t == Jacobian::Line::Shape##X##X0 or \
1020 t == Jacobian::Line::Shape##X##X1 or \
1021 t == Jacobian::Line::Shape##X##X2 or \
1022 t == Jacobian::Line::Shape##X##X3; \
1023 }
1033#undef ISLINESHAPETYPE
1034
1035#define VARISLINESHAPEPARAM(X, Y) (t == Jacobian::Line::Shape##X##Y)
1037 return VARISLINESHAPEPARAM(G0, X0) or VARISLINESHAPEPARAM(D0, X0) or
1041 VARISLINESHAPEPARAM(DV, X0);
1042}
1043
1045 return VARISLINESHAPEPARAM(G0, X1) or VARISLINESHAPEPARAM(D0, X1) or
1049 VARISLINESHAPEPARAM(DV, X1);
1050}
1051
1053 return VARISLINESHAPEPARAM(G0, X2) or VARISLINESHAPEPARAM(D0, X2) or
1057 VARISLINESHAPEPARAM(DV, X2);
1058}
1059
1061 return VARISLINESHAPEPARAM(G0, X3) or VARISLINESHAPEPARAM(D0, X3) or
1065 VARISLINESHAPEPARAM(DV, X3);
1066}
1067#undef VARISLINESHAPEPARAM
1068
1070 const RetrievalQuantity& t) noexcept {
1074}
1075
1082}
1083
1084bool is_line_parameter(const RetrievalQuantity& t) noexcept {
1085 return t == Jacobian::Type::Line;
1086}
1087
1089 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Special::ArrayOfSpeciesTagVMR or j == Jacobian::Line::VMR or is_frequency_parameter(j));});
1090}
1091
1093 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Line::VMR or j == Jacobian::Special::ArrayOfSpeciesTagVMR or is_frequency_parameter(j));});
1094}
1095
1097 ARTS_USER_ERROR_IF (std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_line_parameter(j);}),
1098 "Line specific parameters are not supported while using continuum tags.\nWe do not track what lines are in the continuum.\n")
1099 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Special::ArrayOfSpeciesTagVMR or is_frequency_parameter(j));});
1100}
1101
1103 ARTS_USER_ERROR_IF (std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_line_parameter(j);}),
1104 "Line specific parameters are not supported while\n using the relaxation matrix line mixing routine.\n We do not yet track individual lines in the relaxation matrix calculations.\n")
1105 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or is_frequency_parameter(j));});
1106}
1107
1109 ARTS_USER_ERROR_IF (std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_line_parameter(j);}),
1110 "Line specific parameters are not supported while using Lookup table.\nWe do not track lines in the Lookup.\n")
1111 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Special::ArrayOfSpeciesTagVMR or is_frequency_parameter(j));});
1112}
1113
1115 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return not (j == Jacobian::Type::Sensor);});
1116}
1117
1119 if (rq == Jacobian::Line::VMR) { // Single tag
1120 for (auto& st : ast) {
1121 if ((rq.QuantumIdentity().isotopologue_index == st.spec_ind) or
1122 (rq.QuantumIdentity().Species() == st.Spec() and
1124 st.Isotopologue().isotname == Species::Joker)))
1125 return true;
1126 }
1127 } else {
1128 return rq.Target().species_array_id == ast;
1129 }
1130
1131 return false;
1132}
1133
1134bool species_match(const RetrievalQuantity& rq, const Species::Species species) {
1135 if (rq == Jacobian::Line::VMR and rq.QuantumIdentity().Species() == species)
1136 return true;
1137 return false;
1138}
1139
1141 const Species::IsotopeRecord& ir) {
1142 auto ir2 = rq.QuantumIdentity().Isotopologue();
1143 if (ir.spec == ir2.spec and (ir.isotname == Species::Joker or ir.isotname == ir2.isotname))
1144 return true;
1145 else
1146 return false;
1147}
1148
1150 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return j == Jacobian::Atm::Temperature;});
1151}
1152
1154 const QuantumIdentifier& line_qid) noexcept {
1155 auto p = std::find_if(js.cbegin(), js.cend(), [&line_qid](auto& j) {
1156 return j == Jacobian::Line::VMR
1157 and j.QuantumIdentity().Species() == line_qid.Species()
1158 and j.QuantumIdentity().Isotopologue() == line_qid.Isotopologue();}
1159 );
1160 if (p not_eq js.cend())
1161 return {true, p -> QuantumIdentity()};
1162 else
1163 return {false, line_qid};
1164}
1165
1167 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return j == Jacobian::Line::Center;});
1168}
1169
1171 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_wind_parameter(j);});
1172}
1173
1175 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_frequency_parameter(j);});
1176}
1177
1179 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return j.is_mag();});
1180}
1181
1183 auto p = std::find_if(js.cbegin(), js.cend(), [](auto& j){return j == Jacobian::Atm::Temperature;});
1184 if (p not_eq js.cend())
1185 return p -> Target().perturbation;
1186 else
1187 return std::numeric_limits<Numeric>::quiet_NaN();
1188}
1189
1191 auto p = std::find_if(js.cbegin(), js.cend(), [](auto& j){return is_frequency_parameter(j);});
1192 if (p not_eq js.cend())
1193 return p -> Target().perturbation;
1194 else
1195 return std::numeric_limits<Numeric>::quiet_NaN();
1196}
1197
1199 auto p = std::find_if(js.cbegin(), js.cend(), [](auto& j){return j.is_mag();});
1200 if (p not_eq js.cend())
1201 return p -> Target().perturbation;
1202 else
1203 return std::numeric_limits<Numeric>::quiet_NaN();
1204}
1205
1206namespace Jacobian {
1207std::ostream& operator<<(std::ostream& os, const Target& x) {
1208 os << x.TargetType() << ' ';
1209 switch (toType(x.TargetType())) {
1210 case Type::Atm:
1211 os << x.atm;
1212 break;
1213 case Type::Line:
1214 os << x.line;
1215 break;
1216 case Type::Sensor:
1217 os << x.sensor;
1218 break;
1219 case Type::Special:
1220 os << x.special;
1221 break;
1222 case Type::FINAL:
1223 os << "BAD STATE";
1224 break;
1225 }
1226 if (x.needQuantumIdentity()) os << ' ' << x.qid;
1227 if (x.needArrayOfSpeciesTag()) os << ' ' << x.species_array_id;
1228 if (x.needString()) os << ' ' << x.string_id;
1229 if (std::isnormal(x.perturbation)) os << ' ' << x.perturbation;
1230
1231 return os;
1232}
1233} // namespace Jacobian
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:136
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:190
The global header file for ARTS.
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Definition: artstime.cc:190
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition: check_input.h:149
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1064
Index nrows() const noexcept
Definition: matpackI.h:1075
Index ncols() const noexcept
Definition: matpackI.h:1076
A constant view of a Tensor3.
Definition: matpackIII.h:130
A constant view of a Vector.
Definition: matpackI.h:530
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:554
The MatrixView class.
Definition: matpackI.h:1182
The Matrix class.
Definition: matpackI.h:1283
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 Matrix & TransformationMatrix() const
Definition: jacobian.h:498
bool HasAffine() const
Definition: jacobian.h:495
const String & Mode() const
Returns the mode.
Definition: jacobian.h:393
const String & TransformationFunc() const
Definition: jacobian.h:496
Index nelem() const
Number of elements in the grids.
Definition: jacobian.h:421
const Vector & TFuncParameters() const
Definition: jacobian.h:497
const Vector & OffsetVector() const
Definition: jacobian.h:499
Jacobian::Target & Target()
Get the Jacobian Target.
Definition: jacobian.h:430
const QuantumIdentifier & QuantumIdentity() const
Returns the identity of this Jacobian.
Definition: jacobian.h:444
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 Tensor3View class.
Definition: matpackIII.h:246
The Tensor3 class.
Definition: matpackIII.h:346
The Vector class.
Definition: matpackI.h:922
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.
void dxdvmrscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
Scale factor for conversion of derivatives with respect to VMR.
Definition: jacobian.cc:978
bool species_iso_match(const RetrievalQuantity &rq, const Species::IsotopeRecord &ir)
Returns if the Retrieval quantity is VMR derivative for all the species in the species tags.
Definition: jacobian.cc:1140
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:564
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
Scale factor for conversion between gas species units.
Definition: jacobian.cc:951
#define VARISLINESHAPEPARAM(X, Y)
Definition: jacobian.cc:1035
bool supports_propmat_clearsky(const ArrayOfRetrievalQuantity &js)
Returns if the array supports propagation matrix derivatives.
Definition: jacobian.cc:1114
void diy_from_pos_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstMatrixView diy_dpos, const Index &atmosphere_dim, ConstVectorView rtp_pos)
diy_from_pos_to_rgrids
Definition: jacobian.cc:474
bool do_line_center_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a line center derivative.
Definition: jacobian.cc:1166
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:605
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition: jacobian.cc:847
bool supports_CIA(const ArrayOfRetrievalQuantity &js)
Returns if the array supports CIA derivatives.
Definition: jacobian.cc:1088
void transform_jacobian(Matrix &jacobian, const Vector x, const ArrayOfRetrievalQuantity &jqs)
Applies both functional and affine transformations.
Definition: jacobian.cc:90
bool is_wind_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a wind parameter in propagation matrix calculations.
Definition: jacobian.cc:1001
bool is_lineshape_parameter_X3(const RetrievalQuantity &t) noexcept
Definition: jacobian.cc:1060
bool is_lineshape_parameter_X1(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a X1 derivative.
Definition: jacobian.cc:1044
bool is_lineshape_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0, D0, G2, D2, ETA, FVC, Y, G, DV derivative.
Definition: jacobian.cc:1076
bool supports_relaxation_matrix(const ArrayOfRetrievalQuantity &js)
Returns if the array supports relaxation matrix derivatives.
Definition: jacobian.cc:1102
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1174
bool is_derived_magnetic_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a derived magnetic parameter.
Definition: jacobian.cc:1009
bool supports_continuum(const ArrayOfRetrievalQuantity &js)
Returns if the array supports continuum derivatives.
Definition: jacobian.cc:1096
bool species_match(const RetrievalQuantity &rq, const ArrayOfSpeciesTag &ast)
Returns if the Retrieval quantity is VMR derivative for all the species in the species tags.
Definition: jacobian.cc:1118
bool do_magnetic_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a magnetic derivative.
Definition: jacobian.cc:1178
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1149
#define ISLINESHAPETYPE(X)
Definition: jacobian.cc:1017
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:596
bool supports_hitran_xsec(const ArrayOfRetrievalQuantity &js)
Returns if the array supports HITRAN cross-section derivatives.
Definition: jacobian.cc:1092
bool supports_lookup(const ArrayOfRetrievalQuantity &js)
Returns if the array supports lookup table derivatives.
Definition: jacobian.cc:1108
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:618
const Numeric PI
bool is_lineshape_parameter_X2(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a X2 derivative.
Definition: jacobian.cc:1052
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
Numeric magnetic_field_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the magnetic field perturbation if it exists.
Definition: jacobian.cc:1198
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Maps jacobian data for points along the propagation path, to jacobian retrieval grid data.
Definition: jacobian.cc:312
bool is_nlte_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a NLTE parameter.
Definition: jacobian.cc:1013
bool do_wind_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a wind-based frequency derivative derivative.
Definition: jacobian.cc:1170
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1005
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
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1182
jacobianVMRcheck do_vmr_jacobian(const ArrayOfRetrievalQuantity &js, const QuantumIdentifier &line_qid) noexcept
Returns the required info for VMR Jacobian.
Definition: jacobian.cc:1153
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
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1190
const Numeric NAT_LOG_TEN
void from_dpath_to_dx(MatrixView diy_dx, ConstMatrixView diy_dq, const Numeric &w)
Definition: jacobian.cc:302
void jacobian_type_extrapol(ArrayOfGridPos &gp)
Adopts grid positions to extrapolation used for jacobians.
Definition: jacobian.cc:835
bool is_lineshape_parameter_X0(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a X0 derivative.
Definition: jacobian.cc:1036
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
Definition: jacobian.cc:155
ostream & operator<<(ostream &os, const RetrievalQuantity &ot)
Definition: jacobian.cc:39
void calcBaselineFit(Vector &y_baseline, const Vector &x, const Index &mblock_index, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const RetrievalQuantity &rq, const Index rq_index, const ArrayOfArrayOfIndex &jacobian_indices)
Calculate baseline fit.
Definition: jacobian.cc:871
bool is_lineshape_parameter_bar_linemixing(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0, D0, G2, D2, ETA, FVC derivative.
Definition: jacobian.cc:1069
bool is_line_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is related to the absorption line.
Definition: jacobian.cc:1084
Routines for setting up the jacobian.
bool is_pressure_broadening_D0(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a D0 derivative.
bool is_pressure_broadening_Y(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a Y derivative.
bool is_pressure_broadening_DV(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a DV derivative.
bool is_pressure_broadening_G0(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0 derivative.
bool is_pressure_broadening_FVC(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a FVC derivative.
bool is_pressure_broadening_D2(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a D2 derivative.
bool is_pressure_broadening_G(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G derivative.
bool is_pressure_broadening_G2(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0 derivative.
bool is_pressure_broadening_ETA(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a ETA derivative.
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition: jacobian.h:549
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:541
#define ns
#define min(a, b)
#define dx
#define max(a, b)
Linear algebra functions.
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
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.
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
const Joker joker
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
std::ostream & operator<<(std::ostream &os, const Target &x)
Definition: jacobian.cc:1207
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
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
constexpr std::string_view Joker
Definition: isotopologues.h:13
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
#define w
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:736
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 p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
Header file for special_interp.cc.
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
String string_id
ID for some of the Special types of partial derivatives.
Definition: jacobian.h:133
bool needArrayOfSpeciesTag() const noexcept
Does this type need the ArrayOfSpeciesTag?
Definition: jacobian.h:306
Line line
Type of line quantity.
Definition: jacobian.h:115
ArrayOfSpeciesTag species_array_id
ID for some of the Special types of partial derivatives.
Definition: jacobian.h:130
std::string_view TargetType() const noexcept
Return type as string.
Definition: jacobian.h:198
Sensor sensor
Type of sensor quantity.
Definition: jacobian.h:118
Atm atm
Type of atm quantity.
Definition: jacobian.h:112
bool needString() const noexcept
Does this type need the String?
Definition: jacobian.h:311
QuantumIdentifier qid
ID for the Line types of partial derivatives.
Definition: jacobian.h:127
bool needQuantumIdentity() const noexcept
Does this type need the QuantumIdentifier?
Definition: jacobian.h:301
Special special
Type of special quantity.
Definition: jacobian.h:121
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath_struct.h:33
A logical struct for global quantum numbers with species identifiers.
Species::Species Species() const noexcept
Species::IsotopeRecord Isotopologue() const noexcept
The Sparse class.
Definition: matpackII.h:72
Struct containing all information needed about one isotope.
Definition: isotopologues.h:16
Species spec
Species type as defined in species.h.
Definition: isotopologues.h:18
std::string_view isotname
A custom name that is unique for this Species type.
Definition: isotopologues.h:21
Deals with whether or not we should do a VMR derivative.
Definition: jacobian.h:1271