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