ARTS 2.5.0 (git: 9ee3ac6c)
m_sensor.cc
Go to the documentation of this file.
1/* Copyright (C) 2003-2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3 Stefan Buehler <sbuehler(at)ltu.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
20/*===========================================================================
21 === File description
22 ===========================================================================*/
23
35/*===========================================================================
36 === External declarations
37 ===========================================================================*/
38
39#include <algorithm>
40#include <cmath>
41#include <stdexcept>
42#include <string>
43#include "arts.h"
44#include "auto_md.h"
45#include "check_input.h"
47#include "m_select.h"
48#include "math_funcs.h"
49#include "messages.h"
50#include "ppath.h"
51#include "sensor.h"
52#include "sorting.h"
53#include "special_interp.h"
54#include "xml_io.h"
55
56extern const Numeric PI;
57extern const Numeric NAT_LOG_2;
58extern const Numeric DEG2RAD;
59extern const Numeric RAD2DEG;
60extern const Index GFIELD1_F_GRID;
61extern const Index GFIELD4_FIELD_NAMES;
62extern const Index GFIELD4_F_GRID;
63extern const Index GFIELD4_ZA_GRID;
64extern const Index GFIELD4_AA_GRID;
65extern const Numeric SPEED_OF_LIGHT;
66
67/*===========================================================================
68 === The functions (in alphabetical order)
69 ===========================================================================*/
70
71/* Workspace method: Doxygen documentation will be auto-generated */
73 Matrix& mblock_dlos_grid,
75 Matrix& antenna_dlos,
76 const Index& n_za_grid,
77 const Numeric& fwhm,
78 const Numeric& xwidth_si,
79 const Numeric& dx_si,
80 const Verbosity& verbosity) {
81 if (dx_si > xwidth_si)
82 throw runtime_error("It is demanded that dx_si <= xwidth_si.");
83
84 antenna_dim = 1;
85
86 antenna_dlos.resize(1, 1);
87 antenna_dlos(0, 0) = 0.0;
88
89 antenna_responseGaussian(r, fwhm, xwidth_si, dx_si, 0, verbosity);
90
91 // za grid for response
93 const Index nr = r_za_grid.nelem();
94
95 // Cumulative integral of response (factor /2 skipped, but does not matter)
96 Vector cumtrapz(nr);
97 cumtrapz[0] = 0;
98 for (Index i = 1; i < nr; i++) {
99 cumtrapz[i] = cumtrapz[i - 1] + r.data(0, 0, i - 1, 0) + r.data(0, 0, i, 0);
100 }
101
102 // Equally spaced vector between end points of cumulative sum
103 Vector csp;
104 nlinspace(csp, cumtrapz[0], cumtrapz[nr - 1], n_za_grid);
105
106 // Get mblock_za_grid by interpolation
107 mblock_dlos_grid.resize(n_za_grid, 1);
108 ArrayOfGridPos gp(n_za_grid);
109 gridpos(gp, cumtrapz, csp);
110 Matrix itw(n_za_grid, 2);
111 interpweights(itw, gp);
112 interp(mblock_dlos_grid(joker, 0), itw, r_za_grid, gp);
113}
114
115
116
117/* Workspace method: Doxygen documentation will be auto-generated */
119 Matrix& sensor_los,
120 Matrix& antenna_dlos,
121 Index& antenna_dim,
122 Matrix& mblock_dlos_grid,
123 const Index& atmosphere_dim,
124 const Verbosity& verbosity) {
125 // Sizes
126 const Index nmblock = sensor_pos.nrows();
127 const Index nant = antenna_dlos.nrows();
128
129 // Input checks
130 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
131 chk_if_in_range("antenna_dim", antenna_dim, 1, 2);
132 //
133 if (sensor_pos.ncols() != atmosphere_dim)
134 throw runtime_error(
135 "The number of columns of sensor_pos must be "
136 "equal to the atmospheric dimensionality.");
137 if (atmosphere_dim <= 2 && sensor_los.ncols() != 1)
138 throw runtime_error("For 1D and 2D, sensor_los shall have one column.");
139 if (atmosphere_dim == 3 && sensor_los.ncols() != 2)
140 throw runtime_error("For 3D, sensor_los shall have two columns.");
141 if (sensor_los.nrows() != nmblock) {
142 ostringstream os;
143 os << "The number of rows of sensor_pos and sensor_los must be "
144 << "identical, but sensor_pos has " << nmblock << " rows,\n"
145 << "while sensor_los has " << sensor_los.nrows() << " rows.";
146 throw runtime_error(os.str());
147 }
148 if (antenna_dim == 2 && atmosphere_dim < 3) {
149 throw runtime_error("If *antenna_dim* is 2, *atmosphere_dim* must be 3.");
150 }
151 if (antenna_dlos.empty()) throw runtime_error("*antenna_dlos* is empty.");
152 if (antenna_dlos.ncols() < 1 || antenna_dlos.ncols() > 2)
153 throw runtime_error("*antenna_dlos* must have one or 2 columns.");
154 if (atmosphere_dim < 3 && antenna_dlos.ncols() == 2)
155 throw runtime_error(
156 "*antenna_dlos* can only have two columns for 3D atmosphers.");
157
158 // Claculate new sensor_pos and sensor_los
159 const Matrix pos_copy = sensor_pos;
160 const Matrix los_copy = sensor_los;
161 //
162 sensor_pos.resize(nmblock * nant, pos_copy.ncols());
163 sensor_los.resize(nmblock * nant, los_copy.ncols());
164 //
165 for (Index ib = 0; ib < nmblock; ib++) {
166 for (Index ia = 0; ia < nant; ia++) {
167 const Index i = ib * nant + ia;
168
169 sensor_pos(i, joker) = pos_copy(ib, joker);
170 sensor_los(i, joker) = los_copy(ib, joker);
171
172 sensor_los(i, 0) += antenna_dlos(ia, 0);
173 if (antenna_dlos.ncols() == 2) sensor_los(i, 1) += antenna_dlos(ia, 1);
174 }
175 }
176
177 // Set other variables
178 AntennaOff(antenna_dim, mblock_dlos_grid, verbosity);
179 //
180 antenna_dlos.resize(1, 1);
181 antenna_dlos = 0;
182}
183
184
185
186/* Workspace method: Doxygen documentation will be auto-generated */
187void AntennaOff(Index& antenna_dim,
188 Matrix& mblock_dlos_grid,
189 const Verbosity& verbosity) {
191
192 out2 << " Sets the antenna dimensionality to 1.\n";
193 antenna_dim = 1;
194
195 out2 << " Sets *mblock_dlos_grid* to have one row with value 0.\n";
196 mblock_dlos_grid.resize(1, 1);
197 mblock_dlos_grid = 0;
198}
199
200
201
202/* Workspace method: Doxygen documentation will be auto-generated */
204 const Numeric& fwhm,
205 const Numeric& xwidth_si,
206 const Numeric& dx_si,
207 const Index& do_2d,
208 const Verbosity&) {
209 if (dx_si > xwidth_si)
210 throw runtime_error("It is demanded that dx_si <= xwidth_si.");
211
212 Vector x, y;
213 gaussian_response_autogrid(x, y, 0, fwhm, xwidth_si, dx_si);
214
215 r.set_name("Antenna response");
216
217 r.set_grid_name(0, "Polarisation");
218 r.set_grid(0, {"NaN"});
219
220 r.set_grid_name(1, "Frequency");
221 r.set_grid(1, Vector(1, -999));
222
223 r.set_grid_name(2, "Zenith angle");
224 r.set_grid(2, x);
225
226 // common for 1D and 2D
227 r.set_grid_name(3, "Azimuth angle");
228 const Index n = y.nelem();
229 //
230 if (do_2d) {
231 r.set_grid(3, x);
232 r.data.resize(1, 1, n, n);
233
234 // The code below follows the function *gaussian_response*
235 const Numeric si = fwhm / (2 * sqrt(2 * NAT_LOG_2));
236 const Numeric a = 1 / (si * sqrt(2 * PI));
237
238 for (Index z = 0; z < n; z++) {
239 for (Index b = 0; b < n; b++) {
240 r.data(0, 0, z, b) =
241 a * exp(-0.5 * pow(sqrt(x[z] * x[z] + x[b] * x[b]) / si, 2.0));
242 }
243 }
244 } else {
245 r.set_grid(3, Vector(1, 0));
246 r.data.resize(1, 1, n, 1);
247 r.data(0, 0, joker, 0) = y;
248 }
249}
250
251
252
253/* Workspace method: Doxygen documentation will be auto-generated */
255 const Numeric& leff,
256 const Numeric& xwidth_si,
257 const Numeric& dx_si,
258 const Index& nf,
259 const Numeric& fstart,
260 const Numeric& fstop,
261 const Index& do_2d,
262 const Verbosity& verbosity) {
263 if (dx_si > xwidth_si)
264 throw runtime_error("It is demanded that dx_si <= xwidth_si.");
265
266 // Calculate response for highest frequency, with xwidth_si scaled from
267 // fstart
268 Vector x, y;
269 Numeric fwhm = RAD2DEG * SPEED_OF_LIGHT / (leff * fstop);
271 x, y, 0, fwhm, (fstop / fstart) * xwidth_si, dx_si);
272
273 r.set_name("Antenna response");
274
275 r.set_grid_name(0, "Polarisation");
276 r.set_grid(0, {"NaN"});
277
278 Vector f_grid;
279 VectorNLogSpace(f_grid, nf, fstart, fstop, verbosity);
280 r.set_grid_name(1, "Frequency");
281 r.set_grid(1, f_grid);
282
283 r.set_grid_name(2, "Zenith angle");
284 r.set_grid(2, x);
285
286 // common for 1D and 2D
287 r.set_grid_name(3, "Azimuth angle");
288 const Index n = y.nelem();
289 //
290 if (do_2d) {
291 r.set_grid(3, x);
292 r.data.resize(1, nf, n, n);
293
294 for (Index i = 0; i < nf; i++) {
295 fwhm = RAD2DEG * SPEED_OF_LIGHT / (leff * f_grid[i]);
296 const Numeric si = fwhm / (2 * sqrt(2 * NAT_LOG_2));
297 const Numeric a = 1 / (si * sqrt(2 * PI));
298
299 for (Index z = 0; z < n; z++) {
300 for (Index b = 0; b < n; b++) {
301 r.data(0, i, z, b) =
302 a * exp(-0.5 * pow(sqrt(x[z] * x[z] + x[b] * x[b]) / si, 2.0));
303 }
304 }
305 }
306 } else {
307 r.set_grid(3, Vector(1, 0));
308 r.data.resize(1, nf, n, 1);
309 //
310 r.data(0, nf - 1, joker, 0) = y;
311 //
312 for (Index i = 0; i < nf - 1; i++) {
313 fwhm = RAD2DEG * SPEED_OF_LIGHT / (leff * f_grid[i]);
314 gaussian_response(y, x, 0, fwhm);
315 r.data(0, i, joker, 0) = y;
316 }
317 }
318}
319
320
321
322/* Workspace method: Doxygen documentation will be auto-generated */
324 const Numeric& resolution,
325 const Verbosity&) {
326 r.resize(1);
327 r[0].set_name("Backend channel response function");
328
329 Vector x(2);
330
331 r[0].set_grid_name(0, "Frequency");
332 x[1] = resolution / 2.0;
333 x[0] = -x[1];
334 r[0].set_grid(0, x);
335
336 r[0].data.resize(2);
337 r[0].data[0] = 1 / resolution;
338 r[0].data[1] = r[0].data[0];
339}
340
341
342
343/* Workspace method: Doxygen documentation will be auto-generated */
345 const Vector& fwhm,
346 const Vector& xwidth_si,
347 const Vector& dx_si,
348 const Verbosity&) {
349 if ((fwhm.nelem() != xwidth_si.nelem() && xwidth_si.nelem() != 1) ||
350 (fwhm.nelem() != dx_si.nelem() && dx_si.nelem() != 1)) {
351 std::ostringstream os;
352 os << "*xwidth_si* and *dx_si* must have one element or the same number of\n";
353 os << "elements as *fwhm*.";
354 throw std::runtime_error(os.str());
355 }
356
357 Index nchannels = fwhm.nelem();
358 r.resize(nchannels);
359
360 Vector x, y;
361 Numeric this_xwidth_si = xwidth_si[0];
362 Numeric this_dx_si = dx_si[0];
363 for (Index i = 0; i < nchannels; i++) {
364 if (xwidth_si.nelem() > 1) this_xwidth_si = xwidth_si[i];
365
366 if (dx_si.nelem() > 1) this_dx_si = dx_si[i];
367
368 gaussian_response_autogrid(x, y, 0, fwhm[i], this_xwidth_si, this_dx_si);
369
370 r[i].set_name("Backend channel response function");
371
372 r[i].set_grid_name(0, "Frequency");
373 r[i].set_grid(0, x);
374
375 const Index n = y.nelem();
376 r[i].data.resize(n);
377 for (Index j = 0; j < n; j++) r[i].data[j] = y[j];
378 }
379}
380
381
382
383/* Workspace method: Doxygen documentation will be auto-generated */
384void f_gridFromSensorAMSU( // WS Output:
385 Vector& f_grid,
386 // WS Input:
387 const Vector& lo,
388 const ArrayOfVector& f_backend,
389 const ArrayOfArrayOfGriddedField1& backend_channel_response,
390 // Control Parameters:
391 const Numeric& spacing,
392 const Verbosity& verbosity) {
395
396 // Find out how many channels we have in total:
397 // f_backend is an array of vectors, containing the band frequencies for each Mixer.
398 Index n_chan = 0;
399 for (Index i = 0; i < f_backend.nelem(); ++i)
400 for (Index s = 0; s < f_backend[i].nelem(); ++s) ++n_chan;
401
402 // Checks on input quantities:
403
404 // There must be at least one channel:
405 if (n_chan < 1) {
406 ostringstream os;
407 os << "There must be at least one channel.\n"
408 << "(The vector *lo* must have at least one element.)";
409 throw runtime_error(os.str());
410 }
411
412 // Is number of LOs consistent in all input variables?
413 if ((f_backend.nelem() != lo.nelem()) ||
414 (backend_channel_response.nelem() != lo.nelem())) {
415 ostringstream os;
416 os << "Variables *lo_multi*, *f_backend_multi* and *backend_channel_response_multi*\n"
417 << "must have same number of elements (number of LOs).";
418 throw runtime_error(os.str());
419 }
420
421 // Is number of bands consistent for each LO?
422 for (Index i = 0; i < f_backend.nelem(); ++i)
423 if (f_backend[i].nelem() != backend_channel_response[i].nelem()) {
424 ostringstream os;
425 os << "Variables *f_backend_multi* and *backend_channel_response_multi*\n"
426 << "must have same number of bands for each LO.";
427 throw runtime_error(os.str());
428 }
429
430 // Start the actual work.
431
432 // We construct the necessary input for function find_effective_channel_boundaries,
433 // which will identify channel boundaries, taking care of overlaping channels.
434
435 // A "flat" vector of nominal band frequencies (two for each AMSU channel):
436 Vector f_backend_flat(2 * n_chan);
437
438 // A "flat" list of channel response functions (two for each AMSU channel)
439 ArrayOfGriddedField1 backend_channel_response_flat(2 * n_chan);
440
441 // Counts position inside the flat arrays during construction:
442 Index j = 0;
443
444 for (Index i = 0; i < f_backend.nelem(); ++i)
445 for (Index s = 0; s < f_backend[i].nelem(); ++s) {
446 const GriddedField1& this_grid = backend_channel_response[i][s];
447 const Numeric this_f_backend = f_backend[i][s];
448
449 // Signal sideband:
450 f_backend_flat[j] = this_f_backend;
451 backend_channel_response_flat[j] = this_grid;
452 j++;
453
454 // Image sideband:
455 Numeric offset = this_f_backend - lo[i];
456 Numeric f_image = lo[i] - offset;
457 f_backend_flat[j] = f_image;
458 backend_channel_response_flat[j] = this_grid;
459 j++;
460 }
461
462 // We build up a total list of absolute frequency ranges for
463 // both signal and image sidebands:
464 Vector fmin(2 * n_chan), fmax(2 * n_chan);
465
466 // We have to add some additional margin at the band edges,
467 // otherwise the instrument functions are not happy. Define
468 // this in terms of the grid spacing:
469 Numeric delta = 1 * spacing;
470
471 // Call subfunction to do the actual work of merging overlapping
472 // channels and identifying channel boundaries:
474 fmax,
475 f_backend_flat,
476 backend_channel_response_flat,
477 delta,
478 verbosity);
479
480 // Create f_grid_array. This is an array of Numeric, so that we
481 // can use the STL push_back function.
482 ArrayOfNumeric f_grid_array;
483
484 for (Index i = 0; i < fmin.nelem(); ++i) {
485 // Band width:
486 const Numeric bw = fmax[i] - fmin[i];
487
488 // How many grid intervals do I need?
489 const Numeric npf = ceil(bw / spacing);
490
491 // How many grid points to store? - Number of grid intervals
492 // plus 1.
493 const Index npi = (Index)npf + 1;
494
495 // Create the grid for this band:
496 Vector grid;
497 nlinspace(grid, fmin[i], fmax[i], npi);
498
499 out3 << " Band range " << i << ": " << grid << "\n";
500
501 // Append to f_grid_array:
502 f_grid_array.reserve(f_grid_array.nelem() + npi);
503 for (Index s = 0; s < npi; ++s) f_grid_array.push_back(grid[s]);
504 }
505
506 // Copy result to output vector:
507 f_grid = f_grid_array;
508
509 out2 << " Total number of frequencies in f_grid: " << f_grid.nelem() << "\n";
510 // cout << "f_grid: " << f_grid << "\n";
511}
512
513
514
515/* Workspace method: Doxygen documentation will be auto-generated */
517 Vector& f_grid,
518 // WS Input:
519 const ArrayOfVector& f_backend_multi, // Center frequency for each channel
520 const ArrayOfArrayOfGriddedField1& backend_channel_response_multi,
521 // Control Parameters:
522 const Numeric& spacing,
523 const Vector& verbosityVect,
524 const Verbosity& verbosity) {
527 Index numFpoints = 0;
528 // Calculate the total number of frequency points
529 for (Index idx = 0; idx < backend_channel_response_multi.nelem(); ++idx) {
530 for (Index idy = 0; idy < backend_channel_response_multi[idx].nelem();
531 ++idy) {
532 numFpoints += backend_channel_response_multi[idx][idy].get_grid_size(0);
533 }
534 }
535
536 Index numChan = backend_channel_response_multi.nelem(); // number of channels
537
538 // Checks on input quantities:
539
540 // There must be at least one channel:
541 if (numChan < 1) {
542 ostringstream os;
543 os << "There must be at least one channel.\n"
544 << "(The vector *lo* must have at least one element.)";
545 throw runtime_error(os.str());
546 }
547
548 // Start the actual work.
549 // We construct the necessary input for function find_effective_channel_boundaries,
550 // which will identify channel boundaries, taking care of overlaping channels.
551
552 // A "flat" vector of nominal band frequencies (one for each AMSU channel):
553 Vector f_backend_flat(numChan);
554 // A "flat" list of channel response functions (one for each AMSU channel)
555 ArrayOfGriddedField1 backend_channel_response_nonflat(numChan);
556
557 // Save the correct verbosity value to each passband
558 Vector FminVerbosityVect(numFpoints);
559 Vector FmaxVerbosityVect(numFpoints);
560 Vector VerbosityValVect(numFpoints);
561 Index VerbVectIdx = 0;
562
563 for (Index i = 0; i < f_backend_multi.nelem(); ++i) {
564 for (Index ii = 0; ii < f_backend_multi[i].nelem(); ++ii) {
565 const GriddedField1& this_grid = backend_channel_response_multi[i][ii];
566 const Numeric this_f_backend = f_backend_multi[i][ii];
567 // Signal passband :
568 f_backend_flat[i] = this_f_backend;
569 backend_channel_response_nonflat[i] = this_grid;
570 for (Index idy = 1;
571 idy < backend_channel_response_multi[i][ii].get_grid_size(0);
572 ++idy) {
573 if ((backend_channel_response_multi[i][ii].data[idy - 1] == 0) &&
574 (backend_channel_response_multi[i][ii].data[idy] > 0)) {
575 FminVerbosityVect[VerbVectIdx] =
576 f_backend_multi[i][ii] +
577 backend_channel_response_multi[i][ii].get_numeric_grid(0)[idy];
578 VerbosityValVect[VerbVectIdx] = verbosityVect[i];
579 }
580 if ((backend_channel_response_multi[i][ii].data[idy - 1] > 0) &&
581 (backend_channel_response_multi[i][ii].data[idy] == 0)) {
582 FmaxVerbosityVect[VerbVectIdx] =
583 f_backend_multi[i][ii] +
584 backend_channel_response_multi[i][ii].get_numeric_grid(0)[idy];
585 VerbVectIdx++;
586 }
587 }
588 }
589 }
590 // We build up a total list of absolute frequency ranges for all passbands
591 // Code reused from the function "f_gridFromSensorAMSU"
592 Vector fmin,
593 fmax; // - these variables will be resized, therefore len(1) is enough for now.,
594
595 // We have to add some additional margin at the band edges,
596 // otherwise the instrument functions are not happy. Define
597 // this in terms of the grid spacing:
598 const Numeric delta = 10;
599
600 // Call subfunction to do the actual work of merging overlapping
601 // channels and identifying channel boundaries:
603 fmax,
604 f_backend_flat,
605 backend_channel_response_nonflat,
606 delta,
607 verbosity);
608
609 // Create f_grid_array. This is an array of Numeric, so that we
610 // can use the STL push_back function.
611 ArrayOfNumeric f_grid_array;
612
613 for (Index i = 0; i < fmin.nelem(); ++i) {
614 // Bandwidth:
615 const Numeric bw = fmax[i] - fmin[i];
616 Numeric npf = ceil(bw / spacing); // Set a default value
617
618 // How many grid intervals do I need?
619 Index verbIdx = 0;
620 if (verbosityVect.nelem() > 0) {
621 // find the grid needed for the particular part of passband
622 for (Index ii = 0; ii < VerbVectIdx; ++ii) {
623 if ((FminVerbosityVect[ii] >= fmin[i]) &&
624 (FmaxVerbosityVect[ii] <= fmax[i])) {
625 if (verbIdx == 0) {
626 verbIdx = ii;
627 } else {
628 if (VerbosityValVect[ii] < VerbosityValVect[verbIdx]) {
629 verbIdx = ii;
630 }
631 }
632 }
633 }
634 if (spacing > VerbosityValVect[verbIdx]) {
635 npf = ceil(
636 bw / VerbosityValVect[verbIdx]); // is the default value to coarse?
637 } else {
638 npf = ceil(bw / spacing); // Default value
639 }
640 }
641
642 // How many grid points to store? - Number of grid intervals
643 // plus 1.
644 const Index npi = (Index)npf + 1;
645
646 // What is the actual grid spacing inside the band?
647 const Numeric gs = bw / npf;
648
649 // Create the grid for this band:
650 Vector grid(fmin[i], npi, gs);
651
652 out3 << " Band range " << i << ": " << grid << "\n";
653
654 // Append to f_grid_array:
655 f_grid_array.reserve(f_grid_array.nelem() + npi);
656 for (Index s = 0; s < npi; ++s) f_grid_array.push_back(grid[s]);
657 }
658
659 // Copy result to output vector:
660 f_grid = f_grid_array;
661
662 out2 << " Total number of frequencies in f_grid: " << f_grid.nelem() << "\n";
663}
664
665
666
667/* Workspace method: Doxygen documentation will be auto-generated */
668void f_gridFromSensorHIRS( // WS Output:
669 Vector& f_grid,
670 // WS Input:
671 const Vector& f_backend,
672 const ArrayOfGriddedField1& backend_channel_response,
673 // Control Parameters:
674 const Numeric& spacing,
675 const Verbosity& verbosity) {
678
679 // Check input
680 if (spacing <= 0) {
681 ostringstream os;
682 os << "Expected positive spacing. Found spacing to be: " << spacing << "\n";
683 throw runtime_error(os.str());
684 }
685 // Call subfunction to get channel boundaries. Also does input
686 // consistency checking for us.
687 Vector fmin, fmax;
688
689 // We have to add some additional margin at the band edges,
690 // otherwise the instrument functions are not happy. Define
691 // this in terms of the grid spacing:
692 Numeric delta = 1 * spacing;
693
695 fmin, fmax, f_backend, backend_channel_response, delta, verbosity);
696
697 // Ok, now we just have to create a frequency grid for each of the
698 // fmin/fmax ranges.
699
700 // Create f_grid_array. This is an array of Numeric, so that we
701 // can use the STL push_back function.
702 ArrayOfNumeric f_grid_array;
703
704 for (Index i = 0; i < fmin.nelem(); ++i) {
705 // Band width:
706 const Numeric bw = fmax[i] - fmin[i];
707
708 // How many grid intervals do I need?
709 const Numeric npf = ceil(bw / spacing);
710
711 // How many grid points to store? - Number of grid intervals
712 // plus 1.
713 const Index npi = (Index)npf + 1;
714
715 // Create the grid for this band:
716 Vector grid;
717 nlinspace(grid, fmin[i], fmax[i], npi);
718
719 out3 << " Band range " << i << ": " << grid << "\n";
720
721 // Append to f_grid_array:
722 f_grid_array.reserve(f_grid_array.nelem() + npi);
723 for (Index s = 0; s < npi; ++s) f_grid_array.push_back(grid[s]);
724 }
725
726 // Copy result to output vector:
727 f_grid = f_grid_array;
728
729 out2 << " Total number of frequencies in f_grid: " << f_grid.nelem() << "\n";
730}
731
732
733
734/* Workspace method: Doxygen documentation will be auto-generated */
736 // WS Output:
737 Vector& f_grid,
738 Vector& f_backend,
739 ArrayOfArrayOfIndex& channel2fgrid_indexes,
740 ArrayOfVector& channel2fgrid_weights,
741 // WS Input:
742 const Matrix& mm_back, /* met_mm_backend */
743 // Control Parameters:
744 const Vector& freq_spacing,
745 const ArrayOfIndex& freq_number,
746 const Numeric& freq_merge_threshold,
747 const Verbosity&) {
748 // Some sizes
749 const Index nchannels = mm_back.nrows();
750
751 // Checks of input
752 //
753 chk_met_mm_backend(mm_back);
754 //
755 if (freq_spacing.nelem() != 1 && freq_spacing.nelem() != nchannels)
756 throw std::runtime_error(
757 "Length of *freq_spacing* vector must be either 1 or correspond\n"
758 "to the number of rows in *met_mm_backend*.");
759 //
760 if (freq_number.nelem() != 1 && freq_number.nelem() != nchannels)
761 throw std::runtime_error(
762 "Length of *freq_number* array must be either 1 or correspond\n"
763 "to the number of rows in *met_mm_backend*.");
764 //
765 if (freq_merge_threshold <= 0)
766 throw std::runtime_error("The *freq_merge_threshold* must be > 0.\n");
767 //
768 if (freq_merge_threshold > 100.)
769 throw std::runtime_error(
770 "The *freq_merge_threshold* is only meant to merge frequencies\n"
771 "that are basically identical and only differ slightly due to\n"
772 "numerical inaccuracies. Setting it to >100Hz is not reasonable.");
773
774 ArrayOfNumeric f_grid_unsorted;
775 ArrayOfIndex nf_per_channel(nchannels);
776 ArrayOfIndex index_in_unsorted;
777
778 f_backend.resize(nchannels);
779
780 for (Index ch = 0; ch < nchannels; ch++) {
781 const Numeric lo = mm_back(ch, 0);
782 const Numeric offset1 = mm_back(ch, 1);
783 const Numeric offset2 = mm_back(ch, 2);
784 const Numeric bandwidth = mm_back(ch, 3);
785
786 const Index this_fnumber =
787 (freq_number.nelem() == 1) ? freq_number[0] : freq_number[ch];
788 const Numeric this_spacing =
789 (freq_spacing.nelem() == 1) ? freq_spacing[0] : freq_spacing[ch];
790
791 if (this_spacing <= 0)
792 throw std::runtime_error("*freq_spacing must be > 0.");
793
794 if (this_fnumber == 0) {
795 ostringstream os;
796 os << "*freq_number* must be -1 or greater zero:"
797 << "freq_number[" << ch << "] = " << this_fnumber;
798 std::runtime_error(os.str());
799 }
800
801 // Number of passbands
802 const Index npassb =
803 1 + ((Index)(offset1 > 0)) + (2 * (Index)(offset2 > 0));
804
805 // Number of frequencies per passband
806 Index nfperband = this_fnumber;
807 //
808 if (this_fnumber == -1 ||
809 bandwidth / (Numeric)this_fnumber > this_spacing) {
810 nfperband = (Index)ceil(bandwidth / this_spacing);
811 }
812
813 // Fill the data known so far
814 nf_per_channel[ch] = npassb * nfperband;
815 f_backend[ch] = lo;
816
817 // Distance between each new f_grid value
818 const Numeric df = bandwidth / (Numeric)nfperband;
819
820 // Loop passbands and determine f_grid values
821 for (Index b = 0; b < npassb; b++) {
822 // Centre frequency of passband
823 Numeric fc = lo;
824 if (npassb == 2) {
825 fc += (-1 + 2 * (Numeric)b) * offset1;
826 } else if (npassb == 4) {
827 if (b <= 1) {
828 fc -= offset1;
829 } else {
830 fc += offset1;
831 }
832 if (b == 0 || b == 2) {
833 fc -= offset2;
834 } else {
835 fc += offset2;
836 }
837 }
838
839 // Loop frequencies to add
840 for (Index f_index = 0; f_index < nfperband; f_index++) {
841 const Numeric fnew = fc - bandwidth / 2 + (0.5 + (Numeric)f_index) * df;
842
843 // Does this frequency exist or not?
844 bool found = false;
845 for (size_t f_try = 0; f_try < f_grid_unsorted.size(); f_try++) {
846 if (abs(fnew - f_grid_unsorted[f_try]) < freq_merge_threshold) {
847 found = true;
848 index_in_unsorted.push_back(f_try);
849 break;
850 }
851 }
852 if (!found) {
853 f_grid_unsorted.push_back(fnew);
854 index_in_unsorted.push_back(f_grid_unsorted.size() - 1);
855 }
856 }
857 }
858 }
859
860 // Determine sort order for f_grid
861 const size_t nf = f_grid_unsorted.size();
862 ArrayOfIndex move2index(nf);
863
864 // Create frequency position vector (1...nf)
865 ArrayOfIndex sorted_indices;
866 sorted_indices.resize(nf);
867 for (size_t i = 0; i < nf; i++) {
868 sorted_indices[i] = i;
869 }
870
871 // Sort frequency position vector by frequency
872 std::sort(sorted_indices.begin(),
873 sorted_indices.end(),
874 CmpArrayOfNumeric(f_grid_unsorted));
875
876 // Create vector with indices in target vector
877 for (size_t i = 0; i < nf; i++) {
878 move2index[sorted_indices[i]] = i;
879 }
880
881 // Create f_grid
882 f_grid.resize(nf);
883 //
884 for (size_t f_index = 0; f_index < nf; f_index++) {
885 f_grid[move2index[f_index]] = f_grid_unsorted[f_index];
886 }
887
888 // Create channel2 fgrid variables
889 channel2fgrid_indexes.resize(nchannels);
890 channel2fgrid_weights.resize(nchannels);
891 Index i = 0;
892 for (Index ch = 0; ch < nchannels; ch++) {
893 channel2fgrid_indexes[ch].resize(nf_per_channel[ch]);
894 channel2fgrid_weights[ch].resize(nf_per_channel[ch]);
895 //
896 for (Index j = 0; j < nf_per_channel[ch]; j++) {
897 channel2fgrid_indexes[ch][j] = move2index[index_in_unsorted[i]];
898 channel2fgrid_weights[ch][j] = 1 / (Numeric)nf_per_channel[ch];
899 i++;
900 }
901 }
902}
903
904
905
906/* Workspace method: Doxygen documentation will be auto-generated */
908 const Numeric& spacing,
909 const Numeric& width,
910 const Index& centre,
911 const Verbosity&) {
912 // Create linear, equidistant grid (same in zenith and azimuth)
913 Vector x;
914 Numeric w;
915 if (centre) {
916 w = spacing * ceil(width / spacing);
917 } else {
918 w = spacing * (0.5 + floor(width / spacing));
919 }
920 linspace(x, -w, w, spacing);
921 //
922 const Index l = x.nelem();
923
924 Matrix dlos_try(l * l, 2, 0);
925 Index n_in = 0;
926 const Numeric c = width * width;
927
928 for (Index i = 0; i < l; i++) {
929 const Numeric a = x[i] * x[i];
930
931 for (Index j = 0; j < l; j++) {
932 if (a + x[j] * x[j] <= c) {
933 dlos_try(n_in, 0) = x[i];
934 dlos_try(n_in, 1) = x[j];
935 n_in++;
936 }
937 }
938 }
939
940 mblock_dlos_grid = dlos_try(Range(0, n_in), joker);
941}
942
943
944
945/* Workspace method: Doxygen documentation will be auto-generated */
947 const Numeric& spacing,
948 const Numeric& za_width,
949 const Numeric& aa_width,
950 const Index& centre,
951 const Verbosity&) {
952 // Create za-grid
953 Vector za;
954 Numeric w;
955 if (centre) {
956 w = spacing * ceil(za_width / spacing);
957 } else {
958 w = spacing * (0.5 + floor(za_width / spacing));
959 }
960 linspace(za, -w, w, spacing);
961
962 // Create aa-grid
963 Vector aa;
964 if (centre) {
965 w = spacing * ceil(aa_width / spacing);
966 } else {
967 w = spacing * (0.5 + floor(aa_width / spacing));
968 }
969 linspace(aa, -w, w, spacing);
970
971 const Index nza = za.nelem();
972 const Index naa = aa.nelem();
973
974 mblock_dlos_grid.resize(nza * naa, 2);
975
976 Index n = 0;
977
978 for (Index a = 0; a < naa; a++) {
979 for (Index z = 0; z < nza; z++) {
980 mblock_dlos_grid(n, 0) = za[z];
981 mblock_dlos_grid(n, 1) = aa[a];
982 n++;
983 }
984 }
985}
986
987
988
989/* Workspace method: Doxygen documentation will be auto-generated */
990void sensor_responseAntenna(Sparse& sensor_response,
991 Vector& sensor_response_f,
992 ArrayOfIndex& sensor_response_pol,
993 Matrix& sensor_response_dlos,
994 Matrix& sensor_response_dlos_grid,
995 const Vector& sensor_response_f_grid,
996 const ArrayOfIndex& sensor_response_pol_grid,
997 const Index& atmosphere_dim,
998 const Index& antenna_dim,
999 const Matrix& antenna_dlos,
1000 const GriddedField4& antenna_response,
1001 const Index& sensor_norm,
1002 const String& option_2d,
1003 const Verbosity& verbosity) {
1005
1006 // Basic checks
1007 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1008 chk_if_in_range("antenna_dim", antenna_dim, 1, 2);
1009 chk_if_bool("sensor_norm", sensor_norm);
1010
1011 // Some sizes
1012 const Index nf = sensor_response_f_grid.nelem();
1013 const Index npol = sensor_response_pol_grid.nelem();
1014 const Index nlos = sensor_response_dlos_grid.nrows();
1015 const Index nin = nf * npol * nlos;
1016
1017 // Initialise a output stream for runtime errors and a flag for errors
1018 ostringstream os;
1019 bool error_found = false;
1020
1021 // Check that sensor_response variables are consistent in size
1022 if (sensor_response_f.nelem() != nin) {
1023 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
1024 << "grid variables (sensor_response_f_grid etc.).\n";
1025 error_found = true;
1026 }
1027 if (sensor_response.nrows() != nin) {
1028 os << "The sensor block response matrix *sensor_response* does not have\n"
1029 << "right size compared to the sensor grid variables\n"
1030 << "(sensor_response_f_grid etc.).\n";
1031 error_found = true;
1032 }
1033
1034 // Checks related to antenna dimension
1035 if (antenna_dim == 2 && atmosphere_dim < 3) {
1036 os << "If *antenna_dim* is 2, *atmosphere_dim* must be 3.\n";
1037 error_found = true;
1038 }
1039
1040 // Basic checks of antenna_dlos
1041 if (antenna_dlos.empty()) throw runtime_error("*antenna_dlos* is empty.");
1042 if (antenna_dlos.ncols() < 1 || antenna_dlos.ncols() > 2)
1043 throw runtime_error("*antenna_dlos* must have one or 2 columns.");
1044 if (atmosphere_dim < 3 && antenna_dlos.ncols() == 2)
1045 throw runtime_error(
1046 "*antenna_dlos* can only have two columns for 3D atmosphers.");
1047
1048 // We allow angles in antenna_los to be unsorted
1049
1050 // Checks of antenna_response polarisation dimension
1051 //
1052 const Index lpolgrid =
1053 antenna_response.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
1054 //
1055 if (lpolgrid != 1 && lpolgrid != npol) {
1056 os << "The number of polarisation in *antenna_response* must be 1 or be\n"
1057 << "equal to the number of polarisations used (determined by\n"
1058 << "*stokes_dim* or *instrument_pol*).\n";
1059 error_found = true;
1060 }
1061
1062 // Checks of antenna_response frequency dimension
1063 //
1064 ConstVectorView aresponse_f_grid =
1065 antenna_response.get_numeric_grid(GFIELD4_F_GRID);
1066 //
1067 chk_if_increasing("f_grid of antenna_response", aresponse_f_grid);
1068 //
1069 Numeric f_dlow = 0.0;
1070 Numeric f_dhigh = 0.0;
1071 //
1072 f_dlow = min(sensor_response_f_grid) - aresponse_f_grid[0];
1073 f_dhigh = last(aresponse_f_grid) - max(sensor_response_f_grid);
1074 //
1075 if (aresponse_f_grid.nelem() > 1) {
1076 if (f_dlow < 0) {
1077 os << "The frequency grid of *antenna_response is too narrow. It must\n"
1078 << "cover all considered frequencies (*f_grid*), if the length\n"
1079 << "is > 1. The grid needs to be expanded with " << -f_dlow
1080 << " Hz in\n"
1081 << "the lower end.\n";
1082 error_found = true;
1083 }
1084 if (f_dhigh < 0) {
1085 os << "The frequency grid of *antenna_response is too narrow. It must\n"
1086 << "cover all considered frequencies (*f_grid*), if the length\n"
1087 << "is > 1. The grid needs to be expanded with " << -f_dhigh
1088 << " Hz in\n"
1089 << "the upper end.\n";
1090 error_found = true;
1091 }
1092 }
1093
1094 // Checks of antenna_response za dimension
1095 //
1096 ConstVectorView aresponse_za_grid =
1097 antenna_response.get_numeric_grid(GFIELD4_ZA_GRID);
1098 //
1099 chk_if_increasing("za_grid of *antenna_response*", aresponse_za_grid);
1100 //
1101 if (aresponse_za_grid.nelem() < 2) {
1102 os << "The zenith angle grid of *antenna_response* must have >= 2 values.\n";
1103 error_found = true;
1104 }
1105
1106 // Checks of antenna_response aa dimension
1107 //
1108 ConstVectorView aresponse_aa_grid =
1109 antenna_response.get_numeric_grid(GFIELD4_AA_GRID);
1110 //
1111 if (antenna_dim == 1) {
1112 if (aresponse_aa_grid.nelem() != 1) {
1113 os << "The azimuthal dimension of *antenna_response* must be 1 if\n"
1114 << "*antenna_dim* equals 1.\n";
1115 error_found = true;
1116 }
1117 } else {
1118 chk_if_increasing("aa_grid of antenna_response", aresponse_aa_grid);
1119 //
1120 if (aresponse_za_grid.nelem() < 2) {
1121 os << "The zenith angle grid of *antenna_response* must have >= 2\n"
1122 << "values.\n";
1123 error_found = true;
1124 }
1125 }
1126
1127 // Check of angular grids. These checks differ with antenna_dim
1128 if (antenna_dim == 1) {
1129 if (!(is_increasing(sensor_response_dlos_grid(joker, 0)) ||
1130 is_decreasing(sensor_response_dlos_grid(joker, 0)))) {
1131 os << "For 1D antennas, the zenith angles in *sensor_response_dlos_grid*\n"
1132 << "must be sorted, either in increasing or decreasing order.\n"
1133 << "The original problem is probably found in *mblock_dlos_grid*.\n";
1134 error_found = true;
1135 }
1136
1137 else {
1138 // Check if the za relative grid is outside sensor_response_dlos_grid.
1139 Numeric za_dlow = 0.0;
1140 Numeric za_dhigh = 0.0;
1141 //
1142 za_dlow = antenna_dlos(0, 0) + aresponse_za_grid[0] -
1143 min(sensor_response_dlos_grid(joker, 0));
1144 za_dhigh = max(sensor_response_dlos_grid(joker, 0)) -
1145 (last(antenna_dlos(joker, 0)) + last(aresponse_za_grid));
1146 //
1147 if (za_dlow < 0) {
1148 os << "The WSV zenith angle part of *sensor_response_dlos_grid* is too narrow.\n"
1149 << "It should be expanded with " << -za_dlow
1150 << " deg in the lower end.\n"
1151 << "This change should be probably applied to *mblock_dlos_grid*.\n";
1152 error_found = true;
1153 }
1154 if (za_dhigh < 0) {
1155 os << "The WSV zenith angle part of *sensor_response_dlos_grid* is too narrow.\n"
1156 << "It should be expanded with " << -za_dhigh
1157 << " deg in the upper end.\n"
1158 << "This change should be probably applied to *mblock_dlos_grid*.\n";
1159 error_found = true;
1160 }
1161 }
1162 } else {
1163 // Demands differs between the options and checks are done inside
1164 // sub-functions
1165 }
1166
1167 // If errors where found throw runtime_error with the collected error
1168 // message.
1169 if (error_found) throw runtime_error(os.str());
1170
1171 // And finally check if grids and data size match
1172 antenna_response.checksize_strict();
1173
1174 // Call the core function
1175 //
1176 Sparse hantenna;
1177 //
1178 if (antenna_dim == 1)
1179 antenna1d_matrix(hantenna,
1180 antenna_dim,
1181 antenna_dlos(joker, 0),
1182 antenna_response,
1183 sensor_response_dlos_grid(joker, 0),
1184 sensor_response_f_grid,
1185 npol,
1186 sensor_norm);
1187 else {
1188
1189 if (option_2d == "interp_response" ) {
1191 antenna_dim,
1192 antenna_dlos,
1193 antenna_response,
1194 sensor_response_dlos_grid,
1195 sensor_response_f_grid,
1196 npol);
1197 }
1198 else if (option_2d == "gridded_dlos" ) {
1199 antenna2d_gridded_dlos(hantenna,
1200 antenna_dim,
1201 antenna_dlos,
1202 antenna_response,
1203 sensor_response_dlos_grid,
1204 sensor_response_f_grid,
1205 npol);
1206 }
1207
1208 else {
1209 throw runtime_error( "Unrecognised choice for *option_2d*." );
1210 }
1211 }
1212
1213 // Here we need a temporary sparse that is copy of the sensor_response
1214 // sparse matrix. We need it since the multiplication function can not
1215 // take the same object as both input and output.
1216 Sparse htmp = sensor_response;
1217 sensor_response.resize(hantenna.nrows(), htmp.ncols());
1218 mult(sensor_response, hantenna, htmp);
1219
1220 // Some extra output.
1221 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
1222 << sensor_response.ncols() << "\n";
1223
1224 // Update sensor_response_dlos_grid
1225 sensor_response_dlos_grid = antenna_dlos;
1226
1227 // Set aux variables
1228 sensor_aux_vectors(sensor_response_f,
1229 sensor_response_pol,
1230 sensor_response_dlos,
1231 sensor_response_f_grid,
1232 sensor_response_pol_grid,
1233 sensor_response_dlos_grid);
1234}
1235
1236
1237
1238/* Workspace method: Doxygen documentation will be auto-generated */
1240 Sparse& sensor_response,
1241 Vector& sensor_response_f,
1242 ArrayOfIndex& sensor_response_pol,
1243 Matrix& sensor_response_dlos,
1244 Vector& sensor_response_f_grid,
1245 const ArrayOfIndex& sensor_response_pol_grid,
1246 const Matrix& sensor_response_dlos_grid,
1247 const Vector& f_backend,
1248 const ArrayOfGriddedField1& backend_channel_response,
1249 const Index& sensor_norm,
1250 const Verbosity& verbosity) {
1252
1253 // Some sizes
1254 const Index nf = sensor_response_f_grid.nelem();
1255 const Index npol = sensor_response_pol_grid.nelem();
1256 const Index nlos = sensor_response_dlos_grid.nrows();
1257 const Index nin = nf * npol * nlos;
1258
1259 // Initialise an output stream for runtime errors and a flag for errors
1260 ostringstream os;
1261 bool error_found = false;
1262
1263 // Check that sensor_response variables are consistent in size
1264 if (sensor_response_f.nelem() != nin) {
1265 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
1266 << "grid variables (sensor_response_f_grid etc.).\n";
1267 error_found = true;
1268 }
1269 if (sensor_response.nrows() != nin) {
1270 os << "The sensor block response matrix *sensor_response* does not have\n"
1271 << "right size compared to the sensor grid variables\n"
1272 << "(sensor_response_f_grid etc.).\n";
1273 error_found = true;
1274 }
1275
1276 // We allow f_backend to be unsorted, but must be inside sensor_response_f_grid
1277 if (min(f_backend) < min(sensor_response_f_grid)) {
1278 os << "At least one value in *f_backend* (" << min(f_backend)
1279 << ") below range\ncovered by *sensor_response_f_grid* ("
1280 << min(sensor_response_f_grid) << ").\n";
1281 error_found = true;
1282 }
1283 if (max(f_backend) > max(sensor_response_f_grid)) {
1284 os << "At least one value in *f_backend* (" << max(f_backend)
1285 << ") above range\ncovered by *sensor_response_f_grid* ("
1286 << max(sensor_response_f_grid) << ").\n";
1287 error_found = true;
1288 }
1289
1290 // Check number of columns in backend_channel_response
1291 //
1292 const Index nrp = backend_channel_response.nelem();
1293 //
1294 if (nrp != 1 && nrp != f_backend.nelem()) {
1295 os << "The WSV *backend_channel_response* must have 1 or n elements,\n"
1296 << "where n is the length of *f_backend*.\n";
1297 error_found = true;
1298 }
1299
1300 // If errors where found throw runtime_error with the collected error
1301 // message (before error message gets too long).
1302 if (error_found) throw runtime_error(os.str());
1303
1304 Numeric f_dlow = 0.0;
1305 Numeric f_dhigh = 0.0;
1306
1307 Index freq_full = nrp > 1;
1308 for (Index i = 0; i < f_backend.nelem(); i++) {
1309 const Index irp = i * freq_full;
1310 ConstVectorView bchr_f_grid =
1311 backend_channel_response[irp].get_numeric_grid(GFIELD1_F_GRID);
1312
1313 if (bchr_f_grid.nelem() != backend_channel_response[irp].data.nelem()) {
1314 os << "Mismatch in size of grid and data in element " << i
1315 << "\nof *sideband_response*.\n";
1316 error_found = true;
1317 }
1318
1319 if (!is_increasing(bchr_f_grid)) {
1320 os << "The frequency grid of element " << irp
1321 << " in *backend_channel_response*\nis not strictly increasing.\n";
1322 error_found = true;
1323 }
1324
1325 // Check if the relative grid added to the channel frequencies expands
1326 // outside sensor_response_f_grid.
1327 //
1328 Numeric f1 = f_backend[i] + bchr_f_grid[0] - min(sensor_response_f_grid);
1329 Numeric f2 =
1330 (max(sensor_response_f_grid) - f_backend[i]) - last(bchr_f_grid);
1331 //
1332 f_dlow = min(f_dlow, f1);
1333 f_dhigh = min(f_dhigh, f2);
1334 }
1335
1336 if (f_dlow < 0) {
1337 os << "The WSV *sensor_response_f_grid* is too narrow. It should be\n"
1338 << "expanded with " << -f_dlow << " Hz in the lower end. This change\n"
1339 << "should be applied to either *f_grid* or the sensor part in\n"
1340 << "front of *sensor_responseBackend*.\n";
1341 error_found = true;
1342 }
1343 if (f_dhigh < 0) {
1344 os << "The WSV *sensor_response_f_grid* is too narrow. It should be\n"
1345 << "expanded with " << -f_dhigh << " Hz in the higher end. This change\n"
1346 << "should be applied to either *f_grid* or the sensor part in\n"
1347 << "front of *sensor_responseBackend*.\n";
1348 error_found = true;
1349 }
1350
1351 // If errors where found throw runtime_error with the collected error
1352 // message.
1353 if (error_found) throw runtime_error(os.str());
1354
1355 // Call the core function
1356 //
1357 Sparse hbackend;
1358 //
1359 spectrometer_matrix(hbackend,
1360 f_backend,
1361 backend_channel_response,
1362 sensor_response_f_grid,
1363 npol,
1364 nlos,
1365 sensor_norm);
1366
1367 // Here we need a temporary sparse that is copy of the sensor_response
1368 // sparse matrix. We need it since the multiplication function can not
1369 // take the same object as both input and output.
1370 Sparse htmp = sensor_response;
1371 sensor_response.resize(hbackend.nrows(), htmp.ncols());
1372 mult(sensor_response, hbackend, htmp);
1373
1374 // Some extra output.
1375 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
1376 << sensor_response.ncols() << "\n";
1377
1378 // Update sensor_response_f_grid
1379 sensor_response_f_grid = f_backend;
1380
1381 // Set aux variables
1382 sensor_aux_vectors(sensor_response_f,
1383 sensor_response_pol,
1384 sensor_response_dlos,
1385 sensor_response_f_grid,
1386 sensor_response_pol_grid,
1387 sensor_response_dlos_grid);
1388}
1389
1390
1391
1392/* Workspace method: Doxygen documentation will be auto-generated */
1394 Sparse& sensor_response,
1395 Vector& sensor_response_f,
1396 ArrayOfIndex& sensor_response_pol,
1397 Matrix& sensor_response_dlos,
1398 Vector& sensor_response_f_grid,
1399 const ArrayOfIndex& sensor_response_pol_grid,
1400 const Matrix& sensor_response_dlos_grid,
1401 const Vector& f_backend,
1402 const ArrayOfGriddedField1& backend_channel_response,
1403 const Index& sensor_norm,
1404 const Numeric& df1,
1405 const Numeric& df2,
1406 const Verbosity& verbosity) {
1407 // All needed checks are done in sensor_responseBackend
1408
1409 Sparse H1 = sensor_response, H2 = sensor_response;
1410
1411 // Some needed vectors
1412 Vector f_backend_shifted;
1413 Vector fdummy = sensor_response_f, fdummy_grid = sensor_response_f_grid;
1414
1415 // Cycle 1
1416 f_backend_shifted = f_backend;
1417 f_backend_shifted += df1;
1418 //
1420 fdummy,
1421 sensor_response_pol,
1422 sensor_response_dlos,
1423 fdummy_grid,
1424 sensor_response_pol_grid,
1425 sensor_response_dlos_grid,
1426 f_backend_shifted,
1427 backend_channel_response,
1428 sensor_norm,
1429 verbosity);
1430 // Cycle 2
1431 f_backend_shifted = f_backend;
1432 f_backend_shifted += df2;
1433 //
1435 sensor_response_f,
1436 sensor_response_pol,
1437 sensor_response_dlos,
1438 sensor_response_f_grid,
1439 sensor_response_pol_grid,
1440 sensor_response_dlos_grid,
1441 f_backend_shifted,
1442 backend_channel_response,
1443 sensor_norm,
1444 verbosity);
1445
1446 // Total response
1447 sub(sensor_response, H2, H1);
1448
1449 // sensor_response_f_grid shall be f_backend
1450 sensor_response_f_grid = f_backend;
1451
1452 // Set aux variables
1453 sensor_aux_vectors(sensor_response_f,
1454 sensor_response_pol,
1455 sensor_response_dlos,
1456 sensor_response_f_grid,
1457 sensor_response_pol_grid,
1458 sensor_response_dlos_grid);
1459}
1460
1461
1462
1463/* Workspace method: Doxygen documentation will be auto-generated */
1465 Vector& sensor_response_f,
1466 ArrayOfIndex& sensor_response_pol,
1467 Matrix& sensor_response_dlos,
1468 Matrix& sensor_response_dlos_grid,
1469 const Vector& sensor_response_f_grid,
1470 const ArrayOfIndex& sensor_response_pol_grid,
1471 const Numeric& w1,
1472 const Numeric& w2,
1473 const Verbosity& verbosity) {
1475
1476 if (sensor_response_dlos_grid.nrows() != 2)
1477 throw runtime_error(
1478 "This method requires that the number of observation directions is 2.");
1479
1480 if (sensor_response_pol_grid.nelem() != 1)
1481 throw runtime_error(
1482 "This method handles (so far) only single polarisation cases.");
1483
1484 const Index n = sensor_response_f_grid.nelem();
1485
1486 // Form H matrix representing beam switching
1487 Sparse Hbswitch(n, 2 * n);
1488 Vector hrow(2 * n, 0.0);
1489 //
1490 for (Index i = 0; i < n; i++) {
1491 hrow[i] = w1;
1492 hrow[i + n] = w2;
1493 //
1494 Hbswitch.insert_row(i, hrow);
1495 //
1496 hrow = 0;
1497 }
1498
1499 // Here we need a temporary sparse that is copy of the sensor_response
1500 // sparse matrix. We need it since the multiplication function can not
1501 // take the same object as both input and output.
1502 Sparse Htmp = sensor_response;
1503 sensor_response.resize(Hbswitch.nrows(), Htmp.ncols());
1504 mult(sensor_response, Hbswitch, Htmp);
1505
1506 // Some extra output.
1507 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
1508 << sensor_response.ncols() << "\n";
1509
1510 // Update sensor_response_za_grid
1511 const Vector zaaa = sensor_response_dlos_grid(1, joker);
1512 sensor_response_dlos_grid.resize(1, zaaa.nelem());
1513 sensor_response_dlos_grid(0, joker) = zaaa;
1514
1515 // Set aux variables
1516 sensor_aux_vectors(sensor_response_f,
1517 sensor_response_pol,
1518 sensor_response_dlos,
1519 sensor_response_f_grid,
1520 sensor_response_pol_grid,
1521 sensor_response_dlos_grid);
1522}
1523
1524
1525
1526/* Workspace method: Doxygen documentation will be auto-generated */
1528 Sparse& sensor_response,
1529 Vector& sensor_response_f,
1530 ArrayOfIndex& sensor_response_pol,
1531 Matrix& sensor_response_dlos,
1532 Vector& sensor_response_f_grid,
1533 const ArrayOfIndex& sensor_response_pol_grid,
1534 const Matrix& sensor_response_dlos_grid,
1535 const Verbosity& verbosity) {
1537
1538 if (sensor_response_dlos_grid.nrows() != 1)
1539 throw runtime_error(
1540 "This method requires that the number of observation directions is 1.");
1541
1542 if (sensor_response_pol_grid.nelem() != 1)
1543 throw runtime_error(
1544 "This method handles (so far) only single polarisation cases.");
1545
1546 const Index n = sensor_response_f_grid.nelem();
1547 const Index n2 = n / 2;
1548
1549 if (sensor_response.nrows() != n)
1550 throw runtime_error(
1551 "Assumptions of method are not fulfilled, "
1552 "considering number of rows in *sensor_response* "
1553 "and length of *sensor_response_f_grid*.");
1554
1555 if (!is_multiple(n, 2))
1556 throw runtime_error(
1557 "There is an odd number of total frequencies, "
1558 "which is not consistent with the assumptions of "
1559 "the method.");
1560
1561 // Form H matrix representing frequency switching
1562 Sparse Hbswitch(n2, n);
1563 Vector hrow(n, 0.0);
1564 //
1565 for (Index i = 0; i < n2; i++) {
1566 hrow[i] = -1;
1567 hrow[i + n2] = 1;
1568 //
1569 Hbswitch.insert_row(i, hrow);
1570 //
1571 hrow = 0;
1572 }
1573
1574 // Here we need a temporary sparse that is copy of the sensor_response
1575 // sparse matrix. We need it since the multiplication function can not
1576 // take the same object as both input and output.
1577 Sparse Htmp = sensor_response;
1578 sensor_response.resize(Hbswitch.nrows(), Htmp.ncols());
1579 mult(sensor_response, Hbswitch, Htmp);
1580
1581 // Some extra output.
1582 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
1583 << sensor_response.ncols() << "\n";
1584
1585 // Update sensor_response_f_grid
1586 const Vector f = sensor_response_f_grid;
1587 sensor_response_f_grid.resize(n2);
1588 sensor_response_f_grid = f[Range(n2, n2)];
1589
1590 // Set aux variables
1591 sensor_aux_vectors(sensor_response_f,
1592 sensor_response_pol,
1593 sensor_response_dlos,
1594 sensor_response_f_grid,
1595 sensor_response_pol_grid,
1596 sensor_response_dlos_grid);
1597}
1598
1599
1600
1601/* Workspace method: Doxygen documentation will be auto-generated */
1602void sensor_responseIF2RF( // WS Output:
1603 Vector& sensor_response_f,
1604 Vector& sensor_response_f_grid,
1605 // WS Input:
1606 const Numeric& lo,
1607 const String& sideband_mode,
1608 const Verbosity&) {
1609 // Check that frequencies are not too high. This might be a floating limit.
1610 // For this we use the variable f_lim, given in Hz.
1611 Numeric f_lim = 30e9;
1612 if (max(sensor_response_f_grid) > f_lim)
1613 throw runtime_error("The frequencies seem to already be given in RF.");
1614
1615 // Lower band
1616 if (sideband_mode == "lower") {
1617 sensor_response_f *= -1;
1618 sensor_response_f_grid *= -1;
1619 sensor_response_f += lo;
1620 sensor_response_f_grid += lo;
1621 }
1622
1623 // Upper band
1624 else if (sideband_mode == "upper") {
1625 sensor_response_f += lo;
1626 sensor_response_f_grid += lo;
1627 }
1628
1629 // Unknown option
1630 else {
1631 throw runtime_error(
1632 "Only allowed options for *sideband _mode* are \"lower\" and \"upper\".");
1633 }
1634}
1635
1636
1637
1638/* Workspace method: Doxygen documentation will be auto-generated */
1639void sensor_responseFillFgrid(Sparse& sensor_response,
1640 Vector& sensor_response_f,
1641 ArrayOfIndex& sensor_response_pol,
1642 Matrix& sensor_response_dlos,
1643 Vector& sensor_response_f_grid,
1644 const ArrayOfIndex& sensor_response_pol_grid,
1645 const Matrix& sensor_response_dlos_grid,
1646 const Index& polyorder,
1647 const Index& nfill,
1648 const Verbosity& verbosity) {
1650
1651 // Some sizes
1652 const Index nf = sensor_response_f_grid.nelem();
1653 const Index npol = sensor_response_pol_grid.nelem();
1654 const Index nlos = sensor_response_dlos_grid.nrows();
1655 const Index nin = nf * npol * nlos;
1656
1657 // Initialise a output stream for runtime errors and a flag for errors
1658 ostringstream os;
1659 bool error_found = false;
1660
1661 // Check that sensor_response variables are consistent in size
1662 if (sensor_response_f.nelem() != nin) {
1663 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
1664 << "grid variables (sensor_response_f_grid etc.).\n";
1665 error_found = true;
1666 }
1667 if (sensor_response.nrows() != nin) {
1668 os << "The sensor block response matrix *sensor_response* does not have\n"
1669 << "right size compared to the sensor grid variables\n"
1670 << "(sensor_response_f_grid etc.).\n";
1671 error_found = true;
1672 }
1673
1674 // Check polyorder and nfill
1675 if (polyorder < 2 || polyorder > 7) {
1676 os << "Accepted range for *polyorder* is [3,7].\n";
1677 error_found = true;
1678 }
1679 if (nfill < 1) {
1680 os << "The argument *nfill* must be > 1.\n";
1681 error_found = true;
1682 }
1683
1684 // If errors where found throw runtime_error with the collected error
1685 // message.
1686 if (error_found) throw runtime_error(os.str());
1687
1688 // New frequency grid
1689 //
1690 const Index n1 = nfill + 1;
1691 const Index n2 = nfill + 2;
1692 const Index nnew = (nf - 1) * n1 + 1;
1693 //
1694 Vector fnew(nnew);
1695 //
1696 for (Index i = 0; i < nf - 1; i++) {
1697 Vector fp(n2);
1698 nlinspace(fp, sensor_response_f_grid[i], sensor_response_f_grid[i + 1], n2);
1699 fnew[Range(i * n1, n2)] = fp;
1700 }
1701
1702 const auto lag = Interpolation::LagrangeVector(fnew, sensor_response_f_grid, polyorder);
1703
1704 // Set up H for this part
1705 //
1706 Sparse hpoly(nnew * npol * nlos, nin);
1707 Vector hrow(nin, 0.0);
1708 Index row = 0;
1709 //
1710 for (Index ilos = 0; ilos < nlos; ilos++) {
1711 for (Index iv = 0; iv < nnew; iv++) {
1712 for (Index ip = 0; ip < npol; ip++) {
1713 const Index col0 = ilos * nf * npol;
1714 for (Index i = 0; i < polyorder+1; i++) {
1715 const Numeric w = lag[iv].lx[i];
1716 if (abs(w) > 1e-5) {
1717 hrow[col0 + (lag[iv].pos + i) * npol + ip] = lag[iv].lx[i];
1718 }
1719 }
1720 hpoly.insert_row(row, hrow);
1721 for (Index i = 0; i < polyorder+1; i++) {
1722 hrow[col0 + (lag[iv].pos + i) * npol + ip] = 0;
1723 }
1724 row += 1;
1725 }
1726 }
1727 }
1728
1729 // Here we need a temporary sparse that is copy of the sensor_response
1730 // sparse matrix. We need it since the multiplication function can not
1731 // take the same object as both input and output.
1732 Sparse htmp = sensor_response;
1733 sensor_response.resize(hpoly.nrows(), htmp.ncols());
1734 mult(sensor_response, hpoly, htmp);
1735
1736 // Some extra output.
1737 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
1738 << sensor_response.ncols() << "\n";
1739
1740 // Update sensor_response_f_grid
1741 sensor_response_f_grid = fnew;
1742
1743 // Set aux variables
1744 sensor_aux_vectors(sensor_response_f,
1745 sensor_response_pol,
1746 sensor_response_dlos,
1747 sensor_response_f_grid,
1748 sensor_response_pol_grid,
1749 sensor_response_dlos_grid);
1750}
1751
1752
1753
1754/* Workspace method: Doxygen documentation will be auto-generated */
1755void sensor_responseInit(Sparse& sensor_response,
1756 Vector& sensor_response_f,
1757 ArrayOfIndex& sensor_response_pol,
1758 Matrix& sensor_response_dlos,
1759 Vector& sensor_response_f_grid,
1760 ArrayOfIndex& sensor_response_pol_grid,
1761 Matrix& sensor_response_dlos_grid,
1762 const Vector& f_grid,
1763 const Matrix& mblock_dlos_grid,
1764 const Index& antenna_dim,
1765 const Index& atmosphere_dim,
1766 const Index& stokes_dim,
1767 const Index& sensor_norm,
1768 const Verbosity& verbosity) {
1771
1772 // Check input
1773
1774 // Basic variables
1775 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1776 chk_if_in_range("antenna_dim", antenna_dim, 1, 2);
1777 chk_if_bool("sensor_norm", sensor_norm);
1778
1779 // mblock_dlos_grid
1780 if (mblock_dlos_grid.empty())
1781 throw runtime_error("*mblock_dlos_grid* is empty.");
1782 if (mblock_dlos_grid.ncols() > 2)
1783 throw runtime_error(
1784 "The maximum number of columns in *mblock_dlos_grid* is two.");
1785 if (atmosphere_dim < 3) {
1786 if (mblock_dlos_grid.ncols() != 1)
1787 throw runtime_error(
1788 "For 1D and 2D *mblock_dlos_grid* must have exactly one column.");
1789 if (antenna_dim == 2)
1790 throw runtime_error(
1791 "2D antennas (antenna_dim=2) can only be "
1792 "used with 3D atmospheres.");
1793 }
1794
1795 // Set grid variables
1796 sensor_response_f_grid = f_grid;
1797 sensor_response_dlos_grid = mblock_dlos_grid;
1798 //
1799 sensor_response_pol_grid.resize(stokes_dim);
1800 //
1801 for (Index is = 0; is < stokes_dim; is++) {
1802 sensor_response_pol_grid[is] = is + 1;
1803 }
1804
1805 // Set aux variables
1806 sensor_aux_vectors(sensor_response_f,
1807 sensor_response_pol,
1808 sensor_response_dlos,
1809 sensor_response_f_grid,
1810 sensor_response_pol_grid,
1811 sensor_response_dlos_grid);
1812
1813 //Set response matrix to identity matrix
1814 //
1815 const Index n = sensor_response_f.nelem();
1816 //
1817 out2 << " Initialising *sensor_reponse* as a identity matrix.\n";
1818 out3 << " Size of *sensor_response*: " << n << "x" << n << "\n";
1819 //
1820 sensor_response.resize(n, n);
1821 id_mat(sensor_response);
1822}
1823
1824
1825
1826/* Workspace method: Doxygen documentation will be auto-generated */
1827void sensorOff(Sparse& sensor_response,
1828 Vector& sensor_response_f,
1829 ArrayOfIndex& sensor_response_pol,
1830 Matrix& sensor_response_dlos,
1831 Vector& sensor_response_f_grid,
1832 ArrayOfIndex& sensor_response_pol_grid,
1833 Matrix& sensor_response_dlos_grid,
1834 Matrix& mblock_dlos_grid,
1835 const Index& stokes_dim,
1836 const Vector& f_grid,
1837 const Verbosity& verbosity) {
1838 // Checks are done in sensor_responseInit.
1839 Index antenna_dim;
1840 AntennaOff(antenna_dim, mblock_dlos_grid, verbosity);
1841
1842 // Dummy variables (The method is independent of atmosphere_dim.
1843 // atmosphere_dim used below just for some checks when antenna is active, not
1844 // relevant here. ).
1845 const Index sensor_norm = 1, atmosphere_dim = 1;
1846
1847 sensor_responseInit(sensor_response,
1848 sensor_response_f,
1849 sensor_response_pol,
1850 sensor_response_dlos,
1851 sensor_response_f_grid,
1852 sensor_response_pol_grid,
1853 sensor_response_dlos_grid,
1854 f_grid,
1855 mblock_dlos_grid,
1856 antenna_dim,
1857 atmosphere_dim,
1858 stokes_dim,
1859 sensor_norm,
1860 verbosity);
1861}
1862
1863
1864
1865/* Workspace method: Doxygen documentation will be auto-generated */
1866void sensor_responseMixer(Sparse& sensor_response,
1867 Vector& sensor_response_f,
1868 ArrayOfIndex& sensor_response_pol,
1869 Matrix& sensor_response_dlos,
1870 Vector& sensor_response_f_grid,
1871 const ArrayOfIndex& sensor_response_pol_grid,
1872 const Matrix& sensor_response_dlos_grid,
1873 const Numeric& lo,
1874 const GriddedField1& sideband_response,
1875 const Index& sensor_norm,
1876 const Verbosity& verbosity) {
1878
1879 // Some sizes
1880 const Index nf = sensor_response_f_grid.nelem();
1881 const Index npol = sensor_response_pol_grid.nelem();
1882 const Index nlos = sensor_response_dlos_grid.nrows();
1883 const Index nin = nf * npol * nlos;
1884
1885 // Frequency grid of for sideband response specification
1886 ConstVectorView sbresponse_f_grid =
1887 sideband_response.get_numeric_grid(GFIELD1_F_GRID);
1888
1889 // Initialise a output stream for runtime errors and a flag for errors
1890 ostringstream os;
1891 bool error_found = false;
1892
1893 // Check that sensor_response variables are consistent in size
1894 if (sensor_response_f.nelem() != nin) {
1895 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
1896 << "grid variables (sensor_response_f_grid etc.).\n";
1897 error_found = true;
1898 }
1899 if (sensor_response.nrows() != nin) {
1900 os << "The sensor block response matrix *sensor_response* does not have\n"
1901 << "right size compared to the sensor grid variables\n"
1902 << "(sensor_response_f_grid etc.).\n";
1903 error_found = true;
1904 }
1905
1906 // Check that the lo frequency is within the sensor_response_f_grid
1907 if (lo <= sensor_response_f_grid[0] || lo >= last(sensor_response_f_grid)) {
1908 os << "The given local oscillator frequency is outside the sensor\n"
1909 << "frequency grid. It must be within the *sensor_response_f_grid*.\n";
1910 error_found = true;
1911 }
1912
1913 // Checks of sideband_response, partly in combination with lo
1914 if (sbresponse_f_grid.nelem() != sideband_response.data.nelem()) {
1915 os << "Mismatch in size of grid and data in *sideband_response*.\n";
1916 error_found = true;
1917 }
1918 if (sbresponse_f_grid.nelem() < 2) {
1919 os << "At least two data points must be specified in "
1920 << "*sideband_response*.\n";
1921 error_found = true;
1922 }
1923 if (!is_increasing(sbresponse_f_grid)) {
1924 os << "The frequency grid of *sideband_response* must be strictly\n"
1925 << "increasing.\n";
1926 error_found = true;
1927 }
1928 if (fabs(last(sbresponse_f_grid) + sbresponse_f_grid[0]) > 0) {
1929 os << "The end points of the *sideband_response* frequency grid must be\n"
1930 << "symmetrically placed around 0. That is, the grid shall cover a\n"
1931 << "a range that can be written as [-df,df]. \n";
1932 error_found = true;
1933 }
1934
1935 // Check that response function does not extend outside sensor_response_f_grid
1936 Numeric df_high = lo + last(sbresponse_f_grid) - last(sensor_response_f_grid);
1937 Numeric df_low = sensor_response_f_grid[0] - lo - sbresponse_f_grid[0];
1938 if (df_high > 0 && df_low > 0) {
1939 os << "The *sensor_response_f* grid must be extended by at least\n"
1940 << df_low << " Hz in the lower end and " << df_high << " Hz in the\n"
1941 << "upper end to cover frequency range set by *sideband_response*\n"
1942 << "and *lo*. Or can the frequency grid of *sideband_response* be\n"
1943 << "decreased?";
1944 error_found = true;
1945 } else if (df_high > 0) {
1946 os << "The *sensor_response_f* grid must be extended by at " << df_high
1947 << " Hz\nin the upper end to cover frequency range set by\n"
1948 << "*sideband_response* and *lo*. Or can the frequency grid of\n"
1949 << "*sideband_response* be decreased?";
1950 error_found = true;
1951 } else if (df_low > 0) {
1952 os << "The *sensor_response_f* grid must be extended by at " << df_low
1953 << " Hz\nin the lower end to cover frequency range set by\n"
1954 << "*sideband_response* and *lo*. Or can the frequency grid of\n"
1955 << "*sideband_response* be decreased?";
1956 error_found = true;
1957 }
1958
1959 // If errors where found throw runtime_error with the collected error
1960 // message.
1961 if (error_found) throw runtime_error(os.str());
1962
1963 //Call the core function
1964 //
1965 Sparse hmixer;
1966 Vector f_mixer;
1967 //
1968 mixer_matrix(hmixer,
1969 f_mixer,
1970 lo,
1971 sideband_response,
1972 sensor_response_f_grid,
1973 npol,
1974 nlos,
1975 sensor_norm);
1976
1977 // Here we need a temporary sparse that is copy of the sensor_response
1978 // sparse matrix. We need it since the multiplication function can not
1979 // take the same object as both input and output.
1980 Sparse htmp = sensor_response;
1981 sensor_response.resize(hmixer.nrows(), htmp.ncols());
1982 mult(sensor_response, hmixer, htmp);
1983
1984 // Some extra output.
1985 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
1986 << sensor_response.ncols() << "\n";
1987
1988 // Update sensor_response_f_grid
1989 sensor_response_f_grid = f_mixer;
1990
1991 // Set aux variables
1992 sensor_aux_vectors(sensor_response_f,
1993 sensor_response_pol,
1994 sensor_response_dlos,
1995 sensor_response_f_grid,
1996 sensor_response_pol_grid,
1997 sensor_response_dlos_grid);
1998}
1999
2000
2001
2002/* Workspace method: Doxygen documentation will be auto-generated */
2004 // WS Output:
2005 Index& antenna_dim,
2006 Matrix& mblock_dlos_grid,
2007 Sparse& sensor_response,
2008 Vector& sensor_response_f,
2009 ArrayOfIndex& sensor_response_pol,
2010 Matrix& sensor_response_dlos,
2011 Vector& sensor_response_f_grid,
2012 ArrayOfIndex& sensor_response_pol_grid,
2013 Matrix& sensor_response_dlos_grid,
2014 Index& sensor_norm,
2015 // WS Input:
2016 const Index& atmosphere_dim,
2017 const Index& stokes_dim,
2018 const Vector& f_grid,
2019 const Vector& f_backend,
2020 const ArrayOfArrayOfIndex& channel2fgrid_indexes,
2021 const ArrayOfVector& channel2fgrid_weights,
2022 const String& iy_unit,
2023 const Matrix& antenna_dlos,
2024 const ArrayOfString& mm_pol, /* met_mm_polarisation */
2025 const Vector& mm_ant, /* met_mm_antenna */
2026 // Control Parameters:
2027 const Index& use_antenna,
2028 const Index& mirror_dza,
2029 const Verbosity& verbosity) {
2030 // Input checks
2031
2032 chk_if_bool("use_antenna", use_antenna);
2033 chk_if_bool("mirror_dza", mirror_dza);
2034
2035 if (mm_ant.nelem())
2036 throw std::runtime_error(
2037 "So far inclusion of antenna pattern is NOT supported and\n"
2038 "*met_mm_antenna* must be empty.\n");
2039
2040 if (!(atmosphere_dim == 1 || atmosphere_dim == 3))
2041 throw std::runtime_error(
2042 "This method only supports 1D and 3D atmospheres.");
2043
2044 if (antenna_dlos.empty())
2045 throw std::runtime_error("*antenna_dlos* is empty.");
2046
2047 if (antenna_dlos.ncols() > 2)
2048 throw std::runtime_error(
2049 "The maximum number of columns in *antenna_dlos*"
2050 "is two.");
2051
2052 // Copy, and possibly extend antenna_dlos
2053 Matrix antenna_dlos_local;
2054
2055 // Mirror?
2056 if (mirror_dza) {
2057 if (antenna_dlos.ncols() > 1)
2058 throw std::runtime_error(
2059 "With *mirror_dza* set to true, *antenna_dlos*"
2060 "is only allowed to have a single column.");
2061 if (atmosphere_dim != 3)
2062 throw std::runtime_error("*mirror_dza* only makes sense in 3D.");
2063 // We don't want to duplicate zero!
2064 const Index n = antenna_dlos.nrows();
2065 Index nnew = 0;
2066 for (Index i = 0; i < n; i++) {
2067 if (antenna_dlos(i, 0) != 0) {
2068 nnew += 1;
2069 }
2070 }
2071 antenna_dlos_local.resize(n + nnew, 1);
2072 antenna_dlos_local(Range(0, n), 0) = antenna_dlos(joker, 0);
2073 Index pos = n;
2074 for (Index i = n - 1; i >= 0; i--) {
2075 if (antenna_dlos(i, 0) != 0) {
2076 antenna_dlos_local(pos, 0) = -antenna_dlos(i, 0);
2077 pos += 1;
2078 }
2079 }
2080 } else {
2081 antenna_dlos_local = antenna_dlos;
2082 }
2083
2084 // No normalisation needed here, and set antenna_dim=1 as temporary solution
2085 sensor_norm = 0;
2086 antenna_dim = 1;
2087
2088 // Create sensor response for mixer+backend, matching one viewing angle
2089 Sparse sensor_response_single;
2090 Matrix mblock_dlos_dummy(1, 1);
2091 mblock_dlos_dummy(0, 0) = 0;
2092 sensor_responseInit(sensor_response_single,
2093 sensor_response_f,
2094 sensor_response_pol,
2095 sensor_response_dlos,
2096 sensor_response_f_grid,
2097 sensor_response_pol_grid,
2098 sensor_response_dlos_grid,
2099 f_grid,
2100 mblock_dlos_dummy,
2101 antenna_dim,
2102 atmosphere_dim,
2103 stokes_dim,
2104 sensor_norm,
2105 verbosity);
2106 sensor_responseMixerBackendPrecalcWeights(sensor_response_single,
2107 sensor_response_f,
2108 sensor_response_pol,
2109 sensor_response_dlos,
2110 sensor_response_f_grid,
2111 sensor_response_pol_grid,
2112 sensor_response_dlos_grid,
2113 f_backend,
2114 channel2fgrid_indexes,
2115 channel2fgrid_weights,
2116 verbosity);
2117
2118 // Construct complete sensor_response matrix
2119 const Index num_f = f_grid.nelem();
2120 const Index nchannels = f_backend.nelem();
2121 sensor_response = Sparse(nchannels * antenna_dlos_local.nrows(),
2122 num_f * stokes_dim * antenna_dlos_local.nrows());
2123
2124 sensor_response_pol_grid.resize(1);
2125 sensor_response_pol_grid[0] = 1;
2126
2127 if (stokes_dim > 1) {
2128 // With polarisation
2129 if (mm_pol.nelem() != nchannels) {
2130 ostringstream os;
2131 os << "Length of *met_mm_polarisation* (" << mm_pol.nelem()
2132 << ") must match\n"
2133 << "number of channels in *met_mm_backend* (" << nchannels << ").";
2134 throw std::runtime_error(os.str());
2135 }
2136
2137 Sparse H_pol;
2138 Sparse sensor_response_tmp;
2139
2140 for (Index iza = 0; iza < antenna_dlos_local.nrows(); iza++) {
2141 sensor_response_tmp = Sparse(nchannels, sensor_response_single.ncols());
2143 H_pol, mm_pol, antenna_dlos_local(iza, 0), stokes_dim, iy_unit);
2144 mult(sensor_response_tmp, H_pol, sensor_response_single);
2145 for (Index r = 0; r < sensor_response_tmp.nrows(); r++)
2146 for (Index c = 0; c < sensor_response_tmp.ncols(); c++) {
2147 const Numeric v = sensor_response_tmp(r, c);
2148
2149 if (v != 0.)
2150 sensor_response.rw(iza * nchannels + r,
2151 iza * num_f * stokes_dim + c) = v;
2152 }
2153 }
2154 } else {
2155 // No polarisation
2156 for (Index iza = 0; iza < antenna_dlos_local.nrows(); iza++) {
2157 for (Index r = 0; r < sensor_response_single.nrows(); r++)
2158 for (Index c = 0; c < sensor_response_single.ncols(); c++) {
2159 const Numeric v = sensor_response_single(r, c);
2160
2161 if (v != 0.)
2162 sensor_response.rw(iza * nchannels + r,
2163 iza * num_f * stokes_dim + c) = v;
2164 }
2165 }
2166 }
2167
2168 antenna_dim = 1;
2169 // Setup antenna
2170 if (use_antenna) {
2171 // FIXME: Do something smart here
2172 throw std::runtime_error("The antenna hasn't arrived yet.");
2173 }
2174
2175 // mblock angle grids
2176 mblock_dlos_grid = antenna_dlos_local;
2177
2178 // Set sensor response aux variables
2179 sensor_response_dlos_grid = mblock_dlos_grid;
2180
2181 // Set aux variables
2182 sensor_aux_vectors(sensor_response_f,
2183 sensor_response_pol,
2184 sensor_response_dlos,
2185 sensor_response_f_grid,
2186 sensor_response_pol_grid,
2187 sensor_response_dlos_grid);
2188}
2189
2190
2191
2192/* Workspace method: Doxygen documentation will be auto-generated */
2194 Sparse& sensor_response,
2195 Vector& sensor_response_f,
2196 ArrayOfIndex& sensor_response_pol,
2197 Matrix& sensor_response_dlos,
2198 Vector& sensor_response_f_grid,
2199 const ArrayOfIndex& sensor_response_pol_grid,
2200 const Matrix& sensor_response_dlos_grid,
2201 const Vector& f_backend,
2202 const ArrayOfArrayOfIndex& channel2fgrid_indexes,
2203 const ArrayOfVector& channel2fgrid_weights,
2204 const Verbosity& verbosity) {
2206
2207 // Some sizes
2208 const Index nin_f = sensor_response_f_grid.nelem();
2209 const Index nout_f = f_backend.nelem();
2210 const Index npol = sensor_response_pol_grid.nelem();
2211 const Index nlos = sensor_response_dlos_grid.nrows();
2212 const Index nin = nin_f * npol * nlos;
2213 const Index nout = nout_f * npol * nlos;
2214
2215 // Initialise an output stream for runtime errors and a flag for errors
2216 ostringstream os;
2217 bool error_found = false;
2218
2219 // Check that sensor_response variables are consistent in size
2220 if (sensor_response_f.nelem() != nin) {
2221 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
2222 << "grid variables (sensor_response_f_grid etc.).\n";
2223 error_found = true;
2224 }
2225 if (sensor_response.nrows() != nin) {
2226 os << "The sensor block response matrix *sensor_response* does not have\n"
2227 << "right size compared to the sensor grid variables\n"
2228 << "(sensor_response_f_grid etc.).\n";
2229 error_found = true;
2230 }
2231
2232 // We allow f_backend to be unsorted, but must be inside sensor_response_f_grid
2233 if (nout_f == 0) {
2234 os << "*f_backend* is empty !!!\n";
2235 error_found = true;
2236 }
2237 if (min(f_backend) < min(sensor_response_f_grid)) {
2238 os << "At least one value in *f_backend* (" << min(f_backend)
2239 << ") below range\ncovered by *sensor_response_f_grid* ("
2240 << min(sensor_response_f_grid) << ").\n";
2241 error_found = true;
2242 }
2243 if (max(f_backend) > max(sensor_response_f_grid)) {
2244 os << "At least one value in *f_backend* (" << max(f_backend)
2245 << ") above range\ncovered by *sensor_response_f_grid* ("
2246 << max(sensor_response_f_grid) << ").\n";
2247 error_found = true;
2248 }
2249
2250 // frequency index and weights
2251 if (channel2fgrid_indexes.nelem() != nout_f) {
2252 os << "The first size of *channel2fgrid_indexes* an length of *f_backend* "
2253 << "must be equal.\n";
2254 error_found = true;
2255 }
2256 if (channel2fgrid_weights.nelem() != channel2fgrid_indexes.nelem()) {
2257 os << "Leading sizes of *channel2fgrid_indexes* and "
2258 << "*channel2fgrid_weights* differ.\n";
2259 error_found = true;
2260 }
2261 for (Index i = 0; i < nout_f; i++) {
2262 if (channel2fgrid_indexes[i].nelem() != channel2fgrid_weights[i].nelem()) {
2263 os << "Mismatch in size between *channel2fgrid_indexes* and "
2264 << "*channel2fgrid_weights*, found for array/vector with "
2265 << "index " << i << ".\n";
2266 error_found = true;
2267 }
2268 for (Index j = 0; j < channel2fgrid_indexes[i].nelem(); j++) {
2269 if (channel2fgrid_indexes[i][j] < 0 ||
2270 channel2fgrid_indexes[i][j] >= nin_f) {
2271 os << "At least one value in *channel2fgrid_indexes* is either "
2272 << " < 0 or is too high considering length of "
2273 << "*sensor_response_f_grid*.\n";
2274 error_found = true;
2275 break;
2276 }
2277 }
2278 }
2279
2280 // If errors where found throw runtime_error with the collected error
2281 if (error_found) throw runtime_error(os.str());
2282
2283 // Create response matrix
2284 //
2285 Sparse hmb(nout, nin);
2286 {
2287 // Loop output channels
2288 for (Index ifr = 0; ifr < nout_f; ifr++) {
2289 // The summation vector for 1 polarisation and 1 viewing direction
2290 Vector w1(nin_f, 0.0);
2291 for (Index j = 0; j < channel2fgrid_indexes[ifr].nelem(); j++) {
2292 w1[channel2fgrid_indexes[ifr][j]] = channel2fgrid_weights[ifr][j];
2293 }
2294
2295 // Loop over polarisation and spectra (viewing directions)
2296 // Weights change only with frequency
2297 // (this code is copied from function spectrometer_matrix)
2298 for (Index sp = 0; sp < nlos; sp++) {
2299 for (Index pol = 0; pol < npol; pol++) {
2300 // Distribute the compact weight vector into a complte one
2301 Vector weights_long(nin, 0.0);
2302 weights_long[Range(sp * nin_f * npol + pol, nin_f, npol)] = w1;
2303
2304 // Insert temp_long into H at the correct row
2305 hmb.insert_row(sp * nout_f * npol + ifr * npol + pol, weights_long);
2306 }
2307 }
2308 }
2309 }
2310
2311 // Here we need a temporary sparse that is copy of the sensor_response
2312 // sparse matrix. We need it since the multiplication function can not
2313 // take the same object as both input and output.
2314 Sparse htmp = sensor_response;
2315 sensor_response.resize(hmb.nrows(), htmp.ncols());
2316 mult(sensor_response, hmb, htmp);
2317
2318 // Update sensor_response_f_grid
2319 sensor_response_f_grid = f_backend;
2320
2321 // Set aux variables
2322 sensor_aux_vectors(sensor_response_f,
2323 sensor_response_pol,
2324 sensor_response_dlos,
2325 sensor_response_f_grid,
2326 sensor_response_pol_grid,
2327 sensor_response_dlos_grid);
2328}
2329
2330
2331
2332/* Workspace method: Doxygen documentation will be auto-generated */
2334 Sparse& sensor_response,
2335 Vector& sensor_response_f,
2336 ArrayOfIndex& sensor_response_pol,
2337 Matrix& sensor_response_dlos,
2338 Vector& sensor_response_f_grid,
2339 const ArrayOfIndex& sensor_response_pol_grid,
2340 const Matrix& sensor_response_dlos_grid,
2341 const Vector& lo_multi,
2342 const ArrayOfGriddedField1& sideband_response_multi,
2343 const ArrayOfString& sideband_mode_multi,
2344 const ArrayOfVector& f_backend_multi,
2345 const ArrayOfArrayOfGriddedField1& backend_channel_response_multi,
2346 const Index& sensor_norm,
2347 const Verbosity& verbosity) {
2348 // Some sizes
2349 const Index nf = sensor_response_f_grid.nelem();
2350 const Index npol = sensor_response_pol_grid.nelem();
2351 const Index nlos = sensor_response_dlos_grid.nrows();
2352 const Index nin = nf * npol * nlos;
2353 const Index nlo = lo_multi.nelem();
2354
2355 // Initialise a output stream for runtime errors and a flag for errors
2356 ostringstream os;
2357 bool error_found = false;
2358
2359 // Check that sensor_response variables are consistent in size
2360 if (sensor_response_f.nelem() != nin) {
2361 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
2362 << "grid variables (sensor_response_f_grid etc.).\n";
2363 error_found = true;
2364 }
2365 if (sensor_response.nrows() != nin) {
2366 os << "The sensor block response matrix *sensor_response* does not have\n"
2367 << "right size compared to the sensor grid variables\n"
2368 << "(sensor_response_f_grid etc.).\n";
2369 error_found = true;
2370 }
2371
2372 // Check that response data are consistent with respect to number of
2373 // mixer/reciever chains.
2374 if (sideband_response_multi.nelem() != nlo) {
2375 os << "Inconsistency in length between *lo_mixer* and "
2376 << "*sideband_response_multi*.\n";
2377 error_found = true;
2378 }
2379 if (sideband_mode_multi.nelem() != nlo) {
2380 os << "Inconsistency in length between *lo_mixer* and "
2381 << "*sideband_mode_multi*.\n";
2382 error_found = true;
2383 }
2384 if (f_backend_multi.nelem() != nlo) {
2385 os << "Inconsistency in length between *lo_mixer* and "
2386 << "*f_backend_multi*.\n";
2387 error_found = true;
2388 }
2389 if (backend_channel_response_multi.nelem() != nlo) {
2390 os << "Inconsistency in length between *lo_mixer* and "
2391 << "*backend_channel_response_multi*.\n";
2392 error_found = true;
2393 }
2394
2395 // If errors where found throw runtime_error with the collected error
2396 // message. Data for each mixer and reciever chain are checked below.
2397 if (error_found) throw runtime_error(os.str());
2398
2399 // Variables for data to be appended
2400 Array<Sparse> sr;
2401 ArrayOfVector srfgrid;
2402 ArrayOfIndex cumsumf(nlo + 1, 0);
2403
2404 for (Index ilo = 0; ilo < nlo; ilo++) {
2405 // Copies of variables that will be changed, but must be
2406 // restored for next loop
2407 Sparse sr1 = sensor_response;
2408 Vector srf1 = sensor_response_f;
2409 ArrayOfIndex srpol1 = sensor_response_pol;
2410 Matrix srdlos1 = sensor_response_dlos;
2411 Vector srfgrid1 = sensor_response_f_grid;
2412
2413 // Call single reciever methods. Try/catch for improved error message.
2414 try {
2416 srf1,
2417 srpol1,
2418 srdlos1,
2419 srfgrid1,
2420 sensor_response_pol_grid,
2421 sensor_response_dlos_grid,
2422 lo_multi[ilo],
2423 sideband_response_multi[ilo],
2424 sensor_norm,
2425 verbosity);
2426
2428 srf1, srfgrid1, lo_multi[ilo], sideband_mode_multi[ilo], verbosity);
2429
2431 srf1,
2432 srpol1,
2433 srdlos1,
2434 srfgrid1,
2435 sensor_response_pol_grid,
2436 sensor_response_dlos_grid,
2437 f_backend_multi[ilo],
2438 backend_channel_response_multi[ilo],
2439 sensor_norm,
2440 verbosity);
2441 } catch (const runtime_error& e) {
2442 ostringstream os2;
2443 os2 << "Error when dealing with receiver/mixer chain (1-based index) "
2444 << ilo + 1 << ":\n"
2445 << e.what();
2446 throw runtime_error(os2.str());
2447 }
2448
2449 // Store in temporary arrays
2450 sr.push_back(sr1);
2451 srfgrid.push_back(srfgrid1);
2452 //
2453 cumsumf[ilo + 1] = cumsumf[ilo] + srfgrid1.nelem();
2454 }
2455
2456 // Append data to create sensor_response_f_grid
2457 //
2458 const Index nfnew = cumsumf[nlo];
2459 sensor_response_f_grid.resize(nfnew);
2460 //
2461 for (Index ilo = 0; ilo < nlo; ilo++) {
2462 for (Index i = 0; i < srfgrid[ilo].nelem(); i++) {
2463 sensor_response_f_grid[cumsumf[ilo] + i] = srfgrid[ilo][i];
2464 }
2465 }
2466
2467 // Append data to create total sensor_response
2468 //
2469 const Index ncols = sr[0].ncols();
2470 const Index npolnew = sensor_response_pol_grid.nelem();
2471 const Index nfpolnew = nfnew * npolnew;
2472 //
2473 sensor_response.resize(nlos * nfpolnew, ncols);
2474 //
2475 Vector dummy(ncols, 0.0);
2476 //
2477 for (Index ilo = 0; ilo < nlo; ilo++) {
2478 const Index nfpolthis = (cumsumf[ilo + 1] - cumsumf[ilo]) * npolnew;
2479
2480 ARTS_ASSERT(sr[ilo].nrows() == nlos * nfpolthis);
2481 ARTS_ASSERT(sr[ilo].ncols() == ncols);
2482
2483 for (Index ilos = 0; ilos < nlos; ilos++) {
2484 for (Index i = 0; i < nfpolthis; i++) {
2485 // "Poor mans" transfer of a row from one sparse to another
2486 for (Index ic = 0; ic < ncols; ic++) {
2487 dummy[ic] = sr[ilo](ilos * nfpolthis + i, ic);
2488 }
2489
2490 sensor_response.insert_row(ilos * nfpolnew + cumsumf[ilo] * npolnew + i,
2491 dummy);
2492 }
2493 }
2494 }
2495
2496 // Set aux variables
2497 sensor_aux_vectors(sensor_response_f,
2498 sensor_response_pol,
2499 sensor_response_dlos,
2500 sensor_response_f_grid,
2501 sensor_response_pol_grid,
2502 sensor_response_dlos_grid);
2503}
2504
2505
2506
2507/* Workspace method: Doxygen documentation will be auto-generated */
2509 Vector& sensor_response_f,
2510 ArrayOfIndex& sensor_response_pol,
2511 Matrix& sensor_response_dlos,
2512 ArrayOfIndex& sensor_response_pol_grid,
2513 const Vector& sensor_response_f_grid,
2514 const Matrix& sensor_response_dlos_grid,
2515 const Index& stokes_dim,
2516 const String& iy_unit,
2517 const ArrayOfIndex& instrument_pol,
2518 const Verbosity&) {
2519 // Some sizes
2520 const Index nnew = instrument_pol.nelem();
2521 const Index nf = sensor_response_f_grid.nelem();
2522 const Index npol = sensor_response_pol_grid.nelem();
2523 const Index nlos = sensor_response_dlos_grid.nrows();
2524
2525 // Initialise an output stream for runtime errors and a flag for errors
2526 ostringstream os;
2527 bool error_found = false;
2528
2529 Index nfz = nf * nlos;
2530 Index nin = nfz * npol;
2531
2532 if (sensor_response.nrows() != nin) {
2533 os << "The sensor block response matrix *sensor_response* does not have\n"
2534 << "right size compared to the sensor grid variables\n"
2535 << "(sensor_response_f_grid etc.).\n";
2536 error_found = true;
2537 }
2538
2539 // Check that sensor_response variables are consistent in size
2540 if (sensor_response_f.nelem() != nin) {
2541 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
2542 << "grid variables (sensor_response_f_grid etc.).\n";
2543 error_found = true;
2544 }
2545 if (npol != stokes_dim) {
2546 os << "Number of input polarisation does not match *stokes_dim*.\n";
2547 error_found = true;
2548 }
2549 if (nnew == 0) {
2550 os << "The WSV *instrument_pol* can not be empty.\n";
2551 error_found = true;
2552 }
2553 // If errors where found throw runtime_error with the collected error
2554 // message (before it gets too long)
2555 if (error_found) throw runtime_error(os.str());
2556
2557 // Check polarisation data more in detail
2558 for (Index i = 0; i < npol && !error_found; i++) {
2559 if (sensor_response_pol_grid[i] != i + 1) {
2560 os << "The input polarisations must be I, Q, U and V (up to "
2561 << "stokes_dim). It seems that input data are for other "
2562 << "polarisation components.";
2563 error_found = true;
2564 }
2565 }
2566 for (Index i = 0; i < nnew && !error_found; i++) {
2567 if (instrument_pol[i] < 1 || instrument_pol[i] > 10) {
2568 os << "The elements of *instrument_pol* must be inside the range [1,10].\n";
2569 error_found = true;
2570 }
2571 }
2572 // If errors where found throw runtime_error with the collected error
2573 // message (before it gets too long)
2574 if (error_found) throw runtime_error(os.str());
2575
2576 // If errors where found throw runtime_error with the collected error
2577 // message
2578 if (error_found) throw runtime_error(os.str());
2579
2580 // Normalisation weight to apply
2581 Numeric w = 0.5;
2582 if (iy_unit == "PlanckBT" || iy_unit == "RJBT") {
2583 w = 1.0;
2584 }
2585
2586 // Form H matrix representing polarisation response
2587 //
2588 Sparse Hpol(nfz * nnew, nin);
2589 Vector hrow(nin, 0);
2590 Index row = 0;
2591 //
2592 for (Index i = 0; i < nfz; i++) {
2593 Index col = i * npol;
2594 for (Index in = 0; in < nnew; in++) {
2595 /* Old code, matching older version of stokes2pol
2596 Index p = instrument_pol[in] - 1;
2597 //
2598 for( Index iv=0; iv<pv[p].nelem(); iv++ )
2599 { hrow[col+iv] = pv[p][iv]; }
2600 */
2601 stokes2pol(
2602 hrow[Range(col, stokes_dim)], stokes_dim, instrument_pol[in], w);
2603 //
2604 Hpol.insert_row(row, hrow);
2605 //
2606 hrow = 0;
2607 row += 1;
2608 }
2609 }
2610
2611 // Here we need a temporary sparse that is copy of the sensor_response
2612 // sparse matrix. We need it since the multiplication function can not
2613 // take the same object as both input and output.
2614 Sparse Htmp = sensor_response;
2615 sensor_response.resize(Hpol.nrows(), Htmp.ncols());
2616 mult(sensor_response, Hpol, Htmp);
2617
2618 // Update sensor_response_pol_grid
2619 sensor_response_pol_grid = instrument_pol;
2620
2621 // Set aux variables
2622 sensor_aux_vectors(sensor_response_f,
2623 sensor_response_pol,
2624 sensor_response_dlos,
2625 sensor_response_f_grid,
2626 sensor_response_pol_grid,
2627 sensor_response_dlos_grid);
2628}
2629
2630
2631
2632/* Workspace method: Doxygen documentation will be auto-generated */
2634 const Vector& sensor_response_f_grid,
2635 const ArrayOfIndex& sensor_response_pol_grid,
2636 const Matrix& sensor_response_dlos_grid,
2637 const Index& stokes_dim,
2638 const Vector& stokes_rotation,
2639 const Verbosity&) {
2640 // Basic checks
2641 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2642
2643 // Some sizes
2644 const Index nf = sensor_response_f_grid.nelem();
2645 const Index npol = sensor_response_pol_grid.nelem();
2646 const Index nlos = sensor_response_dlos_grid.nrows();
2647 const Index nin = nf * npol * nlos;
2648
2649 //---------------------------------------------------------------------------
2650 // Initialise a output stream for runtime errors and a flag for errors
2651 ostringstream os;
2652 bool error_found = false;
2653
2654 // Check that sensor_response variables are consistent in size
2655 if (sensor_response.nrows() != nin) {
2656 os << "The sensor block response matrix *sensor_response* does not have\n"
2657 << "right size compared to the sensor grid variables\n"
2658 << "(sensor_response_f_grid etc.).\n";
2659 error_found = true;
2660 }
2661
2662 // Check special stuff for this method
2663 if (stokes_dim < 3) {
2664 os << "To perform a rotation of the Stokes coordinate system,\n"
2665 << "*stokes_dim* must be >= 3.\n";
2666 error_found = true;
2667 }
2668 if (stokes_rotation.nelem() != nlos) {
2669 os << "Incorrect number of angles in *stokes_rotation*. The length\n"
2670 << "of this matrix must match *sensor_response_dlos_grid*.\n";
2671 error_found = true;
2672 }
2673 if (npol != stokes_dim) {
2674 os << "Inconsistency detected. The length of *sensor_response_pol_grid*\n"
2675 << "must be equal to *stokes_dim*, and this is not the case.\n";
2676 error_found = true;
2677 }
2678 for (Index is = 0; is < npol; is++) {
2679 if (sensor_response_pol_grid[is] != is + 1) {
2680 os << "For this method, the values in *sensor_response_pol_grid* must\n"
2681 << "be 1,2...stokes_dim. This is not the case, indicating that\n"
2682 << "some previous sensor part has that the data no longer are\n"
2683 << "Stokes vectors.\n";
2684 error_found = true;
2685 break;
2686 }
2687 }
2688
2689 // If errors where found throw runtime_error with the collected error
2690 // message.
2691 if (error_found) throw runtime_error(os.str());
2692 //---------------------------------------------------------------------------
2693
2694 // Set up complete the H matrix for applying rotation
2695 //
2696 Sparse H(sensor_response.nrows(), sensor_response.ncols());
2697 {
2698 Sparse Hrot(npol, npol); // Mueller matrix for 1 Stokes vec
2699 Vector row(H.ncols(), 0);
2700 Index irow = 0;
2701 //
2702 for (Index ilos = 0; ilos < nlos; ilos++) {
2703 // Rotation matrix for direction of concern
2704 mueller_rotation(Hrot, npol, stokes_rotation[ilos]);
2705
2706 for (Index ifr = 0; ifr < nf; ifr++) {
2707 for (Index ip = 0; ip < npol; ip++) {
2708 // Fill relevant part of row with matching (complete) row
2709 // in Hrot, and instert this row in H
2710 for (Index is = 0; is < npol; is++) {
2711 row[irow + is] = Hrot.ro(ip, is);
2712 }
2713 H.insert_row(irow + ip, row);
2714 // Re-zero row.
2715 for (Index is = 0; is < npol; is++) {
2716 row[irow + is] = 0;
2717 }
2718 }
2719 // Update irow, i.e. jump to next frequency
2720 irow += npol;
2721 }
2722 }
2723 }
2724
2725 // Here we need a temporary sparse that is copy of the sensor_response
2726 // sparse matrix. We need it since the multiplication function can not
2727 // take the same object as both input and output.
2728 Sparse Htmp = sensor_response;
2729 sensor_response.resize(Htmp.nrows(), Htmp.ncols()); //Just in case!
2730 mult(sensor_response, H, Htmp);
2731}
2732
2733
2734
2735/* Workspace method: Doxygen documentation will be auto-generated */
2737 Vector& f_grid,
2738 Index& antenna_dim,
2739 Matrix& mblock_dlos_grid,
2740 Sparse& sensor_response,
2741 Vector& sensor_response_f,
2742 ArrayOfIndex& sensor_response_pol,
2743 Matrix& sensor_response_dlos,
2744 Vector& sensor_response_f_grid,
2745 ArrayOfIndex& sensor_response_pol_grid,
2746 Matrix& sensor_response_dlos_grid,
2747 Index& sensor_norm,
2748 // WS Input:
2749 const Index& atmosphere_dim,
2750 const Index& stokes_dim,
2751 const Matrix& sensor_description_amsu,
2752 // WS Generic Input:
2753 const Numeric& spacing,
2754 const Verbosity& verbosity) {
2755 // Number of instrument channels:
2756 const Index n = sensor_description_amsu.nrows();
2757 const Index m = sensor_description_amsu.ncols();
2758
2759 // FIXME - Oscar Isoz-090413
2760 // add error checks
2761 //
2762
2763 // The meaning of the columns in sensor_description_amsu is:
2764 // LO frequency, channel center offset from LO, channel width, second offset from LO (to be extended if needed);
2765 // (All in Hz.)
2766 //
2767 if (5 > sensor_description_amsu.ncols()) {
2768 ostringstream os;
2769 os << "Input variable sensor_description_amsu must have atleast five columns, but it has "
2770 << sensor_description_amsu.ncols() << ".";
2771 throw runtime_error(os.str());
2772 }
2773
2774 ConstVectorView lo_multi = sensor_description_amsu(Range(joker), 0);
2775 ConstMatrixView offset = sensor_description_amsu(
2776 Range(joker), Range(1, m - 2)); // Remember to ignore column 2..
2777 ConstVectorView verbosityVectIn =
2778 sensor_description_amsu(Range(joker), m - 1);
2779 ConstVectorView width = sensor_description_amsu(Range(joker), m - 2);
2780
2781 //is there any undefined verbosity values in the vector?
2782 //Set the verbosity to one third of the bandwidth to make sure that the passband flanks does't overlap
2783 const Numeric minRatioVerbosityVsFdiff =
2784 10; // To be used when to passbands are closer then one verbosity value
2785 Index totNumPb = 0;
2786 Vector verbosityVect(n);
2787
2788 for (Index idx = 0; idx < n; ++idx) {
2789 if ((verbosityVectIn[idx] == 0) || (verbosityVectIn[idx] > width[idx])) {
2790 verbosityVect[idx] = ((Numeric)width[idx]) / 3;
2791 } else {
2792 verbosityVect[idx] = verbosityVectIn[idx];
2793 }
2794 }
2795
2796 // Create a vector to store the number of passbands (PB) for each channel
2797 ArrayOfIndex numPBpseudo(n); // store values used for calc
2798 ArrayOfIndex numPB(n); // Store the true values
2799 // Find the number of IFs for each channel and calculate the number of passbands
2800 for (Index i = 0; i < n; ++i) {
2801 numPB[i] = 0; // make sure that it is zero
2802 for (Index j = 0; j < (m - 2); ++j) {
2803 if (j != 2) {
2804 if (offset(i, j) > 0) {
2805 numPB[i]++;
2806 }
2807 }
2808 }
2809 numPB[i] = 1 << (int)numPB[i]; // number of passbands= 2^numLO
2810 if (numPB[i] == 1) {
2811 numPBpseudo[i] = 2;
2812 } else {
2813 numPBpseudo[i] = numPB[i];
2814 }
2815
2816 totNumPb += (int)numPB[i];
2817
2818 if (numPB[i] > 4) {
2819 ostringstream os;
2820 os << "This function does currently not support more than 4 passbands per channel"
2821 << numPB[i] << ".";
2822 throw runtime_error(os.str());
2823 }
2824 }
2825
2826 // Find the center frequencies for all sub-channels
2827 // Create one center frequency for each passband
2828 ArrayOfArrayOfGriddedField1 backend_channel_response_multi(n);
2829 ArrayOfVector f_backend_multi(n); // changed !!!
2830 for (Index i = 0; i < n; ++i) {
2831 // Channel frequencies
2832 Vector& f = f_backend_multi[i];
2833 f.resize(1);
2834 f[0] = lo_multi[i] + 0.0 * width[i]; //(offset(i,0)+offset(i,1));
2835
2836 //channel response
2837 const Index numVal = 4;
2838 backend_channel_response_multi[i].resize(1);
2839 GriddedField1& b_resp = backend_channel_response_multi[i][0];
2840 b_resp.set_name("Backend channel response function");
2841 b_resp.resize(numVal * numPBpseudo[i]);
2842 Vector f_range(numVal * numPBpseudo[i]);
2843 Numeric pbOffset = 0;
2844 b_resp.set_grid_name(0, "Frequency");
2845
2846 Numeric slope = 0; // 1900;
2847 // To avoid overlapping passbands in the AMSU-A sensor, reduce the passbands of each channel by a few Hz
2848 for (Index pbOffsetIdx = 0; pbOffsetIdx < numPBpseudo[i];
2849 ++pbOffsetIdx) { // Filter response
2850 slope = width[i] / 100;
2851 f_range[pbOffsetIdx * numVal + 0] = -0.5 * width[i] - 0 * slope;
2852 f_range[pbOffsetIdx * numVal + 1] = -0.5 * width[i] + 1 * slope;
2853 f_range[pbOffsetIdx * numVal + 2] = +0.5 * width[i] - 1 * slope;
2854 f_range[pbOffsetIdx * numVal + 3] = +0.5 * width[i] + 0 * slope;
2855
2856 b_resp.data[pbOffsetIdx * numVal + 0] = 0.0 / (Numeric)numPB[i];
2857 ;
2858 b_resp.data[pbOffsetIdx * numVal + 1] = 1.0 / (Numeric)numPB[i];
2859 b_resp.data[pbOffsetIdx * numVal + 2] = 1.0 / (Numeric)numPB[i];
2860 b_resp.data[pbOffsetIdx * numVal + 3] = 0.0 / (Numeric)numPB[i];
2861
2862 if (numPB[i] == 1) {
2863 if (pbOffsetIdx == 0) {
2864 pbOffset = -0.0 * width[i];
2865 b_resp.data[pbOffsetIdx * numVal + 0] = 0;
2866 b_resp.data[pbOffsetIdx * numVal + 1] = 1.0 / 1;
2867
2868 b_resp.data[pbOffsetIdx * numVal + 2] = 1.0 / 1;
2869 b_resp.data[pbOffsetIdx * numVal + 3] = 1.0 / 1;
2870 f_range[pbOffsetIdx * numVal + 0] = -0.5 * width[i] - 2 * slope;
2871 f_range[pbOffsetIdx * numVal + 1] = -0.5 * width[i] - 1 * slope;
2872 f_range[pbOffsetIdx * numVal + 2] = -0.5 * width[i] + 1 * slope;
2873 f_range[pbOffsetIdx * numVal + 3] = -0.5 * width[i] + 2 * slope;
2874 }
2875 if (pbOffsetIdx == 1) {
2876 pbOffset = 0.0 * width[i]; //just a dummy band
2877 b_resp.data[pbOffsetIdx * numVal + 0] = 1.0 / 1;
2878 b_resp.data[pbOffsetIdx * numVal + 1] = 1.0 / 1;
2879 b_resp.data[pbOffsetIdx * numVal + 2] = 1.0 / 1;
2880 b_resp.data[pbOffsetIdx * numVal + 3] = 0;
2881 f_range[pbOffsetIdx * numVal + 0] = +0.5 * width[i] - 3 * slope;
2882 f_range[pbOffsetIdx * numVal + 1] = +0.5 * width[i] - 2 * slope;
2883 f_range[pbOffsetIdx * numVal + 2] = +0.5 * width[i] - 1 * slope;
2884 f_range[pbOffsetIdx * numVal + 3] = +0.5 * width[i] + 0 * slope - 10;
2885 // without the extra '-10' it will fail due to too narrow backend sensor response.
2886 }
2887 } else if (
2888 numPB[i] ==
2889 2) { // move the passband in frequency to the correct frequency
2890 if (pbOffsetIdx == 0) {
2891 pbOffset = -offset(i, 0);
2892 }
2893 if (pbOffsetIdx == 1) {
2894 pbOffset = offset(i, 0);
2895 }
2896 }
2897 if (numPB[i] == 4) {
2898 if (pbOffsetIdx == 0) {
2899 pbOffset = -offset(i, 0) - offset(i, 1);
2900 }
2901 if (pbOffsetIdx == 1) {
2902 pbOffset = -offset(i, 0) + offset(i, 1);
2903 }
2904 if (pbOffsetIdx == 2) {
2905 pbOffset = offset(i, 0) - offset(i, 1);
2906 }
2907 if (pbOffsetIdx == 3) {
2908 pbOffset = offset(i, 0) + offset(i, 1);
2909 }
2910 }
2911 for (Index iii = 0; iii < numVal; ++iii) {
2912 f_range[pbOffsetIdx * numVal + iii] += 1 * pbOffset;
2913 }
2914 }
2915 // Are any passbands overlapping?
2916 for (Index ii = 2; ii < (f_range.nelem() - 2); ++ii) {
2917 if (((b_resp.data[ii - 1] == 1) && (b_resp.data[ii] == 0) &&
2918 (b_resp.data[ii + 1] == 0) && (b_resp.data[ii + 2] == 1))) {
2919 if ((f_range[ii] >= f_range[ii + 1])) // Overlapping passbands
2920 {
2921 if ((f_range[ii + 2] - f_range[ii - 1]) >
2922 verbosityVectIn[i]) // difference in frequency between passbands
2923 {
2924 f_range[ii + 1] = f_range[ii + 2] - verbosityVectIn[i] / 2;
2925 f_range[ii] = f_range[ii - 1] + verbosityVectIn[i] / 2;
2926 } else {
2927 f_range[ii - 1] = (f_range[ii] + f_range[ii + 2]) / 2 -
2928 2 * verbosityVectIn[i] / minRatioVerbosityVsFdiff;
2929 f_range[ii + 1] = (f_range[ii] + f_range[ii + 2]) / 2 +
2930 verbosityVectIn[i] / minRatioVerbosityVsFdiff;
2931 f_range[ii] =
2932 f_range[ii - 1] + verbosityVectIn[i] / minRatioVerbosityVsFdiff;
2933 f_range[ii + 2] =
2934 f_range[ii + 1] + verbosityVectIn[i] / minRatioVerbosityVsFdiff;
2935 }
2936 }
2937 }
2938 }
2939 b_resp.set_grid(0, f_range);
2940 }
2941
2942 // construct sideband response
2943 ArrayOfGriddedField1 sideband_response_multi(n);
2944 for (Index i = 0; i < n; ++i) {
2945 GriddedField1& r = sideband_response_multi[i];
2946 r.set_name("Sideband response function");
2947 r.resize(numPBpseudo[i]);
2948 Vector f(numPBpseudo[i]);
2949 if (numPB[i] == 1) {
2950 r.data[0] = 0.5;
2951 f[0] = -.0 * width[i];
2952 r.data[1] = 0.5;
2953 f[1] = .0 * width[i];
2954 } else if (numPB[i] == 2) {
2955 r.data[0] = 1. / (Numeric)numPB[i];
2956 r.data[1] = 1. / (Numeric)numPB[i];
2957 f[0] = -1 * offset(i, 0) - 0.5 * width[i];
2958 f[1] = +1 * offset(i, 0) + 0.5 * width[i];
2959 } else if (numPB[i] == 4) {
2960 r.data[0] = 1. / (Numeric)numPB[i];
2961 r.data[1] = 1. / (Numeric)numPB[i];
2962 r.data[2] = 1. / (Numeric)numPB[i];
2963 r.data[3] = 1. / (Numeric)numPB[i];
2964 f[0] = -offset(i, 0) - offset(i, 1) - 0.5 * width[i];
2965 ;
2966 f[1] = -offset(i, 0) + offset(i, 1) - 0.5 * width[i];
2967 ;
2968 f[2] = +offset(i, 0) - offset(i, 1) + 0.5 * width[i];
2969 ;
2970 f[3] = +offset(i, 0) + offset(i, 1) + 0.5 * width[i];
2971 ;
2972 }
2973 r.set_grid_name(0, "Frequency");
2974 r.set_grid(0, f);
2975 }
2976
2977 sensor_norm = 1;
2979 // out
2980 f_grid,
2981 // in
2982 f_backend_multi,
2983 backend_channel_response_multi,
2984 spacing,
2985 verbosityVect,
2986 verbosity);
2987
2988 // do some final work...
2989 AntennaOff(antenna_dim, mblock_dlos_grid, verbosity);
2990
2992 // out
2993 sensor_response,
2994 sensor_response_f,
2995 sensor_response_pol,
2996 sensor_response_dlos,
2997 sensor_response_f_grid,
2998 sensor_response_pol_grid,
2999 sensor_response_dlos_grid,
3000 // in
3001 f_grid,
3002 mblock_dlos_grid,
3003 antenna_dim,
3004 atmosphere_dim,
3005 stokes_dim,
3006 sensor_norm,
3007 verbosity);
3008
3009 Index numLO = lo_multi.nelem();
3010 // Variables for data to be appended
3011 // Based on code from m_sensor->sensor_responseMultiMixerBackend()
3012 Array<Sparse> sr;
3013 ArrayOfVector srfgrid;
3014 ArrayOfIndex cumsumf(numLO + 1, 0);
3015 const Index nlos = sensor_response_dlos_grid.nrows();
3016
3017 // Do this for all channels ....
3018 for (Index idxLO = 0; idxLO < numLO; idxLO++) {
3019 Sparse sr1 = sensor_response;
3020 Vector srf1 = sensor_response_f;
3021 ArrayOfIndex srpol1 = sensor_response_pol;
3022 Matrix srdlos1 = sensor_response_dlos;
3023 Vector srfgrid1 = sensor_response_f_grid;
3024
3026 sr1,
3027 srf1,
3028 srpol1,
3029 srdlos1,
3030 srfgrid1,
3031 //in
3032 sensor_response_pol_grid,
3033 sensor_response_dlos_grid,
3034 f_backend_multi[idxLO],
3035 backend_channel_response_multi[idxLO],
3036 sensor_norm,
3037 verbosity);
3038
3039 // Store in temporary arrays
3040 sr.push_back(sr1);
3041 srfgrid.push_back(srfgrid1);
3042 //
3043 cumsumf[idxLO + 1] = cumsumf[idxLO] + srfgrid1.nelem();
3044 }
3045
3046 // Append data to create sensor_response_f_grid
3047 //
3048 const Index nfnew = cumsumf[numLO];
3049 sensor_response_f_grid.resize(nfnew);
3050 //
3051 for (Index ilo = 0; ilo < numLO; ilo++) {
3052 for (Index i = 0; i < srfgrid[ilo].nelem(); i++) {
3053 sensor_response_f_grid[cumsumf[ilo] + i] = srfgrid[ilo][i];
3054 }
3055 }
3056
3057 const Index ncols = sr[0].ncols();
3058 const Index npolnew = sensor_response_pol_grid.nelem();
3059 const Index nfpolnew = nfnew * npolnew;
3060 //
3061 sensor_response.resize(nlos * nfpolnew, ncols);
3062 //
3063 Vector dummy(ncols, 0.0);
3064 //
3065 for (Index ilo = 0; ilo < numLO; ilo++) {
3066 const Index nfpolthis = (cumsumf[ilo + 1] - cumsumf[ilo]) * npolnew;
3067
3068 ARTS_ASSERT(sr[ilo].nrows() == nlos * nfpolthis);
3069 ARTS_ASSERT(sr[ilo].ncols() == ncols);
3070
3071 for (Index ilos = 0; ilos < nlos; ilos++) {
3072 for (Index i = 0; i < nfpolthis; i++) {
3073 // "Poor mans" transfer of a row from one sparse to another
3074 for (Index ic = 0; ic < ncols; ic++) {
3075 dummy[ic] = sr[ilo](ilos * nfpolthis + i, ic);
3076 }
3077
3078 sensor_response.insert_row(ilos * nfpolnew + cumsumf[ilo] * npolnew + i,
3079 dummy);
3080 }
3081 }
3082 }
3083
3084 sensor_aux_vectors(sensor_response_f,
3085 sensor_response_pol,
3086 sensor_response_dlos,
3087 sensor_response_f_grid,
3088 sensor_response_pol_grid,
3089 sensor_response_dlos_grid);
3090}
3091
3092
3093
3094/* Workspace method: Doxygen documentation will be auto-generated */
3095void sensor_responseSimpleAMSU( // WS Output:
3096 Vector& f_grid,
3097 Index& antenna_dim,
3098 Matrix& mblock_dlos_grid,
3099 Sparse& sensor_response,
3100 Vector& sensor_response_f,
3101 ArrayOfIndex& sensor_response_pol,
3102 Matrix& sensor_response_dlos,
3103 Vector& sensor_response_f_grid,
3104 ArrayOfIndex& sensor_response_pol_grid,
3105 Matrix& sensor_response_dlos_grid,
3106 Index& sensor_norm,
3107 // WS Input:
3108 const Index& atmosphere_dim,
3109 const Index& stokes_dim,
3110 const Matrix& sensor_description_amsu,
3111 // WS Generic Input:
3112 const Numeric& spacing,
3113 const Verbosity& verbosity) {
3114 // Check that sensor_description_amsu has the right dimension:
3115 if (3 != sensor_description_amsu.ncols()) {
3116 ostringstream os;
3117 os << "Input variable sensor_description_amsu must have three columns, but it has "
3118 << sensor_description_amsu.ncols() << ".";
3119 throw runtime_error(os.str());
3120 }
3121
3122 // Number of instrument channels:
3123 const Index n = sensor_description_amsu.nrows();
3124
3125 // The meaning of the columns in sensor_description_amsu is:
3126 // LO frequency, channel center offset from LO, channel width.
3127 // (All in Hz.)
3128 ConstVectorView lo_multi = sensor_description_amsu(Range(joker), 0);
3129 ConstVectorView offset = sensor_description_amsu(Range(joker), 1);
3130 ConstVectorView width = sensor_description_amsu(Range(joker), 2);
3131
3132 // Channel frequencies:
3133 ArrayOfVector f_backend_multi(n);
3134 for (Index i = 0; i < n; ++i) {
3135 Vector& f = f_backend_multi[i];
3136 f.resize(1);
3137 f[0] = lo_multi[i] + offset[i];
3138 }
3139
3140 // Construct channel response
3141 ArrayOfArrayOfGriddedField1 backend_channel_response_multi(n);
3142 for (Index i = 0; i < n; ++i) {
3143 backend_channel_response_multi[i].resize(1);
3144 GriddedField1& r = backend_channel_response_multi[i][0];
3145 r.set_name("Backend channel response function");
3146 r.resize(2);
3147
3148 // Frequency range:
3149 Vector f(2);
3150 f[0] = -0.5 * width[i];
3151 f[1] = +0.5 * width[i];
3152 r.set_grid_name(0, "Frequency");
3153 r.set_grid(0, f);
3154
3155 // Response:
3156 r.data[0] = 1;
3157 r.data[1] = 1;
3158 }
3159
3160 // Construct sideband response:
3161 ArrayOfGriddedField1 sideband_response_multi(n);
3162 for (Index i = 0; i < n; ++i) {
3163 GriddedField1& r = sideband_response_multi[i];
3164 r.set_name("Sideband response function");
3165 r.resize(2);
3166
3167 // Frequency range:
3168 Vector f(2);
3169 f[0] = -(offset[i] + 0.5 * width[i]);
3170 f[1] = +(offset[i] + 0.5 * width[i]);
3171 r.set_grid_name(0, "Frequency");
3172 r.set_grid(0, f);
3173
3174 // Response:
3175 r.data[0] = 0.5;
3176 r.data[1] = 0.5;
3177 }
3178
3179 // Set sideband mode:
3180 ArrayOfString sideband_mode_multi(n, "upper");
3181
3182 // We want to automatically normalize the sensor response data, so set sensor_norm to 1:
3183 sensor_norm = 1;
3184
3185 // Now the rest is just to use some workspace methods:
3186 // ---------------------------------------------------
3187
3188 f_gridFromSensorAMSU(f_grid,
3189 lo_multi,
3190 f_backend_multi,
3191 backend_channel_response_multi,
3192 spacing,
3193 verbosity);
3194
3195 AntennaOff(antenna_dim, mblock_dlos_grid, verbosity);
3196
3197 sensor_responseInit(sensor_response,
3198 sensor_response_f,
3199 sensor_response_pol,
3200 sensor_response_dlos,
3201 sensor_response_f_grid,
3202 sensor_response_pol_grid,
3203 sensor_response_dlos_grid,
3204 f_grid,
3205 mblock_dlos_grid,
3206 antenna_dim,
3207 atmosphere_dim,
3208 stokes_dim,
3209 sensor_norm,
3210 verbosity);
3211
3212 sensor_responseMultiMixerBackend(sensor_response,
3213 sensor_response_f,
3214 sensor_response_pol,
3215 sensor_response_dlos,
3216 sensor_response_f_grid,
3217 sensor_response_pol_grid,
3218 sensor_response_dlos_grid,
3219 lo_multi,
3220 sideband_response_multi,
3221 sideband_mode_multi,
3222 f_backend_multi,
3223 backend_channel_response_multi,
3224 sensor_norm,
3225 verbosity);
3226}
3227
3228
3229/* Workspace method: Doxygen documentation will be auto-generated */
3230void WMRFSelectChannels( // WS Output:
3231 Vector& f_grid,
3232 Sparse& wmrf_weights,
3233 Vector& f_backend,
3234 // WS Input:
3235 const ArrayOfIndex& wmrf_channels,
3236 const Verbosity& verbosity) {
3239
3240 // For error messages:
3241 ostringstream os;
3242
3243 // Some checks of input parameters:
3244
3245 // wmrf_weights must have same number of rows as f_backend, and same
3246 // number of columns as f_grid.
3247 if ((wmrf_weights.nrows() != f_backend.nelem()) ||
3248 (wmrf_weights.ncols() != f_grid.nelem())) {
3249 os << "The WSV *wmrf_weights* must have same number of rows as\n"
3250 << "*f_backend*, and same number of columns as *f_grid*.\n"
3251 << "wmrf_weights.nrows() = " << wmrf_weights.nrows() << "\n"
3252 << "f_backend.nelem() = " << f_backend.nelem() << "\n"
3253 << "wmrf_weights.ncols() = " << wmrf_weights.ncols() << "\n"
3254 << "f_grid.nelem() = " << f_grid.nelem();
3255 throw runtime_error(os.str());
3256 }
3257
3258 // wmrf_channels must be strictly increasing (no repetitions).
3259 chk_if_increasing("wmrf_channels", wmrf_channels);
3260
3261 // All selected channels must be within the original set of
3262 // channels.
3263 if (min(wmrf_channels) < 0) {
3264 os << "Min(wmrf_channels) must be >= 0, but it is " << min(wmrf_channels)
3265 << ".";
3266 }
3267 if (max(wmrf_channels) >= f_backend.nelem()) {
3268 os << "Max(wmrf_channels) must be less than the total number of channels.\n"
3269 << "(We use zero-based indexing!)\n"
3270 << "The actual value you have is " << max(wmrf_channels) << ".";
3271 }
3272
3273 if (wmrf_channels.nelem() == f_backend.nelem()) {
3274 // No channels to be removed, I can return the original grid.
3275 out2 << " Retaining all channels.\n";
3276 } else {
3277 out2 << " Reducing number of channels from " << f_backend.nelem() << " to "
3278 << wmrf_channels.nelem() << ".\n";
3279 }
3280
3281 // Now the real work starts:
3282
3283 // 1. Remove unwanted channels from f_backend:
3284 Select(f_backend, f_backend, wmrf_channels, verbosity);
3285
3286 // 2. Remove unwanted channels from wmrf_weights. (We also have to
3287 // do something about the frequency dimension of wmrf_weights, but
3288 // we'll do that later.)
3289 Select(wmrf_weights, wmrf_weights, wmrf_channels, verbosity);
3290
3291 // 3. Identify, which frequencies are still needed, and which are
3292 // now obsolete. We store the still needed frequencies in an
3293 // ArrayOfIndex.
3294
3295 // Create f_grid_array. We do not store the frequencies themselves,
3296 // but the indices of the frequencies to use.
3297 ArrayOfIndex selection;
3298 // Make sure that selection does not have to be reallocated along
3299 // the way. (This is purely to improve performance a bit.)
3300 selection.reserve(f_grid.nelem());
3301
3302 // Go through f_grid, and check for each frequency whether it is in
3303 // the set of WMRF frequencies for any of the channels.
3304 ARTS_ASSERT(wmrf_weights.nrows() == f_backend.nelem());
3305 ARTS_ASSERT(wmrf_weights.ncols() == f_grid.nelem());
3306 for (Index fi = 0; fi < wmrf_weights.ncols(); ++fi) {
3307 Index i;
3308 for (i = 0; i < wmrf_weights.nrows(); ++i) {
3309 if (wmrf_weights(i, fi) != 0) {
3310 selection.push_back(fi);
3311 break;
3312 }
3313 }
3314 if (i == wmrf_weights.nrows()) {
3315 out3 << " The frequency with index " << fi
3316 << " is not used by any channel.\n";
3317 }
3318 }
3319
3320 if (selection.nelem() == f_grid.nelem()) {
3321 // No frequencies were removed, I can return the original grid.
3322 out2 << " No unnecessary frequencies, leaving f_grid untouched.\n";
3323 } else if (selection.nelem() == 0) {
3324 throw runtime_error("No frequencies found for any selected channels.\n");
3325 } else {
3326 out2 << " Reducing number of frequency grid points from " << f_grid.nelem()
3327 << " to " << selection.nelem() << ".\n";
3328 }
3329
3330 // 4. Select the right frequencies in f_grid:
3331 Select(f_grid, f_grid, selection, verbosity);
3332
3333 // 5. Select the right frequencies in wmrf_weights. This is a bit
3334 // tricky, since Select works on the row dimension. So we have to
3335 // take the transpose.
3336 Sparse wt(wmrf_weights.ncols(), wmrf_weights.nrows());
3337 transpose(wt, wmrf_weights);
3338 Select(wt, wt, selection, verbosity);
3339 wmrf_weights.resize(wt.ncols(), wt.nrows());
3340 transpose(wmrf_weights, wt);
3341}
3342
3343
3344
3345/* Workspace method: Doxygen documentation will be auto-generated */
3346void sensor_responseWMRF( // WS Output:
3347 Sparse& sensor_response,
3348 Vector& sensor_response_f,
3349 ArrayOfIndex& sensor_response_pol,
3350 Matrix& sensor_response_dlos,
3351 Vector& sensor_response_f_grid,
3352 // WS Input:
3353 const ArrayOfIndex& sensor_response_pol_grid,
3354 const Matrix& sensor_response_dlos_grid,
3355 const Sparse& wmrf_weights,
3356 const Vector& f_backend,
3357 const Verbosity& verbosity) {
3359
3360 // Some sizes
3361 const Index nf = sensor_response_f_grid.nelem();
3362 const Index npol = sensor_response_pol_grid.nelem();
3363 const Index nlos = sensor_response_dlos_grid.nrows();
3364 const Index nin = nf * npol * nlos;
3365
3366 // Initialise output stream for runtime errors and a flag for errors
3367 ostringstream os;
3368 bool error_found = false;
3369
3370 // Check that sensor_response variables are consistent in size
3371 if (sensor_response_f.nelem() != nin) {
3372 os << "Inconsistency in size between *sensor_response_f* and the sensor\n"
3373 << "grid variables (sensor_response_f_grid etc.).\n";
3374 error_found = true;
3375 }
3376 if (sensor_response.nrows() != nin) {
3377 os << "The sensor block response matrix *sensor_response* does not have\n"
3378 << "right size compared to the sensor grid variables\n"
3379 << "(sensor_response_f_grid etc.).\n";
3380 error_found = true;
3381 }
3382
3383 if (nin == 0) {
3384 os << "One of f_grid, pol_grid, dlos_grid are empty. Sizes are: (" << nf
3385 << ", " << npol << ", " << nlos << ")"
3386 << "\n";
3387 error_found = true;
3388 }
3389
3390 // Check number of rows in WMRF weight matrix
3391 //
3392 const Index nrw = wmrf_weights.nrows();
3393 //
3394 if (nrw != f_backend.nelem()) {
3395 os << "The WSV *wmrf_weights* must have as many rows\n"
3396 << "as *f_backend* has elements.\n"
3397 << "wmrf_weights.nrows() = " << nrw << "\n"
3398 << "f_backend.nelem() = " << f_backend.nelem() << "\n";
3399 error_found = true;
3400 }
3401
3402 // Check number of columns in WMRF weight matrix
3403 //
3404 const Index ncw = wmrf_weights.ncols();
3405 //
3406 if (ncw != sensor_response_f_grid.nelem()) {
3407 os << "The WSV *wmrf_weights* must have as many columns\n"
3408 << "as *sensor_response_f_grid* has elements.\n"
3409 << "wmrf_weights.ncols() = " << ncw << "\n"
3410 << "sensor_response_f_grid.nelem() = " << sensor_response_f_grid.nelem()
3411 << "\n";
3412 error_found = true;
3413 }
3414
3415 // If errors where found throw runtime_error with the collected error
3416 // message (before error message gets too long).
3417 if (error_found) throw runtime_error(os.str());
3418
3419 // Ok, now the actual work.
3420
3421 // Here we need a temporary sparse that is copy of the sensor_response
3422 // sparse matrix. We need it since the multiplication function can not
3423 // take the same object as both input and output.
3424 Sparse htmp = sensor_response;
3425 sensor_response.resize(wmrf_weights.nrows(), htmp.ncols());
3426 mult(sensor_response, wmrf_weights, htmp);
3427
3428 // Some extra output.
3429 out3 << " Size of *sensor_response*: " << sensor_response.nrows() << "x"
3430 << sensor_response.ncols() << "\n";
3431
3432 // Update sensor_response_f_grid
3433 sensor_response_f_grid = f_backend;
3434
3435 // Set aux variables
3436 sensor_aux_vectors(sensor_response_f,
3437 sensor_response_pol,
3438 sensor_response_dlos,
3439 sensor_response_f_grid,
3440 sensor_response_pol_grid,
3441 sensor_response_dlos_grid);
3442}
3443
3444
3445
3446/* Workspace method: Doxygen documentation will be auto-generated */
3448 Vector& y_f,
3449 const Matrix& iy,
3450 const Index& stokes_dim,
3451 const Vector& f_grid,
3452 const Numeric& df,
3453 const Verbosity& verbosity) {
3454 // Some dummy values
3455 const Index sensor_norm = 1, atmosphere_dim = 1;
3456
3457 // Init sensor reponse
3458 //
3459 Index antenna_dim;
3460 Vector sensor_response_f, sensor_response_f_grid;
3461 Matrix mblock_dlos_grid, sensor_response_dlos_grid, sensor_response_dlos;
3462 ArrayOfIndex sensor_response_pol, sensor_response_pol_grid;
3463 Sparse sensor_response;
3464 //
3465 AntennaOff(antenna_dim, mblock_dlos_grid, verbosity);
3466
3467 sensor_responseInit(sensor_response,
3468 sensor_response_f,
3469 sensor_response_pol,
3470 sensor_response_dlos,
3471 sensor_response_f_grid,
3472 sensor_response_pol_grid,
3473 sensor_response_dlos_grid,
3474 f_grid,
3475 mblock_dlos_grid,
3476 antenna_dim,
3477 atmosphere_dim,
3478 stokes_dim,
3479 sensor_norm,
3480 verbosity);
3481
3482 // Center position of "channels"
3483 Vector f_backend;
3484 linspace(f_backend, f_grid[0] + df / 2, last(f_grid), df);
3485
3486 // Create channel response
3488 backend_channel_responseFlat(r, df, verbosity);
3489
3490 // New sensor response
3491 sensor_responseBackend(sensor_response,
3492 sensor_response_f,
3493 sensor_response_pol,
3494 sensor_response_dlos,
3495 sensor_response_f_grid,
3496 sensor_response_pol_grid,
3497 sensor_response_dlos_grid,
3498 f_backend,
3499 r,
3500 sensor_norm,
3501 verbosity);
3502
3503 // Some sizes
3504 const Index nf = f_grid.nelem();
3505 const Index n = sensor_response.nrows();
3506
3507 // Convert iy to a vector
3508 //
3509 Vector iyb(nf * stokes_dim);
3510 //
3511 for (Index is = 0; is < stokes_dim; is++) {
3512 iyb[Range(is, nf, stokes_dim)] = iy(joker, is);
3513 }
3514
3515 // y and y_f
3516 //
3517 y_f = sensor_response_f;
3518 y.resize(n);
3519 mult(y, sensor_response, iyb);
3520}
3521
3522
3523
3524/* Workspace method: Doxygen documentation will be auto-generated */
3526 Vector& y_f,
3527 ArrayOfIndex& y_pol,
3528 Matrix& y_pos,
3529 Matrix& y_los,
3530 ArrayOfVector& y_aux,
3531 Matrix& y_geo,
3532 Matrix& jacobian,
3533 const Index& stokes_dim,
3534 const Index& jacobian_do,
3535 const Matrix& sensor_pos,
3536 const Matrix& sensor_pol,
3537 const Verbosity&) {
3538 // Some sizes
3539 const Index n1 = y.nelem();
3540 const Index nm = sensor_pol.nrows();
3541 const Index nc = sensor_pol.ncols();
3542 const Index n2 = nm * nc;
3543
3544 // Check consistency of input data
3545 if (y.empty()) throw runtime_error("Input *y* is empty. Use *yCalc*");
3546 if (y_f.nelem() != n1)
3547 throw runtime_error("Lengths of input *y* and *y_f* are inconsistent.");
3548 if (y_pol.nelem() != n1)
3549 throw runtime_error("Lengths of input *y* and *y_pol* are inconsistent.");
3550 if (y_pos.nrows() != n1)
3551 throw runtime_error("Sizes of input *y* and *y_pos* are inconsistent.");
3552 if (y_los.nrows() != n1)
3553 throw runtime_error("Sizes of input *y* and *y_los* are inconsistent.");
3554 if (y_geo.nrows() != n1)
3555 throw runtime_error("Sizes of input *y* and *y_geo* are inconsistent.");
3556 if (jacobian_do) {
3557 if (jacobian.nrows() != n1)
3558 throw runtime_error("Sizes of *y* and *jacobian* are inconsistent.");
3559 }
3560
3561 // Checks associated with the Stokes vector
3562 if (stokes_dim < 3)
3563 throw runtime_error(
3564 "*stokes_dim* must be >= 3 to correctly apply a "
3565 "polarisation rotation.");
3566 if (n1 < stokes_dim)
3567 throw runtime_error("Length of input *y* smaller than *stokes_dim*.");
3568 for (Index i = 0; i < stokes_dim; i++) {
3569 if (y_pol[i] != i + 1)
3570 throw runtime_error(
3571 "*y* must hold Stokes element values. Data in "
3572 "*y_pol* indicates that this is not the case.");
3573 }
3574
3575 // Checks of sensor_pos
3576 if (sensor_pos.nrows() != nm)
3577 throw runtime_error(
3578 "Different number of rows in *sensor_pos* and *sensor_pol*.");
3579 if (n2 * stokes_dim != n1)
3580 throw runtime_error(
3581 "Number of columns in *sensor_pol* not consistent with "
3582 "length of *y* and value of *stokes_dim*.");
3583
3584 // Make copy of all y variables and jacobian
3585 const Vector y1 = y, y_f1 = y_f;
3586 const Matrix y_pos1 = y_pos, y_los1 = y_los, y_geo1 = y_geo;
3587 const ArrayOfIndex y_pol1 = y_pol;
3588 const ArrayOfVector y_aux1 = y_aux;
3589 Matrix jacobian1(0, 0);
3590 if (jacobian_do) {
3591 jacobian1 = jacobian;
3592 }
3593
3594 // Resize the y variables and jacobian
3595 y.resize(n2);
3596 y_f.resize(n2);
3597 y_pol.resize(n2);
3598 y_pos.resize(n2, y_pos1.ncols());
3599 y_los.resize(n2, y_los1.ncols());
3600 y_geo.resize(n2, y_geo1.ncols());
3601 for (Index a = 0; a < y_aux.nelem(); a++) y_aux[a].resize(n2);
3602 if (jacobian_do) {
3603 jacobian.resize(n2, jacobian1.ncols());
3604 }
3605
3606 for (Index r = 0; r < nm; r++) {
3607 for (Index c = 0; c < nc; c++) {
3608 const Index iout = r * nc + c;
3609 const Index iin = iout * stokes_dim;
3610
3611 const Numeric wq = cos(2 * DEG2RAD * sensor_pol(r, c));
3612 const Numeric wu = sin(2 * DEG2RAD * sensor_pol(r, c));
3613
3614 // Extract radiance for polarisation angle of concern
3615 y[iout] = y1[iin] + wq * y1[iin + 1] + wu * y1[iin + 2];
3616
3617 // Same operation for jacobian
3618 if (jacobian_do) {
3619 for (Index q = 0; q < jacobian.ncols(); q++)
3620 jacobian(iout, q) = jacobian1(iin, q) + wq * jacobian1(iin + 1, q) +
3621 wu * jacobian1(iin + 2, q);
3622 }
3623
3624 // Set y_pol
3625 y_pol[iout] = (Index)sensor_pol(r, c);
3626
3627 // For the rest, copy value matching I of in-data
3628 y_f[iout] = y_f1[iin];
3629 y_pos(iout, joker) = y_pos1(iin, joker);
3630 y_los(iout, joker) = y_los1(iin, joker);
3631 y_geo(iout, joker) = y_geo1(iin, joker);
3632 for (Index a = 0; a < y_aux.nelem(); a++) y_aux[a][iout] = y_aux1[a][iin];
3633 }
3634 }
3635}
The global header file for ARTS.
Index nrows
void * data
Index ncols
void Select(T &needles, const T &haystack, const ArrayOfIndex &needleindexes, const Verbosity &verbosity)
WORKSPACE METHOD: Select.
void VectorNLogSpace(Vector &out, const Index &nelem, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLogSpace.
void chk_if_bool(const String &x_name, const Index &x)
chk_if_bool
Definition: check_input.cc:65
void chk_met_mm_backend(const Matrix &mmb)
Check met_mm_backend.
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:86
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:111
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Helper comparison class to sort an array or vector based on an ArrayOfNumeric.
Definition: array.h:331
A constant view of a Matrix.
Definition: matpackI.h:1014
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:430
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Vector.
Definition: matpackI.h:489
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:46
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
void resize(const GriddedField1 &gf)
Make this GriddedField1 the same size as the given one.
void checksize_strict() const final
Strict consistency check.
void set_name(const String &s)
Set name of this gridded field.
const ArrayOfString & get_string_grid(Index i) const
Get a string grid.
void set_grid_name(Index i, const String &s)
Set grid name.
void set_grid(Index i, const Vector &g)
Set a numeric grid.
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
The range class.
Definition: matpackI.h:165
The Sparse class.
Definition: matpackII.h:67
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:91
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Numeric ro(Index r, Index c) const
Read only index operator.
Definition: matpackII.cc:130
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1058
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
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 abs(x)
#define q
#define min(a, b)
#define max(a, b)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2235
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
bool is_multiple(const Index &x, const Index &y)
Checks if an integer is a multiple of another integer.
Definition: logic.cc:65
std::vector< T > linspace(T s, T e, typename std::vector< T >::size_type count) noexcept
const Index GFIELD1_F_GRID
void AntennaOff(Index &antenna_dim, Matrix &mblock_dlos_grid, const Verbosity &verbosity)
WORKSPACE METHOD: AntennaOff.
Definition: m_sensor.cc:187
void f_gridFromSensorAMSU(Vector &f_grid, const Vector &lo, const ArrayOfVector &f_backend, const ArrayOfArrayOfGriddedField1 &backend_channel_response, const Numeric &spacing, const Verbosity &verbosity)
WORKSPACE METHOD: f_gridFromSensorAMSU.
Definition: m_sensor.cc:384
void sensor_responseGenericAMSU(Vector &f_grid, Index &antenna_dim, Matrix &mblock_dlos_grid, Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, Index &sensor_norm, const Index &atmosphere_dim, const Index &stokes_dim, const Matrix &sensor_description_amsu, const Numeric &spacing, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseGenericAMSU.
Definition: m_sensor.cc:2736
void yApplySensorPol(Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, const Index &stokes_dim, const Index &jacobian_do, const Matrix &sensor_pos, const Matrix &sensor_pol, const Verbosity &)
WORKSPACE METHOD: yApplySensorPol.
Definition: m_sensor.cc:3525
void sensor_responseMixer(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Numeric &lo, const GriddedField1 &sideband_response, const Index &sensor_norm, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseMixer.
Definition: m_sensor.cc:1866
void sensor_responseStokesRotation(Sparse &sensor_response, const Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Index &stokes_dim, const Vector &stokes_rotation, const Verbosity &)
WORKSPACE METHOD: sensor_responseStokesRotation.
Definition: m_sensor.cc:2633
void sensorOff(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, Matrix &mblock_dlos_grid, const Index &stokes_dim, const Vector &f_grid, const Verbosity &verbosity)
WORKSPACE METHOD: sensorOff.
Definition: m_sensor.cc:1827
void f_gridMetMM(Vector &f_grid, Vector &f_backend, ArrayOfArrayOfIndex &channel2fgrid_indexes, ArrayOfVector &channel2fgrid_weights, const Matrix &mm_back, const Vector &freq_spacing, const ArrayOfIndex &freq_number, const Numeric &freq_merge_threshold, const Verbosity &)
WORKSPACE METHOD: f_gridMetMM.
Definition: m_sensor.cc:735
const Index GFIELD4_ZA_GRID
void AntennaMultiBeamsToPencilBeams(Matrix &sensor_pos, Matrix &sensor_los, Matrix &antenna_dlos, Index &antenna_dim, Matrix &mblock_dlos_grid, const Index &atmosphere_dim, const Verbosity &verbosity)
WORKSPACE METHOD: AntennaMultiBeamsToPencilBeams.
Definition: m_sensor.cc:118
const Index GFIELD4_FIELD_NAMES
void sensor_responseWMRF(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Sparse &wmrf_weights, const Vector &f_backend, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseWMRF.
Definition: m_sensor.cc:3346
void sensor_responseBackendFrequencySwitching(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Vector &f_backend, const ArrayOfGriddedField1 &backend_channel_response, const Index &sensor_norm, const Numeric &df1, const Numeric &df2, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseBackendFrequencySwitching.
Definition: m_sensor.cc:1393
void backend_channel_responseGaussian(ArrayOfGriddedField1 &r, const Vector &fwhm, const Vector &xwidth_si, const Vector &dx_si, const Verbosity &)
WORKSPACE METHOD: backend_channel_responseGaussian.
Definition: m_sensor.cc:344
void sensor_responseBeamSwitching(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Matrix &sensor_response_dlos_grid, const Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Numeric &w1, const Numeric &w2, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseBeamSwitching.
Definition: m_sensor.cc:1464
void AntennaConstantGaussian1D(Index &antenna_dim, Matrix &mblock_dlos_grid, GriddedField4 &r, Matrix &antenna_dlos, const Index &n_za_grid, const Numeric &fwhm, const Numeric &xwidth_si, const Numeric &dx_si, const Verbosity &verbosity)
WORKSPACE METHOD: AntennaConstantGaussian1D.
Definition: m_sensor.cc:72
void antenna_responseVaryingGaussian(GriddedField4 &r, const Numeric &leff, const Numeric &xwidth_si, const Numeric &dx_si, const Index &nf, const Numeric &fstart, const Numeric &fstop, const Index &do_2d, const Verbosity &verbosity)
WORKSPACE METHOD: antenna_responseVaryingGaussian.
Definition: m_sensor.cc:254
void sensor_responseMetMM(Index &antenna_dim, Matrix &mblock_dlos_grid, Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, Index &sensor_norm, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &f_grid, const Vector &f_backend, const ArrayOfArrayOfIndex &channel2fgrid_indexes, const ArrayOfVector &channel2fgrid_weights, const String &iy_unit, const Matrix &antenna_dlos, const ArrayOfString &mm_pol, const Vector &mm_ant, const Index &use_antenna, const Index &mirror_dza, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseMetMM.
Definition: m_sensor.cc:2003
const Numeric RAD2DEG
void mblock_dlos_gridUniformCircular(Matrix &mblock_dlos_grid, const Numeric &spacing, const Numeric &width, const Index &centre, const Verbosity &)
WORKSPACE METHOD: mblock_dlos_gridUniformCircular.
Definition: m_sensor.cc:907
void antenna_responseGaussian(GriddedField4 &r, const Numeric &fwhm, const Numeric &xwidth_si, const Numeric &dx_si, const Index &do_2d, const Verbosity &)
WORKSPACE METHOD: antenna_responseGaussian.
Definition: m_sensor.cc:203
void sensor_responseFillFgrid(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Index &polyorder, const Index &nfill, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseFillFgrid.
Definition: m_sensor.cc:1639
const Numeric NAT_LOG_2
void sensor_responseSimpleAMSU(Vector &f_grid, Index &antenna_dim, Matrix &mblock_dlos_grid, Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, Index &sensor_norm, const Index &atmosphere_dim, const Index &stokes_dim, const Matrix &sensor_description_amsu, const Numeric &spacing, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseSimpleAMSU.
Definition: m_sensor.cc:3095
void sensor_responseInit(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, const Vector &f_grid, const Matrix &mblock_dlos_grid, const Index &antenna_dim, const Index &atmosphere_dim, const Index &stokes_dim, const Index &sensor_norm, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseInit.
Definition: m_sensor.cc:1755
void mblock_dlos_gridUniformRectangular(Matrix &mblock_dlos_grid, const Numeric &spacing, const Numeric &za_width, const Numeric &aa_width, const Index &centre, const Verbosity &)
WORKSPACE METHOD: mblock_dlos_gridUniformRectangular.
Definition: m_sensor.cc:946
const Index GFIELD4_F_GRID
const Numeric PI
void sensor_responseMixerBackendPrecalcWeights(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Vector &f_backend, const ArrayOfArrayOfIndex &channel2fgrid_indexes, const ArrayOfVector &channel2fgrid_weights, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseMixerBackendPrecalcWeights.
Definition: m_sensor.cc:2193
void sensor_responseMultiMixerBackend(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Vector &lo_multi, const ArrayOfGriddedField1 &sideband_response_multi, const ArrayOfString &sideband_mode_multi, const ArrayOfVector &f_backend_multi, const ArrayOfArrayOfGriddedField1 &backend_channel_response_multi, const Index &sensor_norm, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseMultiMixerBackend.
Definition: m_sensor.cc:2333
const Numeric SPEED_OF_LIGHT
void ySimpleSpectrometer(Vector &y, Vector &y_f, const Matrix &iy, const Index &stokes_dim, const Vector &f_grid, const Numeric &df, const Verbosity &verbosity)
WORKSPACE METHOD: ySimpleSpectrometer.
Definition: m_sensor.cc:3447
void f_gridFromSensorHIRS(Vector &f_grid, const Vector &f_backend, const ArrayOfGriddedField1 &backend_channel_response, const Numeric &spacing, const Verbosity &verbosity)
WORKSPACE METHOD: f_gridFromSensorHIRS.
Definition: m_sensor.cc:668
void backend_channel_responseFlat(ArrayOfGriddedField1 &r, const Numeric &resolution, const Verbosity &)
WORKSPACE METHOD: backend_channel_responseFlat.
Definition: m_sensor.cc:323
void sensor_responseIF2RF(Vector &sensor_response_f, Vector &sensor_response_f_grid, const Numeric &lo, const String &sideband_mode, const Verbosity &)
WORKSPACE METHOD: sensor_responseIF2RF.
Definition: m_sensor.cc:1602
const Numeric DEG2RAD
void sensor_responseAntenna(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Matrix &sensor_response_dlos_grid, const Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Index &atmosphere_dim, const Index &antenna_dim, const Matrix &antenna_dlos, const GriddedField4 &antenna_response, const Index &sensor_norm, const String &option_2d, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseAntenna.
Definition: m_sensor.cc:990
void sensor_responseFrequencySwitching(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseFrequencySwitching.
Definition: m_sensor.cc:1527
void WMRFSelectChannels(Vector &f_grid, Sparse &wmrf_weights, Vector &f_backend, const ArrayOfIndex &wmrf_channels, const Verbosity &verbosity)
WORKSPACE METHOD: WMRFSelectChannels.
Definition: m_sensor.cc:3230
void f_gridFromSensorAMSUgeneric(Vector &f_grid, const ArrayOfVector &f_backend_multi, const ArrayOfArrayOfGriddedField1 &backend_channel_response_multi, const Numeric &spacing, const Vector &verbosityVect, const Verbosity &verbosity)
WORKSPACE METHOD: f_gridFromSensorAMSUgeneric.
Definition: m_sensor.cc:516
const Index GFIELD4_AA_GRID
void sensor_responsePolarisation(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const Index &stokes_dim, const String &iy_unit, const ArrayOfIndex &instrument_pol, const Verbosity &)
WORKSPACE METHOD: sensor_responsePolarisation.
Definition: m_sensor.cc:2508
void sensor_responseBackend(Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Vector &f_backend, const ArrayOfGriddedField1 &backend_channel_response, const Index &sensor_norm, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_responseBackend.
Definition: m_sensor.cc:1239
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
void sub(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse subtraction.
Definition: matpackII.cc:667
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:207
#define CREATE_OUT2
Definition: messages.h:206
Index nelem(const Lines &l)
Number of lines.
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
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
Propagation path structure and functions.
#define v
#define w
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:747
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
void antenna2d_gridded_dlos(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_dlos, const GriddedField4 &antenna_response, ConstMatrixView mblock_dlos, ConstVectorView f_grid, const Index n_pol)
antenna2d_interp_gridded_dlos
Definition: sensor.cc:199
void find_effective_channel_boundaries(Vector &fmin, Vector &fmax, const Vector &f_backend, const ArrayOfGriddedField1 &backend_channel_response, const Numeric &delta, const Verbosity &verbosity)
Calculate channel boundaries from instrument response functions.
Definition: sensor.cc:1165
void met_mm_polarisation_hmatrix(Sparse &H, const ArrayOfString &mm_pol, const Numeric dza, const Index stokes_dim, const String &iy_unit)
Calculate polarisation H-matrix.
Definition: sensor.cc:793
void gaussian_response(Vector &y, const Vector &x, const Numeric &x0, const Numeric &fwhm)
gaussian_response
Definition: sensor.cc:647
void antenna2d_interp_response(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_dlos, const GriddedField4 &antenna_response, ConstMatrixView mblock_dlos, ConstVectorView f_grid, const Index n_pol)
antenna2d_interp_response
Definition: sensor.cc:443
void gaussian_response_autogrid(Vector &x, Vector &y, const Numeric &x0, const Numeric &fwhm, const Numeric &xwidth_si, const Numeric &dx_si)
gaussian_response_autogrid
Definition: sensor.cc:624
void mixer_matrix(Sparse &H, Vector &f_mixer, const Numeric &lo, const GriddedField1 &filter, ConstVectorView f_grid, const Index &n_pol, const Index &n_sp, const Index &do_norm)
mixer_matrix
Definition: sensor.cc:663
void spectrometer_matrix(Sparse &H, ConstVectorView ch_f, const ArrayOfGriddedField1 &ch_response, ConstVectorView sensor_f, const Index &n_pol, const Index &n_sp, const Index &do_norm)
spectrometer_matrix
Definition: sensor.cc:991
void mueller_rotation(Sparse &H, const Index &stokes_dim, const Numeric &rotangle)
mueller_rotation
Definition: sensor.cc:763
void antenna1d_matrix(Sparse &H, const Index &antenna_dim, ConstVectorView antenna_dza, const GriddedField4 &antenna_response, ConstVectorView za_grid, ConstVectorView f_grid, const Index n_pol, const Index do_norm)
antenna1d_matrix
Definition: sensor.cc:58
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Definition: sensor.cc:1057
void sensor_aux_vectors(Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, ConstVectorView sensor_response_f_grid, const ArrayOfIndex &sensor_response_pol_grid, ConstMatrixView sensor_response_dlos_grid)
sensor_aux_vectors
Definition: sensor.cc:954
Sensor modelling functions.
Contains sorting routines.
Header file for special_interp.cc.
This file contains basic functions to handle XML data files.