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