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