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