ARTS  2.2.66
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 
19 
20 
21 /*===========================================================================
22  === File description
23  ===========================================================================*/
24 
37 /*===========================================================================
38  === External declarations
39  ===========================================================================*/
40 
41 #include <cmath>
42 #include <string>
43 #include "absorption.h"
44 #include "arts.h"
45 #include "auto_md.h"
46 #include "check_input.h"
47 #include "math_funcs.h"
48 #include "messages.h"
49 #include "interpolation_poly.h"
50 #include "jacobian.h"
51 #include "physics_funcs.h"
52 #include "rte.h"
53 
54 extern const Numeric PI;
55 
56 extern const String ABSSPECIES_MAINTAG;
57 extern const String FREQUENCY_MAINTAG;
58 extern const String FREQUENCY_SUBTAG_0;
59 extern const String FREQUENCY_SUBTAG_1;
60 extern const String POINTING_MAINTAG;
61 extern const String POINTING_SUBTAG_A;
62 extern const String POINTING_CALCMODE_A;
63 extern const String POINTING_CALCMODE_B;
64 extern const String POLYFIT_MAINTAG;
65 extern const String SINEFIT_MAINTAG;
66 extern const String TEMPERATURE_MAINTAG;
67 extern const String WIND_MAINTAG;
68 
69 
70 
71 /*===========================================================================
72  === The methods, with general methods first followed by the Add/Calc method
73  === pairs for each retrieval quantity.
74  ===========================================================================*/
75 
76 
77 //----------------------------------------------------------------------------
78 // General methods:
79 //----------------------------------------------------------------------------
80 
81 
82 /* Workspace method: Doxygen documentation will be auto-generated */
84  Workspace& ws,
85  Index& jacobian_do,
86  ArrayOfArrayOfIndex& jacobian_indices,
87  Agenda& jacobian_agenda,
88  const ArrayOfRetrievalQuantity& jacobian_quantities,
89  const Matrix& sensor_pos,
90  const Sparse& sensor_response,
91  const Verbosity& verbosity )
92 {
93  // Make sure that the array is not empty
94  if( jacobian_quantities.nelem() == 0 )
95  throw runtime_error(
96  "No retrieval quantities has been added to *jacobian_quantities*." );
97 
98  // Check that sensor_pol and sensor_response has been initialised
99  if( sensor_pos.nrows() == 0 )
100  {
101  ostringstream os;
102  os << "The number of rows in *sensor_pos* is zero, i.e. no measurement\n"
103  << "blocks has been defined. This has to be done before calling\n"
104  << "jacobianClose.";
105  throw runtime_error(os.str());
106  }
107  if( sensor_response.nrows() == 0 )
108  {
109  ostringstream os;
110  os << "The sensor has either to be defined or turned off before calling\n"
111  << "jacobianClose.";
112  throw runtime_error(os.str());
113  }
114 
115  // Loop over retrieval quantities, set JacobianIndices
116  Index ncols = 0;
117  //
118  for( Index it=0; it<jacobian_quantities.nelem(); it++ )
119  {
120  // Store start jacobian index
121  ArrayOfIndex indices(2);
122  indices[0] = ncols;
123 
124  // Count total number of field points, i.e. product of grid lengths
125  Index cols = 1;
126  ArrayOfVector grids = jacobian_quantities[it].Grids();
127  for( Index jt=0; jt<grids.nelem(); jt++ )
128  { cols *= grids[jt].nelem(); }
129 
130  // Store stop index
131  indices[1] = ncols + cols - 1;
132  jacobian_indices.push_back( indices );
133 
134  ncols += cols;
135  }
136 
137  jacobian_agenda.check(ws, verbosity);
138  jacobian_do = 1;
139 }
140 
141 
142 
143 /* Workspace method: Doxygen documentation will be auto-generated */
145  ArrayOfRetrievalQuantity& jacobian_quantities,
146  ArrayOfArrayOfIndex& jacobian_indices,
147  Agenda& jacobian_agenda,
148  const Verbosity& )
149 {
150  jacobian_quantities.resize(0);
151  jacobian_indices.resize(0);
152  jacobian_agenda = Agenda();
153  jacobian_agenda.set_name( "jacobian_agenda" );
154 }
155 
156 
157 
158 /* Workspace method: Doxygen documentation will be auto-generated */
160  Index& jacobian_do,
161  Agenda& jacobian_agenda,
162  ArrayOfRetrievalQuantity& jacobian_quantities,
163  ArrayOfArrayOfIndex& jacobian_indices,
164  const Verbosity& verbosity )
165 {
166  jacobian_do = 0;
167  jacobianInit( jacobian_quantities, jacobian_indices, jacobian_agenda,
168  verbosity );
169 }
170 
171 
172 
173 
174 
175 //----------------------------------------------------------------------------
176 // Absorption species:
177 //----------------------------------------------------------------------------
178 
179 /* Workspace method: Doxygen documentation will be auto-generated */
181  Workspace& ws _U_,
183  Agenda& jacobian_agenda,
184  const Index& atmosphere_dim,
185  const Vector& p_grid,
186  const Vector& lat_grid,
187  const Vector& lon_grid,
188  const Vector& rq_p_grid,
189  const Vector& rq_lat_grid,
190  const Vector& rq_lon_grid,
191  const String& species,
192  const String& method,
193  const String& mode,
194  const Numeric& dx,
195  const Verbosity& verbosity )
196 {
197  CREATE_OUT2;
198  CREATE_OUT3;
199 
200  // Check that this species is not already included in the jacobian.
201  for( Index it=0; it<jq.nelem(); it++ )
202  {
203  if( jq[it].MainTag() == ABSSPECIES_MAINTAG &&
204  jq[it].Subtag() == species )
205  {
206  ostringstream os;
207  os << "The gas species:\n" << species << "\nis already included in "
208  << "*jacobian_quantities*.";
209  throw runtime_error(os.str());
210  }
211  }
212 
213  // Check retrieval grids, here we just check the length of the grids
214  // vs. the atmosphere dimension
215  ArrayOfVector grids(atmosphere_dim);
216  {
217  ostringstream os;
218  if( !check_retrieval_grids( grids, os, p_grid, lat_grid, lon_grid,
219  rq_p_grid, rq_lat_grid, rq_lon_grid,
220  "retrieval pressure grid",
221  "retrieval latitude grid",
222  "retrievallongitude_grid",
223  atmosphere_dim ) )
224  throw runtime_error(os.str());
225  }
226 
227  // Check that method is either "analytical" or "perturbation"
228  bool analytical;
229  if( method == "perturbation" )
230  { analytical = 0; }
231  else if( method == "analytical" )
232  { analytical = 1; }
233  else
234  {
235  ostringstream os;
236  os << "The method for absorption species retrieval can only be "
237  << "\"analytical\"\n or \"perturbation\".";
238  throw runtime_error(os.str());
239  }
240 
241  // Check that mode is either "vmr", "nd" or "rel" with or without prefix log
242  if( mode != "vmr" && mode != "nd" && mode != "rel" && mode != "logrel" )
243  {
244  throw runtime_error( "The retrieval mode can only be \"vmr\", \"nd\" "
245  "\"rel\" or \"logrel\"." );
246  }
247 
248  // If nd, check that not temmperature is retrieved
249  if( mode == "nd" )
250  {
251  for (Index it=0; it<jq.nelem(); it++)
252  {
253  if( jq[it].MainTag() == TEMPERATURE_MAINTAG )
254  {
255  ostringstream os;
256  os <<
257  "Retrieval of temperature and number densities can not be mixed.";
258  throw runtime_error(os.str());
259  }
260  }
261  }
262 
263  // Create the new retrieval quantity
265  rq.MainTag( ABSSPECIES_MAINTAG );
266  rq.Subtag( species );
267  rq.Mode( mode );
268  rq.Analytical( analytical );
269  rq.Perturbation( dx );
270  rq.Grids( grids );
271 
272  // Add it to the *jacobian_quantities*
273  jq.push_back( rq );
274 
275  // Add gas species method to the jacobian agenda
276  if( analytical )
277  {
278  out3 << " Calculations done by semi-analytical expressions.\n";
279  jacobian_agenda.append( "jacobianCalcAbsSpeciesAnalytical", TokVal() );
280  }
281  else
282  {
283  out2 << " Adding absorption species: " << species
284  << " to *jacobian_quantities*\n" << " and *jacobian_agenda*\n";
285  out3 << " Calculations done by perturbation, size " << dx
286  << " " << mode << ".\n";
287 
288  jacobian_agenda.append( "jacobianCalcAbsSpeciesPerturbations", species );
289  }
290 }
291 
292 
293 
294 /* Workspace method: Doxygen documentation will be auto-generated */
296  Matrix& jacobian _U_,
297  const Index& mblock_index _U_,
298  const Vector& iyb _U_,
299  const Vector& yb _U_,
300  const Verbosity& )
301 {
302  /* Nothing to do here for the analytical case, this function just exists
303  to satisfy the required inputs and outputs of the jacobian_agenda */
304 }
305 
306 
307 
308 /* Workspace method: Doxygen documentation will be auto-generated */
310  Workspace& ws,
311  Matrix& jacobian,
312  const Index& mblock_index,
313  const Vector& iyb _U_,
314  const Vector& yb,
315  const Index& atmosphere_dim,
316  const Vector& p_grid,
317  const Vector& lat_grid,
318  const Vector& lon_grid,
319  const Tensor3& t_field,
320  const Tensor3& z_field,
321  const Tensor4& vmr_field,
322  const ArrayOfArrayOfSpeciesTag& abs_species,
323  const Index& cloudbox_on,
324  const Index& stokes_dim,
325  const Vector& f_grid,
326  const Matrix& sensor_pos,
327  const Matrix& sensor_los,
328  const Matrix& transmitter_pos,
329  const Vector& mblock_za_grid,
330  const Vector& mblock_aa_grid,
331  const Index& antenna_dim,
332  const Sparse& sensor_response,
333  const Agenda& iy_main_agenda,
334  const ArrayOfRetrievalQuantity& jacobian_quantities,
335  const ArrayOfArrayOfIndex& jacobian_indices,
336  const String& species,
337  const Verbosity& verbosity)
338 {
339  // Set some useful variables.
341  ArrayOfIndex ji;
342  Index it, pertmode;
343 
344  // Find the retrieval quantity related to this method, i.e. Abs. species -
345  // species. This works since the combined MainTag and Subtag is individual.
346  bool found = false;
347  for( Index n=0; n<jacobian_quantities.nelem() && !found; n++ )
348  {
349  if( jacobian_quantities[n].MainTag() == ABSSPECIES_MAINTAG &&
350  jacobian_quantities[n].Subtag() == species )
351  {
352  found = true;
353  rq = jacobian_quantities[n];
354  ji = jacobian_indices[n];
355  }
356  }
357  if( !found )
358  {
359  ostringstream os;
360  os << "There is no gas species retrieval quantities defined for:\n"
361  << species;
362  throw runtime_error(os.str());
363  }
364 
365  if( rq.Analytical() )
366  {
367  ostringstream os;
368  os << "This WSM handles only perturbation calculations.\n"
369  << "Are you using the method manually?";
370  throw runtime_error(os.str());
371  }
372 
373  // Store the start JacobianIndices and the Grids for this quantity
374  it = ji[0];
375  ArrayOfVector jg = rq.Grids();
376 
377  // Check if a relative pertubation is used or not, this information is needed
378  // by the methods 'perturbation_field_?d'.
379  // Note: both 'vmr' and 'nd' are absolute perturbations
380  if( rq.Mode()=="rel" )
381  pertmode = 0;
382  else
383  pertmode = 1;
384 
385  // For each atmospheric dimension option calculate a ArrayOfGridPos, these
386  // are the base functions for interpolating the perturbations into the
387  // atmospheric grids.
388  ArrayOfGridPos p_gp, lat_gp, lon_gp;
389  Index j_p = jg[0].nelem();
390  Index j_lat = 1;
391  Index j_lon = 1;
392  //
393  get_perturbation_gridpos( p_gp, p_grid, jg[0], true );
394  //
395  if( atmosphere_dim >= 2 )
396  {
397  j_lat = jg[1].nelem();
398  get_perturbation_gridpos( lat_gp, lat_grid, jg[1], false );
399  if( atmosphere_dim == 3 )
400  {
401  j_lon = jg[2].nelem();
402  get_perturbation_gridpos( lon_gp, lon_grid, jg[2], false );
403  }
404  }
405 
406  // Find VMR field for this species.
407  ArrayOfSpeciesTag tags;
408  array_species_tag_from_string( tags, species );
409  Index si = chk_contains( "species", abs_species, tags );
410 
411  // Variables for vmr field perturbation unit conversion
412  Tensor3 nd_field(0,0,0);
413  if( rq.Mode()=="nd" )
414  {
415  nd_field.resize( t_field.npages(), t_field.nrows(), t_field.ncols() );
416  calc_nd_field( nd_field, p_grid, t_field );
417  }
418 
419 
420  // Loop through the retrieval grid and calculate perturbation effect
421  //
422  const Index n1y = sensor_response.nrows();
423  Vector dy( n1y );
424  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
425  //
426  for( Index lon_it=0; lon_it<j_lon; lon_it++ )
427  {
428  for( Index lat_it=0; lat_it<j_lat; lat_it++ )
429  {
430  for (Index p_it=0; p_it<j_p; p_it++)
431  {
432  // Here we calculate the ranges of the perturbation. We want the
433  // perturbation to continue outside the atmospheric grids for the
434  // edge values.
435  Range p_range = Range(0,0);
436  Range lat_range = Range(0,0);
437  Range lon_range = Range(0,0);
438 
439  get_perturbation_range( p_range, p_it, j_p );
440 
441  if( atmosphere_dim>=2 )
442  {
443  get_perturbation_range( lat_range, lat_it, j_lat );
444  if( atmosphere_dim == 3 )
445  {
446  get_perturbation_range( lon_range, lon_it, j_lon );
447  }
448  }
449 
450  // Create VMR field to perturb
451  Tensor4 vmr_p = vmr_field;
452 
453  // If perturbation given in ND convert the vmr-field to ND before
454  // the perturbation is added
455  if( rq.Mode() == "nd" )
456  vmr_p(si,joker,joker,joker) *= nd_field;
457 
458  // Calculate the perturbed field according to atmosphere_dim,
459  // the number of perturbations is the length of the retrieval
460  // grid +2 (for the end points)
461  switch (atmosphere_dim)
462  {
463  case 1:
464  {
465  // Here we perturb a vector
466  perturbation_field_1d( vmr_p(si,joker,lat_it,lon_it),
467  p_gp, jg[0].nelem()+2, p_range,
468  rq.Perturbation(), pertmode );
469  break;
470  }
471  case 2:
472  {
473  // Here we perturb a matrix
474  perturbation_field_2d( vmr_p(si,joker,joker,lon_it),
475  p_gp, lat_gp, jg[0].nelem()+2,
476  jg[1].nelem()+2, p_range, lat_range,
477  rq.Perturbation(), pertmode );
478  break;
479  }
480  case 3:
481  {
482  // Here we need to perturb a tensor3
484  p_gp, lat_gp, lon_gp,
485  jg[0].nelem()+2,
486  jg[1].nelem()+2, jg[2].nelem()+2,
487  p_range, lat_range, lon_range,
488  rq.Perturbation(), pertmode );
489  break;
490  }
491  }
492 
493  // If perturbation given in ND convert back to VMR
494  if (rq.Mode()=="nd")
495  vmr_p(si,joker,joker,joker) /= nd_field;
496 
497  // Calculate the perturbed spectrum
498  //
499  Vector iybp;
500  ArrayOfVector dummy3;
501  ArrayOfMatrix dummy4;
502  //
503  iyb_calc( ws, iybp, dummy3, dummy4, mblock_index,
504  atmosphere_dim, t_field, z_field, vmr_p, cloudbox_on,
505  stokes_dim, f_grid, sensor_pos, sensor_los,
506  transmitter_pos, mblock_za_grid,
507  mblock_aa_grid, antenna_dim, iy_main_agenda,
509  ArrayOfArrayOfIndex(), ArrayOfString(), verbosity );
510  //
511  mult( dy, sensor_response, iybp );
512 
513  // Difference spectrum
514  for( Index i=0; i<n1y; i++ )
515  { dy[i] = ( dy[i]- yb[i] ) / rq.Perturbation(); }
516 
517  // Put into jacobian
518  jacobian(rowind,it) = dy;
519 
520  // Result from next loop shall go into next column of J
521  it++;
522  }
523  }
524  }
525 }
526 
527 
528 
529 
530 
531 //----------------------------------------------------------------------------
532 // Frequency shift
533 //----------------------------------------------------------------------------
534 
535 /* Workspace method: Doxygen documentation will be auto-generated */
537  Workspace& ws _U_,
538  ArrayOfRetrievalQuantity& jacobian_quantities,
539  Agenda& jacobian_agenda,
540  const Vector& f_grid,
541  const Matrix& sensor_pos,
542  const Vector& sensor_time,
543  const Index& poly_order,
544  const Numeric& df,
545  const Verbosity& )
546 {
547  // Check that poly_order is -1 or positive
548  if( poly_order < -1 )
549  throw runtime_error(
550  "The polynomial order has to be positive or -1 for gitter." );
551 
552  // Check that this jacobian type is not already included.
553  for( Index it=0; it<jacobian_quantities.nelem(); it++ )
554  {
555  if (jacobian_quantities[it].MainTag()== FREQUENCY_MAINTAG &&
556  jacobian_quantities[it].Subtag() == FREQUENCY_SUBTAG_0 )
557  {
558  ostringstream os;
559  os << "Fit of frequency shift is already included in\n"
560  << "*jacobian_quantities*.";
561  throw runtime_error(os.str());
562  }
563  }
564 
565  // Checks of df
566  if( df <= 0 )
567  throw runtime_error( "The argument *df* must be > 0." );
568  if( df > 1e6 )
569  throw runtime_error( "The argument *df* is not allowed to exceed 1 MHz." );
570  const Index nf = f_grid.nelem();
571  const Numeric maxdf = f_grid[nf-1] - f_grid[nf-2];
572  if( df > maxdf )
573  {
574  ostringstream os;
575  os << "The value of *df* is too big with respect to spacing of "
576  << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
577  << "between the two last elements of *f_grid*.\n"
578  << "This spacing is : " <<maxdf/1e3 << " kHz\n"
579  << "The value of df is: " << df/1e3 << " kHz";
580  throw runtime_error(os.str());
581  }
582 
583  // Check that sensor_time is consistent with sensor_pos
584  if( sensor_time.nelem() != sensor_pos.nrows() )
585  {
586  ostringstream os;
587  os << "The WSV *sensor_time* must be defined for every "
588  << "measurement block.\n";
589  throw runtime_error(os.str());
590  }
591 
592  // Do not allow that *poly_order* is not too large compared to *sensor_time*
593  if( poly_order > sensor_time.nelem()-1 )
594  { throw runtime_error(
595  "The polynomial order can not be >= length of *sensor_time*." ); }
596 
597  // Create the new retrieval quantity
599  rq.MainTag( FREQUENCY_MAINTAG );
600  rq.Subtag( FREQUENCY_SUBTAG_0 );
601  rq.Mode( "" );
602  rq.Analytical( 0 );
603  rq.Perturbation( df );
604 
605  // To store the value or the polynomial order, create a vector with length
606  // poly_order+1, in case of gitter set the size of the grid vector to be the
607  // number of measurement blocks, all elements set to -1.
608  Vector grid(0,poly_order+1,1);
609  if( poly_order == -1 )
610  {
611  grid.resize(sensor_pos.nrows());
612  grid = -1.0;
613  }
614  ArrayOfVector grids(1,grid);
615  rq.Grids(grids);
616 
617  // Add it to the *jacobian_quantities*
618  jacobian_quantities.push_back( rq );
619 
620  // Add corresponding calculation method to the jacobian agenda
621  jacobian_agenda.append( "jacobianCalcFreqShift", "" );
622 }
623 
624 
625 
626 /* Workspace method: Doxygen documentation will be auto-generated */
628  Matrix& jacobian,
629  const Index& mblock_index,
630  const Vector& iyb,
631  const Vector& yb,
632  const Index& stokes_dim,
633  const Vector& f_grid,
634  const Matrix& sensor_los,
635  const Vector& mblock_za_grid,
636  const Vector& mblock_aa_grid,
637  const Index& antenna_dim,
638  const Sparse& sensor_response,
639  const Vector& sensor_time,
640  const ArrayOfRetrievalQuantity& jacobian_quantities,
641  const ArrayOfArrayOfIndex& jacobian_indices,
642  const Verbosity& )
643 {
644  // Set some useful (and needed) variables.
646  ArrayOfIndex ji;
647 
648  // Find the retrieval quantity related to this method.
649  // This works since the combined MainTag and Subtag is individual.
650  bool found = false;
651  for( Index n=0; n<jacobian_quantities.nelem() && !found; n++ )
652  {
653  if( jacobian_quantities[n].MainTag() == FREQUENCY_MAINTAG &&
654  jacobian_quantities[n].Subtag() == FREQUENCY_SUBTAG_0 )
655  {
656  found = true;
657  rq = jacobian_quantities[n];
658  ji = jacobian_indices[n];
659  }
660  }
661  if( !found )
662  {
663  throw runtime_error(
664  "There is no such frequency retrieval quantity defined.\n" );
665  }
666 
667  // Check that sensor_response is consistent with yb and iyb
668  //
669  if( sensor_response.nrows() != yb.nelem() )
670  throw runtime_error(
671  "Mismatch in size between *sensor_response* and *yb*." );
672  if( sensor_response.ncols() != iyb.nelem() )
673  throw runtime_error(
674  "Mismatch in size between *sensor_response* and *iyb*." );
675 
676  // Get disturbed (part of) y
677  //
678  const Index n1y = sensor_response.nrows();
679  Vector dy( n1y );
680  {
681  const Index nf2 = f_grid.nelem();
682  const Index nza2 = mblock_za_grid.nelem();
683  Index naa2 = mblock_aa_grid.nelem();
684  if( antenna_dim == 1 )
685  { naa2 = 1; }
686  const Index niyb = nf2 * nza2 * naa2 * stokes_dim;
687 
688  // Interpolation weights
689  //
690  const Index porder = 3;
691  //
692  ArrayOfGridPosPoly gp( nf2 );
693  Matrix itw( nf2, porder+1) ;
694  Vector fg_new = f_grid, iyb2(niyb);
695  //
696  fg_new += rq.Perturbation();
697  gridpos_poly( gp, f_grid, fg_new, porder, 1.0 );
698  interpweights( itw, gp );
699 
700  // Do interpolation
701  for( Index iza=0; iza<nza2; iza++ )
702  {
703  for( Index iaa=0; iaa<naa2; iaa++ )
704  {
705  const Index row0 =( iza*naa2 + iaa ) * nf2 * stokes_dim;
706 
707  for( Index is=0; is<stokes_dim; is++ )
708  {
709  interp( iyb2[Range(row0+is,nf2,stokes_dim)], itw,
710  iyb[Range(row0+is,nf2,stokes_dim)], gp );
711  }
712  }
713  }
714 
715  // Determine difference
716  //
717  mult( dy, sensor_response, iyb2 );
718  //
719  for( Index i=0; i<n1y; i++ )
720  { dy[i] = ( dy[i]- yb[i] ) / rq.Perturbation(); }
721  }
722 
723  //--- Create jacobians ---
724 
725  const Index lg = rq.Grids()[0].nelem();
726  const Index it = ji[0];
727  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
728  const Index row0 = rowind.get_start();
729 
730  // Handle gitter seperately
731  if( rq.Grids()[0][0] == -1 ) // Not all values are set here,
732  { // but should already have been
733  assert( lg == sensor_los.nrows() ); // set to 0
734  assert( rq.Grids()[0][mblock_index] == -1 );
735  jacobian(rowind,it+mblock_index) = dy;
736  }
737 
738  // Polynomial representation
739  else
740  {
741  Vector w;
742  for( Index c=0; c<lg; c++ )
743  {
744  assert( Numeric(c) == rq.Grids()[0][c] );
745  //
746  polynomial_basis_func( w, sensor_time, c );
747  //
748  for( Index i=0; i<n1y; i++ )
749  { jacobian(row0+i,it+c) = w[mblock_index] * dy[i]; }
750  }
751  }
752 }
753 
754 
755 
756 
757 //----------------------------------------------------------------------------
758 // Frequency stretch
759 //----------------------------------------------------------------------------
760 
761 /* Workspace method: Doxygen documentation will be auto-generated */
763  Workspace& ws _U_,
764  ArrayOfRetrievalQuantity& jacobian_quantities,
765  Agenda& jacobian_agenda,
766  const Vector& f_grid,
767  const Matrix& sensor_pos,
768  const Vector& sensor_time,
769  const Index& poly_order,
770  const Numeric& df,
771  const Verbosity& )
772 {
773  // Check that poly_order is -1 or positive
774  if( poly_order < -1 )
775  throw runtime_error(
776  "The polynomial order has to be positive or -1 for gitter." );
777 
778  // Check that this jacobian type is not already included.
779  for( Index it=0; it<jacobian_quantities.nelem(); it++ )
780  {
781  if (jacobian_quantities[it].MainTag()== FREQUENCY_MAINTAG &&
782  jacobian_quantities[it].Subtag() == FREQUENCY_SUBTAG_1 )
783  {
784  ostringstream os;
785  os << "Fit of frequency stretch is already included in\n"
786  << "*jacobian_quantities*.";
787  throw runtime_error(os.str());
788  }
789  }
790 
791  // Checks of df
792  if( df <= 0 )
793  throw runtime_error( "The argument *df* must be > 0." );
794  if( df > 1e6 )
795  throw runtime_error( "The argument *df* is not allowed to exceed 1 MHz." );
796  const Index nf = f_grid.nelem();
797  const Numeric maxdf = f_grid[nf-1] - f_grid[nf-2];
798  if( df > maxdf )
799  {
800  ostringstream os;
801  os << "The value of *df* is too big with respect to spacing of "
802  << "*f_grid*. The maximum\nallowed value of *df* is the spacing "
803  << "between the two last elements of *f_grid*.\n"
804  << "This spacing is : " <<maxdf/1e3 << " kHz\n"
805  << "The value of df is: " << df/1e3 << " kHz";
806  throw runtime_error(os.str());
807  }
808 
809  // Check that sensor_time is consistent with sensor_pos
810  if( sensor_time.nelem() != sensor_pos.nrows() )
811  {
812  ostringstream os;
813  os << "The WSV *sensor_time* must be defined for every "
814  << "measurement block.\n";
815  throw runtime_error(os.str());
816  }
817 
818  // Do not allow that *poly_order* is not too large compared to *sensor_time*
819  if( poly_order > sensor_time.nelem()-1 )
820  { throw runtime_error(
821  "The polynomial order can not be >= length of *sensor_time*." ); }
822 
823  // Create the new retrieval quantity
825  rq.MainTag( FREQUENCY_MAINTAG );
826  rq.Subtag( FREQUENCY_SUBTAG_1 );
827  rq.Mode( "" );
828  rq.Analytical( 0 );
829  rq.Perturbation( df );
830 
831  // To store the value or the polynomial order, create a vector with length
832  // poly_order+1, in case of gitter set the size of the grid vector to be the
833  // number of measurement blocks, all elements set to -1.
834  Vector grid(0,poly_order+1,1);
835  if( poly_order == -1 )
836  {
837  grid.resize(sensor_pos.nrows());
838  grid = -1.0;
839  }
840  ArrayOfVector grids(1,grid);
841  rq.Grids(grids);
842 
843  // Add it to the *jacobian_quantities*
844  jacobian_quantities.push_back( rq );
845 
846  // Add corresponding calculation method to the jacobian agenda
847  jacobian_agenda.append( "jacobianCalcFreqStretch", "" );
848 }
849 
850 
851 
852 /* Workspace method: Doxygen documentation will be auto-generated */
854  Matrix& jacobian,
855  const Index& mblock_index,
856  const Vector& iyb,
857  const Vector& yb,
858  const Index& stokes_dim,
859  const Vector& f_grid,
860  const Matrix& sensor_los,
861  const Vector& mblock_za_grid,
862  const Vector& mblock_aa_grid,
863  const Index& antenna_dim,
864  const Sparse& sensor_response,
865  const ArrayOfIndex& sensor_response_pol_grid,
866  const Vector& sensor_response_f_grid,
867  const Vector& sensor_response_za_grid,
868  const Vector& sensor_time,
869  const ArrayOfRetrievalQuantity& jacobian_quantities,
870  const ArrayOfArrayOfIndex& jacobian_indices,
871  const Verbosity& )
872 {
873  // The code here is close to identical to the one for Shift. The main
874  // difference is that dy is weighted with poly_order 1 basis function.
875 
876  // Set some useful (and needed) variables.
878  ArrayOfIndex ji;
879 
880  // Find the retrieval quantity related to this method.
881  // This works since the combined MainTag and Subtag is individual.
882  bool found = false;
883  for( Index n=0; n<jacobian_quantities.nelem() && !found; n++ )
884  {
885  if( jacobian_quantities[n].MainTag() == FREQUENCY_MAINTAG &&
886  jacobian_quantities[n].Subtag() == FREQUENCY_SUBTAG_1 )
887  {
888  found = true;
889  rq = jacobian_quantities[n];
890  ji = jacobian_indices[n];
891  }
892  }
893  if( !found )
894  {
895  throw runtime_error(
896  "There is no such frequency retrieval quantity defined.\n" );
897  }
898 
899  // Check that sensor_response is consistent with yb and iyb
900  //
901  if( sensor_response.nrows() != yb.nelem() )
902  throw runtime_error(
903  "Mismatch in size between *sensor_response* and *yb*." );
904  if( sensor_response.ncols() != iyb.nelem() )
905  throw runtime_error(
906  "Mismatch in size between *sensor_response* and *iyb*." );
907 
908  // Get disturbed (part of) y
909  //
910  const Index n1y = sensor_response.nrows();
911  Vector dy( n1y );
912  {
913  const Index nf2 = f_grid.nelem();
914  const Index nza2 = mblock_za_grid.nelem();
915  Index naa2 = mblock_aa_grid.nelem();
916  if( antenna_dim == 1 )
917  { naa2 = 1; }
918  const Index niyb = nf2 * nza2 * naa2 * stokes_dim;
919 
920  // Interpolation weights
921  //
922  const Index porder = 3;
923  //
924  ArrayOfGridPosPoly gp( nf2 );
925  Matrix itw( nf2, porder+1) ;
926  Vector fg_new = f_grid, iyb2(niyb);
927  //
928  fg_new += rq.Perturbation();
929  gridpos_poly( gp, f_grid, fg_new, porder, 1.0 );
930  interpweights( itw, gp );
931 
932  // Do interpolation
933  for( Index iza=0; iza<nza2; iza++ )
934  {
935  for( Index iaa=0; iaa<naa2; iaa++ )
936  {
937  const Index row0 =( iza*naa2 + iaa ) * nf2 * stokes_dim;
938 
939  for( Index is=0; is<stokes_dim; is++ )
940  {
941  interp( iyb2[Range(row0+is,nf2,stokes_dim)], itw,
942  iyb[Range(row0+is,nf2,stokes_dim)], gp );
943  }
944  }
945  }
946 
947  // Determine difference
948  //
949  mult( dy, sensor_response, iyb2 );
950  //
951  for( Index i=0; i<n1y; i++ )
952  { dy[i] = ( dy[i]- yb[i] ) / rq.Perturbation(); }
953 
954  // dy above corresponds now to shift. Convert to stretch:
955  //
956  Vector w;
957  polynomial_basis_func( w, sensor_response_f_grid, 1 );
958  //
959  const Index nf = sensor_response_f_grid.nelem();
960  const Index npol = sensor_response_pol_grid.nelem();
961  const Index nza = sensor_response_za_grid.nelem();
962  //
963  for( Index l=0; l<nza; l++ )
964  {
965  for( Index f=0; f<nf; f++ )
966  {
967  const Index row1 = (l*nf + f)*npol;
968  for( Index p=0; p<npol; p++ )
969  { dy[row1+p] *= w[f]; }
970  }
971  }
972  }
973 
974  //--- Create jacobians ---
975 
976  const Index lg = rq.Grids()[0].nelem();
977  const Index it = ji[0];
978  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
979  const Index row0 = rowind.get_start();
980 
981  // Handle gitter seperately
982  if( rq.Grids()[0][0] == -1 ) // Not all values are set here,
983  { // but should already have been
984  assert( lg == sensor_los.nrows() ); // set to 0
985  assert( rq.Grids()[0][mblock_index] == -1 );
986  jacobian(rowind,it+mblock_index) = dy;
987  }
988 
989  // Polynomial representation
990  else
991  {
992  Vector w;
993  for( Index c=0; c<lg; c++ )
994  {
995  assert( Numeric(c) == rq.Grids()[0][c] );
996  //
997  polynomial_basis_func( w, sensor_time, c );
998  //
999  for( Index i=0; i<n1y; i++ )
1000  { jacobian(row0+i,it+c) = w[mblock_index] * dy[i]; }
1001  }
1002  }
1003 }
1004 
1005 
1006 
1007 
1008 
1009 //----------------------------------------------------------------------------
1010 // Pointing:
1011 //----------------------------------------------------------------------------
1012 
1013 /* Workspace method: Doxygen documentation will be auto-generated */
1015  Workspace& ws _U_,
1016  ArrayOfRetrievalQuantity& jacobian_quantities,
1017  Agenda& jacobian_agenda,
1018  const Matrix& sensor_pos,
1019  const Vector& sensor_time,
1020  const Index& poly_order,
1021  const String& calcmode,
1022  const Numeric& dza,
1023  const Verbosity& )
1024 {
1025  // Check that poly_order is -1 or positive
1026  if( poly_order < -1 )
1027  throw runtime_error(
1028  "The polynomial order has to be positive or -1 for gitter." );
1029 
1030  // Check that this jacobian type is not already included.
1031  for( Index it=0; it<jacobian_quantities.nelem(); it++ )
1032  {
1033  if (jacobian_quantities[it].MainTag()== POINTING_MAINTAG &&
1034  jacobian_quantities[it].Subtag() == POINTING_SUBTAG_A )
1035  {
1036  ostringstream os;
1037  os << "Fit of zenith angle pointing off-set is already included in\n"
1038  << "*jacobian_quantities*.";
1039  throw runtime_error(os.str());
1040  }
1041  }
1042 
1043  // Checks of dza
1044  if( dza <= 0 )
1045  throw runtime_error( "The argument *dza* must be > 0." );
1046  if( dza > 0.1 )
1047  throw runtime_error(
1048  "The argument *dza* is not allowed to exceed 0.1 deg." );
1049 
1050  // Check that sensor_time is consistent with sensor_pos
1051  if( sensor_time.nelem() != sensor_pos.nrows() )
1052  {
1053  ostringstream os;
1054  os << "The WSV *sensor_time* must be defined for every "
1055  << "measurement block.\n";
1056  throw runtime_error(os.str());
1057  }
1058 
1059  // Do not allow that *poly_order* is not too large compared to *sensor_time*
1060  if( poly_order > sensor_time.nelem()-1 )
1061  { throw runtime_error(
1062  "The polynomial order can not be >= length of *sensor_time*." ); }
1063 
1064  // Create the new retrieval quantity
1065  RetrievalQuantity rq;
1066  rq.MainTag( POINTING_MAINTAG );
1067  rq.Subtag( POINTING_SUBTAG_A );
1068  rq.Analytical( 0 );
1069  rq.Perturbation( dza );
1070 
1071 
1072  // To store the value or the polynomial order, create a vector with length
1073  // poly_order+1, in case of gitter set the size of the grid vector to be the
1074  // number of measurement blocks, all elements set to -1.
1075  Vector grid(0,poly_order+1,1);
1076  if( poly_order == -1 )
1077  {
1078  grid.resize(sensor_pos.nrows());
1079  grid = -1.0;
1080  }
1081  ArrayOfVector grids(1,grid);
1082  rq.Grids(grids);
1083 
1084  if( calcmode == "recalc" )
1085  {
1086  rq.Mode( POINTING_CALCMODE_A );
1087  jacobian_agenda.append( "jacobianCalcPointingZaRecalc", "" );
1088  }
1089  else if( calcmode == "interp" )
1090  {
1091  rq.Mode( POINTING_CALCMODE_B );
1092  jacobian_agenda.append( "jacobianCalcPointingZaInterp", "" );
1093  }
1094  else
1095  throw runtime_error(
1096  "Possible choices for *calcmode* are \"recalc\" and \"interp\"." );
1097 
1098  // Add it to the *jacobian_quantities*
1099  jacobian_quantities.push_back( rq );
1100 }
1101 
1102 
1103 
1104 /* Workspace method: Doxygen documentation will be auto-generated */
1106  Matrix& jacobian,
1107  const Index& mblock_index,
1108  const Vector& iyb,
1109  const Vector& yb _U_,
1110  const Index& stokes_dim,
1111  const Vector& f_grid,
1112  const Matrix& sensor_los,
1113  const Vector& mblock_za_grid,
1114  const Vector& mblock_aa_grid,
1115  const Index& antenna_dim,
1116  const Sparse& sensor_response,
1117  const Vector& sensor_time,
1118  const ArrayOfRetrievalQuantity& jacobian_quantities,
1119  const ArrayOfArrayOfIndex& jacobian_indices,
1120  const Verbosity& )
1121 {
1122  if( mblock_za_grid.nelem() < 2 )
1123  throw runtime_error( "The method demands that *mblock_za_grid* has a "
1124  "length of > 1." );
1125 
1126  // Set some useful variables.
1127  RetrievalQuantity rq;
1128  ArrayOfIndex ji;
1129 
1130  // Find the retrieval quantity related to this method.
1131  // This works since the combined MainTag and Subtag is individual.
1132  bool found = false;
1133  for( Index n=0; n<jacobian_quantities.nelem() && !found; n++ )
1134  {
1135  if( jacobian_quantities[n].MainTag() == POINTING_MAINTAG &&
1136  jacobian_quantities[n].Subtag() == POINTING_SUBTAG_A &&
1137  jacobian_quantities[n].Mode() == POINTING_CALCMODE_B )
1138  {
1139  found = true;
1140  rq = jacobian_quantities[n];
1141  ji = jacobian_indices[n];
1142  }
1143  }
1144  if( !found )
1145  { throw runtime_error(
1146  "There is no such pointing retrieval quantity defined.\n" );
1147  }
1148 
1149 
1150  // Get "dy", by inter/extra-polation of existing iyb
1151  //
1152  const Index n1y = sensor_response.nrows();
1153  Vector dy( n1y );
1154  {
1155  // Sizes
1156  const Index nf = f_grid.nelem();
1157  const Index nza = mblock_za_grid.nelem();
1158  Index naa = mblock_aa_grid.nelem();
1159  if( antenna_dim == 1 )
1160  { naa = 1; }
1161 
1162  // Shifted zenith angles
1163  Vector za1 = mblock_za_grid; za1 -= rq.Perturbation();
1164  Vector za2 = mblock_za_grid; za2 += rq.Perturbation();
1165 
1166  // Find interpolation weights
1167  ArrayOfGridPos gp1(nza), gp2(nza);
1168  gridpos( gp1, mblock_za_grid, za1, 1e6 ); // Note huge extrapolation!
1169  gridpos( gp2, mblock_za_grid, za2, 1e6 ); // Note huge extrapolation!
1170  Matrix itw1(nza,2), itw2(nza,2);
1171  interpweights( itw1, gp1 );
1172  interpweights( itw2, gp2 );
1173 
1174  // Make interpolation (for all azimuth angles, frequencies and Stokes)
1175  //
1176  Vector iyb1(iyb.nelem()), iyb2(iyb.nelem());
1177  //
1178  for( Index iaa=0; iaa<naa; iaa++ )
1179  {
1180  for( Index iv=0; iv<nf; iv++ )
1181  {
1182  for( Index is=0; is<stokes_dim; is++ )
1183  {
1184  const Range r( iaa*nza*nf*stokes_dim+iv*stokes_dim+is,
1185  nza, nf*stokes_dim );
1186  interp( iyb1[r], itw1, iyb[r], gp1 );
1187  interp( iyb2[r], itw2, iyb[r], gp2 );
1188  }
1189  }
1190  }
1191 
1192  // Apply sensor and take difference
1193  //
1194  Vector y1(n1y), y2(n1y);
1195  mult( y1, sensor_response, iyb1 );
1196  mult( y2, sensor_response, iyb2 );
1197  //
1198  for( Index i=0; i<n1y; i++ )
1199  { dy[i] = ( y2[i]- y1[i] ) / ( 2* rq.Perturbation() ); }
1200  }
1201 
1202  //--- Create jacobians ---
1203 
1204  const Index lg = rq.Grids()[0].nelem();
1205  const Index it = ji[0];
1206  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
1207  const Index row0 = rowind.get_start();
1208 
1209  // Handle gitter seperately
1210  if( rq.Grids()[0][0] == -1 ) // Not all values are set here,
1211  { // but should already have been
1212  assert( lg == sensor_los.nrows() ); // set to 0
1213  assert( rq.Grids()[0][mblock_index] == -1 );
1214  jacobian(rowind,it+mblock_index) = dy;
1215  }
1216 
1217  // Polynomial representation
1218  else
1219  {
1220  Vector w;
1221  for( Index c=0; c<lg; c++ )
1222  {
1223  assert( Numeric(c) == rq.Grids()[0][c] );
1224  //
1225  polynomial_basis_func( w, sensor_time, c );
1226  //
1227  for( Index i=0; i<n1y; i++ )
1228  { jacobian(row0+i,it+c) = w[mblock_index] * dy[i]; }
1229  }
1230  }
1231 }
1232 
1233 
1234 
1235 /* Workspace method: Doxygen documentation will be auto-generated */
1237  Workspace& ws,
1238  Matrix& jacobian,
1239  const Index& mblock_index,
1240  const Vector& iyb _U_,
1241  const Vector& yb,
1242  const Index& atmosphere_dim,
1243  const Tensor3& t_field,
1244  const Tensor3& z_field,
1245  const Tensor4& vmr_field,
1246  const Index& cloudbox_on,
1247  const Index& stokes_dim,
1248  const Vector& f_grid,
1249  const Matrix& sensor_pos,
1250  const Matrix& sensor_los,
1251  const Matrix& transmitter_pos,
1252  const Vector& mblock_za_grid,
1253  const Vector& mblock_aa_grid,
1254  const Index& antenna_dim,
1255  const Sparse& sensor_response,
1256  const Vector& sensor_time,
1257  const Agenda& iy_main_agenda,
1258  const ArrayOfRetrievalQuantity& jacobian_quantities,
1259  const ArrayOfArrayOfIndex& jacobian_indices,
1260  const Verbosity& verbosity )
1261 {
1262  // Set some useful variables.
1263  RetrievalQuantity rq;
1264  ArrayOfIndex ji;
1265 
1266  // Find the retrieval quantity related to this method.
1267  // This works since the combined MainTag and Subtag is individual.
1268  bool found = false;
1269  for( Index n=0; n<jacobian_quantities.nelem() && !found; n++ )
1270  {
1271  if( jacobian_quantities[n].MainTag() == POINTING_MAINTAG &&
1272  jacobian_quantities[n].Subtag() == POINTING_SUBTAG_A &&
1273  jacobian_quantities[n].Mode() == POINTING_CALCMODE_A )
1274  {
1275  found = true;
1276  rq = jacobian_quantities[n];
1277  ji = jacobian_indices[n];
1278  }
1279  }
1280  if( !found )
1281  { throw runtime_error(
1282  "There is no such pointing retrieval quantity defined.\n" );
1283  }
1284 
1285 
1286  // Get "dy", by calling iyb_calc with shifted sensor_los.
1287  //
1288  const Index n1y = sensor_response.nrows();
1289  Vector dy( n1y );
1290  {
1291  Vector iyb2;
1292  Matrix los = sensor_los;
1293  ArrayOfVector iyb_aux;
1294  ArrayOfMatrix diyb_dx;
1295 
1296  los(joker,0) += rq.Perturbation();
1297 
1298  iyb_calc( ws, iyb2, iyb_aux, diyb_dx, mblock_index,
1299  atmosphere_dim,
1300  t_field, z_field, vmr_field, cloudbox_on, stokes_dim,
1301  f_grid, sensor_pos, los, transmitter_pos, mblock_za_grid,
1302  mblock_aa_grid, antenna_dim, iy_main_agenda,
1304  ArrayOfString(), verbosity );
1305 
1306  // Apply sensor and take difference
1307  //
1308  mult( dy, sensor_response, iyb2 );
1309  //
1310  for( Index i=0; i<n1y; i++ )
1311  { dy[i] = ( dy[i]- yb[i] ) / rq.Perturbation(); }
1312  }
1313 
1314  //--- Create jacobians ---
1315 
1316  const Index lg = rq.Grids()[0].nelem();
1317  const Index it = ji[0];
1318  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
1319  const Index row0 = rowind.get_start();
1320 
1321  // Handle gitter seperately
1322  if( rq.Grids()[0][0] == -1 ) // Not all values are set here,
1323  { // but should already have been
1324  assert( lg == sensor_los.nrows() ); // set to 0
1325  assert( rq.Grids()[0][mblock_index] == -1 );
1326  jacobian(rowind,it+mblock_index) = dy;
1327  }
1328 
1329  // Polynomial representation
1330  else
1331  {
1332  Vector w;
1333  for( Index c=0; c<lg; c++ )
1334  {
1335  assert( Numeric(c) == rq.Grids()[0][c] );
1336  //
1337  polynomial_basis_func( w, sensor_time, c );
1338  //
1339  for( Index i=0; i<n1y; i++ )
1340  { jacobian(row0+i,it+c) = w[mblock_index] * dy[i]; }
1341  }
1342  }
1343 }
1344 
1345 
1346 
1347 
1348 
1349 //----------------------------------------------------------------------------
1350 // Polynomial baseline fits:
1351 //----------------------------------------------------------------------------
1352 
1353 /* Workspace method: Doxygen documentation will be auto-generated */
1355  Workspace& ws _U_,
1357  Agenda& jacobian_agenda,
1358  const ArrayOfIndex& sensor_response_pol_grid,
1359  const Vector& sensor_response_za_grid,
1360  const Matrix& sensor_pos,
1361  const Index& poly_order,
1362  const Index& no_pol_variation,
1363  const Index& no_los_variation,
1364  const Index& no_mblock_variation,
1365  const Verbosity& )
1366 {
1367  // Check that poly_order is >= 0
1368  if( poly_order < 0 )
1369  throw runtime_error( "The polynomial order has to be >= 0.");
1370 
1371  // Check that polyfit is not already included in the jacobian.
1372  for( Index it=0; it<jq.nelem(); it++ )
1373  {
1374  if( jq[it].MainTag() == POLYFIT_MAINTAG )
1375  {
1376  ostringstream os;
1377  os << "Sinusoidal baseline fit is already included in\n"
1378  << "*jacobian_quantities*.";
1379  throw runtime_error(os.str());
1380  }
1381  }
1382 
1383  // "Grids"
1384  //
1385  // Grid dimensions correspond here to
1386  // 1: polynomial order
1387  // 2: polarisation
1388  // 3: viewing direction
1389  // 4: measurement block
1390  //
1391  ArrayOfVector grids(4);
1392  //
1393  if( no_pol_variation )
1394  grids[1] = Vector(1,1);
1395  else
1396  grids[1] = Vector(0,sensor_response_pol_grid.nelem(),1);
1397  if( no_los_variation )
1398  grids[2] = Vector(1,1);
1399  else
1400  grids[2] = Vector(0,sensor_response_za_grid.nelem(),1);
1401  if( no_mblock_variation )
1402  grids[3] = Vector(1,1);
1403  else
1404  grids[3] = Vector(0,sensor_pos.nrows(),1);
1405 
1406  // Create the new retrieval quantity
1407  RetrievalQuantity rq;
1408  rq.MainTag( POLYFIT_MAINTAG );
1409  rq.Mode( "" );
1410  rq.Analytical( 0 );
1411  rq.Perturbation( 0 );
1412 
1413  // Each polynomial coeff. is treated as a retrieval quantity
1414  //
1415  for( Index i=0; i<=poly_order; i++ )
1416  {
1417  ostringstream sstr;
1418  sstr << "Coefficient " << i;
1419  rq.Subtag( sstr.str() );
1420 
1421  // Grid is a scalar, use polynomial coeff.
1422  grids[0] = Vector(1,(Numeric)i);
1423  rq.Grids( grids );
1424 
1425  // Add it to the *jacobian_quantities*
1426  jq.push_back( rq );
1427 
1428  // Add pointing method to the jacobian agenda
1429  jacobian_agenda.append( "jacobianCalcPolyfit", i );
1430  }
1431 }
1432 
1433 
1434 
1435 /* Workspace method: Doxygen documentation will be auto-generated */
1437  Matrix& jacobian,
1438  const Index& mblock_index,
1439  const Vector& iyb _U_,
1440  const Vector& yb _U_,
1441  const Sparse& sensor_response,
1442  const ArrayOfIndex& sensor_response_pol_grid,
1443  const Vector& sensor_response_f_grid,
1444  const Vector& sensor_response_za_grid,
1445  const ArrayOfRetrievalQuantity& jacobian_quantities,
1446  const ArrayOfArrayOfIndex& jacobian_indices,
1447  const Index& poly_coeff,
1448  const Verbosity& )
1449 {
1450  // Find the retrieval quantity related to this method
1451  RetrievalQuantity rq;
1452  ArrayOfIndex ji;
1453  bool found = false;
1454  Index iq;
1455  ostringstream sstr;
1456  sstr << "Coefficient " << poly_coeff;
1457  for( iq=0; iq<jacobian_quantities.nelem() && !found; iq++ )
1458  {
1459  if( jacobian_quantities[iq].MainTag() == POLYFIT_MAINTAG &&
1460  jacobian_quantities[iq].Subtag() == sstr.str() )
1461  {
1462  found = true;
1463  break;
1464  }
1465  }
1466  if( !found )
1467  {
1468  throw runtime_error( "There is no Polyfit jacobian defined, in general "
1469  "or for the selected polynomial coefficient.\n");
1470  }
1471 
1472  // Size and check of sensor_response
1473  //
1474  const Index nf = sensor_response_f_grid.nelem();
1475  const Index npol = sensor_response_pol_grid.nelem();
1476  const Index nza = sensor_response_za_grid.nelem();
1477 
1478  // Make a vector with values to distribute over *jacobian*
1479  //
1480  Vector w;
1481  //
1482  polynomial_basis_func( w, sensor_response_f_grid, poly_coeff );
1483 
1484  // Fill J
1485  //
1486  ArrayOfVector jg = jacobian_quantities[iq].Grids();
1487  const Index n1 = jg[1].nelem();
1488  const Index n2 = jg[2].nelem();
1489  const Index n3 = jg[3].nelem();
1490  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
1491  const Index row4 = rowind.get_start();
1492  Index col4 = jacobian_indices[iq][0];
1493 
1494  if( n3 > 1 )
1495  { col4 += mblock_index*n2*n1; }
1496 
1497  for( Index l=0; l<nza; l++ )
1498  {
1499  const Index row3 = row4 + l*nf*npol;
1500  const Index col3 = col4 + l*n1;
1501 
1502  for( Index f=0; f<nf; f++ )
1503  {
1504  const Index row2 = row3 + f*npol;
1505 
1506  for( Index p=0; p<npol; p++ )
1507  {
1508  Index col1 = col3;
1509  if( n1 > 1 )
1510  { col1 += p; }
1511 
1512  jacobian(row2+p,col1) = w[f];
1513  }
1514  }
1515  }
1516 }
1517 
1518 
1519 
1520 
1521 
1522 //----------------------------------------------------------------------------
1523 // Sinusoidal baseline fits:
1524 //----------------------------------------------------------------------------
1525 
1526 /* Workspace method: Doxygen documentation will be auto-generated */
1528  Workspace& ws _U_,
1530  Agenda& jacobian_agenda,
1531  const ArrayOfIndex& sensor_response_pol_grid,
1532  const Vector& sensor_response_za_grid,
1533  const Matrix& sensor_pos,
1534  const Vector& period_lengths,
1535  const Index& no_pol_variation,
1536  const Index& no_los_variation,
1537  const Index& no_mblock_variation,
1538  const Verbosity& )
1539 {
1540  const Index np = period_lengths.nelem();
1541 
1542  // Check that poly_order is >= 0
1543  if( np == 0 )
1544  throw runtime_error( "No sinusoidal periods has benn given.");
1545 
1546  // Check that polyfit is not already included in the jacobian.
1547  for( Index it=0; it<jq.nelem(); it++ )
1548  {
1549  if( jq[it].MainTag() == SINEFIT_MAINTAG )
1550  {
1551  ostringstream os;
1552  os << "Polynomial baseline fit is already included in\n"
1553  << "*jacobian_quantities*.";
1554  throw runtime_error(os.str());
1555  }
1556  }
1557 
1558  // "Grids"
1559  //
1560  // Grid dimensions correspond here to
1561  // 1: polynomial order
1562  // 2: polarisation
1563  // 3: viewing direction
1564  // 4: measurement block
1565  //
1566  ArrayOfVector grids(4);
1567  //
1568  if( no_pol_variation )
1569  grids[1] = Vector(1,1);
1570  else
1571  grids[1] = Vector(0,sensor_response_pol_grid.nelem(),1);
1572  if( no_los_variation )
1573  grids[2] = Vector(1,1);
1574  else
1575  grids[2] = Vector(0,sensor_response_za_grid.nelem(),1);
1576  if( no_mblock_variation )
1577  grids[3] = Vector(1,1);
1578  else
1579  grids[3] = Vector(0,sensor_pos.nrows(),1);
1580 
1581  // Create the new retrieval quantity
1582  RetrievalQuantity rq;
1583  rq.MainTag( SINEFIT_MAINTAG );
1584  rq.Mode( "" );
1585  rq.Analytical( 0 );
1586  rq.Perturbation( 0 );
1587 
1588  // Each sinefit coeff. pair is treated as a retrieval quantity
1589  //
1590  for( Index i=0; i<np; i++ )
1591  {
1592  ostringstream sstr;
1593  sstr << "Period " << i;
1594  rq.Subtag( sstr.str() );
1595 
1596  // "Grid" has length 2, set to period length
1597  grids[0] = Vector( 2, period_lengths[i] );
1598  rq.Grids( grids );
1599 
1600  // Add it to the *jacobian_quantities*
1601  jq.push_back( rq );
1602 
1603  // Add pointing method to the jacobian agenda
1604  jacobian_agenda.append( "jacobianCalcSinefit", i );
1605  }
1606 }
1607 
1608 
1609 
1610 /* Workspace method: Doxygen documentation will be auto-generated */
1612  Matrix& jacobian,
1613  const Index& mblock_index,
1614  const Vector& iyb _U_,
1615  const Vector& yb _U_,
1616  const Sparse& sensor_response,
1617  const ArrayOfIndex& sensor_response_pol_grid,
1618  const Vector& sensor_response_f_grid,
1619  const Vector& sensor_response_za_grid,
1620  const ArrayOfRetrievalQuantity& jacobian_quantities,
1621  const ArrayOfArrayOfIndex& jacobian_indices,
1622  const Index& period_index,
1623  const Verbosity& )
1624 {
1625  // Find the retrieval quantity related to this method
1626  RetrievalQuantity rq;
1627  ArrayOfIndex ji;
1628  bool found = false;
1629  Index iq;
1630  ostringstream sstr;
1631  sstr << "Period " << period_index;
1632  for( iq=0; iq<jacobian_quantities.nelem() && !found; iq++ )
1633  {
1634  if( jacobian_quantities[iq].MainTag() == SINEFIT_MAINTAG &&
1635  jacobian_quantities[iq].Subtag() == sstr.str() )
1636  {
1637  found = true;
1638  break;
1639  }
1640  }
1641  if( !found )
1642  {
1643  throw runtime_error( "There is no Sinefit jacobian defined, in general "
1644  "or for the selected period length.\n");
1645  }
1646 
1647  // Size and check of sensor_response
1648  //
1649  const Index nf = sensor_response_f_grid.nelem();
1650  const Index npol = sensor_response_pol_grid.nelem();
1651  const Index nza = sensor_response_za_grid.nelem();
1652 
1653  // Make vectors with values to distribute over *jacobian*
1654  //
1655  // (period length stored in grid 0)
1656  ArrayOfVector jg = jacobian_quantities[iq].Grids();
1657  //
1658  Vector s(nf), c(nf);
1659  //
1660  for( Index f=0; f<nf; f++ )
1661  {
1662  Numeric a = (sensor_response_f_grid[f]-sensor_response_f_grid[0]) *
1663  2 * PI / jg[0][0];
1664  s[f] = sin( a );
1665  c[f] = cos( a );
1666  }
1667 
1668 
1669  // Fill J
1670  //
1671  const Index n1 = jg[1].nelem();
1672  const Index n2 = jg[2].nelem();
1673  const Index n3 = jg[3].nelem();
1674  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
1675  const Index row4 = rowind.get_start();
1676  Index col4 = jacobian_indices[iq][0];
1677 
1678  if( n3 > 1 )
1679  { col4 += mblock_index*n2*n1*2; }
1680 
1681  for( Index l=0; l<nza; l++ )
1682  {
1683  const Index row3 = row4 + l*nf*npol;
1684  const Index col3 = col4 + l*n1*2;
1685 
1686  for( Index f=0; f<nf; f++ )
1687  {
1688  const Index row2 = row3 + f*npol;
1689 
1690  for( Index p=0; p<npol; p++ )
1691  {
1692  Index col1 = col3;
1693  if( n1 > 1 )
1694  { col1 += p*2; }
1695 
1696  jacobian(row2+p,col1) = s[f];
1697  jacobian(row2+p,col1+1) = c[f];
1698  }
1699  }
1700  }
1701 }
1702 
1703 
1704 
1705 
1706 
1707 //----------------------------------------------------------------------------
1708 // Temperatures (atmospheric):
1709 //----------------------------------------------------------------------------
1710 
1711 /* Workspace method: Doxygen documentation will be auto-generated */
1713  Workspace& ws _U_,
1715  Agenda& jacobian_agenda,
1716  const Index& atmosphere_dim,
1717  const Vector& p_grid,
1718  const Vector& lat_grid,
1719  const Vector& lon_grid,
1720  const Vector& rq_p_grid,
1721  const Vector& rq_lat_grid,
1722  const Vector& rq_lon_grid,
1723  const String& hse,
1724  const String& method,
1725  const Numeric& dx,
1726  const Verbosity& verbosity )
1727 {
1728  CREATE_OUT3;
1729 
1730  // Check that temperature is not already included in the jacobian.
1731  // We only check the main tag.
1732  for (Index it=0; it<jq.nelem(); it++)
1733  {
1734  if( jq[it].MainTag() == TEMPERATURE_MAINTAG )
1735  {
1736  ostringstream os;
1737  os << "Temperature is already included in *jacobian_quantities*.";
1738  throw runtime_error(os.str());
1739  }
1740  }
1741 
1742  // Check that no number density retrieval has been added
1743  for (Index it=0; it<jq.nelem(); it++)
1744  {
1745  if( jq[it].MainTag() == ABSSPECIES_MAINTAG && jq[it].Mode() == "nd" )
1746  {
1747  ostringstream os;
1748  os <<
1749  "Retrieval of temperature and number densities can not be mixed.";
1750  throw runtime_error(os.str());
1751  }
1752  }
1753 
1754  // Check retrieval grids, here we just check the length of the grids
1755  // vs. the atmosphere dimension
1756  ArrayOfVector grids(atmosphere_dim);
1757  {
1758  ostringstream os;
1759  if( !check_retrieval_grids( grids, os, p_grid, lat_grid, lon_grid,
1760  rq_p_grid, rq_lat_grid, rq_lon_grid,
1761  "retrieval pressure grid",
1762  "retrieval latitude grid",
1763  "retrievallongitude_grid",
1764  atmosphere_dim ) )
1765  throw runtime_error(os.str());
1766  }
1767 
1768  // Check that method is either "analytic" or "perturbation"
1769  bool analytical;
1770  if( method == "perturbation" )
1771  { analytical = 0; }
1772  else if( method == "analytical" )
1773  { analytical = 1; }
1774  else
1775  {
1776  ostringstream os;
1777  os << "The method for atmospheric temperature retrieval can only be "
1778  << "\"analytical\"\n or \"perturbation\".";
1779  throw runtime_error(os.str());
1780  }
1781 
1782  // Set subtag
1783  String subtag;
1784  if( hse == "on" )
1785  { subtag = "HSE on"; }
1786  else if( hse == "off" )
1787  { subtag = "HSE off"; }
1788  else
1789  {
1790  ostringstream os;
1791  os << "The keyword for hydrostatic equilibrium can only be set to\n"
1792  << "\"on\" or \"off\"\n";
1793  throw runtime_error(os.str());
1794  }
1795 
1796  // Create the new retrieval quantity
1797  RetrievalQuantity rq;
1798  rq.MainTag( TEMPERATURE_MAINTAG );
1799  rq.Subtag( subtag );
1800  rq.Mode( "abs" );
1801  rq.Analytical( analytical );
1802  rq.Perturbation( dx );
1803  rq.Grids( grids );
1804 
1805  // Add it to the *jacobian_quantities*
1806  jq.push_back( rq );
1807 
1808  if( analytical )
1809  {
1810  out3 << " Calculations done by semi-analytical expression.\n";
1811  jacobian_agenda.append( "jacobianCalcTemperatureAnalytical", TokVal() );
1812  }
1813  else
1814  {
1815  out3 << " Calculations done by perturbations, of size " << dx << ".\n";
1816 
1817  jacobian_agenda.append( "jacobianCalcTemperaturePerturbations", "" );
1818  }
1819 }
1820 
1821 
1822 
1823 /* Workspace method: Doxygen documentation will be auto-generated */
1825  Matrix& jacobian _U_,
1826  const Index& mblock_index _U_,
1827  const Vector& iyb _U_,
1828  const Vector& yb _U_,
1829  const Verbosity& )
1830 {
1831  /* Nothing to do here for the analytical case, this function just exists
1832  to satisfy the required inputs and outputs of the jacobian_agenda */
1833 }
1834 
1835 
1836 
1837 
1838 /* Workspace method: Doxygen documentation will be auto-generated */
1840  Workspace& ws,
1841  Matrix& jacobian,
1842  const Index& mblock_index,
1843  const Vector& iyb _U_,
1844  const Vector& yb,
1845  const Index& atmosphere_dim,
1846  const Vector& p_grid,
1847  const Vector& lat_grid,
1848  const Vector& lon_grid,
1849  const Vector& lat_true,
1850  const Vector& lon_true,
1851  const Tensor3& t_field,
1852  const Tensor3& z_field,
1853  const Tensor4& vmr_field,
1854  const ArrayOfArrayOfSpeciesTag& abs_species,
1855  const Vector& refellipsoid,
1856  const Matrix& z_surface,
1857  const Index& cloudbox_on,
1858  const Index& stokes_dim,
1859  const Vector& f_grid,
1860  const Matrix& sensor_pos,
1861  const Matrix& sensor_los,
1862  const Matrix& transmitter_pos,
1863  const Vector& mblock_za_grid,
1864  const Vector& mblock_aa_grid,
1865  const Index& antenna_dim,
1866  const Sparse& sensor_response,
1867  const Agenda& iy_main_agenda,
1868  const Agenda& g0_agenda,
1869  const Numeric& molarmass_dry_air,
1870  const Numeric& p_hse,
1871  const Numeric& z_hse_accuracy,
1872  const ArrayOfRetrievalQuantity& jacobian_quantities,
1873  const ArrayOfArrayOfIndex& jacobian_indices,
1874  const Verbosity& verbosity )
1875 {
1876  // Set some useful variables.
1877  RetrievalQuantity rq;
1878  ArrayOfIndex ji;
1879  Index it;
1880 
1881  // Find the retrieval quantity related to this method, i.e. Temperature.
1882  // For temperature only the main tag is checked.
1883  bool found = false;
1884  for( Index n=0; n<jacobian_quantities.nelem() && !found; n++ )
1885  {
1886  if( jacobian_quantities[n].MainTag() == TEMPERATURE_MAINTAG )
1887  {
1888  found = true;
1889  rq = jacobian_quantities[n];
1890  ji = jacobian_indices[n];
1891  }
1892  }
1893  if( !found )
1894  {
1895  ostringstream os;
1896  os << "There is no temperature retrieval quantities defined.\n";
1897  throw runtime_error(os.str());
1898  }
1899 
1900  if( rq.Analytical() )
1901  {
1902  ostringstream os;
1903  os << "This WSM handles only perturbation calculations.\n"
1904  << "Are you using the method manually?";
1905  throw runtime_error(os.str());
1906  }
1907 
1908  // Store the start JacobianIndices and the Grids for this quantity
1909  it = ji[0];
1910  ArrayOfVector jg = rq.Grids();
1911 
1912  // "Perturbation mode". 1 means absolute perturbations
1913  const Index pertmode = 1;
1914 
1915  // For each atmospheric dimension option calculate a ArrayOfGridPos,
1916  // these will be used to interpolate a perturbation into the atmospheric
1917  // grids.
1918  ArrayOfGridPos p_gp, lat_gp, lon_gp;
1919  Index j_p = jg[0].nelem();
1920  Index j_lat = 1;
1921  Index j_lon = 1;
1922  //
1923  get_perturbation_gridpos( p_gp, p_grid, jg[0], true );
1924  //
1925  if( atmosphere_dim >= 2 )
1926  {
1927  j_lat = jg[1].nelem();
1928  get_perturbation_gridpos( lat_gp, lat_grid, jg[1], false );
1929  if( atmosphere_dim == 3 )
1930  {
1931  j_lon = jg[2].nelem();
1932  get_perturbation_gridpos( lon_gp, lon_grid, jg[2], false );
1933  }
1934  }
1935 
1936  // Local copy of z_field.
1937  Tensor3 z = z_field;
1938 
1939  // Loop through the retrieval grid and calculate perturbation effect
1940  //
1941  const Index n1y = sensor_response.nrows();
1942  Vector dy( n1y );
1943  const Range rowind = get_rowindex_for_mblock( sensor_response, mblock_index );
1944  //
1945  for( Index lon_it=0; lon_it<j_lon; lon_it++ )
1946  {
1947  for( Index lat_it=0; lat_it<j_lat; lat_it++ )
1948  {
1949  for( Index p_it=0; p_it<j_p; p_it++ )
1950  {
1951  // Perturbed temperature field
1952  Tensor3 t_p = t_field;
1953 
1954  // Here we calculate the ranges of the perturbation. We want the
1955  // perturbation to continue outside the atmospheric grids for the
1956  // edge values.
1957  Range p_range = Range(0,0);
1958  Range lat_range = Range(0,0);
1959  Range lon_range = Range(0,0);
1960  get_perturbation_range( p_range, p_it, j_p );
1961  if( atmosphere_dim >= 2 )
1962  {
1963  get_perturbation_range( lat_range, lat_it, j_lat );
1964  if( atmosphere_dim == 3 )
1965  {
1966  get_perturbation_range( lon_range, lon_it, j_lon );
1967  }
1968  }
1969 
1970  // Calculate the perturbed field according to atmosphere_dim,
1971  // the number of perturbations is the length of the retrieval
1972  // grid +2 (for the end points)
1973  switch (atmosphere_dim)
1974  {
1975  case 1:
1976  {
1977  // Here we perturb a vector
1978  perturbation_field_1d( t_p(joker,lat_it,lon_it),
1979  p_gp, jg[0].nelem()+2, p_range,
1980  rq.Perturbation(), pertmode );
1981  break;
1982  }
1983  case 2:
1984  {
1985  // Here we perturb a matrix
1986  perturbation_field_2d( t_p(joker,joker,lon_it),
1987  p_gp, lat_gp, jg[0].nelem()+2,
1988  jg[1].nelem()+2, p_range, lat_range,
1989  rq.Perturbation(), pertmode );
1990  break;
1991  }
1992  case 3:
1993  {
1994  // Here we need to perturb a tensor3
1995  perturbation_field_3d( t_p(joker,joker,joker), p_gp,
1996  lat_gp, lon_gp, jg[0].nelem()+2,
1997  jg[1].nelem()+2, jg[2].nelem()+2,
1998  p_range, lat_range, lon_range,
1999  rq.Perturbation(), pertmode );
2000  break;
2001  }
2002  }
2003 
2004  // Apply HSE, if selected
2005  if( rq.Subtag() == "HSE on" )
2006  {
2007  z_fieldFromHSE( ws,z, atmosphere_dim, p_grid, lat_grid,
2008  lon_grid, lat_true, lon_true, abs_species,
2009  t_p, vmr_field, refellipsoid, z_surface, 1,
2010  g0_agenda, molarmass_dry_air,
2011  p_hse, z_hse_accuracy, verbosity );
2012  }
2013 
2014  // Calculate the perturbed spectrum
2015  Vector iybp;
2016  ArrayOfVector dummy3;
2017  ArrayOfMatrix dummy4;
2018  //
2019  iyb_calc( ws, iybp, dummy3, dummy4, mblock_index,
2020  atmosphere_dim, t_p, z, vmr_field, cloudbox_on,
2021  stokes_dim, f_grid, sensor_pos, sensor_los,
2022  transmitter_pos, mblock_za_grid, mblock_aa_grid,
2023  antenna_dim, iy_main_agenda,
2025  ArrayOfArrayOfIndex(), ArrayOfString(), verbosity );
2026  //
2027  mult( dy, sensor_response, iybp );
2028 
2029  // Difference spectrum
2030  for( Index i=0; i<n1y; i++ )
2031  { dy[i] = ( dy[i]- yb[i] ) / rq.Perturbation(); }
2032 
2033  // Put into jacobian
2034  jacobian(rowind,it) = dy;
2035 
2036  // Result from next loop shall go into next column of J
2037  it++;
2038  }
2039  }
2040  }
2041 }
2042 
2043 
2044 
2045 
2046 //----------------------------------------------------------------------------
2047 // Winds:
2048 //----------------------------------------------------------------------------
2049 
2050 /* Workspace method: Doxygen documentation will be auto-generated */
2052  Workspace& ws _U_,
2054  Agenda& jacobian_agenda,
2055  const Index& atmosphere_dim,
2056  const Vector& p_grid,
2057  const Vector& lat_grid,
2058  const Vector& lon_grid,
2059  const Vector& rq_p_grid,
2060  const Vector& rq_lat_grid,
2061  const Vector& rq_lon_grid,
2062  const String& component,
2063  const Verbosity& verbosity )
2064 {
2065  CREATE_OUT2;
2066  CREATE_OUT3;
2067 
2068  // Check that this species is not already included in the jacobian.
2069  for( Index it=0; it<jq.nelem(); it++ )
2070  {
2071  if( jq[it].MainTag() == WIND_MAINTAG &&
2072  jq[it].Subtag() == component )
2073  {
2074  ostringstream os;
2075  os << "The wind component:\n" << component << "\nis already included "
2076  << "in *jacobian_quantities*.";
2077  throw runtime_error(os.str());
2078  }
2079  }
2080 
2081  // Check retrieval grids, here we just check the length of the grids
2082  // vs. the atmosphere dimension
2083  ArrayOfVector grids(atmosphere_dim);
2084  {
2085  ostringstream os;
2086  if( !check_retrieval_grids( grids, os, p_grid, lat_grid, lon_grid,
2087  rq_p_grid, rq_lat_grid, rq_lon_grid,
2088  "retrieval pressure grid",
2089  "retrieval latitude grid",
2090  "retrievallongitude_grid",
2091  atmosphere_dim ) )
2092  throw runtime_error(os.str());
2093  }
2094 
2095  // Check that component is either "u", "v" or "w"
2096  if( component != "u" && component != "v" && component != "w" )
2097  {
2098  throw runtime_error(
2099  "The selection for *component* can only be \"u\", \"u\" or \"w\"." );
2100  }
2101 
2102  // Create the new retrieval quantity
2103  RetrievalQuantity rq;
2104  rq.MainTag( WIND_MAINTAG );
2105  rq.Subtag( component );
2106  rq.Analytical( 1 );
2107  rq.Grids( grids );
2108 
2109  // Add it to the *jacobian_quantities*
2110  jq.push_back( rq );
2111 
2112  // Add gas species method to the jacobian agenda
2113  jacobian_agenda.append( "jacobianCalcWindAnalytical", TokVal() );
2114 }
2115 
2116 
2117 
2118 /* Workspace method: Doxygen documentation will be auto-generated */
2120  Matrix& jacobian _U_,
2121  const Index& mblock_index _U_,
2122  const Vector& iyb _U_,
2123  const Vector& yb _U_,
2124  const Verbosity& )
2125 {
2126  /* Nothing to do here for the analytical case, this function just exists
2127  to satisfy the required inputs and outputs of the jacobian_agenda */
2128 }
2129 
2130 
2131 
2132 
2133 
2134 
2135 
2136 
2137 
2138 
2139 
2140 //----------------------------------------------------------------------------
2141 // Old:
2142 //----------------------------------------------------------------------------
2143 
2144 
2145 // /* Workspace method: Doxygen documentation will be auto-generated */
2146 // void jacobianAddParticle(// WS Output:
2147 // ArrayOfRetrievalQuantity& jq,
2148 // Agenda& jacobian_agenda,
2149 // // WS Input:
2150 // const Matrix& jac,
2151 // const Index& atmosphere_dim,
2152 // const Vector& p_grid,
2153 // const Vector& lat_grid,
2154 // const Vector& lon_grid,
2155 // const Tensor4& pnd_field,
2156 // const Tensor5& pnd_perturb,
2157 // const ArrayOfIndex& cloudbox_limits,
2158 // // WS Generic Input:
2159 // const Vector& rq_p_grid,
2160 // const Vector& rq_lat_grid,
2161 // const Vector& rq_lon_grid,
2162 // const Verbosity& verbosity)
2163 // {
2164 // throw runtime_error("Particle jacobians not yet handled correctly.");
2165 
2166 // // Check that the jacobian matrix is empty. Otherwise it is either
2167 // // not initialised or it is closed.
2168 // if( jac.nrows()!=0 && jac.ncols()!=0 )
2169 // {
2170 // ostringstream os;
2171 // os << "The Jacobian matrix is not initialised correctly or closed.\n"
2172 // << "New retrieval quantities can not be added at this point.";
2173 // throw runtime_error(os.str());
2174 // }
2175 
2176 // // Check that pnd_perturb is consistent with pnd_field
2177 // if( pnd_perturb.nbooks()!=pnd_field.nbooks() ||
2178 // pnd_perturb.npages()!=pnd_field.npages() ||
2179 // pnd_perturb.nrows()!=pnd_field.nrows() ||
2180 // pnd_perturb.ncols()!=pnd_field.ncols() )
2181 // {
2182 // ostringstream os;
2183 // os << "The perturbation field *pnd_field_perturb* is not consistent with"
2184 // << "*pnd_field*,\none or several dimensions do not match.";
2185 // throw runtime_error(os.str());
2186 // }
2187 
2188 // // Check that particles are not already included in the jacobian.
2189 // for( Index it=0; it<jq.nelem(); it++ )
2190 // {
2191 // if( jq[it].MainTag()=="Particles" )
2192 // {
2193 // ostringstream os;
2194 // os << "The particles number densities are already included in "
2195 // << "*jacobian_quantities*.";
2196 // throw runtime_error(os.str());
2197 // }
2198 // }
2199 
2200 // // Particle Jacobian only defined for 1D and 3D atmosphere. check the
2201 // // retrieval grids, here we just check the length of the grids vs. the
2202 // // atmosphere dimension
2203 // if (atmosphere_dim==2)
2204 // {
2205 // ostringstream os;
2206 // os << "Atmosphere dimension not equal to 1 or 3. "
2207 // << "Jacobians for particle number\n"
2208 // << "density only available for 1D and 3D atmosphere.";
2209 // throw runtime_error(os.str());
2210 // }
2211 
2212 // ArrayOfVector grids(atmosphere_dim);
2213 // // The retrieval grids should only consists of gridpoints within
2214 // // the cloudbox. Setup local atmospheric fields inside the cloudbox
2215 // {
2216 // Vector p_cbox = p_grid;
2217 // Vector lat_cbox = lat_grid;
2218 // Vector lon_cbox = lon_grid;
2219 // switch (atmosphere_dim)
2220 // {
2221 // case 3:
2222 // {
2223 // lon_cbox = lon_grid[Range(cloudbox_limits[4],
2224 // cloudbox_limits[5]-cloudbox_limits[4]+1)];
2225 // }
2226 // case 2:
2227 // {
2228 // lat_cbox = lat_grid[Range(cloudbox_limits[2],
2229 // cloudbox_limits[3]-cloudbox_limits[2]+1)];
2230 // }
2231 // case 1:
2232 // {
2233 // p_cbox = p_grid[Range(cloudbox_limits[0],
2234 // cloudbox_limits[1]-cloudbox_limits[0]+1)];
2235 // }
2236 // }
2237 // ostringstream os;
2238 // if( !check_retrieval_grids( grids, os, p_cbox, lat_cbox, lon_cbox,
2239 // rq_p_grid, rq_lat_grid, rq_lon_grid,
2240 // // FIXMEOLE: These strings have to replaced later with the proper
2241 // // names from the WSM documentation in methods.cc
2242 // "rq_p_grid", "rq_lat_grid", "rq_lon_grid", atmosphere_dim ))
2243 // throw runtime_error(os.str());
2244 // }
2245 
2246 // // Common part for all particle variables
2247 // RetrievalQuantity rq;
2248 // rq.MainTag("Particles");
2249 // rq.Grids(grids);
2250 // rq.Analytical(0);
2251 // rq.Perturbation(-999.999);
2252 // rq.Mode("Fields *mode* and *perturbation* are not defined");
2253 
2254 // // Set info for each particle variable
2255 // for( Index ipt=0; ipt<pnd_perturb.nshelves(); ipt++ )
2256 // {
2257 // out2 << " Adding particle variable " << ipt +1
2258 // << " to *jacobian_quantities / agenda*.\n";
2259 
2260 // ostringstream os;
2261 // os << "Variable " << ipt+1;
2262 // rq.Subtag(os.str());
2263 
2264 // jq.push_back(rq);
2265 // }
2266 
2267 // // Add gas species method to the jacobian agenda
2268 // String methodname = "jacobianCalcParticle";
2269 // String kwv = "";
2270 // jacobian_agenda.append (methodname, kwv);
2271 // }
2272 
2273 
2274 
2275 
2276 
2277 
2278 
2279 
2280 
2281 
2282 
2283 
2284 // /* Workspace method: Doxygen documentation will be auto-generated */
2285 // void jacobianCalcParticle(
2286 // Workspace& ws,
2287 // // WS Output:
2288 // Matrix& jacobian,
2289 // // WS Input:
2290 // const Vector& y,
2291 // const ArrayOfRetrievalQuantity& jq,
2292 // const ArrayOfArrayOfIndex& jacobian_indices,
2293 // const Tensor5& pnd_field_perturb,
2294 // const Agenda& jacobian_particle_update_agenda,
2295 // const Agenda& ppath_step_agenda,
2296 // const Agenda& rte_agenda,
2297 // const Agenda& iy_space_agenda,
2298 // const Agenda& surface_rtprop_agenda,
2299 // const Agenda& iy_cloudbox_agenda,
2300 // const Index& atmosphere_dim,
2301 // const Vector& p_grid,
2302 // const Vector& lat_grid,
2303 // const Vector& lon_grid,
2304 // const Tensor3& z_field,
2305 // const Tensor3& t_field,
2306 // const Tensor4& vmr_field,
2307 // const Vector& refellipsoid,
2308 // const Matrix& z_surface,
2309 // const Index& cloudbox_on,
2310 // const ArrayOfIndex& cloudbox_limits,
2311 // const Tensor4& pnd_field,
2312 // const Sparse& sensor_response,
2313 // const Matrix& sensor_pos,
2314 // const Matrix& sensor_los,
2315 // const Vector& f_grid,
2316 // const Index& stokes_dim,
2317 // const Index& antenna_dim,
2318 // const Vector& mblock_za_grid,
2319 // const Vector& mblock_aa_grid,
2320 // const Verbosity& verbosity)
2321 // {
2322 // // Set some useful (and needed) variables.
2323 // Index n_jq = jq.nelem();
2324 // RetrievalQuantity rq;
2325 // ArrayOfIndex ji;
2326 
2327 // // Setup local atmospheric fields inside the cloudbox
2328 // Vector p_cbox = p_grid;
2329 // Vector lat_cbox = lat_grid;
2330 // Vector lon_cbox = lon_grid;
2331 // switch (atmosphere_dim)
2332 // {
2333 // case 3:
2334 // {
2335 // lon_cbox = lon_grid[Range(cloudbox_limits[4],
2336 // cloudbox_limits[5]-cloudbox_limits[4]+1)];
2337 // }
2338 // case 2:
2339 // {
2340 // lat_cbox = lat_grid[Range(cloudbox_limits[2],
2341 // cloudbox_limits[3]-cloudbox_limits[2]+1)];
2342 // }
2343 // case 1:
2344 // {
2345 // p_cbox = p_grid[Range(cloudbox_limits[0],
2346 // cloudbox_limits[1]-cloudbox_limits[0]+1)];
2347 // }
2348 // }
2349 
2350 
2351 // // Variables to handle and store perturbations
2352 // Vector yp;
2353 // Tensor4 pnd_p, base_pert = pnd_field;
2354 
2355 
2356 // // Loop particle variables (indexed by *ipt*, where *ipt* is zero based)
2357 // //
2358 // Index ipt = -1;
2359 // bool not_ready = true;
2360 
2361 // while( not_ready )
2362 // {
2363 // // Step *ipt*
2364 // ipt++;
2365 
2366 // // Define sub-tag string
2367 // ostringstream os;
2368 // os << "Variable " << ipt+1;
2369 
2370 // // Find the retrieval quantity related to this particle type
2371 // //
2372 // bool found = false;
2373 // //
2374 // for( Index n=0; n<n_jq; n++ )
2375 // {
2376 // if( jq[n].MainTag()=="Particles" && jq[n].Subtag()== os.str() )
2377 // {
2378 // found = true;
2379 // rq = jq[n];
2380 // ji = jacobian_indices[n];
2381 // n = n_jq; // To jump out of for-loop
2382 // }
2383 // }
2384 
2385 // // At least one particle type must be found
2386 // assert( !( ipt==0 && !found ) );
2387 
2388 // // Ready or something to do?
2389 // if( !found )
2390 // {
2391 // not_ready = false;
2392 // }
2393 // else
2394 // {
2395 // // Counters for report string
2396 // Index it = 0;
2397 // Index nit = ji[1] -ji[0] + 1;
2398 
2399 // // Counter for column in *jacobian*
2400 // Index icol = ji[0];
2401 
2402 // // Retrieval grid positions
2403 // ArrayOfVector jg = rq.Grids();
2404 // ArrayOfGridPos p_gp, lat_gp, lon_gp;
2405 // Index j_p = jg[0].nelem();
2406 // Index j_lat = 1;
2407 // Index j_lon = 1;
2408 // get_perturbation_gridpos( p_gp, p_cbox, jg[0], true );
2409 // if (atmosphere_dim==3)
2410 // {
2411 // j_lat = jg[1].nelem();
2412 // get_perturbation_gridpos( lat_gp, lat_cbox, jg[1], false );
2413 
2414 // j_lon = jg[2].nelem();
2415 // get_perturbation_gridpos( lon_gp, lon_cbox, jg[2], false );
2416 // }
2417 
2418 // // Give verbose output
2419 // out1 << " Calculating retrieval quantity:" << rq << "\n";
2420 
2421 
2422 // // Loop through the retrieval grid and calculate perturbation effect
2423 // for (Index lon_it=0; lon_it<j_lon; lon_it++)
2424 // {
2425 // for (Index lat_it=0; lat_it<j_lat; lat_it++)
2426 // {
2427 // for (Index p_it=0; p_it<j_p; p_it++)
2428 // {
2429 // // Update the perturbation field
2430 // pnd_p =
2431 // pnd_field_perturb( ipt, joker, joker, joker, joker);
2432 
2433 // it++;
2434 // out1 << " Calculating perturbed spectra no. " << it
2435 // << " of " << nit << "\n";
2436 
2437 // // Here we calculate the ranges of the perturbation.
2438 // // We want the perturbation to continue outside the
2439 // // atmospheric grids for the edge values.
2440 // Range p_range = Range(0,0);
2441 // Range lat_range = Range(0,0);
2442 // Range lon_range = Range(0,0);
2443 // get_perturbation_range( p_range, p_it, j_p );
2444 // if (atmosphere_dim==3)
2445 // {
2446 // get_perturbation_range( lat_range, lat_it, j_lat);
2447 // get_perturbation_range( lon_range, lon_it, j_lon);
2448 // }
2449 
2450 // // Make empty copy of pnd_pert for base functions
2451 // base_pert *= 0;
2452 
2453 // // Calculate the perturbed field according to atm_dim,
2454 // switch (atmosphere_dim)
2455 // {
2456 // case 1:
2457 // {
2458 // for( Index typ_it=0; typ_it<pnd_field.nbooks();
2459 // typ_it++ )
2460 // {
2461 // perturbation_field_1d(
2462 // base_pert(typ_it,joker,lat_it,lon_it),
2463 // p_gp, jg[0].nelem()+2, p_range, 1.0, 1 );
2464 // }
2465 // break;
2466 // }
2467 // case 3:
2468 // {
2469 // for( Index typ_it=0; typ_it<pnd_field.nrows();
2470 // typ_it++ )
2471 // {
2472 // perturbation_field_3d(
2473 // base_pert(typ_it,joker,joker,joker),
2474 // p_gp, lat_gp, lon_gp, jg[0].nelem()+2,
2475 // jg[1].nelem()+2, jg[2].nelem()+2,
2476 // p_range, lat_range, lon_range, 1.0, 1);
2477 // }
2478 // break;
2479 // }
2480 // }
2481 
2482 // // Now add the weighted perturbation field to the
2483 // // reference field and recalculate the scattered field
2484 // pnd_p *= base_pert;
2485 // pnd_p += pnd_field;
2486 // jacobian_particle_update_agendaExecute( ws, pnd_p,
2487 // jacobian_particle_update_agenda );
2488 
2489 // // Calculate the perturbed spectrum
2490 // yCalc( ws, yp, ppath_step_agenda, rte_agenda,
2491 // iy_space_agenda, surface_rtprop_agenda,
2492 // iy_cloudbox_agenda, atmosphere_dim,
2493 // p_grid, lat_grid, lon_grid, z_field, t_field,
2494 // vmr_field, refellipsoid, z_surface, cloudbox_on,
2495 // cloudbox_limits, sensor_response, sensor_pos,
2496 // sensor_los, f_grid, stokes_dim, antenna_dim,
2497 // mblock_za_grid, mblock_aa_grid);
2498 
2499 // // Add dy as column in jacobian. Note that we just return
2500 // // the difference between the two spectra.
2501 // for( Index y_it=0; y_it<yp.nelem(); y_it++ )
2502 // {
2503 // jacobian(y_it,icol) = yp[y_it]-y[y_it];
2504 // }
2505 
2506 // // Step *icol*
2507 // icol++;
2508 // }
2509 // }
2510 // }
2511 // }
2512 // }
2513 // }
2514 
2515 
2516 
2517 
2518 
2519 
2520 
2521 
2522 
2523 
2524 
2525 
2526 
2527 
const String FREQUENCY_MAINTAG
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void jacobianAddAbsSpecies(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 String &species, const String &method, const String &mode, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddAbsSpecies.
Definition: m_jacobian.cc:180
Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:196
const ArrayOfVector & Grids() const
Grids.
Definition: jacobian.h:91
void perturbation_field_1d(VectorView field, const ArrayOfGridPos &p_gp, const Index &p_pert_n, const Range &p_range, const Numeric &size, const Index &method)
Calculate the 1D perturbation for a relative perturbation.
Definition: jacobian.cc:662
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:1014
const String POINTING_CALCMODE_B
void jacobianAddFreqShift(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Matrix &sensor_pos, const Vector &sensor_time, const Index &poly_order, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqShift.
Definition: m_jacobian.cc:536
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
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 Vector &sensor_response_za_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Index &poly_coeff, const Verbosity &)
WORKSPACE METHOD: jacobianCalcPolyfit.
Definition: m_jacobian.cc:1436
const Index & Analytical() const
Boolean to make analytical calculations (if possible).
Definition: jacobian.h:85
void jacobianCalcPointingZaRecalc(Workspace &ws, Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &atmosphere_dim, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_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 Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const Vector &sensor_time, const Agenda &iy_main_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianCalcPointingZaRecalc.
Definition: m_jacobian.cc:1236
void check(Workspace &ws, const Verbosity &verbosity)
Checks consistency of an agenda.
Definition: agenda_class.cc:84
Declarations having to do with the four output streams.
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, const Index &mblock_index, const Index &atmosphere_dim, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, const Index &cloudbox_on, const Index &stokes_dim, ConstVectorView f_grid, ConstMatrixView sensor_pos, ConstMatrixView sensor_los, ConstMatrixView transmitter_pos, ConstVectorView mblock_za_grid, ConstVectorView mblock_aa_grid, const Index &antenna_dim, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
iyb_calc
Definition: rte.cc:2101
Declarations required for the calculation of jacobians.
The Vector class.
Definition: matpackI.h:556
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void jacobianClose(Workspace &ws, Index &jacobian_do, ArrayOfArrayOfIndex &jacobian_indices, Agenda &jacobian_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Matrix &sensor_pos, const Sparse &sensor_response, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianClose.
Definition: m_jacobian.cc:83
The Sparse class.
Definition: matpackII.h:55
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
const String POINTING_SUBTAG_A
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:109
const String ABSSPECIES_MAINTAG
void get_perturbation_range(Range &range, const Index &index, const Index &length)
Get range for perturbation.
Definition: jacobian.cc:632
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:62
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
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:396
void jacobianCalcAbsSpeciesPerturbations(Workspace &ws, Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, 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 Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const Agenda &iy_main_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const String &species, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianCalcAbsSpeciesPerturbations.
Definition: m_jacobian.cc:309
void jacobianAddFreqStretch(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Matrix &sensor_pos, const Vector &sensor_time, const Index &poly_order, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqStretch.
Definition: m_jacobian.cc:762
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:146
void perturbation_field_3d(Tensor3View field, const ArrayOfGridPos &p_gp, const ArrayOfGridPos &lat_gp, const ArrayOfGridPos &lon_gp, const Index &p_pert_n, const Index &lat_pert_n, const Index &lon_pert_n, const Range &p_range, const Range &lat_range, const Range &lon_range, const Numeric &size, const Index &method)
Calculate the 3D perturbation for a relative perturbation.
Definition: jacobian.cc:758
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
void jacobianCalcPointingZaInterp(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_los, const Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const Vector &sensor_time, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &)
WORKSPACE METHOD: jacobianCalcPointingZaInterp.
Definition: m_jacobian.cc:1105
Contains the data for one retrieval quantity.
Definition: jacobian.h:45
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
The implementation for String, the ARTS string class.
Definition: mystring.h:63
const String SINEFIT_MAINTAG
const String POINTING_CALCMODE_A
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition: check_input.h:164
void jacobianAddTemperature(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 String &hse, const String &method, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddTemperature.
Definition: m_jacobian.cc:1712
void jacobianInit(ArrayOfRetrievalQuantity &jacobian_quantities, ArrayOfArrayOfIndex &jacobian_indices, Agenda &jacobian_agenda, const Verbosity &)
WORKSPACE METHOD: jacobianInit.
Definition: m_jacobian.cc:144
void jacobianCalcTemperatureAnalytical(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Verbosity &)
WORKSPACE METHOD: jacobianCalcTemperatureAnalytical.
Definition: m_jacobian.cc:1824
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
const String & Mode() const
Calculation mode.
Definition: jacobian.h:82
const String FREQUENCY_SUBTAG_0
void jacobianCalcTemperaturePerturbations(Workspace &ws, Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &refellipsoid, const Matrix &z_surface, 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 Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const Agenda &iy_main_agenda, const Agenda &g0_agenda, const Numeric &molarmass_dry_air, const Numeric &p_hse, const Numeric &z_hse_accuracy, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianCalcTemperaturePerturbations.
Definition: m_jacobian.cc:1839
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:56
const String & MainTag() const
Main tag.
Definition: jacobian.h:76
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
const String POINTING_MAINTAG
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
#define dx
Definition: continua.cc:21928
The Matrix class.
Definition: matpackI.h:788
Declarations required for the calculation of absorption coefficients.
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
Array< String > ArrayOfString
An array of Strings.
Definition: mystring.h:321
void jacobianCalcAbsSpeciesAnalytical(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Verbosity &)
WORKSPACE METHOD: jacobianCalcAbsSpeciesAnalytical.
Definition: m_jacobian.cc:295
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
void jacobianOff(Index &jacobian_do, Agenda &jacobian_agenda, ArrayOfRetrievalQuantity &jacobian_quantities, ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianOff.
Definition: m_jacobian.cc:159
void jacobianCalcWindAnalytical(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Verbosity &)
WORKSPACE METHOD: jacobianCalcWindAnalytical.
Definition: m_jacobian.cc:2119
void jacobianAddWind(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 String &component, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddWind.
Definition: m_jacobian.cc:2051
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
Header file for interpolation_poly.cc.
void set_name(const String &nname)
Set agenda name.
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
void get_perturbation_gridpos(ArrayOfGridPos &gp, const Vector &atm_grid, const Vector &jac_grid, const bool &is_pressure)
Calculate array of GridPos for perturbation interpolation.
Definition: jacobian.cc:521
Array< ArrayOfIndex > ArrayOfArrayOfIndex
Definition: array.h:45
void jacobianAddPolyfit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_za_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:1354
const String POLYFIT_MAINTAG
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
get_rowindex_for_mblock
Definition: rte.cc:1960
void calc_nd_field(Tensor3View &nd, const VectorView &p, const Tensor3View &t)
Calculate the number density field.
Definition: jacobian.cc:347
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition: jacobian.cc:805
Workspace class.
Definition: workspace_ng.h:47
#define _U_
Definition: config.h:167
void jacobianCalcFreqShift(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_los, const Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const Vector &sensor_time, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &)
WORKSPACE METHOD: jacobianCalcFreqShift.
Definition: m_jacobian.cc:627
#define CREATE_OUT3
Definition: messages.h:214
const String FREQUENCY_SUBTAG_1
void jacobianCalcFreqStretch(Matrix &jacobian, const Index &mblock_index, const Vector &iyb, const Vector &yb, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_los, const Vector &mblock_za_grid, const Vector &mblock_aa_grid, const Index &antenna_dim, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Vector &sensor_response_za_grid, const Vector &sensor_time, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Verbosity &)
WORKSPACE METHOD: jacobianCalcFreqStretch.
Definition: m_jacobian.cc:853
void perturbation_field_2d(MatrixView field, const ArrayOfGridPos &p_gp, const ArrayOfGridPos &lat_gp, const Index &p_pert_n, const Index &lat_pert_n, const Range &p_range, const Range &lat_range, const Numeric &size, const Index &method)
Calculate the 2D perturbation for a relative perturbation.
Definition: jacobian.cc:707
const String TEMPERATURE_MAINTAG
void z_fieldFromHSE(Workspace &ws, Tensor3 &z_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &t_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Matrix &z_surface, const Index &atmfields_checked, const Agenda &g0_agenda, const Numeric &molarmass_dry_air, const Numeric &p_hse, const Numeric &z_hse_accuracy, const Verbosity &verbosity)
WORKSPACE METHOD: z_fieldFromHSE.
const String & Subtag() const
Subtag.
Definition: jacobian.h:79
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define CREATE_OUT2
Definition: messages.h:213
const String WIND_MAINTAG
void array_species_tag_from_string(ArrayOfSpeciesTag &tags, const String &names)
Converts a String to ArrayOfSpeciesTag.
This stores arbitrary token values and remembers the type.
Definition: token.h:33
const Numeric PI
const Numeric & Perturbation() const
Size of perturbation used for perturbation calculations.
Definition: jacobian.h:88
void jacobianAddSinefit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_za_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:1527
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
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 Vector &sensor_response_za_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const Index &period_index, const Verbosity &)
WORKSPACE METHOD: jacobianCalcSinefit.
Definition: m_jacobian.cc:1611
Declaration of functions in rte.cc.
void append(const String &methodname, const TokVal &keywordvalue)
Appends methods to an agenda.
Definition: agenda_class.cc:59