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