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