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