ARTS  2.2.66
jacobian.cc
Go to the documentation of this file.
1 /* Copyright (C) 2004-2012 Mattias Ekstrom <ekstrom@rss.chalmers.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
26 #include "arts.h"
27 #include "auto_md.h"
28 #include "check_input.h"
29 #include "jacobian.h"
30 #include "special_interp.h"
31 #include "physics_funcs.h"
32 
33 extern const String ABSSPECIES_MAINTAG;
34 extern const String TEMPERATURE_MAINTAG;
35 extern const String WIND_MAINTAG;
36 
37 
38 ostream& operator << (ostream& os, const RetrievalQuantity& ot)
39 {
40  return os << "\n Main tag = " << ot.MainTag()
41  << "\n Sub tag = " << ot.Subtag()
42  << "\n Mode = " << ot.Mode()
43  << "\n Analytical = " << ot.Analytical();
44 }
45 
46 
47 
48 
49 /*===========================================================================
50  === Help sub-functions to handle analytical jacobians (in alphabetical order)
51  ===========================================================================*/
52 
53 
55 
69 // Small help functions, to make the code below cleaner
71  MatrixView diy_dx,
72  ConstMatrixView diy_dq,
73  const Numeric& w )
74 {
75  for( Index irow=0; irow<diy_dx.nrows(); irow++ )
76  {
77  for( Index icol=0; icol<diy_dx.ncols(); icol++ )
78  { diy_dx(irow,icol) += w * diy_dq(irow,icol); }
79  }
80 }
82 {
83  for( Index i=0; i<gp.nelem(); i++ )
84  {
85  gp[i].idx = 0;
86  gp[i].fd[0] = 0;
87  gp[i].fd[1] = 1;
88  }
89 }
90 // Kept as a local function as long as just used here.
91 // We trust here gridpos, and "extrapolation points" are identified simply
92 // by a fd outside [0,1].
94 {
95  for( Index i=0; i<gp.nelem(); i++ )
96  {
97  if( gp[i].fd[0] < 0 )
98  {
99  gp[i].fd[0] = 0;
100  gp[i].fd[1] = 1;
101  }
102  else if( gp[i].fd[0] > 1 )
103  {
104  gp[i].fd[0] = 1;
105  gp[i].fd[1] = 0;
106  }
107  }
108 }
109 //
111  Tensor3View diy_dx,
112  const RetrievalQuantity& jacobian_quantity,
113  ConstTensor3View diy_dpath,
114  const Index& atmosphere_dim,
115  const Ppath& ppath,
116  ConstVectorView ppath_p )
117 {
118  // We want here an extrapolation to infinity ->
119  // extremly high extrapolation factor
120  const Numeric extpolfac = 1.0e99;
121 
122  if( ppath.np > 1 ) // Otherwise nothing to do here
123  {
124  // Pressure
125  Index nr1 = jacobian_quantity.Grids()[0].nelem();
126  ArrayOfGridPos gp_p(ppath.np);
127  if( nr1 > 1 )
128  {
129  p2gridpos( gp_p, jacobian_quantity.Grids()[0], ppath_p, extpolfac );
130  add_extrap( gp_p );
131  }
132  else
133  { gp4length1grid( gp_p ); }
134 
135  // Latitude
136  Index nr2 = 1;
137  ArrayOfGridPos gp_lat;
138  if( atmosphere_dim > 1 )
139  {
140  gp_lat.resize(ppath.np);
141  nr2 = jacobian_quantity.Grids()[1].nelem();
142  if( nr2 > 1 )
143  {
144  gridpos( gp_lat, jacobian_quantity.Grids()[1],
145  ppath.pos(joker,1), extpolfac );
146  add_extrap( gp_lat );
147  }
148  else
149  { gp4length1grid( gp_lat ); }
150  }
151 
152  // Longitude
153  ArrayOfGridPos gp_lon;
154  if( atmosphere_dim > 2 )
155  {
156  gp_lon.resize(ppath.np);
157  if( jacobian_quantity.Grids()[2].nelem() > 1 )
158  {
159  gridpos( gp_lon, jacobian_quantity.Grids()[2],
160  ppath.pos(joker,2), extpolfac );
161  add_extrap( gp_lon );
162  }
163  else
164  { gp4length1grid( gp_lon ); }
165  }
166 
167  //- 1D
168  if( atmosphere_dim == 1 )
169  {
170  for( Index ip=0; ip<ppath.np; ip++ )
171  {
172  if( gp_p[ip].fd[1] > 0 )
173  {
174  from_dpath_to_dx( diy_dx(gp_p[ip].idx,joker,joker),
175  diy_dpath(ip,joker,joker), gp_p[ip].fd[1] );
176  }
177  if( gp_p[ip].fd[0] > 0 )
178  {
179  from_dpath_to_dx( diy_dx(gp_p[ip].idx+1,joker,joker),
180  diy_dpath(ip,joker,joker), gp_p[ip].fd[0] );
181  }
182  }
183  }
184 
185  //- 2D
186  else if( atmosphere_dim == 2 )
187  {
188  for( Index ip=0; ip<ppath.np; ip++ )
189  {
190  Index ix = nr1*gp_lat[ip].idx + gp_p[ip].idx;
191  // Low lat, low p
192  if( gp_lat[ip].fd[1]>0 && gp_p[ip].fd[1]>0 )
193  from_dpath_to_dx( diy_dx(ix,joker,joker),
194  diy_dpath(ip,joker,joker),
195  gp_lat[ip].fd[1]*gp_p[ip].fd[1] );
196  // Low lat, high p
197  if( gp_lat[ip].fd[1]>0 && gp_p[ip].fd[0]>0 )
198  from_dpath_to_dx( diy_dx(ix+1,joker,joker),
199  diy_dpath(ip,joker,joker),
200  gp_lat[ip].fd[1]*gp_p[ip].fd[0] );
201  // High lat, low p
202  if( gp_lat[ip].fd[0]>0 && gp_p[ip].fd[1]>0 )
203  from_dpath_to_dx( diy_dx(ix+nr1,joker,joker),
204  diy_dpath(ip,joker,joker),
205  gp_lat[ip].fd[0]*gp_p[ip].fd[1] );
206  // High lat, high p
207  if( gp_lat[ip].fd[0]>0 && gp_p[ip].fd[0]>0 )
208  from_dpath_to_dx( diy_dx(ix+nr1+1,joker,joker),
209  diy_dpath(ip,joker,joker),
210  gp_lat[ip].fd[0]*gp_p[ip].fd[0] );
211  }
212  }
213 
214  //- 3D
215  else if( atmosphere_dim == 3 )
216  {
217  for( Index ip=0; ip<ppath.np; ip++ )
218  {
219  Index ix = nr2*nr1*gp_lon[ip].idx +
220  nr1*gp_lat[ip].idx + gp_p[ip].idx;
221  // Low lon, low lat, low p
222  if( gp_lon[ip].fd[1]>0 && gp_lat[ip].fd[1]>0 && gp_p[ip].fd[1]>0 )
223  from_dpath_to_dx( diy_dx(ix,joker,joker),
224  diy_dpath(ip,joker,joker),
225  gp_lon[ip].fd[1]*gp_lat[ip].fd[1]*gp_p[ip].fd[1]);
226  // Low lon, low lat, high p
227  if( gp_lon[ip].fd[1]>0 && gp_lat[ip].fd[1]>0 && gp_p[ip].fd[0]>0 )
228  from_dpath_to_dx( diy_dx(ix+1,joker,joker),
229  diy_dpath(ip,joker,joker),
230  gp_lon[ip].fd[1]*gp_lat[ip].fd[1]*gp_p[ip].fd[0]);
231  // Low lon, high lat, low p
232  if( gp_lon[ip].fd[1]>0 && gp_lat[ip].fd[0]>0 && gp_p[ip].fd[1]>0 )
233  from_dpath_to_dx( diy_dx(ix+nr1,joker,joker),
234  diy_dpath(ip,joker,joker),
235  gp_lon[ip].fd[1]*gp_lat[ip].fd[0]*gp_p[ip].fd[1]);
236  // Low lon, high lat, high p
237  if( gp_lon[ip].fd[1]>0 && gp_lat[ip].fd[0]>0 && gp_p[ip].fd[0]>0 )
238  from_dpath_to_dx( diy_dx(ix+nr1+1,joker,joker),
239  diy_dpath(ip,joker,joker),
240  gp_lon[ip].fd[1]*gp_lat[ip].fd[0]*gp_p[ip].fd[0]);
241 
242  // Increase *ix* (to be valid for high lon level)
243  ix += nr2*nr1;
244 
245  // High lon, low lat, low p
246  if( gp_lon[ip].fd[0]>0 && gp_lat[ip].fd[1]>0 && gp_p[ip].fd[1]>0 )
247  from_dpath_to_dx( diy_dx(ix,joker,joker),
248  diy_dpath(ip,joker,joker),
249  gp_lon[ip].fd[0]*gp_lat[ip].fd[1]*gp_p[ip].fd[1]);
250  // High lon, low lat, high p
251  if( gp_lon[ip].fd[0]>0 && gp_lat[ip].fd[1]>0 && gp_p[ip].fd[0]>0 )
252  from_dpath_to_dx( diy_dx(ix+1,joker,joker),
253  diy_dpath(ip,joker,joker),
254  gp_lon[ip].fd[0]*gp_lat[ip].fd[1]*gp_p[ip].fd[0]);
255  // High lon, high lat, low p
256  if( gp_lon[ip].fd[0]>0 && gp_lat[ip].fd[0]>0 && gp_p[ip].fd[1]>0 )
257  from_dpath_to_dx( diy_dx(ix+nr1,joker,joker),
258  diy_dpath(ip,joker,joker),
259  gp_lon[ip].fd[0]*gp_lat[ip].fd[0]*gp_p[ip].fd[1]);
260  // High lon, high lat, high p
261  if( gp_lon[ip].fd[0]>0 && gp_lat[ip].fd[0]>0 && gp_p[ip].fd[0]>0 )
262  from_dpath_to_dx( diy_dx(ix+nr1+1,joker,joker),
263  diy_dpath(ip,joker,joker),
264  gp_lon[ip].fd[0]*gp_lat[ip].fd[0]*gp_p[ip].fd[0]);
265  }
266  }
267  }
268 }
269 
270 
271 
273 
293  ArrayOfIndex& abs_species_i,
294  ArrayOfIndex& is_t,
295  ArrayOfIndex& wind_i,
296  const ArrayOfRetrievalQuantity& jacobian_quantities,
297  const ArrayOfArrayOfSpeciesTag& abs_species )
298 {
299 
301  //
302  if( jacobian_quantities[iq].MainTag() == TEMPERATURE_MAINTAG )
303  { is_t[iq] = 1; }
304  else
305  { is_t[iq] = 0; }
306  //
307  if( jacobian_quantities[iq].MainTag() == ABSSPECIES_MAINTAG )
308  {
309  ArrayOfSpeciesTag atag;
310  array_species_tag_from_string( atag, jacobian_quantities[iq].Subtag() );
311  abs_species_i[iq] = chk_contains( "abs_species", abs_species, atag );
312  }
313  else
314  { abs_species_i[iq] = -1; }
315  //
316  if( jacobian_quantities[iq].MainTag() == WIND_MAINTAG )
317  {
318  // Map u, v and w to 1, 2 and 3, respectively
319  char c = jacobian_quantities[iq].Subtag()[0];
320  wind_i[iq] = Index( c ) - 116;
321  }
322  else
323  { wind_i[iq] = 0; }
324  )
325 }
326 
327 
328 
329 
330 
331 /*===========================================================================
332  === Other functions, in alphabetical order
333  ===========================================================================*/
334 
336 
348  const VectorView& p,
349  const Tensor3View& t)
350 {
351  assert( nd.npages()==t.npages() );
352  assert( nd.nrows()==t.nrows() );
353  assert( nd.ncols()==t.ncols() );
354  assert( nd.npages()==p.nelem() );
355 
356  for (Index p_it=0; p_it<nd.npages(); p_it++)
357  {
358  for (Index lat_it=0; lat_it<nd.nrows(); lat_it++)
359  {
360  for (Index lon_it=0; lon_it<nd.ncols(); lon_it++)
361  {
362  nd(p_it,lat_it,lon_it) = number_density( p[p_it], t(p_it,lat_it,lon_it));
363  }
364  }
365  }
366 }
367 
368 
369 
371 
397  ostringstream& os,
398  const Vector& p_grid,
399  const Vector& lat_grid,
400  const Vector& lon_grid,
401  const Vector& p_retr,
402  const Vector& lat_retr,
403  const Vector& lon_retr,
404  const String& p_retr_name,
405  const String& lat_retr_name,
406  const String& lon_retr_name,
407  const Index& dim)
408 {
409  if ( p_retr.nelem()==0 )
410  {
411  os << "The grid vector *" << p_retr_name << "* is empty,"
412  << " at least one pressure level\n"
413  << "should be specified.";
414  return false;
415  }
416  else if( !is_decreasing( p_retr ) )
417  {
418  os << "The pressure grid vector *" << p_retr_name << "* is not a\n"
419  << "strictly decreasing vector, which is required.";
420  return false;
421  }
422  else if ( log(p_retr[0])> 1.5*log(p_grid[0])-0.5*log(p_grid[1]) ||
423  log(p_retr[p_retr.nelem()-1])<1.5*log(p_grid[p_grid.nelem()-1])-
424  0.5*log(p_grid[p_grid.nelem()-2]))
425  {
426  os << "The grid vector *" << p_retr_name << "* is not covered by the\n"
427  << "corresponding atmospheric grid.";
428  return false;
429  }
430  else
431  {
432  // Pressure grid ok, add it to grids
433  grids[0]=p_retr;
434  }
435 
436  if (dim>=2)
437  {
438  // If 2D and 3D atmosphere, check latitude grid
439  if ( lat_retr.nelem()==0 )
440  {
441  os << "The grid vector *" << lat_retr_name << "* is empty,"
442  << " at least one latitude\n"
443  << "should be specified for a 2D/3D atmosphere.";
444  return false;
445  }
446  else if( !is_increasing( lat_retr ) )
447  {
448  os << "The latitude grid vector *" << lat_retr_name << "* is not a\n"
449  << "strictly increasing vector, which is required.";
450  return false;
451  }
452  else if ( lat_retr[0]<1.5*lat_grid[0]-0.5*lat_grid[1] ||
453  lat_retr[lat_retr.nelem()-1]>1.5*lat_grid[lat_grid.nelem()-1]-
454  0.5*lat_grid[lat_grid.nelem()-2] )
455  {
456  os << "The grid vector *" << lat_retr_name << "* is not covered by the\n"
457  << "corresponding atmospheric grid.";
458  return false;
459  }
460  else
461  {
462  // Latitude grid ok, add it to grids
463  grids[1]=lat_retr;
464  }
465  if (dim==3)
466  {
467  // For 3D atmosphere check longitude grid
468  if ( lon_retr.nelem()==0 )
469  {
470  os << "The grid vector *" << lon_retr_name << "* is empty,"
471  << " at least one longitude\n"
472  << "should be specified for a 3D atmosphere.";
473  return false;
474  }
475  else if( !is_increasing( lon_retr ) )
476  {
477  os << "The longitude grid vector *" << lon_retr_name << "* is not a\n"
478  << "strictly increasing vector, which is required.";
479  return false;
480  }
481  else if ( lon_retr[0]<1.5*lon_grid[0]-0.5*lon_grid[1] ||
482  lon_retr[lon_retr.nelem()-1]>1.5*lon_grid[lon_grid.nelem()-1]-
483  0.5*lon_grid[lon_grid.nelem()-2] )
484  {
485  os << "The grid vector *" << lon_retr_name << "* is not covered by the\n"
486  << "corresponding atmospheric grid.";
487  return false;
488  }
489  else
490  {
491  // Longitude grid ok, add it to grids
492  grids[2]=lon_retr;
493  }
494  }
495  }
496  return true;
497 }
498 
499 
500 
502 
522  const Vector& atm_grid,
523  const Vector& jac_grid,
524  const bool& is_pressure)
525 {
526  Index nj = jac_grid.nelem();
527  Index na = atm_grid.nelem();
528  Vector pert(nj+2);
529 
530  // Create perturbation grid, with extension outside the atmospheric grid
531  if ( is_pressure )
532  {
533  pert[0] = atm_grid[0]*10.0;
534  pert[nj+1] = atm_grid[na-1]*0.1;
535  }
536  else
537  {
538  pert[0] = atm_grid[0]-1.0;
539  pert[nj+1] = atm_grid[na-1]+1.0;
540  }
541  pert[Range(1,nj)] = jac_grid;
542 
543  // Calculate the gridpos
544  gp.resize(na);
545  if( is_pressure ){
546  p2gridpos( gp, pert, atm_grid);
547  }
548  else
549  {
550  gridpos( gp, pert, atm_grid);
551  }
552 }
553 
554 
555 
557 
581  const Vector& pert_grid,
582  const Vector& atm_limit)
583 {
584  limit.resize(2);
585 // Index np = pert_grid.nelem()-1;
586  Index na = atm_limit.nelem()-1;
587 
588  // If the field is ordered in decreasing order set the
589  // increment factor to -1
590  Numeric inc = 1;
591  if (is_decreasing(pert_grid))
592  inc = -1;
593 
594  // Check that the pert_grid is encompassing atm_limit
595 // assert( inc*pert_grid[0] < inc*atm_limit[0] &&
596 // inc*pert_grid[np] > inc*atm_limit[na]);
597 
598  // Find first limit, check if following value is above lower limit
599  limit[0]=0;
600  while (inc*pert_grid[limit[0]+1] < inc*atm_limit[0])
601  {
602  limit[0]++;
603  }
604 
605  // Find last limit, check if previous value is below upper limit
606  limit[1]=pert_grid.nelem();
607  while (inc*pert_grid[limit[1]-1] > inc*atm_limit[na])
608  {
609  limit[1]--;
610  }
611  // Check that the limits are ok
612  assert(limit[1]>limit[0]);
613 
614 }
615 
616 
617 
619 
633  const Index& index,
634  const Index& length)
635 {
636  if (index==0)
637  range = Range(index,2);
638  else if (index==length-1)
639  range = Range(index+1,2);
640  else
641  range = Range(index+1,1);
642 
643 }
644 
645 
646 
648 
663  const ArrayOfGridPos& p_gp,
664  const Index& p_pert_n,
665  const Range& p_range,
666  const Numeric& size,
667  const Index& method)
668 {
669  // Here we only perturb a vector
670  Vector pert(field.nelem());
671  Matrix itw(p_gp.nelem(),2);
672  interpweights(itw,p_gp);
673  // For relative pert_field should be 1.0 and for absolute 0.0
674  Vector pert_field(p_pert_n,1.0-(Numeric)method);
675  pert_field[p_range] += size;
676  interp( pert, itw, pert_field, p_gp);
677  if (method==0)
678  {
679  field *= pert;
680  }
681  else
682  {
683  field += pert;
684  }
685 }
686 
687 
688 
690 
708  const ArrayOfGridPos& p_gp,
709  const ArrayOfGridPos& lat_gp,
710  const Index& p_pert_n,
711  const Index& lat_pert_n,
712  const Range& p_range,
713  const Range& lat_range,
714  const Numeric& size,
715  const Index& method)
716 {
717  // Here we perturb a matrix
718  Matrix pert(field.nrows(),field.ncols());
719  Tensor3 itw(p_gp.nelem(),lat_gp.nelem(),4);
720  interpweights(itw,p_gp,lat_gp);
721  // Init pert_field to 1.0 for relative and 0.0 for absolute
722  Matrix pert_field(p_pert_n,lat_pert_n,1.0-(Numeric)method);
723  pert_field(p_range,lat_range) += size;
724  interp( pert, itw, pert_field, p_gp, lat_gp);
725  if (method==0)
726  {
727  field *= pert;
728  }
729  else
730  {
731  field += pert;
732  }
733 }
734 
735 
736 
738 
759  const ArrayOfGridPos& p_gp,
760  const ArrayOfGridPos& lat_gp,
761  const ArrayOfGridPos& lon_gp,
762  const Index& p_pert_n,
763  const Index& lat_pert_n,
764  const Index& lon_pert_n,
765  const Range& p_range,
766  const Range& lat_range,
767  const Range& lon_range,
768  const Numeric& size,
769  const Index& method)
770 {
771  // Here we need to perturb a tensor3
772  Tensor3 pert(field.npages(),field.nrows(),field.ncols());
773  Tensor4 itw(p_gp.nelem(),lat_gp.nelem(),lon_gp.nelem(),8);
774  interpweights(itw,p_gp,lat_gp,lon_gp);
775  // Init pert_field to 1.0 for relative and 0.0 for absolute
776  Tensor3 pert_field(p_pert_n,lat_pert_n,lon_pert_n,1.0-(Numeric)method);
777  pert_field(p_range,lat_range,lon_range) += size;
778  interp( pert, itw, pert_field, p_gp, lat_gp, lon_gp);
779  if (method==0)
780  {
781  field *= pert;
782  }
783  else
784  {
785  field += pert;
786  }
787 }
788 
789 
790 
792 
806  Vector& b,
807  const Vector& x,
808  const Index& poly_coeff )
809 {
810  const Index l = x.nelem();
811 
812  assert( l > poly_coeff );
813 
814  if( b.nelem() != l )
815  b.resize( l );
816 
817  if( poly_coeff == 0 )
818  { b = 1.0; }
819  else
820  {
821  const Numeric xmin = min( x );
822  const Numeric dx = 0.5 * ( max( x ) - xmin );
823  //
824  for( Index i=0; i<l; i++ )
825  {
826  b[i] = ( x[i] - xmin ) / dx - 1.0;
827  b[i] = pow( b[i], int(poly_coeff) );
828  }
829  //
830  b -= mean( b );
831  }
832 }
833 
834 
835 
837 
854  Numeric& x,
855  const String& unit,
856  const Numeric& vmr,
857  const Numeric& p,
858  const Numeric& t )
859 {
860  if( unit == "rel" || unit == "logrel" )
861  { x = 1; }
862  else if( unit == "vmr" )
863  { x = 1 / vmr; }
864  else if( unit == "nd" )
865  { x = 1 / ( vmr * number_density( p, t ) ); }
866  else
867  {
868  throw runtime_error( "Allowed options for gas species jacobians are "
869  "\"rel\", \"vmr\", \"nd\" and \"logrel\"." );
870  }
871 }
872 
873 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void get_perturbation_limit(ArrayOfIndex &limit, const Vector &pert_grid, const Vector &atm_limit)
Get limits for perturbation of a box.
Definition: jacobian.cc:580
The VectorView class.
Definition: matpackI.h:372
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
const String TEMPERATURE_MAINTAG
Index nelem() const
Number of elements.
Definition: array.h:176
const Index & Analytical() const
Boolean to make analytical calculations (if possible).
Definition: jacobian.h:85
void gp4length1grid(ArrayOfGridPos &gp)
Definition: jacobian.cc:81
Declarations required for the calculation of jacobians.
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:275
void from_dpath_to_dx(MatrixView diy_dx, ConstMatrixView diy_dq, const Numeric &w)
diy_from_path_to_rgrids
Definition: jacobian.cc:70
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
const String ABSSPECIES_MAINTAG
void get_perturbation_range(Range &range, const Index &index, const Index &length)
Get range for perturbation.
Definition: jacobian.cc:632
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
Matrix pos
Definition: ppath.h:67
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
p2gridpos
const String WIND_MAINTAG
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
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:114
Numeric number_density(const Numeric &p, const Numeric &t)
number_density
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Contains the data for one retrieval quantity.
Definition: jacobian.h:45
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
void get_pointers_for_analytical_jacobians(ArrayOfIndex &abs_species_i, ArrayOfIndex &is_t, ArrayOfIndex &wind_i, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:292
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
ostream & operator<<(ostream &os, const RetrievalQuantity &ot)
Output operator for RetrievalQuantity.
Definition: jacobian.cc:38
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
vmrunitscf
Definition: jacobian.cc:853
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
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
const String & Mode() const
Calculation mode.
Definition: jacobian.h:82
The Tensor3View class.
Definition: matpackIII.h:232
#define max(a, b)
Definition: continua.cc:20461
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 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
void add_extrap(ArrayOfGridPos &gp)
Definition: jacobian.cc:93
Header file for special_interp.cc.
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
#define min(a, b)
Definition: continua.cc:20460
A constant view of a Tensor3.
Definition: matpackIII.h:139
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
A constant view of a Vector.
Definition: matpackI.h:292
Index np
Definition: ppath.h:61
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Definition: jacobian.cc:110
Numeric mean(const ConstVectorView &x)
Mean function, vector version.
Definition: matpackI.cc:1972
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
A constant view of a Matrix.
Definition: matpackI.h:596
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:324
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
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
const String & Subtag() const
Subtag.
Definition: jacobian.h:79
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void array_species_tag_from_string(ArrayOfSpeciesTag &tags, const String &names)
Converts a String to ArrayOfSpeciesTag.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832