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