ARTS  2.2.66
sensor.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2012 Mattias Ekström <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 
29 /*===========================================================================
30  === External declarations
31  ===========================================================================*/
32 
33 #include <cmath>
34 #include <list>
35 #include <stdexcept>
36 #include "arts.h"
37 #include "logic.h"
38 #include "make_vector.h"
39 #include "matpackI.h"
40 #include "matpackII.h"
41 #include "messages.h"
42 #include "sorting.h"
43 #include "sensor.h"
44 
45 extern const Numeric PI;
46 extern const Numeric NAT_LOG_2;
47 extern const Index GFIELD1_F_GRID;
48 extern const Index GFIELD4_FIELD_NAMES;
49 extern const Index GFIELD4_F_GRID;
50 extern const Index GFIELD4_ZA_GRID;
51 extern const Index GFIELD4_AA_GRID;
52 
53 
54 
55 /*===========================================================================
56  === The functions (in alphabetical order)
57  ===========================================================================*/
58 
60 
80  Sparse& H,
81 #ifndef NDEBUG
82  const Index& antenna_dim,
83 #else
84  const Index& antenna_dim _U_,
85 #endif
86  ConstMatrixView antenna_los,
87  const GriddedField4& antenna_response,
88  ConstVectorView za_grid,
89  ConstVectorView f_grid,
90  const Index n_pol,
91  const Index do_norm )
92 {
93  // Number of input za and frequency angles
94  const Index n_za = za_grid.nelem();
95  const Index n_f = f_grid.nelem();
96 
97  // Calculate number of antenna beams
98  const Index n_ant = antenna_los.nrows();
99 
100  // Asserts for variables beside antenna_response
101  assert( antenna_dim == 1 );
102  assert( antenna_los.ncols() == antenna_dim );
103  assert( n_za >= 2 );
104  assert( n_pol >= 1 );
105  assert( do_norm >= 0 && do_norm <= 1 );
106 
107  // Extract antenna_response grids
108  const Index n_ar_pol =
109  antenna_response.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
110  ConstVectorView aresponse_f_grid =
111  antenna_response.get_numeric_grid(GFIELD4_F_GRID);
112  ConstVectorView aresponse_za_grid =
113  antenna_response.get_numeric_grid(GFIELD4_ZA_GRID);
114  DEBUG_ONLY( const Index n_ar_aa =
115  antenna_response.get_numeric_grid(GFIELD4_AA_GRID).nelem(); )
116 
117  //
118  const Index n_ar_f = aresponse_f_grid.nelem();
119  const Index n_ar_za = aresponse_za_grid.nelem();
120  const Index pol_step = n_ar_pol > 1;
121 
122  // Asserts for antenna_response
123  assert( n_ar_pol == 1 || n_ar_pol == n_pol );
124  assert( n_ar_f );
125  assert( n_ar_za > 1 );
126  assert( n_ar_aa == 1 );
127 
128  // If response data extend outside za_grid is checked in
129  // sensor_integration_vector
130 
131  // Some size(s)
132  const Index nfpol = n_f * n_pol;
133 
134  // Resize H
135  H.resize( n_ant*nfpol, n_za*nfpol );
136 
137  // Storage vectors for response weights
138  Vector hrow( H.ncols(), 0.0 );
139  Vector hza( n_za, 0.0 );
140 
141  // Antenna response to apply (possibly obtained by frequency interpolation)
142  Vector aresponse( n_ar_za, 0.0 );
143 
144 
145  // Antenna beam loop
146  for( Index ia=0; ia<n_ant; ia++ )
147  {
148  Vector shifted_aresponse_za_grid = aresponse_za_grid;
149  shifted_aresponse_za_grid += antenna_los(ia,0);
150 
151  // Order of loops assumes that the antenna response more often
152  // changes with frequency than for polarisation
153 
154  // Frequency loop
155  for( Index f=0; f<n_f; f++ )
156  {
157 
158  // Polarisation loop
159  for( Index ip=0; ip<n_pol; ip++ )
160  {
161  // Determine antenna pattern to apply
162  //
163  // Interpolation needed only if response has a frequency grid
164  // New antenna for each loop of response changes with polarisation
165  //
166  Index new_antenna = 1;
167  //
168  if( n_ar_f > 1 )
169  {
170  // Interpolation (do this in "green way")
171  ArrayOfGridPos gp_f( 1 ), gp_za(n_ar_za);
172  gridpos( gp_f, aresponse_f_grid, Vector(1,f_grid[f]) );
173  gridpos( gp_za, aresponse_za_grid, aresponse_za_grid );
174  Tensor3 itw( 1, n_ar_za, 4 );
175  interpweights( itw, gp_f, gp_za );
176  Matrix aresponse_matrix(1,n_ar_za);
177  interp( aresponse_matrix, itw,
178  antenna_response.data(ip,joker,joker,0), gp_f, gp_za );
179  aresponse = aresponse_matrix(0,joker);
180  }
181  else if( pol_step ) // Response changes with polarisation
182  {
183  aresponse = antenna_response.data(ip,0,joker,0);
184  }
185  else if( f == 0 ) // Same response for all f and polarisations
186  {
187  aresponse = antenna_response.data(0,0,joker,0);
188  }
189  else
190  {
191  new_antenna = 0;
192  }
193 
194  // Calculate response weights
195  if( new_antenna )
196  {
197  sensor_integration_vector( hza, aresponse,
198  shifted_aresponse_za_grid,
199  za_grid );
200  // Normalisation?
201  if( do_norm )
202  { hza /= hza.sum(); }
203  }
204 
205  // Put weights into H
206  //
207  const Index ii = f*n_pol + ip;
208  //
209  hrow[ Range(ii,n_za,nfpol) ] = hza;
210  //
211  H.insert_row( ia*nfpol+ii, hrow );
212  //
213  hrow = 0;
214  }
215  }
216  }
217 }
218 
219 
220 
222 
248  Sparse& H,
249 #ifndef NDEBUG
250  const Index& antenna_dim,
251 #else
252  const Index& antenna_dim _U_,
253 #endif
254  ConstMatrixView antenna_los,
255  const GriddedField4& antenna_response,
256  ConstVectorView za_grid,
257  ConstVectorView aa_grid,
258  ConstVectorView f_grid,
259  const Index n_pol,
260  const Index do_norm )
261 {
262  // Sizes
263  const Index n_f = f_grid.nelem();
264  const Index nfpol = n_f * n_pol;
265  const Index n_aa = aa_grid.nelem();
266  const Index n_za = za_grid.nelem();
267  const Index n_ant = antenna_los.nrows();
268  const Index n_ar_pol = antenna_response.data.nbooks();
269  const Index n_ar_f = antenna_response.data.npages();
270  const Index n_ar_za = antenna_response.data.nrows();
271 
272  // Asserts for variables beside antenna_response (not done in antenna1d_matrix)
273  assert( antenna_dim == 2 );
274  assert( n_aa >= 2 );
275  assert( do_norm >= 0 && do_norm <= 1 );
276 
277  // Make copy of antenna response suitable as input to antenna1d_matrix
278  //
279  GriddedField4 aresponse = antenna_response;
280  //
281  ConstVectorView response_aa_grid =
282  antenna_response.get_numeric_grid(GFIELD4_AA_GRID);
283  //
284  aresponse.resize( n_ar_pol, n_ar_f, n_ar_za, 1 );
285  aresponse.set_grid( GFIELD4_AA_GRID, Vector(1,0) );
286 
287  // Resize H
288  H.resize( n_ant*nfpol, n_aa*n_za*nfpol );
289 
290  // Loop antenna_los
291  for( Index il=0; il<n_ant; il++ )
292  {
293 
294  // Set up an ArrayOfVector that can hold all data for one antenna_los
295  ArrayOfVector hrows(nfpol);
296  for( Index row=0; row<nfpol; row++ )
297  {
298  hrows[row].resize(n_aa*n_za*nfpol);
299  hrows[row] = 0; // To get correct value for aa_grid
300  } // points outside response aa grid
301 
302  // Loop azimuth angles
303  for( Index ia=0; ia<n_aa; ia++ )
304  {
305  const Numeric aa_point = aa_grid[ia] - antenna_los(il,1);
306 
307  if( aa_point >= response_aa_grid[0] &&
308  aa_point <= last(response_aa_grid) )
309  {
310  // Interpolate antenna patterns to aa_grid[ia]
311  // Use grid position function to find weights
312  //
313  ArrayOfGridPos gp( 1 );
314  gridpos( gp, response_aa_grid, Vector(1,aa_point) );
315  //
316  for( Index i4=0; i4<n_ar_pol; i4++ )
317  {
318  for( Index i3=0; i3<n_ar_f; i3++ )
319  {
320  for( Index i2=0; i2<n_ar_za; i2++ )
321  {
322  aresponse.data(i4,i3,i2,0) =
323  gp[0].fd[1] * antenna_response.data(i4,i3,i2,gp[0].idx) +
324  gp[0].fd[0] * antenna_response.data(i4,i3,i2,gp[0].idx+1);
325  }
326  }
327  }
328 
329  // Find the aa width for present angle
330  //
331  // Lower and upper end of "coverage" for present aa angle
332  Numeric aa_low = response_aa_grid[0];
333  if( ia > 0 )
334  {
335  const Numeric aam = antenna_los(il,1) +
336  ( aa_grid[ia] + aa_grid[ia-1] ) / 2.0;
337  if( aam > aa_low )
338  { aa_low = aam; };
339  }
340  Numeric aa_high = last(response_aa_grid);
341  if( ia < n_aa-1 )
342  {
343  const Numeric aam = antenna_los(il,1) +
344  ( aa_grid[ia+1] + aa_grid[ia] ) / 2.0;
345  if( aam < aa_high )
346  { aa_high = aam; };
347  }
348  //
349  const Numeric aa_width = aa_high - aa_low;
350 
351  // Do 1D calculations
352  //
353  Sparse Hza;
354  //
355  antenna1d_matrix( Hza, 1, Matrix(1,1,antenna_los(il,0)),
356  aresponse, za_grid, f_grid, n_pol, 0 );
357 
358  for( Index row=0; row<nfpol; row++ )
359  {
360  for( Index iz=0; iz<n_za; iz++ )
361  {
362  for( Index i=0; i<nfpol; i++ )
363  {
364  hrows[row][(iz*n_aa+ia)*nfpol+i] =
365  aa_width * Hza(row,iz*nfpol+i);
366  }
367  }
368  }
369  } // if-statement
370  } // aa loop
371 
372  // Move results to H
373  for( Index row=0; row<nfpol; row++ )
374  {
375  if( do_norm )
376  {
377  hrows[row] /= hrows[row].sum();
378  }
379  H.insert_row( il*nfpol+row, hrows[row] );
380  }
381 
382  } // antenna_los loop
383 }
384 
385 
386 
388 
413  Vector& x,
414  Vector& y,
415  const Numeric& x0,
416  const Numeric& fwhm,
417  const Numeric& xwidth_si,
418  const Numeric& dx_si )
419 {
420  assert( dx_si <= xwidth_si );
421 
422  const Numeric si = fwhm / ( 2 * sqrt( 2 * NAT_LOG_2 ) );
423 
424  // Number of points needed to enure that spacing is max dx_si
425  const Index n = (Index) floor(2*xwidth_si/dx_si) + 1;
426 
427  // Grid for response
428  const Numeric dd = si * xwidth_si;
429  nlinspace( x, -dd, dd, n );
430 
431  // Calculate response
432  gaussian_response( y, x, x0, fwhm );
433 }
434 
435 
436 
438 
453  Vector& y,
454  const Vector& x,
455  const Numeric& x0,
456  const Numeric& fwhm )
457 {
458  const Numeric si = fwhm / ( 2 * sqrt( 2 * NAT_LOG_2 ) );
459  const Numeric a = 1 / ( si * sqrt( 2 * PI ) );
460  const Index n = x.nelem();
461 
462  y.resize( n );
463  //
464  for( Index i=0; i<n; i++ )
465  y[i] = a * exp( -0.5 * pow((x[i]-x0)/si,2.0) );
466 }
467 
468 
469 
471 
497  Sparse& H,
498  Vector& f_mixer,
499  const Numeric& lo,
500  const GriddedField1& filter,
501  ConstVectorView f_grid,
502  const Index& n_pol,
503  const Index& n_sp,
504  const Index& do_norm )
505 {
506  // Frequency grid of for sideband response specification
507  ConstVectorView filter_grid = filter.get_numeric_grid(GFIELD1_F_GRID);
508 
509  DEBUG_ONLY( const Index nrp = filter.data.nelem(); )
510 
511  // Asserts
512  assert( lo > f_grid[0] );
513  assert( lo < last(f_grid) );
514  assert( filter_grid.nelem() == nrp );
515  assert( fabs(last(filter_grid)+filter_grid[0]) < 1e3 );
516  // If response data extend outside f_grid is checked in sensor_summation_vector
517 
518  // Find indices in f_grid where f_grid is just below and above the
519  // lo frequency.
520  /*
521  Index i_low = 0, i_high = f_grid.nelem()-1, i_mean;
522  while( i_high-i_low > 1 )
523  {
524  i_mean = (Index) (i_high+i_low)/2;
525  if (f_grid[i_mean]<lo)
526  {
527  i_low = i_mean;
528  }
529  else
530  {
531  i_high = i_mean;
532  }
533  }
534  if (f_grid[i_high]==lo)
535  {
536  i_high++;
537  }
538  const Numeric lim_low = max( lo-f_grid[i_low], f_grid[i_high]-lo );
539  */
540 
541  // Determine IF limits for new frequency grid
542  const Numeric lim_low = 0;
543  const Numeric lim_high = -filter_grid[0];
544 
545  // Convert sidebands to IF and use list to make a unique sorted
546  // vector, this sorted vector is f_mixer.
547  list<Numeric> l_mixer;
548  for( Index i=0; i<f_grid.nelem(); i++ )
549  {
550  if( fabs(f_grid[i]-lo)>=lim_low && fabs(f_grid[i]-lo)<=lim_high )
551  {
552  l_mixer.push_back(fabs(f_grid[i]-lo));
553  }
554  }
555  l_mixer.push_back(lim_high); // Not necessarily a point in f_grid
556  l_mixer.sort();
557  l_mixer.unique();
558  f_mixer.resize((Index) l_mixer.size());
559  Index e=0;
560  for (list<Numeric>::iterator li=l_mixer.begin(); li != l_mixer.end(); li++)
561  {
562  f_mixer[e] = *li;
563  e++;
564  }
565 
566  // Resize H
567  H.resize( f_mixer.nelem()*n_pol*n_sp, f_grid.nelem()*n_pol*n_sp );
568 
569  // Calculate the sensor summation vector and insert the values in the
570  // final matrix taking number of polarisations and zenith angles into
571  // account.
572  Vector row_temp( f_grid.nelem() );
573  Vector row_final( f_grid.nelem()*n_pol*n_sp );
574  //
575  Vector if_grid = f_grid;
576  if_grid -= lo;
577  //
578  for( Index i=0; i<f_mixer.nelem(); i++ )
579  {
580  sensor_summation_vector( row_temp, filter.data, filter_grid,
581  if_grid, f_mixer[i], -f_mixer[i] );
582 
583  // Normalise if flag is set
584  if (do_norm)
585  row_temp /= row_temp.sum();
586 
587  // Loop over number of polarisations
588  for (Index p=0; p<n_pol; p++)
589  {
590  // Loop over number of zenith angles/antennas
591  for (Index a=0; a<n_sp; a++)
592  {
593  // Distribute elements of row_temp to row_final.
594  row_final = 0.0;
595  row_final[Range(a*f_grid.nelem()*n_pol+p,f_grid.nelem(),n_pol)]
596  = row_temp;
597  H.insert_row(a*f_mixer.nelem()*n_pol+p+i*n_pol,row_final);
598  }
599  }
600  }
601 }
602 
603 
604 
606 
628  Vector& sensor_response_f,
629  ArrayOfIndex& sensor_response_pol,
630  Vector& sensor_response_za,
631  Vector& sensor_response_aa,
632  ConstVectorView sensor_response_f_grid,
633  const ArrayOfIndex& sensor_response_pol_grid,
634  ConstVectorView sensor_response_za_grid,
635  ConstVectorView sensor_response_aa_grid,
636  const Index za_aa_independent )
637 {
638  // Sizes
639  const Index nf = sensor_response_f_grid.nelem();
640  const Index npol = sensor_response_pol_grid.nelem();
641  const Index nza = sensor_response_za_grid.nelem();
642  Index naa = sensor_response_aa_grid.nelem();
643  Index empty_aa = 0;
644  //
645  if( naa == 0 )
646  {
647  empty_aa = 1;
648  naa = 1;
649  }
650  if( !za_aa_independent )
651  { naa = 1; }
652  //
653  const Index n = nf * npol * nza * naa;
654 
655  // Allocate
656  sensor_response_f.resize( n );
657  sensor_response_pol.resize( n );
658  sensor_response_za.resize( n );
659  if( empty_aa )
660  { sensor_response_aa.resize( 0 ); }
661  else
662  { sensor_response_aa.resize( n ); }
663 
664  // Fill
665  for( Index iaa=0; iaa<naa; iaa++ )
666  {
667  const Index i1 = iaa * nza * nf * npol;
668  //
669  for( Index iza=0; iza<nza; iza++ )
670  {
671  const Index i2 = i1 + iza * nf * npol;
672  //
673  for( Index ifr=0; ifr<nf; ifr++ )
674  {
675  const Index i3 = i2 + ifr * npol;
676  //
677  for( Index ip=0; ip<npol; ip++ )
678  {
679  const Index i = i3 + ip;
680  //
681  sensor_response_f[i] = sensor_response_f_grid[ifr];
682  sensor_response_pol[i] = sensor_response_pol_grid[ip];
683  sensor_response_za[i] = sensor_response_za_grid[iza];
684  if( !empty_aa )
685  {
686  if( za_aa_independent )
687  sensor_response_aa[i] = sensor_response_aa_grid[iaa];
688  else
689  sensor_response_aa[i] = sensor_response_aa_grid[iza];
690  }
691  }
692  }
693  }
694  }
695 }
696 
697 
698 
700 
724  VectorView h,
725  ConstVectorView f,
726  ConstVectorView x_f_in,
727  ConstVectorView x_g_in )
728 {
729  // Basic sizes
730  const Index nf = x_f_in.nelem();
731  const Index ng = x_g_in.nelem();
732 
733  // Asserts
734  assert( h.nelem() == ng );
735  assert( f.nelem() == nf );
736  assert( is_increasing( x_f_in ) );
737  assert( is_increasing( x_g_in ) || is_decreasing( x_g_in ) );
738  // More asserts below
739 
740  // End points of x_f
741  Numeric xfmin = x_f_in[0];
742  Numeric xfmax = x_f_in[nf-1];
743 
744  // Handle possibly reversed x_g.
745  Vector x_g;
746  Index xg_reversed = 0;
747  //
748  if( is_decreasing( x_g_in ) )
749  {
750  x_g = x_g_in[Range(ng-1,ng,-1)];
751  xg_reversed = 1;
752  }
753  else
754  {
755  x_g = x_g_in;
756  }
757  //
758  assert( x_g[0] <= xfmin );
759  assert( x_g[ng-1] >= xfmax );
760 
761  // Normalise grids so x_f covers [0,1]
762  const Numeric df = xfmax - xfmin;
763  Vector x_f(nf);
764  //
765  for( Index i=0; i<nf; i++ )
766  { x_f[i] = ( x_f_in[i] - xfmin ) / df; }
767  for( Index i=0; i<ng; i++ )
768  { x_g[i] = ( x_g[i] - xfmin ) / df; }
769  xfmin = 0;
770  xfmax = 1;
771  // To test without normalisation, comment out lines above and use:
772  //const Numeric df = 1;
773  //const Vector x_f = x_f_in;
774 
775  // Create a reference grid vector, x_ref that containing the values
776  // of x_f and x_g strictly sorted. Only g points inside the f range
777  // are of concern.
778  list<Numeric> l_x;
779  for( Index it=0; it<nf; it++ )
780  { l_x.push_back(x_f[it]); }
781  for (Index it=0; it<ng; it++)
782  {
783  if( x_g[it] > xfmin && x_g[it] < xfmax )
784  { l_x.push_back(x_g[it]); }
785  }
786  l_x.sort();
787  l_x.unique();
788  //
789  Vector x_ref(l_x.size());
790  Index e = 0;
791  for( list<Numeric>::iterator li=l_x.begin(); li != l_x.end(); li++ )
792  {
793  x_ref[e] = *li;
794  e++;
795  }
796 
797  // Initiate output vector, with equal size as x_g, with zeros.
798  // Start calculations
799  h = 0.0;
800  Index i_f = 0;
801  Index i_g = 0;
802  Numeric dx,a0,b0,c0,a1,b1,c1,x3,x2,x1;
803  //
804  for( Index i=0; i<x_ref.nelem()-1; i++ )
805  {
806  // Find for what index in x_g (which is the same as for h) and f
807  // calculation corresponds to
808  while( x_g[i_g+1] <= x_ref[i] )
809  { i_g++; }
810  while( x_f[i_f+1] <= x_ref[i] )
811  { i_f++; }
812 
813  // If x_ref[i] is out of x_f's range then that part of the integral is 0,
814  // and no calculations should be done
815  if( x_ref[i] >= xfmin && x_ref[i] < xfmax )
816  {
817  // Product of steps in x_f and x_g
818  dx = ( x_f[i_f+1] - x_f[i_f] ) * ( x_g[i_g+1] - x_g[i_g] );
819 
820  // Calculate a, b and c coefficients; h[i]=ax^3+bx^2+cx
821  a0 = ( f[i_f] - f[i_f+1] ) / 3.0;
822  b0 = (-f[i_f] * ( x_g[i_g+1] + x_f[i_f+1] ) +
823  f[i_f+1] * (x_g[i_g+1]+x_f[i_f]) ) / 2.0;
824  c0 = x_g[i_g+1] * ( f[i_f] * x_f[i_f+1] - f[i_f+1] * x_f[i_f] );
825 
826  a1 = -a0;
827  b1 = ( f[i_f] * ( x_g[i_g] + x_f[i_f+1] ) -
828  f[i_f+1] * ( x_g[i_g] + x_f[i_f]) ) / 2.0;
829  c1 = x_g[i_g] * ( -f[i_f] * x_f[i_f+1] + f[i_f+1] * x_f[i_f] );
830 
831  x1 = x_ref[i+1]-x_ref[i];
832  // Simple way, but sensitive to numerical problems:
833  //x2 = pow(x_ref[i+1],2) - pow(x_ref[i],2);
834  //x3 = pow(x_ref[i+1],3) - pow(x_ref[i],3);
835  // The same but a numerically better way:
836  x2 = x1 * ( 2*x_ref[i] + x1 );
837  x3 = x1 * ( 3*x_ref[i]*(x_ref[i]+x1) + x1*x1 );
838 
839  // Calculate h[i] and h[i+1] increment
840  // (df-factor to compensate for normalisation)
841  h[i_g] += df * ( a0*x3 + b0*x2 + c0*x1 ) / dx;
842  h[i_g+1] += df * ( a1*x3 + b1*x2 + c1*x1 ) / dx;
843  }
844  }
845 
846  // Flip back if x_g was decreasing
847  if( xg_reversed )
848  {
849  Vector tmp = h[Range(ng-1,ng,-1)]; // Flip order
850  h = tmp;
851  }
852 
853  // The expressions are sensitive to numerical issues if two points in x_ref
854  // are very close compared to the values in x_ref. A test trying to warn when
855  // this happens:
856  if( min(f) >= 0 && min(h) < -1e-15 )
857  throw runtime_error( "Significant negative response value obtained, "
858  "despite sensor reponse is strictly positive. This "
859  "indicates numerical problems. Is there any very "
860  "small spacing of the sensor response grid?"
861  "Please, send a report to Patrick if you see this!" );
862 }
863 
865 
887  VectorView h,
888  ConstVectorView f,
889  ConstVectorView x_f,
890  ConstVectorView x_g_in )
891 {
892  // Basic sizes
893  const Index nf = x_f.nelem();
894  const Index ng = x_g_in.nelem();
895 
896  // Asserts
897  assert( h.nelem() == ng );
898  assert( f.nelem() == nf );
899  assert( is_increasing( x_f ) );
900  assert( is_increasing( x_g_in ) || is_decreasing( x_g_in ) );
901  // More asserts below
902 
903  // End points of x_f
904  const Numeric xfmin = x_f[0];
905  const Numeric xfmax = x_f[nf-1];
906 
907  // Handle possibly reversed x_g.
908  Vector x_g;
909  Index xg_reversed = 0;
910  //
911  if( is_decreasing( x_g_in ) )
912  {
913  x_g = x_g_in[Range(ng-1,ng,-1)];
914  xg_reversed = 1;
915  }
916  else
917  {
918  x_g = x_g_in;
919  }
920  //
921  assert( x_g[0] <= xfmin );
922  assert( x_g[ng-1] >= xfmax );
923 
924  // Create a reference grid vector, x_ref that containing the values
925  // of x_f and x_g strictly sorted. Only g points inside the f range
926  // are of concern.
927  list<Numeric> l_x;
928  for( Index it=0; it<nf; it++ )
929  l_x.push_back(x_f[it]);
930  for( Index it=0; it<ng; it++ )
931  {
932  if( x_g[it] > xfmin && x_g[it] < xfmax )
933  { l_x.push_back(x_g[it]); }
934  }
935  l_x.sort();
936  l_x.unique();
937  //
938  Vector x_ref(l_x.size());
939  Index e = 0;
940  for( list<Numeric>::iterator li=l_x.begin(); li != l_x.end(); li++ )
941  {
942  x_ref[e] = *li;
943  e++;
944  }
945 
946  // Initiate output vector, with equal size as x_g, with zeros.
947  h = 0.0;
948 
949  // Start calculations
950  //
951  Index i_f = 0;
952  Index i_g = 0;
953  //
954  for( Index i=0; i<x_ref.nelem()-1; i++ )
955  {
956  // Find for what index in x_g (which is the same as for h) and f
957  // calculation corresponds to
958  while( x_g[i_g+1] <= x_ref[i] )
959  { i_g++; }
960  while( x_f[i_f+1] <= x_ref[i] )
961  { i_f++; }
962 
963  // If x_ref[i] is out of x_f's range then that part of the integral
964  // is 0, and no calculations should be done
965  if( x_ref[i] >= xfmin && x_ref[i] < xfmax )
966  {
967  // ( x_k+1 - x_k )/2
968  const Numeric dxk = 0.5 * ( x_ref[i+1] - x_ref[i] );
969 
970  // ( x_i+1 - x_i )
971  const Numeric dxi = x_g[i_g+1] - x_g[i_g];
972 
973  // The ratio of the two quantities above
974  const Numeric r = dxk / dxi;
975 
976  // Values of f at x_ref[i] and x_ref[i+1]
977  const Numeric dx = x_f[i_f+1] - x_f[i_f];
978  const Numeric f0 = ( f[i_f] * (x_f[i_f+1]-x_ref[i]) +
979  f[i_f+1] * (x_ref[i]-x_f[i_f]) ) / dx;
980  const Numeric f1 = ( f[i_f] * (x_f[i_f+1]-x_ref[i+1]) +
981  f[i_f+1] * (x_ref[i+1]-x_f[i_f]) ) / dx;
982 
983  //cout << "f0/f1 = " << f0-x_ref[i] << " / " << f1-x_ref[i+1] << endl;
984 
985  // Add increaments to h
986  h[i_g] += r * ( f0 * ( x_g[i_g+1] - x_ref[i] ) +
987  f1 * ( x_g[i_g+1] - x_ref[i+1] ) );
988  h[i_g+1] += r * ( f0 * ( x_ref[i] - x_g[i_g] ) +
989  f1 * ( x_ref[i+1] - x_g[i_g] ) );
990  }
991  }
992 
993  // Flip back if x_g was decreasing
994  if( xg_reversed )
995  {
996  Vector tmp = h[Range(ng-1,ng,-1)]; // Flip order
997  h = tmp;
998  }
999 }
1000 
1001 
1002 
1004 
1031  VectorView h,
1032  ConstVectorView f,
1033  ConstVectorView x_f,
1034  ConstVectorView x_g,
1035  const Numeric x1,
1036  const Numeric x2 )
1037 {
1038  // Asserts
1039  assert( h.nelem() == x_g.nelem() );
1040  assert( f.nelem() == x_f.nelem() );
1041  assert( x_g[0] <= x_f[0] );
1042  assert( last(x_g) >= last(x_f) );
1043  assert( x1 >= x_f[0] );
1044  assert( x2 >= x_f[0] );
1045  assert( x1 <= last(x_f) );
1046  assert( x2 <= last(x_f) );
1047 
1048  // Determine grid positions for point 1 (both with respect to f and g grids)
1049  // and interpolate response function.
1050  ArrayOfGridPos gp1g(1), gp1f(1);
1051  gridpos( gp1g, x_g, x1 );
1052  gridpos( gp1f, x_f, x1 );
1053  Matrix itw1(1,2);
1054  interpweights( itw1, gp1f );
1055  Numeric f1;
1056  interp( f1, itw1, f, gp1f );
1057 
1058  // Same for point 2
1059  ArrayOfGridPos gp2g(1), gp2f(1);
1060  gridpos( gp2g, x_g, x2 );
1061  gridpos( gp2f, x_f, x2 );
1062  Matrix itw2(1,2);
1063  interpweights( itw2, gp2f );
1064  Numeric f2;
1065  interp( f2, itw2, f, gp2f );
1066 
1067  //Initialise h at zero and store calculated weighting components
1068  h = 0.0;
1069  h[gp1g[0].idx] += f1 * gp1g[0].fd[1];
1070  h[gp1g[0].idx+1] += f1 * gp1g[0].fd[0];
1071  h[gp2g[0].idx] += f2 * gp2g[0].fd[1];
1072  h[gp2g[0].idx+1] += f2 * gp2g[0].fd[0];
1073 }
1074 
1075 
1076 
1078 
1097  Sparse& H,
1098  ConstVectorView ch_f,
1099  const ArrayOfGriddedField1& ch_response,
1100  ConstVectorView sensor_f,
1101  const Index& n_pol,
1102  const Index& n_sp,
1103  const Index& do_norm )
1104 {
1105  // Check if matrix has one frequency column or one for every channel
1106  // frequency
1107  //
1108  assert( ch_response.nelem()==1 || ch_response.nelem()==ch_f.nelem() );
1109  //
1110  Index freq_full = ch_response.nelem() > 1;
1111 
1112  // If response data extend outside sensor_f is checked in
1113  // sensor_integration_vector
1114 
1115  // Reisze H
1116  //
1117  const Index nin_f = sensor_f.nelem();
1118  const Index nout_f = ch_f.nelem();
1119  const Index nin = n_sp * nin_f * n_pol;
1120  const Index nout = n_sp * nout_f * n_pol;
1121  //
1122  H.resize( nout, nin );
1123 
1124  // Calculate the sensor integration vector and put values in the temporary
1125  // vector, then copy vector to the transfer matrix
1126  //
1127  Vector ch_response_f;
1128  Vector weights( nin_f );
1129  Vector weights_long( nin, 0.0 );
1130  //
1131  for( Index ifr=0; ifr<nout_f; ifr++ )
1132  {
1133  const Index irp = ifr * freq_full;
1134 
1135  //The spectrometer response is shifted for each centre frequency step
1136  ch_response_f = ch_response[irp].get_numeric_grid(GFIELD1_F_GRID);
1137  ch_response_f += ch_f[ifr];
1138 
1139  // Call sensor_integration_vector and store it in the temp vector
1140  sensor_integration_vector( weights, ch_response[irp].data,
1141  ch_response_f, sensor_f );
1142 
1143  // Normalise if flag is set
1144  if( do_norm )
1145  weights /= weights.sum();
1146 
1147  // Loop over polarisation and spectra (viewing directions)
1148  // Weights change only with frequency
1149  for( Index sp=0; sp<n_sp; sp++ )
1150  {
1151  for( Index pol=0; pol<n_pol; pol++ )
1152  {
1153  // Distribute the compact weight vector into weight_long
1154  weights_long[Range(sp*nin_f*n_pol+pol,nin_f,n_pol)] = weights;
1155 
1156  // Insert temp_long into H at the correct row
1157  H.insert_row( sp*nout_f*n_pol + ifr*n_pol + pol, weights_long );
1158 
1159  // Reset weight_long to zero.
1160  weights_long = 0.0;
1161  }
1162  }
1163  }
1164 }
1165 
1166 
1167 
1169 
1197  ArrayOfVector& s2p,
1198  const Numeric& nv )
1199 {
1200  s2p.resize(10);
1201 
1202  s2p[0] = MakeVector( 1 ); // I
1203  s2p[1] = MakeVector( 0, 1 ); // Q
1204  s2p[2] = MakeVector( 0, 0, 1 ); // U
1205  s2p[3] = MakeVector( 0, 0, 0, 1 ); // V
1206  s2p[4] = MakeVector( nv, nv ); // Iv
1207  s2p[5] = MakeVector( nv, -nv ); // Ih
1208  s2p[6] = MakeVector( nv, 0, nv ); // I+45
1209  s2p[7] = MakeVector( nv, 0, -nv ); // I-45
1210  s2p[8] = MakeVector( nv, 0, 0, nv ); // Irhc
1211  s2p[9] = MakeVector( nv, 0, 0, -nv ); // Ilhc
1212 }
1213 
1214 
1215 
1216 
1217 // Functions by Stefan, needed for HIRS:
1218 
1219 
1221 
1250  Vector& fmax,
1251  Index i,
1252  Index j)
1253 {
1254  const Index nf = fmin.nelem();
1255  assert(fmax.nelem()==nf);
1256  assert(i>=0 && i<nf);
1257  assert(j>=0 && j<nf);
1258  assert(fmin[i]<=fmin[j]);
1259  assert(i<j);
1260 
1261  // There are three cases to consider:
1262  // a) The two channels are separate: fmax[i] < fmin[j]
1263  // b) They overlap: fmax[i] >= fmin[j]
1264  // c) j is inside i: fmax[i] > fmax[j]
1265 
1266  // In the easiest case (a), we do not have to do anything.
1267  if (fmax[i] >= fmin[j])
1268  {
1269  // We are in case (b) or (c), so we know that we have to combine
1270  // the channels. The new minimum frequency is fmin[i]. The new
1271  // maximum frequency is the larger one of the two channels we
1272  // are combining:
1273  if (fmax[j] > fmax[i])
1274  fmax[i] = fmax[j];
1275 
1276  // We now have to kick out element j from both vectors.
1277 
1278  // Number of elements behind j:
1279  Index n_behind = nf-1 - j;
1280 
1281  Vector dummy = fmin;
1282  fmin.resize(nf-1);
1283  fmin[Range(0,j)] = dummy[Range(0,j)];
1284  if (n_behind > 0)
1285  fmin[Range(j,n_behind)] = dummy[Range(j+1,n_behind)];
1286 
1287  dummy = fmax;
1288  fmax.resize(nf-1);
1289  fmax[Range(0,j)] = dummy[Range(0,j)];
1290  if (n_behind > 0)
1291  fmax[Range(j,n_behind)] = dummy[Range(j+1,n_behind)];
1292 
1293  return true;
1294  }
1295 
1296  return false;
1297 }
1298 
1299 
1301 
1326  Vector& fmin,
1327  Vector& fmax,
1328  // Input:
1329  const Vector& f_backend,
1330  const ArrayOfGriddedField1& backend_channel_response,
1331  const Numeric& delta,
1332  const Verbosity& verbosity)
1333 {
1334  CREATE_OUT2;
1335 
1336  // How many channels in total:
1337  const Index n_chan = f_backend.nelem();
1338 
1339  // Checks on input quantities:
1340 
1341  // There must be at least one channel.
1342  if (n_chan < 1)
1343  {
1344  ostringstream os;
1345  os << "There must be at least one channel.\n"
1346  << "(The vector *f_backend* must have at least one element.)";
1347  throw runtime_error(os.str());
1348  }
1349 
1350  // There must be a response function for each channel.
1351  if (n_chan != backend_channel_response.nelem())
1352  {
1353  ostringstream os;
1354  os << "Variables *f_backend_multi* and *backend_channel_response_multi*\n"
1355  << "must have same number of bands for each LO.";
1356  throw runtime_error(os.str());
1357  }
1358 
1359  // Frequency grids for response functions must be strictly increasing.
1360  for (Index i=0; i<n_chan; ++i)
1361  {
1362  // Frequency grid for this response function:
1363  const Vector& backend_f_grid = backend_channel_response[i].get_numeric_grid(0);
1364 
1365  if ( !is_increasing(backend_f_grid) )
1366  {
1367  ostringstream os;
1368  os << "The frequency grid for the backend channel response of\n"
1369  << "channel " << i << " is not strictly increasing.\n";
1370  os << "It is: " << backend_f_grid << "\n";
1371  throw runtime_error( os.str() );
1372  }
1373  }
1374 
1375 
1376  // Start the actual work.
1377 
1378  out2 << " Original channel characteristics:\n"
1379  << " min nominal max (all in Hz):\n";
1380 
1381  // count the number of passbands as defined by segments of filtershapes, backend_channel_response.data = [0,>0,>0,0]
1382  // Borders between passbands are identified as [...0,0...]
1383 
1384  // Get a list of original channel boundaries:
1385  Index numPB = 0;
1386  for(Index idx=0;idx<n_chan;++idx)
1387  {
1388  const Vector& backend_filter = backend_channel_response[idx].data;
1389  if(backend_filter.nelem()>2)
1390  { // only run this code when there is more then two elements in the backend
1391  for(Index idy =1;idy < backend_filter.nelem();++idy)
1392  {
1393  if((backend_filter[idy] >0)&& (backend_filter[idy-1]==0))
1394  {
1395  numPB ++; // Two consecutive zeros gives the border between two passbands
1396  }
1397  }
1398  }
1399  else
1400  {
1401  numPB ++;
1402  }
1403  }
1404 
1405  if (!numPB)
1406  {
1407  throw std::runtime_error("No passbands found.\n"
1408  "*backend_channel_response* must be zero around the passbands.\n"
1409  "backend_channel_response.data = [0, >0, >0, 0]\n"
1410  "Borders between passbands are identified as [...0,0...]");
1411  }
1412 
1413  Vector fmin_pb(numPB);
1414  Vector fmax_pb(numPB);
1415  Index pbIdx = 0;
1416 
1417  for(Index idx=0;idx<n_chan;++idx)
1418  {
1419  // Some handy shortcuts:
1420  //
1421  // We have to find the first and last frequency where the
1422  // response is actually different from 0. (No point in making
1423  // calculations for frequencies where the response is 0.)
1424  // Index j=0;
1425  // while (backend_response[j] <= 0) ++j;
1426  // Numeric bf_min = backend_f_grid[j];
1427 
1428  // j=nf-1;
1429  // while (backend_response[j] <= 0) --j;
1430  // Numeric bf_max = backend_f_grid[j];
1431  //
1432  // No, aparently the sensor part want values also where the
1433  // response is zero. So we simply take the grid boundaries here.
1434  //
1435  // We need to add a bit of extra margin at both sides,
1436  // otherwise there is a numerical problem in the sensor WSMs.
1437  //
1438  // PE 081003: The accuracy for me (double on 64 bit machine) appears to
1439  // be about 3 Hz. Select 1 MHz to have a clear margin. Hopefully OK
1440  // for other machines.
1441  //
1442  // SAB 2010-04-14: The approach with a constant delta does not seem to work
1443  // well in practice. What I do now is that I add a delta corresponding to a
1444  // fraction of the grid spacing. But that is done outside of this function.
1445  // So we set delta = 0 here.
1446  //
1447  // SAB 2010-05-03: Now we pass delta as a parameter (with a default value of 0).
1448  //
1449  // Isoz 2013-05-21: Added methods to ignore areas between passbands
1450  //
1451  const Vector& backend_f_grid = backend_channel_response[idx].get_numeric_grid(0);
1452  const Vector& backend_filter = backend_channel_response[idx].data;
1453  if(backend_filter.nelem()>=4)// Is the passband frequency response given explicitly ? e.g. [0,>0,>0,0]
1454  {
1455  for(Index idy =1;idy < backend_filter.nelem();++idy)
1456  {
1457  if(idy==1)
1458  {
1459  fmin_pb[pbIdx] = f_backend[idx]+backend_f_grid[0];
1460  }
1461  else if((backend_filter[idy] >0) && (backend_filter[idy-1]==0))
1462  {
1463  fmin_pb[pbIdx]= f_backend[idx] + backend_f_grid[idy-1]-delta;
1464  }
1465  if((backend_filter[idy] ==0) && (backend_filter[idy-1]>0))
1466  {
1467  fmax_pb[pbIdx]= f_backend[idx] + backend_f_grid[idy]+delta;
1468  out2 << " " << "fmin_pb "<< fmin_pb[pbIdx]
1469  << " " << "f_backend" <<f_backend[idx]
1470  << " " << "fmax_pb "<<fmax_pb[pbIdx]
1471  << " " << "diff " << fmax_pb[pbIdx]-fmin_pb[pbIdx]
1472  << "\n";
1473  pbIdx++;
1474  }
1475  }
1476  fmax_pb[pbIdx-1] = f_backend[idx]+last(backend_f_grid);
1477  }
1478  else // Or are the passbands given implicitly - such as the default for AMSUA and MHS
1479  {
1480  fmin_pb[pbIdx]= f_backend[idx] + backend_f_grid[0]-delta ; //delta;
1481  fmax_pb[pbIdx]= f_backend[idx] + backend_f_grid[backend_f_grid.nelem()-1]+delta;
1482  out2 << " " << "fmin_pb "<< fmin_pb[pbIdx]
1483  << " " << "f_backend" <<f_backend[pbIdx]
1484  << " " << "fmax_pb "<<fmax_pb[pbIdx]
1485  << "\n";
1486  pbIdx++;
1487  }
1488  }
1489 
1490  // The problem is that channels may be overlapping. In that case, we
1491  // want to create a frequency grid that covers their entire range,
1492  // but we do not want to duplicate any frequencies.
1493 
1494  // To make matters worse, one or even several channels may be
1495  // completely inside another very broad channel.
1496 
1497  // Sort channels by frequency:
1498  // Caveat: A channel may be higher in
1499  // characteristic frequency f_backend, but also wider, so that it
1500  // has a lower minimum frequency fmin_orig. (This is the case for
1501  // some HIRS channels.) We sort by the minimum frequency here, not
1502  // by f_backend. This is necessary for function
1503  // test_and_merge_two_channels to work correctly.
1504  ArrayOfIndex isorted;
1505  get_sorted_indexes(isorted,fmin_pb);
1506 
1507  fmin.resize(numPB);
1508  fmax.resize(numPB);
1509  out2 <<" resize numPb "<<numPB<<"\n";
1510  for(Index idx = 0;idx<numPB;++idx)
1511  {
1512  fmin[idx] = fmin_pb[isorted[idx]];
1513  fmax[idx] = fmax_pb[isorted[idx]];
1514  }
1515  // We will be testing pairs of channels, and combine them if
1516  // possible. We have to test always only against the direct
1517  // neighbour. If that has no overlap, higher channels can not have
1518  // any either, due to the sorting by fmin.
1519  //
1520  // Note that fmin.nelem() changes, as the loop is
1521  // iterated. Nevertheless this is the correct stop condition.
1522  for (Index i=0; i<fmin.nelem()-1; ++i)
1523  {
1524  bool continue_checking = true;
1525  // The "i<fmin.nelem()" condition below is necessary, since
1526  // fmin.nelem() can decrease while the loop is executed, due to
1527  // merging.
1528  while (continue_checking && i<fmin.nelem()-1)
1529  {
1530  continue_checking =
1531  test_and_merge_two_channels(fmin, fmax, i, i+1);
1532 
1533  // Function returns true if merging has taken place.
1534  // In this case, we have to check again.
1535 
1536  }
1537  }
1538 
1539  out2 << " New channel characteristics:\n"
1540  << " min max (all in Hz):\n";
1541  for (Index i=0; i<fmin.nelem(); ++i)
1542  out2 << " " << fmin[i] << " " << fmax[i] << "\n";
1543 
1544 }
1545 
1546 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
The VectorView class.
Definition: matpackI.h:372
void gaussian_response(Vector &y, const Vector &x, const Numeric &x0, const Numeric &fwhm)
gaussian_response
Definition: sensor.cc:452
void gaussian_response_autogrid(Vector &x, Vector &y, const Numeric &x0, const Numeric &fwhm, const Numeric &xwidth_si, const Numeric &dx_si)
gaussian_response_autogrid
Definition: sensor.cc:412
void spectrometer_matrix(Sparse &H, ConstVectorView ch_f, const ArrayOfGriddedField1 &ch_response, ConstVectorView sensor_f, const Index &n_pol, const Index &n_sp, const Index &do_norm)
spectrometer_matrix
Definition: sensor.cc:1096
Index nelem() const
Number of elements.
Definition: array.h:176
const Index GFIELD4_FIELD_NAMES
Declarations having to do with the four output streams.
ConstVectorView get_numeric_grid(Index i) const
Get a numeric grid.
The Vector class.
Definition: matpackI.h:556
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void antenna1d_matrix(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_los, const GriddedField4 &antenna_response, ConstVectorView za_grid, ConstVectorView f_grid, const Index n_pol, const Index do_norm)
antenna1d_matrix
Definition: sensor.cc:79
const ArrayOfString & get_string_grid(Index i) const
Get a string grid.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:275
The Sparse class.
Definition: matpackII.h:55
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:183
The range class.
Definition: matpackI.h:148
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:62
void sensor_aux_vectors(Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Vector &sensor_response_za, Vector &sensor_response_aa, ConstVectorView sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, ConstVectorView sensor_response_za_grid, ConstVectorView sensor_response_aa_grid, const Index za_aa_independent)
sensor_aux_vectors
Definition: sensor.cc:627
void sensor_integration_vector2(VectorView h, ConstVectorView f, ConstVectorView x_f, ConstVectorView x_g_in)
sensor_integration_vector
Definition: sensor.cc:886
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:75
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Contains sorting routines.
void sensor_integration_vector(VectorView h, ConstVectorView f, ConstVectorView x_f_in, ConstVectorView x_g_in)
sensor_integration_vector
Definition: sensor.cc:723
Sensor modelling functions.
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
#define DEBUG_ONLY(...)
Definition: debug.h:37
Header file for sparse matrices.
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:186
#define b0
Definition: continua.cc:21471
void set_grid(Index i, const Vector &g)
Set a numeric grid.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:261
bool test_and_merge_two_channels(Vector &fmin, Vector &fmax, Index i, Index j)
Test if two instrument channels overlap, and if so, merge them.
Definition: sensor.cc:1249
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
const Index GFIELD4_ZA_GRID
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 mixer_matrix(Sparse &H, Vector &f_mixer, const Numeric &lo, const GriddedField1 &filter, ConstVectorView f_grid, const Index &n_pol, const Index &n_sp, const Index &do_norm)
mixer_matrix
Definition: sensor.cc:496
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:79
void sensor_summation_vector(VectorView h, ConstVectorView f, ConstVectorView x_f, ConstVectorView x_g, const Numeric x1, const Numeric x2)
sensor_summation_vector
Definition: sensor.cc:1030
const Numeric NAT_LOG_2
Header file for logic.cc.
const Index GFIELD1_F_GRID
The class MakeVector is a special kind of Vector that can be initialized explicitly from one or more ...
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 Vector.
Definition: matpackI.h:292
A constant view of a Matrix.
Definition: matpackI.h:596
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
#define _U_
Definition: config.h:167
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:469
void antenna2d_simplified(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_los, const GriddedField4 &antenna_response, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView f_grid, const Index n_pol, const Index do_norm)
antenna2d_simplified
Definition: sensor.cc:247
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:324
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:264
#define CREATE_OUT2
Definition: messages.h:213
const Numeric PI
void stokes2pol(ArrayOfVector &s2p, const Numeric &nv)
stokes2pol
Definition: sensor.cc:1196
const Index GFIELD4_AA_GRID
const Index GFIELD4_F_GRID
void find_effective_channel_boundaries(Vector &fmin, Vector &fmax, const Vector &f_backend, const ArrayOfGriddedField1 &backend_channel_response, const Numeric &delta, const Verbosity &verbosity)
Calculate channel boundaries from instrument response functions.
Definition: sensor.cc:1325