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