ARTS  2.4.0(git:4fb77825)
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 "sensor.h"
34 #include <cmath>
35 #include <list>
36 #include <stdexcept>
37 #include "arts.h"
38 #include "logic.h"
39 #include "matpackI.h"
40 #include "matpackII.h"
41 #include "messages.h"
42 #include "sorting.h"
43 
44 extern const Numeric PI;
45 extern const Numeric NAT_LOG_2;
46 extern const Numeric DEG2RAD;
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  === The functions, besides the core integration and sum functions that are
55  === placed together at the end
56  ===========================================================================*/
57 
59 #ifndef NDEBUG
60  const Index& antenna_dim,
61 #else
62  const Index& antenna_dim _U_,
63 #endif
64  ConstVectorView antenna_dza,
68  const Index n_pol,
69  const Index do_norm) {
70  // Number of input pencil beam directions and frequency
71  const Index n_za = za_grid.nelem();
72  const Index n_f = f_grid.nelem();
73 
74  // Calculate number of antenna beams
75  const Index n_ant = antenna_dza.nelem();
76 
77  // Asserts for variables beside antenna_response
78  assert(antenna_dim == 1);
79  assert(n_za >= 2);
80  assert(n_pol >= 1);
81  assert(do_norm >= 0 && do_norm <= 1);
82 
83  // Extract antenna_response grids
84  const Index n_ar_pol =
85  antenna_response.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
86  ConstVectorView aresponse_f_grid =
87  antenna_response.get_numeric_grid(GFIELD4_F_GRID);
88  ConstVectorView aresponse_za_grid =
89  antenna_response.get_numeric_grid(GFIELD4_ZA_GRID);
90  DEBUG_ONLY(const Index n_ar_aa =
91  antenna_response.get_numeric_grid(GFIELD4_AA_GRID).nelem();)
92 
93  //
94  const Index n_ar_f = aresponse_f_grid.nelem();
95  const Index n_ar_za = aresponse_za_grid.nelem();
96  const Index pol_step = n_ar_pol > 1;
97 
98  // Asserts for antenna_response
99  assert(n_ar_pol == 1 || n_ar_pol >= n_pol);
100  assert(n_ar_f);
101  assert(n_ar_za > 1);
102  assert(n_ar_aa == 1);
103 
104  // If response data extend outside za_grid is checked in
105  // integration_func_by_vecmult
106 
107  // Some size(s)
108  const Index nfpol = n_f * n_pol;
109 
110  // Resize H
111  H.resize(n_ant * nfpol, n_za * nfpol);
112 
113  // Storage vectors for response weights
114  Vector hrow(H.ncols(), 0.0);
115  Vector hza(n_za, 0.0);
116 
117  // Antenna response to apply (possibly obtained by frequency interpolation)
118  Vector aresponse(n_ar_za, 0.0);
119 
120  // Antenna beam loop
121  for (Index ia = 0; ia < n_ant; ia++) {
122  Vector shifted_aresponse_za_grid = aresponse_za_grid;
123  shifted_aresponse_za_grid += antenna_dza[ia];
124 
125  // Order of loops assumes that the antenna response more often
126  // changes with frequency than for polarisation
127 
128  // Frequency loop
129  for (Index f = 0; f < n_f; f++) {
130  // Polarisation loop
131  for (Index ip = 0; ip < n_pol; ip++) {
132  // Determine antenna pattern to apply
133  //
134  // Interpolation needed only if response has a frequency grid
135  //
136  Index new_antenna = 1;
137  //
138  if (n_ar_f == 1) // No frequency variation
139  {
140  if (pol_step) // Polarisation variation, update always needed
141  {
142  aresponse = antenna_response.data(ip, 0, joker, 0);
143  } else if (f == 0 && ip == 0) // Set fully constant pattern
144  {
145  aresponse = antenna_response.data(0, 0, joker, 0);
146  } else // The one set just above can be reused
147  {
148  new_antenna = 0;
149  }
150  } else {
151  if (ip == 0 || pol_step) // Interpolation required
152  {
153  // Interpolation (do this in "green way")
154  ArrayOfGridPos gp_f(1), gp_za(n_ar_za);
155  gridpos(gp_f, aresponse_f_grid, Vector(1, f_grid[f]));
156  gridpos(gp_za, aresponse_za_grid, aresponse_za_grid);
157  Tensor3 itw(1, n_ar_za, 4);
158  interpweights(itw, gp_f, gp_za);
159  Matrix aresponse_matrix(1, n_ar_za);
160  interp(aresponse_matrix,
161  itw,
162  antenna_response.data(ip, joker, joker, 0),
163  gp_f,
164  gp_za);
165  aresponse = aresponse_matrix(0, joker);
166  } else // Reuse pattern for ip==0
167  {
168  new_antenna = 0;
169  }
170  }
171 
172  // Calculate response weights
173  if (new_antenna) {
175  hza, aresponse, shifted_aresponse_za_grid, za_grid);
176 
177  // Normalisation?
178  if (do_norm) {
179  hza /= hza.sum();
180  }
181  }
182 
183  // Put weights into H
184  //
185  const Index ii = f * n_pol + ip;
186  //
187  hrow[Range(ii, n_za, nfpol)] = hza;
188  //
189  H.insert_row(ia * nfpol + ii, hrow);
190  //
191  hrow = 0;
192  }
193  }
194  }
195 }
196 
197 
198 
200 #ifndef NDEBUG
201  const Index& antenna_dim,
202 #else
203  const Index& antenna_dim _U_,
204 #endif
207  ConstMatrixView mblock_dlos,
209  const Index n_pol)
210 {
211  // Number of input pencil beam directions and frequency
212  const Index n_dlos = mblock_dlos.nrows();
213  const Index n_f = f_grid.nelem();
214 
215  // Decompose mblock_dlos into za and aa grids, including checking
216  if( mblock_dlos.ncols() != 2 )
217  throw runtime_error("For the gridded_dlos option, *mblock_dlos_grid* "
218  "must have two columns.");
219 
220 
221  Index nza = 1;
222  for(Index i=0; i<n_dlos-1 && mblock_dlos(i+1,0) > mblock_dlos(i,0); i++ ) {
223  nza++;
224  }
225  if( nza < 2 )
226  throw runtime_error("For the gridded_dlos option, the number of za angles "
227  "(among dlos directions) must be >= 2.");
228  if( n_dlos % nza > 0 )
229  throw runtime_error("For the gridded_dlos option, the number of dlos angles "
230  "must be a product of two integers.");
231  const Index naa = n_dlos / nza;
232  const Vector za_grid = mblock_dlos(Range(0,nza),0);
233  const Vector aa_grid = mblock_dlos(Range(0,naa,nza),1);
234  for(Index i=0; i<n_dlos; i++ ) {
235  if(i>=nza && abs(mblock_dlos(i,0)-mblock_dlos(i-nza,0)) > 1e-6 ) {
236  ostringstream os;
237  os << "Zenith angle of dlos " << i << " (0-based) differs to zenith "
238  << "angle of dlos " << i-nza << ", while they need to be equal "
239  << "to form rectangular grid.";
240  throw std::runtime_error(os.str());
241  }
242  if(abs(mblock_dlos(i,1)-aa_grid[i/nza]) > 1e-6) {
243  ostringstream os;
244  os << "Azimuth angle of dlos " << i << " (0-based) differs to azimuth "
245  << "angle " << (i/nza)*nza << ", while they need to be equal "
246  << "to form rectangular grid.";
247  throw std::runtime_error(os.str());
248  }
249  }
250 
251  // Calculate number of antenna beams
252  const Index n_ant = antenna_dlos.nrows();
253 
254  // Asserts for variables beside antenna_response
255  assert(antenna_dim == 2);
256  assert(n_dlos >= 1);
257  assert(n_pol >= 1);
258 
259  // Extract antenna_response grids
260  const Index n_ar_pol =
261  antenna_response.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
262  ConstVectorView aresponse_f_grid =
263  antenna_response.get_numeric_grid(GFIELD4_F_GRID);
264  ConstVectorView aresponse_za_grid =
265  antenna_response.get_numeric_grid(GFIELD4_ZA_GRID);
266  ConstVectorView aresponse_aa_grid =
267  antenna_response.get_numeric_grid(GFIELD4_AA_GRID);
268  //
269  const Index n_ar_f = aresponse_f_grid.nelem();
270  const Index n_ar_za = aresponse_za_grid.nelem();
271  const Index n_ar_aa = aresponse_aa_grid.nelem();
272  const Index pol_step = n_ar_pol > 1;
273 
274  // Asserts for antenna_response
275  assert(n_ar_pol == 1 || n_ar_pol >= n_pol);
276  assert(n_ar_f);
277  assert(n_ar_za > 1);
278  assert(n_ar_aa > 1);
279 
280  // Some size(s)
281  const Index nfpol = n_f * n_pol;
282 
283  // Resize H
284  H.resize(n_ant * nfpol, n_dlos * nfpol);
285 
286  // Storage vectors for response weights
287  Vector hrow(H.ncols(), 0.0);
288  Vector hza(n_dlos, 0.0);
289 
290  // Antenna response to apply (possibly obtained by frequency interpolation)
291  Matrix aresponse(n_ar_za, n_ar_aa, 0.0);
292 
293  // If you find a bug or change something, likely also change other 2D antenna
294  // function(s) as they are similar
295 
296  // Antenna beam loop
297  for (Index ia = 0; ia < n_ant; ia++) {
298  // Order of loops assumes that the antenna response more often
299  // changes with frequency than for polarisation
300 
301  // Frequency loop
302  for (Index f = 0; f < n_f; f++) {
303  // Polarisation loop
304  for (Index ip = 0; ip < n_pol; ip++) {
305  // Determine antenna pattern to apply
306  //
307  // Interpolation needed only if response has a frequency grid
308  //
309  Index new_antenna = 1;
310  //
311  if (n_ar_f == 1) // No frequency variation
312  {
313  if (pol_step) // Polarisation variation, update always needed
314  {
315  aresponse = antenna_response.data(ip, 0, joker, joker);
316  } else if (f == 0 && ip == 0) // Set fully constant pattern
317  {
318  aresponse = antenna_response.data(0, 0, joker, joker);
319  } else // The one set just above can be reused
320  {
321  new_antenna = 0;
322  }
323  } else {
324  if (ip == 0 || pol_step) {
325  // Interpolation (do this in "green way")
326  ArrayOfGridPos gp_f(1), gp_za(n_ar_za), gp_aa(n_ar_aa);
327  gridpos(gp_f, aresponse_f_grid, Vector(1, f_grid[f]));
328  gridpos(gp_za, aresponse_za_grid, aresponse_za_grid);
329  gridpos(gp_aa, aresponse_aa_grid, aresponse_aa_grid);
330  Tensor4 itw(1, n_ar_za, n_ar_aa, 8);
331  interpweights(itw, gp_f, gp_za, gp_aa);
332  Tensor3 aresponse_matrix(1, n_ar_za, n_ar_aa);
333  interp(aresponse_matrix,
334  itw,
335  antenna_response.data(ip, joker, joker, joker),
336  gp_f,
337  gp_za,
338  gp_aa);
339  aresponse = aresponse_matrix(0, joker, joker);
340  } else // Reuse pattern for ip==0
341  {
342  new_antenna = 0;
343  }
344  }
345 
346  // Calculate response weights, by using grid positions and "itw"
347  if (new_antenna) {
348 
349  // za grid positions
350  Vector zas = aresponse_za_grid;
351  zas += antenna_dlos(ia, 0);
352  if( zas[0] < za_grid[0] ) {
353  ostringstream os;
354  os << "The zenith angle grid in *mblock_dlos_grid* is too narrow. "
355  << "It must be extended downwards with at least "
356  << za_grid[0]-zas[0] << " deg.";
357  throw std::runtime_error(os.str());
358  }
359  if( zas[n_ar_za-1] > za_grid[nza-1] ) {
360  ostringstream os;
361  os << "The zenith angle grid in *mblock_dlos_grid* is too narrow. "
362  << "It must be extended upwards with at least "
363  << zas[n_ar_za-1] - za_grid[nza-1] << " deg.";
364  throw std::runtime_error(os.str());
365  }
366  ArrayOfGridPos gp_za(n_ar_za);
367  gridpos(gp_za, za_grid, zas);
368 
369  // aa grid positions
370  Vector aas = aresponse_aa_grid;
371  if (antenna_dlos.ncols() > 1) { aas += antenna_dlos(ia, 1); }
372  if( aas[0] < aa_grid[0] ) {
373  ostringstream os;
374  os << "The azimuth angle grid in *mblock_dlos_grid* is too narrow. "
375  << "It must be extended downwards with at least "
376  << aa_grid[0]-aas[0] << " deg.";
377  throw std::runtime_error(os.str());
378  }
379  if( aas[n_ar_aa-1] > aa_grid[naa-1] ) {
380  ostringstream os;
381  os << "The azimuth angle grid in *mblock_dlos_grid* is too narrow. "
382  << "It must be extended upwards with at least "
383  << aas[n_ar_aa-1] - aa_grid[naa-1] << " deg.";
384  throw std::runtime_error(os.str());
385  }
386  ArrayOfGridPos gp_aa(n_ar_aa);
387  gridpos(gp_aa, aa_grid, aas);
388 
389 
390  // Derive interpolation weights
391  Tensor3 itw(n_ar_za, n_ar_za, 4);
392  interpweights(itw, gp_za, gp_aa);
393 
394  // Convert iwt to weights for H
395  //
396  hza = 0; // Note that values in hza must be accumulated
397  //
398  for (Index iaa = 0; iaa < n_ar_aa; iaa++) {
399  const Index a = gp_aa[iaa].idx;
400  const Index b = a + 1;
401 
402  for (Index iza = 0; iza < n_ar_za; iza++) {
403 
404  const Index z = gp_za[iza].idx;
405  const Index x = z + 1;
406 
407  if( itw(iza,iaa,0) > 1e-9 ) {
408  hza[a*nza+z] += aresponse(iza,iaa) * itw(iza,iaa,0);
409  }
410  if( itw(iza,iaa,1) > 1e-9 ) {
411  hza[b*nza+z] += aresponse(iza,iaa) * itw(iza,iaa,1);
412  }
413  if( itw(iza,iaa,2) > 1e-9 ) {
414  hza[a*nza+x] += aresponse(iza,iaa) * itw(iza,iaa,2);
415  }
416  if( itw(iza,iaa,3) > 1e-9 ) {
417  hza[b*nza+x] += aresponse(iza,iaa) * itw(iza,iaa,3);
418  }
419  }
420  }
421 
422  // For 2D antennas we always normalise
423  hza /= hza.sum();
424  }
425 
426  // Put weights into H
427  //
428  const Index ii = f * n_pol + ip;
429  //
430  hrow[Range(ii, n_dlos, nfpol)] = hza;
431  //
432  H.insert_row(ia * nfpol + ii, hrow);
433  //
434  hrow = 0;
435  }
436  }
437  }
438 }
439 
440 
441 
443 #ifndef NDEBUG
444  const Index& antenna_dim,
445 #else
446  const Index& antenna_dim _U_,
447 #endif
450  ConstMatrixView mblock_dlos,
452  const Index n_pol)
453 {
454  // Number of input pencil beam directions and frequency
455  const Index n_dlos = mblock_dlos.nrows();
456  const Index n_f = f_grid.nelem();
457 
458  // Calculate number of antenna beams
459  const Index n_ant = antenna_dlos.nrows();
460 
461  // Asserts for variables beside antenna_response
462  assert(antenna_dim == 2);
463  assert(n_dlos >= 1);
464  assert(n_pol >= 1);
465 
466  // Extract antenna_response grids
467  const Index n_ar_pol =
468  antenna_response.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
469  ConstVectorView aresponse_f_grid =
470  antenna_response.get_numeric_grid(GFIELD4_F_GRID);
471  ConstVectorView aresponse_za_grid =
472  antenna_response.get_numeric_grid(GFIELD4_ZA_GRID);
473  ConstVectorView aresponse_aa_grid =
474  antenna_response.get_numeric_grid(GFIELD4_AA_GRID);
475  //
476  const Index n_ar_f = aresponse_f_grid.nelem();
477  const Index n_ar_za = aresponse_za_grid.nelem();
478  const Index n_ar_aa = aresponse_aa_grid.nelem();
479  const Index pol_step = n_ar_pol > 1;
480 
481  // Asserts for antenna_response
482  assert(n_ar_pol == 1 || n_ar_pol >= n_pol);
483  assert(n_ar_f);
484  assert(n_ar_za > 1);
485  assert(n_ar_aa > 1);
486 
487  // Some size(s)
488  const Index nfpol = n_f * n_pol;
489 
490  // Resize H
491  H.resize(n_ant * nfpol, n_dlos * nfpol);
492 
493  // Storage vectors for response weights
494  Vector hrow(H.ncols(), 0.0);
495  Vector hza(n_dlos, 0.0);
496 
497  // Antenna response to apply (possibly obtained by frequency interpolation)
498  Matrix aresponse(n_ar_za, n_ar_aa, 0.0);
499 
500  // If you find a bug or change something, likely also change other 2D antenna
501  // function(s) as they are similar
502 
503  // Antenna beam loop
504  for (Index ia = 0; ia < n_ant; ia++) {
505  // Order of loops assumes that the antenna response more often
506  // changes with frequency than for polarisation
507 
508  // Frequency loop
509  for (Index f = 0; f < n_f; f++) {
510  // Polarisation loop
511  for (Index ip = 0; ip < n_pol; ip++) {
512  // Determine antenna pattern to apply
513  //
514  // Interpolation needed only if response has a frequency grid
515  //
516  Index new_antenna = 1;
517  //
518  if (n_ar_f == 1) // No frequency variation
519  {
520  if (pol_step) // Polarisation variation, update always needed
521  {
522  aresponse = antenna_response.data(ip, 0, joker, joker);
523  } else if (f == 0 && ip == 0) // Set fully constant pattern
524  {
525  aresponse = antenna_response.data(0, 0, joker, joker);
526  } else // The one set just above can be reused
527  {
528  new_antenna = 0;
529  }
530  } else {
531  if (ip == 0 || pol_step) {
532  // Interpolation (do this in "green way")
533  ArrayOfGridPos gp_f(1), gp_za(n_ar_za), gp_aa(n_ar_aa);
534  gridpos(gp_f, aresponse_f_grid, Vector(1, f_grid[f]));
535  gridpos(gp_za, aresponse_za_grid, aresponse_za_grid);
536  gridpos(gp_aa, aresponse_aa_grid, aresponse_aa_grid);
537  Tensor4 itw(1, n_ar_za, n_ar_aa, 8);
538  interpweights(itw, gp_f, gp_za, gp_aa);
539  Tensor3 aresponse_matrix(1, n_ar_za, n_ar_aa);
540  interp(aresponse_matrix,
541  itw,
542  antenna_response.data(ip, joker, joker, joker),
543  gp_f,
544  gp_za,
545  gp_aa);
546  aresponse = aresponse_matrix(0, joker, joker);
547  } else // Reuse pattern for ip==0
548  {
549  new_antenna = 0;
550  }
551  }
552 
553  // Calculate response weights
554  if (new_antenna) {
555  for (Index l = 0; l < n_dlos; l++) {
556  const Numeric za = mblock_dlos(l, 0) - antenna_dlos(ia, 0);
557  Numeric aa = 0.0;
558  if (mblock_dlos.ncols() > 1) {
559  aa += mblock_dlos(l, 1);
560  }
561  if (antenna_dlos.ncols() > 1) {
562  aa -= antenna_dlos(ia, 1);
563  }
564 
565  // The response is zero if mblock_dlos is outside of
566  // antennna pattern
567  if (za < aresponse_za_grid[0] ||
568  za > aresponse_za_grid[n_ar_za - 1] ||
569  aa < aresponse_aa_grid[0] ||
570  aa > aresponse_aa_grid[n_ar_aa - 1]) {
571  hza[l] = 0;
572  }
573  // Otherwise we make an (blue) interpolation
574  else {
575  ArrayOfGridPos gp_za(1), gp_aa(1);
576  gridpos(gp_za, aresponse_za_grid, Vector(1, za));
577  gridpos(gp_aa, aresponse_aa_grid, Vector(1, aa));
578  Matrix itw(1, 4);
579  interpweights(itw, gp_za, gp_aa);
580  Vector value(1);
581  interp(value, itw, aresponse, gp_za, gp_aa);
582  hza[l] = value[0];
583  }
584  }
585 
586  // For 2D antennas we always normalise
587  hza /= hza.sum();
588  }
589 
590  // Put weights into H
591  //
592  const Index ii = f * n_pol + ip;
593  //
594  hrow[Range(ii, n_dlos, nfpol)] = hza;
595  //
596  H.insert_row(ia * nfpol + ii, hrow);
597  //
598  hrow = 0;
599  }
600  }
601  }
602 }
603 
604 
605 
607  Vector& y,
608  const Numeric& x0,
609  const Numeric& fwhm,
610  const Numeric& xwidth_si,
611  const Numeric& dx_si) {
612  assert(dx_si <= xwidth_si);
613 
614  const Numeric si = fwhm / (2 * sqrt(2 * NAT_LOG_2));
615 
616  // Number of points needed to enure that spacing is max dx_si
617  const Index n = (Index)floor(2 * xwidth_si / dx_si) + 1;
618 
619  // Grid for response
620  const Numeric dd = si * xwidth_si;
621  nlinspace(x, -dd, dd, n);
622 
623  // Calculate response
624  gaussian_response(y, x, x0, fwhm);
625 }
626 
627 
628 
630  const Vector& x,
631  const Numeric& x0,
632  const Numeric& fwhm) {
633  const Numeric si = fwhm / (2 * sqrt(2 * NAT_LOG_2));
634  const Numeric a = 1 / (si * sqrt(2 * PI));
635  const Index n = x.nelem();
636 
637  y.resize(n);
638  //
639  for (Index i = 0; i < n; i++)
640  y[i] = a * exp(-0.5 * pow((x[i] - x0) / si, 2.0));
641 }
642 
643 
644 
646  Vector& f_mixer,
647  const Numeric& lo,
648  const GriddedField1& filter,
650  const Index& n_pol,
651  const Index& n_sp,
652  const Index& do_norm) {
653  // Frequency grid of for sideband response specification
654  ConstVectorView filter_grid = filter.get_numeric_grid(GFIELD1_F_GRID);
655 
656  DEBUG_ONLY(const Index nrp = filter.data.nelem();)
657 
658  // Asserts
659  assert(lo > f_grid[0]);
660  assert(lo < last(f_grid));
661  assert(filter_grid.nelem() == nrp);
662  assert(fabs(last(filter_grid) + filter_grid[0]) < 1e3);
663  // If response data extend outside f_grid is checked in summation_by_vecmult
664 
665  // Find indices in f_grid where f_grid is just below and above the
666  // lo frequency.
667  /*
668  Index i_low = 0, i_high = f_grid.nelem()-1, i_mean;
669  while( i_high-i_low > 1 )
670  {
671  i_mean = (Index) (i_high+i_low)/2;
672  if (f_grid[i_mean]<lo)
673  {
674  i_low = i_mean;
675  }
676  else
677  {
678  i_high = i_mean;
679  }
680  }
681  if (f_grid[i_high]==lo)
682  {
683  i_high++;
684  }
685  const Numeric lim_low = max( lo-f_grid[i_low], f_grid[i_high]-lo );
686  */
687 
688  // Determine IF limits for new frequency grid
689  const Numeric lim_low = 0;
690  const Numeric lim_high = -filter_grid[0];
691 
692  // Convert sidebands to IF and use list to make a unique sorted
693  // vector, this sorted vector is f_mixer.
694  list<Numeric> l_mixer;
695  for (Index i = 0; i < f_grid.nelem(); i++) {
696  if (fabs(f_grid[i] - lo) >= lim_low && fabs(f_grid[i] - lo) <= lim_high) {
697  l_mixer.push_back(fabs(f_grid[i] - lo));
698  }
699  }
700  l_mixer.push_back(lim_high); // Not necessarily a point in f_grid
701  l_mixer.sort();
702  l_mixer.unique();
703  f_mixer.resize((Index)l_mixer.size());
704  Index e = 0;
705  for (list<Numeric>::iterator li = l_mixer.begin(); li != l_mixer.end();
706  li++) {
707  f_mixer[e] = *li;
708  e++;
709  }
710 
711  // Resize H
712  H.resize(f_mixer.nelem() * n_pol * n_sp, f_grid.nelem() * n_pol * n_sp);
713 
714  // Calculate the sensor summation vector and insert the values in the
715  // final matrix taking number of polarisations and zenith angles into
716  // account.
717  Vector row_temp(f_grid.nelem());
718  Vector row_final(f_grid.nelem() * n_pol * n_sp);
719  //
720  Vector if_grid = f_grid;
721  if_grid -= lo;
722  //
723  for (Index i = 0; i < f_mixer.nelem(); i++) {
725  row_temp, filter.data, filter_grid, if_grid, f_mixer[i], -f_mixer[i]);
726 
727  // Normalise if flag is set
728  if (do_norm) row_temp /= row_temp.sum();
729 
730  // Loop over number of polarisations
731  for (Index p = 0; p < n_pol; p++) {
732  // Loop over number of zenith angles/antennas
733  for (Index a = 0; a < n_sp; a++) {
734  // Distribute elements of row_temp to row_final.
735  row_final = 0.0;
736  row_final[Range(
737  a * f_grid.nelem() * n_pol + p, f_grid.nelem(), n_pol)] = row_temp;
738  H.insert_row(a * f_mixer.nelem() * n_pol + p + i * n_pol, row_final);
739  }
740  }
741  }
742 }
743 
744 
746  const Index& stokes_dim,
747  const Numeric& rotangle) {
748  assert(stokes_dim > 1);
749  assert(H.nrows() == stokes_dim);
750  assert(H.ncols() == stokes_dim);
751  assert(H(0, 1) == 0);
752  assert(H(1, 0) == 0);
753 
754  H.rw(0, 0) = 1;
755  const Numeric a = cos(2 * DEG2RAD * rotangle);
756  H.rw(1, 1) = a;
757  if (stokes_dim > 2) {
758  assert(H(2, 0) == 0);
759  assert(H(0, 2) == 0);
760 
761  const Numeric b = sin(2 * DEG2RAD * rotangle);
762  H.rw(1, 2) = b;
763  H.rw(2, 1) = -b;
764  H.rw(2, 2) = a;
765  if (stokes_dim > 3) {
766  // More values should be checked, but to save time we just assert one
767  assert(H(2, 3) == 0);
768  H.rw(3, 3) = 1;
769  }
770  }
771 }
772 
773 
774 
776  const ArrayOfString& mm_pol,
777  const Numeric dza,
778  const Index stokes_dim,
779  const String& iy_unit) {
780  assert(stokes_dim > 1);
781 
782  // Set "Stokes vector weights" according to iy_unit
783  Numeric w = 0.5;
784  if (iy_unit == "PlanckBT" || iy_unit == "RJBT") {
785  w = 1.0;
786  }
787 
788  // Identify (basic) polarisation response and possible sensor specific
789  // "rotation"
790  const Index nch = mm_pol.nelem(); // Number of channels
791  ArrayOfString pol(nch);
792  ArrayOfString rot(nch);
793  for (Index i = 0; i < nch; i++) {
794  if (mm_pol[i] == "AMSU-H") {
795  rot[i] = "AMSU";
796  pol[i] = "H";
797  } else if (mm_pol[i] == "AMSU-V") {
798  rot[i] = "AMSU";
799  pol[i] = "V";
800  } else if (mm_pol[i] == "ISMAR-H") {
801  rot[i] = "ISMAR";
802  pol[i] = "H";
803  } else if (mm_pol[i] == "ISMAR-V") {
804  rot[i] = "ISMAR";
805  pol[i] = "V";
806  } else if (mm_pol[i] == "MARSS-H") {
807  rot[i] = "MARSS";
808  pol[i] = "H";
809  } else if (mm_pol[i] == "MARSS-V") {
810  rot[i] = "MARSS";
811  pol[i] = "V";
812  } else if (mm_pol[i] == "H" || mm_pol[i] == "V" || mm_pol[i] == "LHC" ||
813  mm_pol[i] == "RHC") {
814  rot[i] = "none";
815  pol[i] = mm_pol[i];
816  } else {
817  ostringstream os;
818  os << "Unknown polarisation " << mm_pol[i];
819  throw std::runtime_error(os.str());
820  }
821  }
822 
823  // Vectors representing standard cases of sensor polarisation response
824  /*
825  ArrayOfVector pv;
826  stokes2pol( pv, w );
827  */
828 
829  // Complete H, for all channels
830  H = Sparse(nch, nch * stokes_dim);
831 
832  for (Index i = 0; i < nch; i++) {
833  /*
834  // See stokes2pol for index order used in pv
835  Index ipv = -1;
836  if( pol[i] == "V" )
837  { ipv = 4; }
838  else if( pol[i] == "H" )
839  { ipv = 5; }
840  else if( pol[i] == "LHC" ) // Left hand circular
841  { ipv = 8; }
842  else if( pol[i] == "RHC" ) // Right hand circular
843  { ipv = 9; }
844  else
845  { assert( 0 ); }
846  */
847  // See instrument_pol for index order
848  Index ipol = -1;
849  if (pol[i] == "V") {
850  ipol = 5;
851  } else if (pol[i] == "H") {
852  ipol = 6;
853  } else if (pol[i] == "LHC") // Left hand circular
854  {
855  ipol = 9;
856  } else if (pol[i] == "RHC") // Right hand circular
857  {
858  ipol = 10;
859  } else {
860  assert(0);
861  }
862 
863  /*
864  // Maybe this error messages should be mofified
865  if( pv[ipv].nelem() > stokes_dim )
866  {
867  ostringstream os;
868  os << "You have selected a channel polarisation that is not "
869  << "covered by present value of *stokes_dim* (the later has to be "
870  << "increased).";
871  throw runtime_error(os.str());
872  }*/
873 
874  // No rotation, just plane polarisation response
875  if (rot[i] == "none") {
876  // Here we just need to fill the row H
877  Vector hrow(nch * stokes_dim, 0.0);
878  /* Old code, matching older version of stokes2pol:
879  hrow[Range(i*stokes_dim,pv[ipv].nelem())] = pv[ipv];
880  */
881  stokes2pol(hrow[Range(i * stokes_dim, stokes_dim)], stokes_dim, ipol, w);
882  H.insert_row(i, hrow);
883  }
884 
885  // Rotation + pol-response
886  else {
887  // Rotation part
889  if (rot[i] == "AMSU") {
890  // No idea about the sign. Not important if U=0,
891  // but matter for U != 0.
892  mueller_rotation(Hrot, stokes_dim, abs(dza));
893  } else if (rot[i] == "ISMAR") {
894  // No rotation at -53 (= forward direction). But no idea about the
895  // sign, as for AMSU above
896  mueller_rotation(Hrot, stokes_dim, dza + 50);
897  } else if (rot[i] == "MARSS") {
898  // MARSS special, as 48 deg between different polarisation (not 90,
899  // as for ISMAR. This is best interpretation of information
900  // from Stuart. Should be double-checked with him at some point.
901  if (pol[i] == "H") {
902  mueller_rotation(Hrot, stokes_dim, dza + 42);
903  } else {
904  mueller_rotation(Hrot, stokes_dim, dza);
905  }
906  } else {
907  assert(0);
908  }
909 
910  // H-matrix matching polarization
911  Sparse Hpol(1, stokes_dim);
912  { /* Old code, matching older version of stokes2pol
913  Vector hrow( stokes_dim, 0.0 );
914  hrow[Range(0,pv[ipv].nelem())] = pv[ipv];
915  */
916  Vector hrow(stokes_dim);
917  stokes2pol(hrow, stokes_dim, ipol, w);
918  Hpol.insert_row(0, hrow);
919  }
920 
921  // H for the individual channel
922  Sparse Hc(1, stokes_dim);
923  mult(Hc, Hpol, Hrot);
924 
925  // Put Hc into H
926  Vector hrow(nch * stokes_dim, 0.0);
927  const Index i0 = i * stokes_dim;
928  for (Index s = 0; s < stokes_dim; s++) {
929  hrow[i0 + s] = Hc(0, s);
930  }
931  H.insert_row(i, hrow);
932  }
933  }
934 }
935 
936 
937 
944  // Sizes
945  const Index nf = sensor_response_f_grid.nelem();
946  const Index npol = sensor_response_pol_grid.nelem();
947  const Index nlos = sensor_response_dlos_grid.nrows();
948  const Index n = nf * npol * nlos;
949 
950  // Allocate
951  sensor_response_f.resize(n);
952  sensor_response_pol.resize(n);
954 
955  // Fill
956  for (Index ilos = 0; ilos < nlos; ilos++) {
957  const Index i2 = ilos * nf * npol;
958  //
959  for (Index ifr = 0; ifr < nf; ifr++) {
960  const Index i3 = i2 + ifr * npol;
961  //
962  for (Index ip = 0; ip < npol; ip++) {
963  const Index i = i3 + ip;
964  //
968  }
969  }
970  }
971 }
972 
973 
974 
976  ConstVectorView ch_f,
977  const ArrayOfGriddedField1& ch_response,
978  ConstVectorView sensor_f,
979  const Index& n_pol,
980  const Index& n_sp,
981  const Index& do_norm) {
982  // Check if matrix has one frequency column or one for every channel
983  // frequency
984  //
985  assert(ch_response.nelem() == 1 || ch_response.nelem() == ch_f.nelem());
986  //
987  Index freq_full = ch_response.nelem() > 1;
988 
989  // If response data extend outside sensor_f is checked in
990  // integration_func_by_vecmult
991 
992  // Reisze H
993  //
994  const Index nin_f = sensor_f.nelem();
995  const Index nout_f = ch_f.nelem();
996  const Index nin = n_sp * nin_f * n_pol;
997  const Index nout = n_sp * nout_f * n_pol;
998  //
999  H.resize(nout, nin);
1000 
1001  // Calculate the sensor integration vector and put values in the temporary
1002  // vector, then copy vector to the transfer matrix
1003  //
1004  Vector ch_response_f;
1005  Vector weights(nin_f);
1006  Vector weights_long(nin, 0.0);
1007  //
1008  for (Index ifr = 0; ifr < nout_f; ifr++) {
1009  const Index irp = ifr * freq_full;
1010 
1011  //The spectrometer response is shifted for each centre frequency step
1012  ch_response_f = ch_response[irp].get_numeric_grid(GFIELD1_F_GRID);
1013  ch_response_f += ch_f[ifr];
1014 
1015  // Call *integration_func_by_vecmult* and store it in the temp vector
1017  weights, ch_response[irp].data, ch_response_f, sensor_f);
1018 
1019  // Normalise if flag is set
1020  if (do_norm) weights /= weights.sum();
1021 
1022  // Loop over polarisation and spectra (viewing directions)
1023  // Weights change only with frequency
1024  for (Index sp = 0; sp < n_sp; sp++) {
1025  for (Index pol = 0; pol < n_pol; pol++) {
1026  // Distribute the compact weight vector into weight_long
1027  weights_long[Range(sp * nin_f * n_pol + pol, nin_f, n_pol)] = weights;
1028 
1029  // Insert temp_long into H at the correct row
1030  H.insert_row(sp * nout_f * n_pol + ifr * n_pol + pol, weights_long);
1031 
1032  // Reset weight_long to zero.
1033  weights_long = 0.0;
1034  }
1035  }
1036  }
1037 }
1038 
1039 
1040 
1042  const Index& stokes_dim,
1043  const Index& ipol_1based,
1044  const Numeric nv) {
1045  assert(w.nelem() == stokes_dim);
1046 
1047  if (ipol_1based < 1 || ipol_1based > 10)
1048  throw runtime_error("Valid polarization indices are 1 to 10 (1-based).");
1049 
1050  ArrayOfVector s2p(10);
1051  //
1052  s2p[0] = {1}; // I
1053  s2p[1] = {0, 1}; // Q
1054  s2p[2] = {0, 0, 1}; // U
1055  s2p[3] = {0, 0, 0, 1}; // V
1056  s2p[4] = {nv, nv}; // Iv
1057  s2p[5] = {nv, -nv}; // Ih
1058  s2p[6] = {nv, 0, nv}; // I+45
1059  s2p[7] = {nv, 0, -nv}; // I-45
1060  s2p[8] = {nv, 0, 0, nv}; // Ilhc
1061  s2p[9] = {nv, 0, 0, -nv}; // Irhc
1062 
1063  const Index l = s2p[ipol_1based - 1].nelem();
1064  if (l > stokes_dim) {
1065  ostringstream os;
1066  os << "You have selected polarization with 1-based index: " << ipol_1based
1067  << endl
1068  << "but this polarization demands stokes_dim >= " << l << endl
1069  << "while the actual values of stokes_dim is " << stokes_dim;
1070  throw std::runtime_error(os.str());
1071  }
1072 
1073  w[Range(0, l)] = s2p[ipol_1based - 1];
1074  if (l < stokes_dim) {
1075  w[Range(l, stokes_dim - l)] = 0;
1076  }
1077 }
1078 
1079 // Functions by Stefan, needed for HIRS:
1080 
1082 
1111  const Index nf = fmin.nelem();
1112  assert(fmax.nelem() == nf);
1113  assert(i >= 0 && i < nf);
1114  assert(j >= 0 && j < nf);
1115  assert(fmin[i] <= fmin[j]);
1116  assert(i < j);
1117 
1118  // There are three cases to consider:
1119  // a) The two channels are separate: fmax[i] < fmin[j]
1120  // b) They overlap: fmax[i] >= fmin[j]
1121  // c) j is inside i: fmax[i] > fmax[j]
1122 
1123  // In the easiest case (a), we do not have to do anything.
1124  if (fmax[i] >= fmin[j]) {
1125  // We are in case (b) or (c), so we know that we have to combine
1126  // the channels. The new minimum frequency is fmin[i]. The new
1127  // maximum frequency is the larger one of the two channels we
1128  // are combining:
1129  if (fmax[j] > fmax[i]) fmax[i] = fmax[j];
1130 
1131  // We now have to kick out element j from both vectors.
1132 
1133  // Number of elements behind j:
1134  Index n_behind = nf - 1 - j;
1135 
1136  Vector dummy = fmin;
1137  fmin.resize(nf - 1);
1138  fmin[Range(0, j)] = dummy[Range(0, j)];
1139  if (n_behind > 0) fmin[Range(j, n_behind)] = dummy[Range(j + 1, n_behind)];
1140 
1141  dummy = fmax;
1142  fmax.resize(nf - 1);
1143  fmax[Range(0, j)] = dummy[Range(0, j)];
1144  if (n_behind > 0) fmax[Range(j, n_behind)] = dummy[Range(j + 1, n_behind)];
1145 
1146  return true;
1147  }
1148 
1149  return false;
1150 }
1151 
1153  Vector& fmin,
1154  Vector& fmax,
1155  // Input:
1156  const Vector& f_backend,
1158  const Numeric& delta,
1159  const Verbosity& verbosity) {
1160  CREATE_OUT2;
1161 
1162  // How many channels in total:
1163  const Index n_chan = f_backend.nelem();
1164 
1165  // Checks on input quantities:
1166 
1167  // There must be at least one channel.
1168  if (n_chan < 1) {
1169  ostringstream os;
1170  os << "There must be at least one channel.\n"
1171  << "(The vector *f_backend* must have at least one element.)";
1172  throw runtime_error(os.str());
1173  }
1174 
1175  // There must be a response function for each channel.
1176  if (n_chan != backend_channel_response.nelem()) {
1177  ostringstream os;
1178  os << "Variables *f_backend_multi* and *backend_channel_response_multi*\n"
1179  << "must have same number of bands for each LO.";
1180  throw runtime_error(os.str());
1181  }
1182 
1183  // Frequency grids for response functions must be strictly increasing.
1184  for (Index i = 0; i < n_chan; ++i) {
1185  // Frequency grid for this response function:
1186  const Vector& backend_f_grid =
1187  backend_channel_response[i].get_numeric_grid(0);
1188 
1189  if (!is_increasing(backend_f_grid)) {
1190  ostringstream os;
1191  os << "The frequency grid for the backend channel response of\n"
1192  << "channel " << i << " is not strictly increasing.\n";
1193  os << "It is: " << backend_f_grid << "\n";
1194  throw runtime_error(os.str());
1195  }
1196  }
1197 
1198  // Start the actual work.
1199 
1200  out2 << " Original channel characteristics:\n"
1201  << " min nominal max (all in Hz):\n";
1202 
1203  // count the number of passbands as defined by segments of filtershapes, backend_channel_response.data = [0,>0,>0,0]
1204  // Borders between passbands are identified as [...0,0...]
1205 
1206  // Get a list of original channel boundaries:
1207  Index numPB = 0;
1208  for (Index idx = 0; idx < n_chan; ++idx) {
1209  const Vector& backend_filter = backend_channel_response[idx].data;
1210  if (backend_filter.nelem() >
1211  2) { // only run this code when there is more then two elements in the backend
1212  for (Index idy = 1; idy < backend_filter.nelem(); ++idy) {
1213  if ((backend_filter[idy] > 0) && (backend_filter[idy - 1] == 0)) {
1214  numPB++; // Two consecutive zeros gives the border between two passbands
1215  }
1216  }
1217  } else {
1218  numPB++;
1219  }
1220  }
1221 
1222  if (!numPB) {
1223  throw std::runtime_error(
1224  "No passbands found.\n"
1225  "*backend_channel_response* must be zero around the passbands.\n"
1226  "backend_channel_response.data = [0, >0, >0, 0]\n"
1227  "Borders between passbands are identified as [...0,0...]");
1228  }
1229 
1230  Vector fmin_pb(numPB);
1231  Vector fmax_pb(numPB);
1232  Index pbIdx = 0;
1233 
1234  for (Index idx = 0; idx < n_chan; ++idx) {
1235  // Some handy shortcuts:
1236  //
1237  // We have to find the first and last frequency where the
1238  // response is actually different from 0. (No point in making
1239  // calculations for frequencies where the response is 0.)
1240  // Index j=0;
1241  // while (backend_response[j] <= 0) ++j;
1242  // Numeric bf_min = backend_f_grid[j];
1243 
1244  // j=nf-1;
1245  // while (backend_response[j] <= 0) --j;
1246  // Numeric bf_max = backend_f_grid[j];
1247  //
1248  // No, aparently the sensor part want values also where the
1249  // response is zero. So we simply take the grid boundaries here.
1250  //
1251  // We need to add a bit of extra margin at both sides,
1252  // otherwise there is a numerical problem in the sensor WSMs.
1253  //
1254  // PE 081003: The accuracy for me (double on 64 bit machine) appears to
1255  // be about 3 Hz. Select 1 MHz to have a clear margin. Hopefully OK
1256  // for other machines.
1257  //
1258  // SAB 2010-04-14: The approach with a constant delta does not seem to work
1259  // well in practice. What I do now is that I add a delta corresponding to a
1260  // fraction of the grid spacing. But that is done outside of this function.
1261  // So we set delta = 0 here.
1262  //
1263  // SAB 2010-05-03: Now we pass delta as a parameter (with a default value of 0).
1264  //
1265  // Isoz 2013-05-21: Added methods to ignore areas between passbands
1266  //
1267  const Vector& backend_f_grid =
1268  backend_channel_response[idx].get_numeric_grid(0);
1269  const Vector& backend_filter = backend_channel_response[idx].data;
1270  if (backend_filter.nelem() >=
1271  4) // Is the passband frequency response given explicitly ? e.g. [0,>0,>0,0]
1272  {
1273  for (Index idy = 1; idy < backend_filter.nelem(); ++idy) {
1274  if (idy == 1) {
1275  fmin_pb[pbIdx] = f_backend[idx] + backend_f_grid[0];
1276  } else if ((backend_filter[idy] > 0) &&
1277  (backend_filter[idy - 1] == 0)) {
1278  fmin_pb[pbIdx] = f_backend[idx] + backend_f_grid[idy - 1] - delta;
1279  }
1280  if ((backend_filter[idy] == 0) && (backend_filter[idy - 1] > 0)) {
1281  fmax_pb[pbIdx] = f_backend[idx] + backend_f_grid[idy] + delta;
1282  out2 << " "
1283  << "fmin_pb " << fmin_pb[pbIdx] << " "
1284  << "f_backend" << f_backend[idx] << " "
1285  << "fmax_pb " << fmax_pb[pbIdx] << " "
1286  << "diff " << fmax_pb[pbIdx] - fmin_pb[pbIdx] << "\n";
1287  pbIdx++;
1288  }
1289  }
1290  fmax_pb[pbIdx - 1] = f_backend[idx] + last(backend_f_grid);
1291  } else // Or are the passbands given implicitly - such as the default for AMSUA and MHS
1292  {
1293  fmin_pb[pbIdx] = f_backend[idx] + backend_f_grid[0] - delta; //delta;
1294  fmax_pb[pbIdx] = f_backend[idx] +
1295  backend_f_grid[backend_f_grid.nelem() - 1] +
1296  delta;
1297  out2 << " "
1298  << "fmin_pb " << fmin_pb[pbIdx] << " "
1299  << "f_backend" << f_backend[pbIdx] << " "
1300  << "fmax_pb " << fmax_pb[pbIdx] << "\n";
1301  pbIdx++;
1302  }
1303  }
1304 
1305  // The problem is that channels may be overlapping. In that case, we
1306  // want to create a frequency grid that covers their entire range,
1307  // but we do not want to duplicate any frequencies.
1308 
1309  // To make matters worse, one or even several channels may be
1310  // completely inside another very broad channel.
1311 
1312  // Sort channels by frequency:
1313  // Caveat: A channel may be higher in
1314  // characteristic frequency f_backend, but also wider, so that it
1315  // has a lower minimum frequency fmin_orig. (This is the case for
1316  // some HIRS channels.) We sort by the minimum frequency here, not
1317  // by f_backend. This is necessary for function
1318  // test_and_merge_two_channels to work correctly.
1319  ArrayOfIndex isorted;
1320  get_sorted_indexes(isorted, fmin_pb);
1321 
1322  fmin.resize(numPB);
1323  fmax.resize(numPB);
1324  out2 << " resize numPb " << numPB << "\n";
1325  for (Index idx = 0; idx < numPB; ++idx) {
1326  fmin[idx] = fmin_pb[isorted[idx]];
1327  fmax[idx] = fmax_pb[isorted[idx]];
1328  }
1329  // We will be testing pairs of channels, and combine them if
1330  // possible. We have to test always only against the direct
1331  // neighbour. If that has no overlap, higher channels can not have
1332  // any either, due to the sorting by fmin.
1333  //
1334  // Note that fmin.nelem() changes, as the loop is
1335  // iterated. Nevertheless this is the correct stop condition.
1336  for (Index i = 0; i < fmin.nelem() - 1; ++i) {
1337  bool continue_checking = true;
1338  // The "i<fmin.nelem()" condition below is necessary, since
1339  // fmin.nelem() can decrease while the loop is executed, due to
1340  // merging.
1341  while (continue_checking && i < fmin.nelem() - 1) {
1342  continue_checking = test_and_merge_two_channels(fmin, fmax, i, i + 1);
1343 
1344  // Function returns true if merging has taken place.
1345  // In this case, we have to check again.
1346  }
1347  }
1348 
1349  out2 << " New channel characteristics:\n"
1350  << " min max (all in Hz):\n";
1351  for (Index i = 0; i < fmin.nelem(); ++i)
1352  out2 << " " << fmin[i] << " " << fmax[i] << "\n";
1353 }
1354 
1355 
1356 
1357 /*===========================================================================
1358  === Core integration and sum functions:
1359  ===========================================================================*/
1360 
1362  ConstVectorView f,
1363  ConstVectorView x_f_in,
1364  ConstVectorView x_g_in) {
1365  // Basic sizes
1366  const Index nf = x_f_in.nelem();
1367  const Index ng = x_g_in.nelem();
1368 
1369  // Asserts
1370  assert(h.nelem() == ng);
1371  assert(f.nelem() == nf);
1372  assert(is_increasing(x_f_in));
1373  assert(is_increasing(x_g_in) || is_decreasing(x_g_in));
1374  // More asserts below
1375 
1376  // End points of x_f
1377  Numeric xfmin = x_f_in[0];
1378  Numeric xfmax = x_f_in[nf - 1];
1379 
1380  // Handle possibly reversed x_g.
1381  Vector x_g;
1382  Index xg_reversed = 0;
1383  //
1384  if (is_decreasing(x_g_in)) {
1385  x_g = x_g_in[Range(ng - 1, ng, -1)];
1386  xg_reversed = 1;
1387  } else {
1388  x_g = x_g_in;
1389  }
1390  //
1391  assert(x_g[0] <= xfmin);
1392  assert(x_g[ng - 1] >= xfmax);
1393 
1394  // Normalise grids so x_f covers [0,1]
1395  const Numeric df = xfmax - xfmin;
1396  Vector x_f(nf);
1397  //
1398  for (Index i = 0; i < nf; i++) {
1399  x_f[i] = (x_f_in[i] - xfmin) / df;
1400  }
1401  for (Index i = 0; i < ng; i++) {
1402  x_g[i] = (x_g[i] - xfmin) / df;
1403  }
1404  xfmin = 0;
1405  xfmax = 1;
1406  // To test without normalisation, comment out lines above and use:
1407  //const Numeric df = 1;
1408  //const Vector x_f = x_f_in;
1409 
1410  // Create a reference grid vector, x_ref that containing the values
1411  // of x_f and x_g strictly sorted. Only g points inside the f range
1412  // are of concern.
1413  list<Numeric> l_x;
1414  for (Index it = 0; it < nf; it++) {
1415  l_x.push_back(x_f[it]);
1416  }
1417  for (Index it = 0; it < ng; it++) {
1418  if (x_g[it] > xfmin && x_g[it] < xfmax) {
1419  l_x.push_back(x_g[it]);
1420  }
1421  }
1422  l_x.sort();
1423  l_x.unique();
1424  //
1425  Vector x_ref(l_x.size());
1426  Index e = 0;
1427  for (list<Numeric>::iterator li = l_x.begin(); li != l_x.end(); li++) {
1428  x_ref[e] = *li;
1429  e++;
1430  }
1431 
1432  // Initiate output vector, with equal size as x_g, with zeros.
1433  // Start calculations
1434  h = 0.0;
1435  Index i_f = 0;
1436  Index i_g = 0;
1437  Numeric dx, a0, b0, c0, a1, b1, c1, x3, x2, x1;
1438  //
1439  for (Index i = 0; i < x_ref.nelem() - 1; i++) {
1440  // Find for what index in x_g (which is the same as for h) and f
1441  // calculation corresponds to
1442  while (x_g[i_g + 1] <= x_ref[i]) {
1443  i_g++;
1444  }
1445  while (x_f[i_f + 1] <= x_ref[i]) {
1446  i_f++;
1447  }
1448 
1449  // If x_ref[i] is out of x_f's range then that part of the integral is 0,
1450  // and no calculations should be done
1451  if (x_ref[i] >= xfmin && x_ref[i] < xfmax) {
1452  // Product of steps in x_f and x_g
1453  dx = (x_f[i_f + 1] - x_f[i_f]) * (x_g[i_g + 1] - x_g[i_g]);
1454 
1455  // Calculate a, b and c coefficients; h[i]=ax^3+bx^2+cx
1456  a0 = (f[i_f] - f[i_f + 1]) / 3.0;
1457  b0 = (-f[i_f] * (x_g[i_g + 1] + x_f[i_f + 1]) +
1458  f[i_f + 1] * (x_g[i_g + 1] + x_f[i_f])) /
1459  2.0;
1460  c0 = x_g[i_g + 1] * (f[i_f] * x_f[i_f + 1] - f[i_f + 1] * x_f[i_f]);
1461 
1462  a1 = -a0;
1463  b1 = (f[i_f] * (x_g[i_g] + x_f[i_f + 1]) -
1464  f[i_f + 1] * (x_g[i_g] + x_f[i_f])) /
1465  2.0;
1466  c1 = x_g[i_g] * (-f[i_f] * x_f[i_f + 1] + f[i_f + 1] * x_f[i_f]);
1467 
1468  x1 = x_ref[i + 1] - x_ref[i];
1469  // Simple way, but sensitive to numerical problems:
1470  //x2 = pow(x_ref[i+1],2) - pow(x_ref[i],2);
1471  //x3 = pow(x_ref[i+1],3) - pow(x_ref[i],3);
1472  // The same but a numerically better way:
1473  x2 = x1 * (2 * x_ref[i] + x1);
1474  x3 = x1 * (3 * x_ref[i] * (x_ref[i] + x1) + x1 * x1);
1475 
1476  // Calculate h[i] and h[i+1] increment
1477  // (df-factor to compensate for normalisation)
1478  h[i_g] += df * (a0 * x3 + b0 * x2 + c0 * x1) / dx;
1479  h[i_g + 1] += df * (a1 * x3 + b1 * x2 + c1 * x1) / dx;
1480  }
1481  }
1482 
1483  // Flip back if x_g was decreasing
1484  if (xg_reversed) {
1485  Vector tmp = h[Range(ng - 1, ng, -1)]; // Flip order
1486  h = tmp;
1487  }
1488 
1489  // The expressions are sensitive to numerical issues if two points in x_ref
1490  // are very close compared to the values in x_ref. A test trying to warn when
1491  // this happens:
1492  if (min(f) >= 0 && min(h) < -1e-15)
1493  throw runtime_error(
1494  "Significant negative response value obtained, "
1495  "despite sensor reponse is strictly positive. This "
1496  "indicates numerical problems. Is there any very "
1497  "small spacing of the sensor response grid?"
1498  "Please, send a report to Patrick if you see this!");
1499 }
1500 
1501 
1502 
1504  ConstVectorView x_g_in,
1505  const Numeric& limit1,
1506  const Numeric& limit2) {
1507  // Basic sizes
1508  const Index ng = x_g_in.nelem();
1509 
1510  // Asserts
1511  assert(ng > 1);
1512  assert(h.nelem() == ng);
1513  assert(limit1 <= limit2);
1514 
1515  // Handle possibly reversed x_g.
1516  Vector x_g;
1517  Index xg_reversed = 0;
1518  //
1519  if (is_decreasing(x_g_in)) {
1520  x_g = x_g_in[Range(ng - 1, ng, -1)];
1521  xg_reversed = 1;
1522  } else {
1523  x_g = x_g_in;
1524  }
1525  //
1526  assert(x_g[0] <= limit1);
1527  assert(x_g[ng - 1] >= limit2);
1528 
1529  // Handle extreme cases
1530  // Bin has zero width
1531  if (limit1 == limit2) {
1532  h = 0.0;
1533  return;
1534  }
1535  // Bin covers complete x_g:
1536  else if (limit1 == x_g[0] && limit2 == x_g[ng - 1]) {
1537  h[0] = (x_g[1] - x_g[0]) / 2.0;
1538  for (Index i = 1; i < ng - 1; i++) {
1539  h[i] = (x_g[i + 1] - x_g[i - 1]) / 2.0;
1540  }
1541  h[ng - 1] = (x_g[ng - 1] - x_g[ng - 2]) / 2.0;
1542  return;
1543  }
1544 
1545  // The general case
1546  Numeric x1 = 0,
1547  x2 =
1548  0; // End points of bin, inside basis range of present grid point
1549  for (Index i = 0; i < ng; i++) {
1550  bool inside = false;
1551 
1552  if (i == 0) {
1553  if (limit1 < x_g[1]) {
1554  inside = true;
1555  x1 = limit1;
1556  x2 = min(limit2, x_g[1]);
1557  }
1558  } else if (i == ng - 1) {
1559  if (limit2 > x_g[ng - 2]) {
1560  inside = true;
1561  x1 = max(limit1, x_g[ng - 2]);
1562  x2 = limit2;
1563  }
1564  } else {
1565  if ((limit1 < x_g[i + 1] && limit2 > x_g[i - 1]) ||
1566  (limit2 > x_g[i - 1] && limit1 < x_g[i + 1])) {
1567  inside = true;
1568  x1 = max(limit1, x_g[i - 1]);
1569  x2 = min(limit2, x_g[i + 1]);
1570  }
1571  }
1572 
1573  // h is zero if no overlap between bin and basis range of point
1574  if (!inside) {
1575  h[i] = 0.0;
1576  } else {
1577  // Lower part
1578  if (x1 < x_g[i]) {
1579  const Numeric r = 1.0 / (x_g[i] - x_g[i - 1]);
1580  const Numeric y1 = r * (x1 - x_g[i - 1]);
1581  const Numeric dx = min(x2, x_g[i]) - x1;
1582  const Numeric y2 = y1 + r * dx;
1583  h[i] = 0.5 * dx * (y1 + y2);
1584  } else {
1585  h[i] = 0.0;
1586  }
1587 
1588  // Upper part
1589  if (x2 > x_g[i]) {
1590  const Numeric r = 1.0 / (x_g[i + 1] - x_g[i]);
1591  const Numeric y2 = r * (x_g[i + 1] - x2);
1592  const Numeric dx = x2 - max(x1, x_g[i]);
1593  const Numeric y1 = y2 + r * dx;
1594  h[i] += 0.5 * dx * (y1 + y2);
1595  }
1596  }
1597  }
1598 
1599  // Flip back if x_g was decreasing
1600  if (xg_reversed) {
1601  Vector tmp = h[Range(ng - 1, ng, -1)]; // Flip order
1602  h = tmp;
1603  }
1604 }
1605 
1606 
1608  ConstVectorView f,
1609  ConstVectorView x_f,
1610  ConstVectorView x_g,
1611  const Numeric x1,
1612  const Numeric x2) {
1613  // Asserts
1614  assert(h.nelem() == x_g.nelem());
1615  assert(f.nelem() == x_f.nelem());
1616  assert(x_g[0] <= x_f[0]);
1617  assert(last(x_g) >= last(x_f));
1618  assert(x1 >= x_f[0]);
1619  assert(x2 >= x_f[0]);
1620  assert(x1 <= last(x_f));
1621  assert(x2 <= last(x_f));
1622 
1623  // Determine grid positions for point 1 (both with respect to f and g grids)
1624  // and interpolate response function.
1625  ArrayOfGridPos gp1g(1), gp1f(1);
1626  gridpos(gp1g, x_g, x1);
1627  gridpos(gp1f, x_f, x1);
1628  Matrix itw1(1, 2);
1629  interpweights(itw1, gp1f);
1630  Numeric f1;
1631  interp(f1, itw1, f, gp1f);
1632 
1633  // Same for point 2
1634  ArrayOfGridPos gp2g(1), gp2f(1);
1635  gridpos(gp2g, x_g, x2);
1636  gridpos(gp2f, x_f, x2);
1637  Matrix itw2(1, 2);
1638  interpweights(itw2, gp2f);
1639  Numeric f2;
1640  interp(f2, itw2, f, gp2f);
1641 
1642  //Initialise h at zero and store calculated weighting components
1643  h = 0.0;
1644  h[gp1g[0].idx] += f1 * gp1g[0].fd[1];
1645  h[gp1g[0].idx + 1] += f1 * gp1g[0].fd[0];
1646  h[gp2g[0].idx] += f2 * gp2g[0].fd[1];
1647  h[gp2g[0].idx + 1] += f2 * gp2g[0].fd[0];
1648 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Group::Sparse
Sparse Sparse
Definition: autoarts.h:99
ARTS::Var::sensor_response_dlos_grid
Matrix sensor_response_dlos_grid(Workspace &ws) noexcept
Definition: autoarts.h:6356
ARTS::Var::antenna_dlos
Matrix antenna_dlos(Workspace &ws) noexcept
Definition: autoarts.h:2354
find_effective_channel_boundaries
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:1152
gaussian_response_autogrid
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:606
ARTS::Var::lo
Numeric lo(Workspace &ws) noexcept
Definition: autoarts.h:4039
Sparse::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
b0
#define b0
Definition: legacy_continua.cc:21714
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
antenna2d_gridded_dlos
void antenna2d_gridded_dlos(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_dlos, const GriddedField4 &antenna_response, ConstMatrixView mblock_dlos, ConstVectorView f_grid, const Index n_pol)
antenna2d_interp_gridded_dlos
Definition: sensor.cc:199
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
spectrometer_matrix
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:975
integration_bin_by_vecmult
void integration_bin_by_vecmult(VectorView h, ConstVectorView x_g_in, const Numeric &limit1, const Numeric &limit2)
integration_bin_by_vecmult
Definition: sensor.cc:1503
ARTS::Var::backend_channel_response
ArrayOfGriddedField1 backend_channel_response(Workspace &ws) noexcept
Definition: autoarts.h:2548
antenna1d_matrix
void antenna1d_matrix(Sparse &H, const Index &antenna_dim, ConstVectorView antenna_dza, const GriddedField4 &antenna_response, ConstVectorView za_grid, ConstVectorView f_grid, const Index n_pol, const Index do_norm)
antenna1d_matrix
Definition: sensor.cc:58
antenna2d_interp_response
void antenna2d_interp_response(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_dlos, const GriddedField4 &antenna_response, ConstMatrixView mblock_dlos, ConstVectorView f_grid, const Index n_pol)
antenna2d_interp_response
Definition: sensor.cc:442
joker
const Joker joker
last
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:165
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
ARTS::Var::antenna_response
GriddedField4 antenna_response(Workspace &ws) noexcept
Definition: autoarts.h:2396
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
Sparse
The Sparse class.
Definition: matpackII.h:60
gaussian_response
void gaussian_response(Vector &y, const Vector &x, const Numeric &x0, const Numeric &fwhm)
gaussian_response
Definition: sensor.cc:629
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
DEBUG_ONLY
#define DEBUG_ONLY(...)
Definition: debug.h:36
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
ARTS::Var::antenna_dim
Index antenna_dim(Workspace &ws) noexcept
Definition: autoarts.h:2327
data
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Definition: arts_api_classes.cc:232
GriddedField1::data
Vector data
Definition: gridded_fields.h:234
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
ARTS::Var::iy_unit
String iy_unit(Workspace &ws) noexcept
Definition: autoarts.h:3856
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
x0
#define x0
Definition: lineshapemodel.h:413
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
Sparse::insert_row
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
matpackI.h
Implementation of Matrix, Vector, and such stuff.
Array
This can be used to make arrays out of anything.
Definition: array.h:108
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
is_decreasing
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
get_sorted_indexes
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:73
ARTS::Var::sensor_response_f
Vector sensor_response_f(Workspace &ws) noexcept
Definition: autoarts.h:6376
sensor.h
Sensor modelling functions.
met_mm_polarisation_hmatrix
void met_mm_polarisation_hmatrix(Sparse &H, const ArrayOfString &mm_pol, const Numeric dza, const Index stokes_dim, const String &iy_unit)
Calculate polarisation H-matrix.
Definition: sensor.cc:775
GriddedField::get_numeric_grid
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
Definition: gridded_fields.cc:87
min
#define min(a, b)
Definition: legacy_continua.cc:20628
test_and_merge_two_channels
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:1110
messages.h
Declarations having to do with the four output streams.
GFIELD4_ZA_GRID
const Index GFIELD4_ZA_GRID
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
_U_
#define _U_
Definition: config.h:183
my_basic_string< char >
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
summation_by_vecmult
void summation_by_vecmult(VectorView h, ConstVectorView f, ConstVectorView x_f, ConstVectorView x_g, const Numeric x1, const Numeric x2)
summation_by_vecmult
Definition: sensor.cc:1607
ConstVectorView::sum
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:53
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
GriddedField4
Definition: gridded_fields.h:340
ARTS::Var::f_backend
Vector f_backend(Workspace &ws) noexcept
Definition: autoarts.h:3418
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
max
#define max(a, b)
Definition: legacy_continua.cc:20629
ARTS::Var::sensor_response_f_grid
Vector sensor_response_f_grid(Workspace &ws) noexcept
Definition: autoarts.h:6392
a1
#define a1
Definition: complex.h:56
ARTS::Var::za_grid
Vector za_grid(Workspace &ws) noexcept
Definition: autoarts.h:7771
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
GFIELD1_F_GRID
const Index GFIELD1_F_GRID
integration_func_by_vecmult
void integration_func_by_vecmult(VectorView h, ConstVectorView f, ConstVectorView x_f_in, ConstVectorView x_g_in)
integration_func_by_vecmult
Definition: sensor.cc:1361
Sparse::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
Range
The range class.
Definition: matpackI.h:160
ARTS::Var::sensor_response_pol
ArrayOfIndex sensor_response_pol(Workspace &ws) noexcept
Definition: autoarts.h:6411
mixer_matrix
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:645
ARTS::Var::sensor_response_pol_grid
ArrayOfIndex sensor_response_pol_grid(Workspace &ws) noexcept
Definition: autoarts.h:6428
is_increasing
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
logic.h
Header file for logic.cc.
Sparse::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
sensor_aux_vectors
void sensor_aux_vectors(Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, ConstVectorView sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, ConstMatrixView sensor_response_dlos_grid)
sensor_aux_vectors
Definition: sensor.cc:938
stokes2pol
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Definition: sensor.cc:1041
Sparse::rw
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:91
x1
#define x1
Definition: lineshapemodel.h:414
GriddedField1
Definition: gridded_fields.h:196
matpackII.h
Header file for sparse matrices.
dx
#define dx
Definition: legacy_continua.cc:22180
NAT_LOG_2
const Numeric NAT_LOG_2
GFIELD4_AA_GRID
const Index GFIELD4_AA_GRID
PI
const Numeric PI
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ARTS::Var::sensor_response_dlos
Matrix sensor_response_dlos(Workspace &ws) noexcept
Definition: autoarts.h:6340
GFIELD4_F_GRID
const Index GFIELD4_F_GRID
Vector
The Vector class.
Definition: matpackI.h:860
ARTS::Var::aa_grid
Vector aa_grid(Workspace &ws) noexcept
Definition: autoarts.h:1717
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
sorting.h
Contains sorting routines.
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
DEG2RAD
const Numeric DEG2RAD
x2
#define x2
Definition: lineshapemodel.h:415
b1
#define b1
Definition: complex.h:57
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
arts.h
The global header file for ARTS.
mueller_rotation
void mueller_rotation(Sparse &H, const Index &stokes_dim, const Numeric &rotangle)
mueller_rotation
Definition: sensor.cc:745
GFIELD4_FIELD_NAMES
const Index GFIELD4_FIELD_NAMES
x3
#define x3
Definition: lineshapemodel.h:416