33#include "matpack_math.h"
64 const Index& atmosphere_dim,
67 const Index nmblock = sensor_pos.nrows();
68 const Index nant = antenna_dlos.nrows();
74 if (sensor_pos.ncols() != atmosphere_dim)
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) {
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());
89 if (antenna_dim == 2 && atmosphere_dim < 3) {
90 throw runtime_error(
"If *antenna_dim* is 2, *atmosphere_dim* must be 3.");
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)
97 "*antenna_dlos* can only have two columns for 3D atmosphers.");
100 const Matrix pos_copy = sensor_pos;
101 const Matrix los_copy = sensor_los;
103 sensor_pos.resize(nmblock * nant, pos_copy.ncols());
104 sensor_los.resize(nmblock * nant, los_copy.ncols());
106 for (Index ib = 0; ib < nmblock; ib++) {
107 for (Index ia = 0; ia < nant; ia++) {
108 const Index i = ib * nant + ia;
110 sensor_pos(i, joker) = pos_copy(ib, joker);
111 sensor_los(i, joker) = los_copy(ib, joker);
113 sensor_los(i, 0) += antenna_dlos(ia, 0);
114 if (antenna_dlos.ncols() == 2) sensor_los(i, 1) += antenna_dlos(ia, 1);
119 AntennaOff(antenna_dim, mblock_dlos, verbosity);
121 antenna_dlos.resize(1, 1);
133 out2 <<
" Sets the antenna dimensionality to 1.\n";
136 out2 <<
" Sets *mblock_dlos* to have one row with value 0.\n";
137 mblock_dlos.resize(1, 1);
145 const Vector& f_points,
147 const Numeric& grid_width,
148 const Index& grid_npoints,
151 const Index nf = f_points.nelem();
155 "*f_points* and *fwhm* must have the same length.");
159 const Numeric gwidth = grid_width > 0 ? grid_width : 2 * max(fwhm);
163 nlinspace(grid, -gwidth/2, gwidth/2, grid_npoints);
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
");
173 r.set_grid_name(3, "Azimuth angle
");
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.
");
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;
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.
");
196 VectorGaussian(y, grid, 0, -1.0, fwhm[i], verbosity);
197 r.data(0, i, joker, 0) = y;
204 /* Workspace method: Doxygen documentation will be auto-generated */
205void antenna_responseGaussianConstant(GriddedField4& r,
207 const Numeric& grid_width,
208 const Index& grid_npoints,
210 const Verbosity& verbosity) {
211 ARTS_USER_ERROR_IF (fwhm <= 0, "*fwhm* must be > 0.
");
213 antenna_responseGaussian(r,
224/* Workspace method: Doxygen documentation will be auto-generated */
225void antenna_responseGaussianEffectiveSize(GriddedField4& r,
227 const Numeric& grid_width,
228 const Index& grid_npoints,
230 const Numeric& fstart,
231 const Numeric& fstop,
233 const Verbosity& verbosity) {
234 ARTS_USER_ERROR_IF (grid_npoints < 2, "*grid_npoints* must be > 1.
");
236 Numeric gwidth = grid_width;
238 gwidth = 2.0 * RAD2DEG * SPEED_OF_LIGHT / (leff * fstart);
243 nlinspace(grid, -gwidth/2, gwidth/2, grid_npoints);
246 r.set_name("Antenna response
");
247 r.set_grid_name(0, "Polarisation
");
248 r.set_grid(0, {"NaN
"});
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
");
255 r.set_grid_name(3, "Azimuth angle
");
259 r.data.resize(1, nf, grid_npoints, grid_npoints);
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;
270 r.set_grid(3, Vector(1, 0));
271 r.data.resize(1, nf, grid_npoints, 1);
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;
283/* Workspace method: Doxygen documentation will be auto-generated */
284void backend_channel_responseFlat(ArrayOfGriddedField1& r,
285 const Numeric& resolution,
288 r[0].set_name("Backend channel response function
");
292 r[0].set_grid_name(0, "Frequency
");
293 x[1] = resolution / 2.0;
298 r[0].data[0] = 1 / resolution;
299 r[0].data[1] = r[0].data[0];
304/* Workspace method: Doxygen documentation will be auto-generated */
305void backend_channel_responseGaussian(ArrayOfGriddedField1& r,
306 const Vector& f_backend,
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.
");
320 for (Index i=0; i<nf; ++i) {
321 ARTS_USER_ERROR_IF (fwhm[i] <= 0, "All values in *fwhm* must be >= 0.
");
323 const Numeric gwidth = grid_width > 0 ? grid_width : 2 * fwhm[i];
325 nlinspace(grid, -gwidth/2, gwidth/2, grid_npoints);
328 VectorGaussian(y, grid, 0, -1.0, fwhm[i], verbosity);
330 r[i].set_name("Backend channel response function
");
331 r[i].set_grid_name(0, "Frequency
");
332 r[i].set_grid(0, grid);
339/* Workspace method: Doxygen documentation will be auto-generated */
340void backend_channel_responseGaussianConstant(ArrayOfGriddedField1& r,
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.
");
347 backend_channel_responseGaussian(r,
357/* Workspace method: Doxygen documentation will be auto-generated */
358void f_gridFromSensorAMSU( // WS Output:
362 const ArrayOfVector& f_backend,
363 const ArrayOfArrayOfGriddedField1& backend_channel_response,
364 // Control Parameters:
365 const Numeric& spacing,
366 const Verbosity& verbosity) {
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.
373 for (Index i = 0; i < f_backend.nelem(); ++i)
374 for (Index s = 0; s < f_backend[i].nelem(); ++s) ++n_chan;
376 // Checks on input quantities:
378 // There must be at least one channel:
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());
386 // Is number of LOs consistent in all input variables?
387 if ((f_backend.nelem() != lo.nelem()) ||
388 (backend_channel_response.nelem() != lo.nelem())) {
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());
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()) {
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());
404 // Start the actual work.
406 // We construct the necessary input for function find_effective_channel_boundaries,
407 // which will identify channel boundaries, taking care of overlaping channels.
409 // A "flat
" vector of nominal band frequencies (two for each AMSU channel):
410 Vector f_backend_flat(2 * n_chan);
412 // A "flat" list of channel response functions (two for each AMSU channel)
413 ArrayOfGriddedField1 backend_channel_response_flat(2 * n_chan);
415 // Counts position inside the flat arrays during construction:
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];
424 f_backend_flat[j] = this_f_backend;
425 backend_channel_response_flat[j] = this_grid;
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;
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);
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;
445 // Call subfunction to do the actual work of merging overlapping
446 // channels and identifying channel boundaries:
447 find_effective_channel_boundaries(fmin,
450 backend_channel_response_flat,
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;
458 for (Index i = 0; i < fmin.nelem(); ++i) {
460 const Numeric bw = fmax[i] - fmin[i];
462 // How many grid intervals do I need?
463 const Numeric npf = ceil(bw / spacing);
465 // How many grid points to store? - Number of grid intervals
467 const Index npi = (Index)npf + 1;
469 // Create the grid for this band:
471 nlinspace(grid, fmin[i], fmax[i], npi);
473 out3 << " Band range
" << i << ":
" << grid << "\n
";
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]);
480 // Copy result to output vector:
481 f_grid = f_grid_array;
483 out2 << " Total number of frequencies in f_grid:
" << f_grid.nelem() << "\n
";
488/* Workspace method: Doxygen documentation will be auto-generated */
489void f_gridFromSensorAMSUgeneric( // WS Output:
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) {
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();
505 numFpoints += backend_channel_response_multi[idx][idy].get_grid_size(0);
509 Index numChan = backend_channel_response_multi.nelem(); // number of channels
511 // Checks on input quantities:
513 // There must be at least one channel:
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());
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.
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);
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;
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];
541 f_backend_flat[i] = this_f_backend;
542 backend_channel_response_nonflat[i] = this_grid;
544 idy < backend_channel_response_multi[i][ii].get_grid_size(0);
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];
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];
563 // We build up a total list of absolute frequency ranges for all passbands
566 fmax; // - these variables will be resized, therefore len(1) is enough for now.,
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;
573 // Call subfunction to do the actual work of merging overlapping
574 // channels and identifying channel boundaries:
575 find_effective_channel_boundaries(fmin,
578 backend_channel_response_nonflat,
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;
586 for (Index i = 0; i < fmin.nelem(); ++i) {
588 const Numeric bw = fmax[i] - fmin[i];
589 Numeric npf = ceil(bw / spacing); // Set a default value
591 // How many grid intervals do I need?
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])) {
601 if (VerbosityValVect[ii] < VerbosityValVect[verbIdx]) {
607 if (spacing > VerbosityValVect[verbIdx]) {
609 bw / VerbosityValVect[verbIdx]); // is the default value to coarse?
611 npf = ceil(bw / spacing); // Default value
615 // How many grid points to store? - Number of grid intervals
617 const Index npi = (Index)npf + 1;
619 // What is the actual grid spacing inside the band?
620 const Numeric gs = bw / npf;
622 // Create the grid for this band:
623 Vector grid=uniform_grid(fmin[i], npi, gs);
625 out3 << " Band range
" << i << ":
" << grid << "\n
";
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]);
632 // Copy result to output vector:
633 f_grid = f_grid_array;
635 out2 << " Total number of frequencies in f_grid:
" << f_grid.nelem() << "\n
";
640/* Workspace method: Doxygen documentation will be auto-generated */
641void f_gridFromSensorHIRS( // WS Output:
644 const Vector& f_backend,
645 const ArrayOfGriddedField1& backend_channel_response,
646 // Control Parameters:
647 const Numeric& spacing,
648 const Verbosity& verbosity) {
655 os << "Expected positive spacing. Found spacing to be:
" << spacing << "\n
";
656 throw runtime_error(os.str());
658 // Call subfunction to get channel boundaries. Also does input
659 // consistency checking for us.
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;
667 find_effective_channel_boundaries(
668 fmin, fmax, f_backend, backend_channel_response, delta, verbosity);
670 // Ok, now we just have to create a frequency grid for each of the
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;
677 for (Index i = 0; i < fmin.nelem(); ++i) {
679 const Numeric bw = fmax[i] - fmin[i];
681 // How many grid intervals do I need?
682 const Numeric npf = ceil(bw / spacing);
684 // How many grid points to store? - Number of grid intervals
686 const Index npi = (Index)npf + 1;
688 // Create the grid for this band:
690 nlinspace(grid, fmin[i], fmax[i], npi);
692 out3 << " Band range
" << i << ":
" << grid << "\n
";
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]);
699 // Copy result to output vector:
700 f_grid = f_grid_array;
702 out2 << " Total number of frequencies in f_grid:
" << f_grid.nelem() << "\n
";
707/* Workspace method: Doxygen documentation will be auto-generated */
712 ArrayOfArrayOfIndex& channel2fgrid_indexes,
713 ArrayOfVector& channel2fgrid_weights,
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,
722 const Index nchannels = mm_back.nrows();
726 chk_met_mm_backend(mm_back);
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*.
");
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*.
");
738 if (freq_merge_threshold <= 0)
739 throw std::runtime_error("The *freq_merge_threshold* must be > 0.\n
");
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.
");
747 ArrayOfNumeric f_grid_unsorted;
748 ArrayOfIndex nf_per_channel(nchannels);
749 ArrayOfIndex index_in_unsorted;
751 f_backend.resize(nchannels);
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);
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];
764 if (this_spacing <= 0)
765 throw std::runtime_error("*freq_spacing must be > 0.
");
767 if (this_fnumber == 0) {
769 os << "*freq_number* must be -1 or greater zero:
"
770 << "freq_number[
" << ch << "] =
" << this_fnumber;
771 std::runtime_error(os.str());
774 // Number of passbands
776 1 + ((Index)(offset1 > 0)) + (2 * (Index)(offset2 > 0));
778 // Number of frequencies per passband
779 Index nfperband = this_fnumber;
781 if (this_fnumber == -1 ||
782 bandwidth / (Numeric)this_fnumber > this_spacing) {
783 nfperband = (Index)ceil(bandwidth / this_spacing);
786 // Fill the data known so far
787 nf_per_channel[ch] = npassb * nfperband;
790 // Distance between each new f_grid value
791 const Numeric df = bandwidth / (Numeric)nfperband;
793 // Loop passbands and determine f_grid values
794 for (Index b = 0; b < npassb; b++) {
795 // Centre frequency of passband
798 fc += (-1 + 2 * (Numeric)b) * offset1;
799 } else if (npassb == 4) {
805 if (b == 0 || b == 2) {
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;
816 // Does this frequency exist or not?
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) {
821 index_in_unsorted.push_back(f_try);
826 f_grid_unsorted.push_back(fnew);
827 index_in_unsorted.push_back(f_grid_unsorted.size() - 1);
833 // Determine sort order for f_grid
834 const size_t nf = f_grid_unsorted.size();
835 ArrayOfIndex move2index(nf);
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;
844 // Sort frequency position vector by frequency
845 std::sort(sorted_indices.begin(),
846 sorted_indices.end(),
847 CmpArrayOfNumeric(f_grid_unsorted));
849 // Create vector with indices in target vector
850 for (size_t i = 0; i < nf; i++) {
851 move2index[sorted_indices[i]] = i;
857 for (size_t f_index = 0; f_index < nf; f_index++) {
858 f_grid[move2index[f_index]] = f_grid_unsorted[f_index];
861 // Create channel2 fgrid variables
862 channel2fgrid_indexes.resize(nchannels);
863 channel2fgrid_weights.resize(nchannels);
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]);
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];
879/* Workspace method: Doxygen documentation will be auto-generated */
880void mblock_dlosFrom1dAntenna(Matrix& mblock_dlos,
881 const GriddedField4& antenna_response,
882 const Index& npoints,
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.
");
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();
892 // Cumulative integral of response (factor /2 skipped, but does not matter)
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);
900 // Equally spaced vector between end points of cumulative sum
902 nlinspace(csp, cumtrapz[0], cumtrapz[nr - 1], npoints);
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);
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) {
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);
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;
944 // Initialise a output stream for runtime errors and a flag for errors
946 bool error_found = false;
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
";
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
";
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
";
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)
973 "*antenna_dlos* can only have two columns
for 3D atmosphers.
");
975 // We allow angles in antenna_los to be unsorted
977 // Checks of antenna_response polarisation dimension
979 const Index lpolgrid =
980 antenna_response.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
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
";
989 // Checks of antenna_response frequency dimension
991 ConstVectorView aresponse_f_grid =
992 antenna_response.get_numeric_grid(GFIELD4_F_GRID);
994 chk_if_increasing("f_grid of antenna_response
", aresponse_f_grid);
996 Numeric f_dlow = 0.0;
997 Numeric f_dhigh = 0.0;
999 f_dlow = min(sensor_response_f_grid) - aresponse_f_grid[0];
1000 f_dhigh = last(aresponse_f_grid) - max(sensor_response_f_grid);
1002 if (aresponse_f_grid.nelem() > 1) {
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
1008 << "the lower end.\n
";
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
1016 << "the upper end.\n
";
1021 // Checks of antenna_response za dimension
1023 ConstVectorView aresponse_za_grid =
1024 antenna_response.get_numeric_grid(GFIELD4_ZA_GRID);
1026 chk_if_increasing("za_grid of *antenna_response*
", aresponse_za_grid);
1028 if (aresponse_za_grid.nelem() < 2) {
1029 os << "The zenith angle grid of *antenna_response* must have >= 2 values.\n
";
1033 // Checks of antenna_response aa dimension
1035 ConstVectorView aresponse_aa_grid =
1036 antenna_response.get_numeric_grid(GFIELD4_AA_GRID);
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
";
1045 chk_if_increasing("aa_grid of antenna_response
", aresponse_aa_grid);
1047 if (aresponse_za_grid.nelem() < 2) {
1048 os << "The zenith angle grid of *antenna_response* must have >= 2\n
"
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
";
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;
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));
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
";
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
";
1090 // Demands differs between the options and checks are done inside
1094 // If errors where found throw runtime_error with the collected error
1096 if (error_found) throw runtime_error(os.str());
1098 // And finally check if grids and data size match
1099 antenna_response.checksize_strict();
1101 // Call the core function
1105 if (antenna_dim == 1)
1106 antenna1d_matrix(hantenna,
1108 antenna_dlos(joker, 0),
1110 sensor_response_dlos_grid(joker, 0),
1111 sensor_response_f_grid,
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,
1123 sensor_response_dlos_grid,
1125 sensor_response_f_grid,
1128 else if (option_2d == "gridded_dlos
" ) {
1129 antenna2d_gridded_dlos(hantenna,
1133 sensor_response_dlos_grid,
1134 sensor_response_f_grid,
1139 throw runtime_error( "Unrecognised choice
for *option_2d*.
" );
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);
1150 // Some extra output.
1151 out3 << " Size of *sensor_response*:
" << sensor_response.nrows() << "x
"
1152 << sensor_response.ncols() << "\n
";
1154 // Update sensor_response_dlos_grid
1155 sensor_response_dlos_grid = antenna_dlos;
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);
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) {
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;
1189 // Initialise an output stream for runtime errors and a flag for errors
1191 bool error_found = false;
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
";
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
";
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
";
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
";
1220 // Check number of columns in backend_channel_response
1222 const Index nrp = backend_channel_response.nelem();
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
";
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());
1234 Numeric f_dlow = 0.0;
1235 Numeric f_dhigh = 0.0;
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);
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
";
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
";
1255 // Check if the relative grid added to the channel frequencies expands
1256 // outside sensor_response_f_grid.
1258 Numeric f1 = f_backend[i] + bchr_f_grid[0] - min(sensor_response_f_grid);
1260 (max(sensor_response_f_grid) - f_backend[i]) - last(bchr_f_grid);
1262 f_dlow = min(f_dlow, f1);
1263 f_dhigh = min(f_dhigh, f2);
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
"
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
"
1281 // If errors where found throw runtime_error with the collected error
1283 if (error_found) throw runtime_error(os.str());
1285 // Call the core function
1289 spectrometer_matrix(hbackend,
1291 backend_channel_response,
1292 sensor_response_f_grid,
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);
1304 // Some extra output.
1305 out3 << " Size of *sensor_response*:
" << sensor_response.nrows() << "x
"
1306 << sensor_response.ncols() << "\n
";
1308 // Update sensor_response_f_grid
1309 sensor_response_f_grid = f_backend;
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);
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,
1336 const Verbosity& verbosity) {
1337 // All needed checks are done in sensor_responseBackend
1339 Sparse H1 = sensor_response, H2 = sensor_response;
1341 // Some needed vectors
1342 Vector f_backend_shifted;
1343 Vector fdummy = sensor_response_f, fdummy_grid = sensor_response_f_grid;
1346 f_backend_shifted = f_backend;
1347 f_backend_shifted += df1;
1349 sensor_responseBackend(H1,
1351 sensor_response_pol,
1352 sensor_response_dlos,
1354 sensor_response_pol_grid,
1355 sensor_response_dlos_grid,
1357 backend_channel_response,
1361 f_backend_shifted = f_backend;
1362 f_backend_shifted += df2;
1364 sensor_responseBackend(H2,
1366 sensor_response_pol,
1367 sensor_response_dlos,
1368 sensor_response_f_grid,
1369 sensor_response_pol_grid,
1370 sensor_response_dlos_grid,
1372 backend_channel_response,
1377 sub(sensor_response, H2, H1);
1379 // sensor_response_f_grid shall be f_backend
1380 sensor_response_f_grid = f_backend;
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);
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,
1403 const Verbosity& verbosity) {
1406 if (sensor_response_dlos_grid.nrows() != 2)
1407 throw runtime_error(
1408 "This method
requires that the number of observation directions is 2.
");
1410 if (sensor_response_pol_grid.nelem() != 1)
1411 throw runtime_error(
1412 "This method handles (so far) only single polarisation cases.
");
1414 const Index n = sensor_response_f_grid.nelem();
1416 // Form H matrix representing beam switching
1417 Sparse Hbswitch(n, 2 * n);
1418 Vector hrow(2 * n, 0.0);
1420 for (Index i = 0; i < n; i++) {
1424 Hbswitch.insert_row(i, hrow);
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);
1436 // Some extra output.
1437 out3 << " Size of *sensor_response*:
" << sensor_response.nrows() << "x
"
1438 << sensor_response.ncols() << "\n
";
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;
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);
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) {
1468 if (sensor_response_dlos_grid.nrows() != 1)
1469 throw runtime_error(
1470 "This method
requires that the number of observation directions is 1.
");
1472 if (sensor_response_pol_grid.nelem() != 1)
1473 throw runtime_error(
1474 "This method handles (so far) only single polarisation cases.
");
1476 const Index n = sensor_response_f_grid.nelem();
1477 const Index n2 = n / 2;
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*.
");
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
"
1491 // Form H matrix representing frequency switching
1492 Sparse Hbswitch(n2, n);
1493 Vector hrow(n, 0.0);
1495 for (Index i = 0; i < n2; i++) {
1499 Hbswitch.insert_row(i, hrow);
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);
1511 // Some extra output.
1512 out3 << " Size of *sensor_response*:
" << sensor_response.nrows() << "x
"
1513 << sensor_response.ncols() << "\n
";
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)];
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);
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,
1537 const String& sideband_mode,
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.
");
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;
1554 else if (sideband_mode == "upper
") {
1555 sensor_response_f += lo;
1556 sensor_response_f_grid += lo;
1561 throw runtime_error(
1562 "Only allowed options
for *sideband _mode* are \
"lower\" and \"upper\".");
1570 Vector& sensor_response_f,
1572 Matrix& sensor_response_dlos,
1573 Vector& sensor_response_f_grid,
1575 const Matrix& sensor_response_dlos_grid,
1576 const Index& polyorder,
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;
1589 bool error_found =
false;
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";
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";
1605 if (polyorder < 2 || polyorder > 7) {
1606 os <<
"Accepted range for *polyorder* is [3,7].\n";
1610 os <<
"The argument *nfill* must be > 1.\n";
1616 if (error_found)
throw runtime_error(os.str());
1620 const Index n1 = nfill + 1;
1621 const Index n2 = nfill + 2;
1622 const Index nnew = (nf - 1) * n1 + 1;
1626 for (Index i = 0; i < nf - 1; i++) {
1628 nlinspace(fp, sensor_response_f_grid[i], sensor_response_f_grid[i + 1], n2);
1629 fnew[Range(i * n1, n2)] = fp;
1632 const auto lag = my_interp::lagrange_interpolation_list<LagrangeInterpolation>(fnew, sensor_response_f_grid, polyorder);
1636 Sparse hpoly(nnew * npol * nlos, nin);
1637 Vector hrow(nin, 0.0);
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];
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;
1662 Sparse htmp = sensor_response;
1663 sensor_response.resize(hpoly.nrows(), htmp.ncols());
1664 mult(sensor_response, hpoly, htmp);
1667 out3 <<
" Size of *sensor_response*: " << sensor_response.nrows() <<
"x"
1668 << sensor_response.ncols() <<
"\n";
1671 sensor_response_f_grid = fnew;
1675 sensor_response_pol,
1676 sensor_response_dlos,
1677 sensor_response_f_grid,
1678 sensor_response_pol_grid,
1679 sensor_response_dlos_grid);
1686 Vector& sensor_response_f,
1688 Matrix& sensor_response_dlos,
1689 Vector& sensor_response_f_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,
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.");
1726 sensor_response_f_grid = f_grid;
1727 sensor_response_dlos_grid = mblock_dlos;
1729 sensor_response_pol_grid.resize(stokes_dim);
1731 for (Index is = 0; is < stokes_dim; is++) {
1732 sensor_response_pol_grid[is] = is + 1;
1737 sensor_response_pol,
1738 sensor_response_dlos,
1739 sensor_response_f_grid,
1740 sensor_response_pol_grid,
1741 sensor_response_dlos_grid);
1745 const Index n = sensor_response_f.nelem();
1747 out2 <<
" Initialising *sensor_reponse* as a identity matrix.\n";
1748 out3 <<
" Size of *sensor_response*: " << n <<
"x" << n <<
"\n";
1750 sensor_response.resize(n, n);
1751 id_mat(sensor_response);
1758 Vector& sensor_response_f,
1760 Matrix& sensor_response_dlos,
1761 Vector& sensor_response_f_grid,
1763 Matrix& sensor_response_dlos_grid,
1764 Matrix& mblock_dlos,
1765 const Index& stokes_dim,
1766 const Vector& f_grid,
1770 AntennaOff(antenna_dim, mblock_dlos, verbosity);
1775 const Index sensor_norm = 1, atmosphere_dim = 1;
1779 sensor_response_pol,
1780 sensor_response_dlos,
1781 sensor_response_f_grid,
1782 sensor_response_pol_grid,
1783 sensor_response_dlos_grid,
1797 Vector& sensor_response_f,
1799 Matrix& sensor_response_dlos,
1800 Vector& sensor_response_f_grid,
1802 const Matrix& sensor_response_dlos_grid,
1805 const Index& sensor_norm,
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;
1816 ConstVectorView sbresponse_f_grid =
1821 bool error_found =
false;
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";
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";
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";
1844 if (sbresponse_f_grid.nelem() != sideband_response.
data.nelem()) {
1845 os <<
"Mismatch in size of grid and data in *sideband_response*.\n";
1848 if (sbresponse_f_grid.nelem() < 2) {
1849 os <<
"At least two data points must be specified in "
1850 <<
"*sideband_response*.\n";
1853 if (!is_increasing(sbresponse_f_grid)) {
1854 os <<
"The frequency grid of *sideband_response* must be strictly\n"
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";
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"
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?";
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?";
1891 if (error_found)
throw runtime_error(os.str());
1902 sensor_response_f_grid,
1910 Sparse htmp = sensor_response;
1911 sensor_response.resize(hmixer.nrows(), htmp.ncols());
1912 mult(sensor_response, hmixer, htmp);
1915 out3 <<
" Size of *sensor_response*: " << sensor_response.nrows() <<
"x"
1916 << sensor_response.ncols() <<
"\n";
1919 sensor_response_f_grid = f_mixer;
1923 sensor_response_pol,
1924 sensor_response_dlos,
1925 sensor_response_f_grid,
1926 sensor_response_pol_grid,
1927 sensor_response_dlos_grid);
1936 Matrix& mblock_dlos,
1937 Sparse& sensor_response,
1938 Vector& sensor_response_f,
1940 Matrix& sensor_response_dlos,
1941 Vector& sensor_response_f_grid,
1943 Matrix& sensor_response_dlos_grid,
1946 const Index& atmosphere_dim,
1947 const Index& stokes_dim,
1948 const Vector& f_grid,
1949 const Vector& f_backend,
1951 const ArrayOfVector& channel2fgrid_weights,
1953 const Matrix& antenna_dlos,
1955 const Vector& mm_ant,
1957 const Index& use_antenna,
1958 const Index& mirror_dza,
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");
1970 if (!(atmosphere_dim == 1 || atmosphere_dim == 3))
1971 throw std::runtime_error(
1972 "This method only supports 1D and 3D atmospheres.");
1974 if (antenna_dlos.empty())
1975 throw std::runtime_error(
"*antenna_dlos* is empty.");
1977 if (antenna_dlos.ncols() > 2)
1978 throw std::runtime_error(
1979 "The maximum number of columns in *antenna_dlos*"
1983 Matrix antenna_dlos_local;
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.");
1994 const Index n = antenna_dlos.nrows();
1996 for (Index i = 0; i < n; i++) {
1997 if (antenna_dlos(i, 0) != 0) {
2001 antenna_dlos_local.resize(n + nnew, 1);
2002 antenna_dlos_local(Range(0, n), 0) = antenna_dlos(joker, 0);
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);
2011 antenna_dlos_local = antenna_dlos;
2019 Sparse sensor_response_single;
2020 Matrix mblock_dlos_dummy(1, 1);
2021 mblock_dlos_dummy(0, 0) = 0;
2024 sensor_response_pol,
2025 sensor_response_dlos,
2026 sensor_response_f_grid,
2027 sensor_response_pol_grid,
2028 sensor_response_dlos_grid,
2038 sensor_response_pol,
2039 sensor_response_dlos,
2040 sensor_response_f_grid,
2041 sensor_response_pol_grid,
2042 sensor_response_dlos_grid,
2044 channel2fgrid_indexes,
2045 channel2fgrid_weights,
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());
2054 sensor_response_pol_grid.resize(1);
2055 sensor_response_pol_grid[0] = 1;
2057 if (stokes_dim > 1) {
2059 if (mm_pol.
nelem() != nchannels) {
2061 os <<
"Length of *met_mm_polarisation* (" << mm_pol.
nelem()
2063 <<
"number of channels in *met_mm_backend* (" << nchannels <<
").";
2064 throw std::runtime_error(os.str());
2068 Sparse sensor_response_tmp;
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);
2080 sensor_response.rw(iza * nchannels + r,
2081 iza * num_f * stokes_dim +
c) =
v;
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);
2092 sensor_response.rw(iza * nchannels + r,
2093 iza * num_f * stokes_dim +
c) =
v;
2102 throw std::runtime_error(
"The antenna hasn't arrived yet.");
2106 mblock_dlos = antenna_dlos_local;
2109 sensor_response_dlos_grid = mblock_dlos;
2113 sensor_response_pol,
2114 sensor_response_dlos,
2115 sensor_response_f_grid,
2116 sensor_response_pol_grid,
2117 sensor_response_dlos_grid);
2124 Sparse& sensor_response,
2125 Vector& sensor_response_f,
2127 Matrix& sensor_response_dlos,
2128 Vector& sensor_response_f_grid,
2130 const Matrix& sensor_response_dlos_grid,
2131 const Vector& f_backend,
2133 const ArrayOfVector& channel2fgrid_weights,
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;
2147 bool error_found =
false;
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";
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";
2164 os <<
"*f_backend* is empty !!!\n";
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";
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";
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";
2186 if (channel2fgrid_weights.nelem() != channel2fgrid_indexes.
nelem()) {
2187 os <<
"Leading sizes of *channel2fgrid_indexes* and "
2188 <<
"*channel2fgrid_weights* differ.\n";
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";
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";
2211 if (error_found)
throw runtime_error(os.str());
2215 Sparse hmb(nout, nin);
2218 for (Index ifr = 0; ifr < nout_f; ifr++) {
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];
2228 for (Index sp = 0; sp < nlos; sp++) {
2229 for (Index pol = 0; pol < npol; pol++) {
2231 Vector weights_long(nin, 0.0);
2232 weights_long[Range(sp * nin_f * npol + pol, nin_f, npol)] = w1;
2235 hmb.insert_row(sp * nout_f * npol + ifr * npol + pol, weights_long);
2244 Sparse htmp = sensor_response;
2245 sensor_response.resize(hmb.nrows(), htmp.ncols());
2246 mult(sensor_response, hmb, htmp);
2249 sensor_response_f_grid = f_backend;
2253 sensor_response_pol,
2254 sensor_response_dlos,
2255 sensor_response_f_grid,
2256 sensor_response_pol_grid,
2257 sensor_response_dlos_grid);
2264 Sparse& sensor_response,
2265 Vector& sensor_response_f,
2267 Matrix& sensor_response_dlos,
2268 Vector& sensor_response_f_grid,
2270 const Matrix& sensor_response_dlos_grid,
2271 const Vector& lo_multi,
2274 const ArrayOfVector& f_backend_multi,
2276 const Index& sensor_norm,
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();
2287 bool error_found =
false;
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";
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";
2304 if (sideband_response_multi.
nelem() != nlo) {
2305 os <<
"Inconsistency in length between *lo_mixer* and "
2306 <<
"*sideband_response_multi*.\n";
2309 if (sideband_mode_multi.
nelem() != nlo) {
2310 os <<
"Inconsistency in length between *lo_mixer* and "
2311 <<
"*sideband_mode_multi*.\n";
2314 if (f_backend_multi.nelem() != nlo) {
2315 os <<
"Inconsistency in length between *lo_mixer* and "
2316 <<
"*f_backend_multi*.\n";
2319 if (backend_channel_response_multi.
nelem() != nlo) {
2320 os <<
"Inconsistency in length between *lo_mixer* and "
2321 <<
"*backend_channel_response_multi*.\n";
2327 if (error_found)
throw runtime_error(os.str());
2331 ArrayOfVector srfgrid;
2334 for (Index ilo = 0; ilo < nlo; ilo++) {
2337 Sparse sr1 = sensor_response;
2338 Vector srf1 = sensor_response_f;
2340 Matrix srdlos1 = sensor_response_dlos;
2341 Vector srfgrid1 = sensor_response_f_grid;
2350 sensor_response_pol_grid,
2351 sensor_response_dlos_grid,
2353 sideband_response_multi[ilo],
2358 srf1, srfgrid1, lo_multi[ilo], sideband_mode_multi[ilo], verbosity);
2365 sensor_response_pol_grid,
2366 sensor_response_dlos_grid,
2367 f_backend_multi[ilo],
2368 backend_channel_response_multi[ilo],
2371 }
catch (
const runtime_error& e) {
2373 os2 <<
"Error when dealing with receiver/mixer chain (1-based index) "
2376 throw runtime_error(os2.str());
2381 srfgrid.push_back(srfgrid1);
2383 cumsumf[ilo + 1] = cumsumf[ilo] + srfgrid1.
nelem();
2388 const Index nfnew = cumsumf[nlo];
2389 sensor_response_f_grid.resize(nfnew);
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];
2399 const Index ncols = sr[0].ncols();
2400 const Index npolnew = sensor_response_pol_grid.
nelem();
2401 const Index nfpolnew = nfnew * npolnew;
2403 sensor_response.resize(nlos * nfpolnew, ncols);
2405 Vector dummy(ncols, 0.0);
2407 for (Index ilo = 0; ilo < nlo; ilo++) {
2408 const Index nfpolthis = (cumsumf[ilo + 1] - cumsumf[ilo]) * npolnew;
2413 for (Index ilos = 0; ilos < nlos; ilos++) {
2414 for (Index i = 0; i < nfpolthis; i++) {
2416 for (Index ic = 0; ic < ncols; ic++) {
2417 dummy[ic] = sr[ilo](ilos * nfpolthis + i, ic);
2420 sensor_response.insert_row(ilos * nfpolnew + cumsumf[ilo] * npolnew + i,
2428 sensor_response_pol,
2429 sensor_response_dlos,
2430 sensor_response_f_grid,
2431 sensor_response_pol_grid,
2432 sensor_response_dlos_grid);
2439 Vector& sensor_response_f,
2441 Matrix& sensor_response_dlos,
2443 const Vector& sensor_response_f_grid,
2444 const Matrix& sensor_response_dlos_grid,
2445 const Index& stokes_dim,
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();
2457 bool error_found =
false;
2459 Index nfz = nf * nlos;
2460 Index nin = nfz * npol;
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";
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";
2475 if (npol != stokes_dim) {
2476 os <<
"Number of input polarisation does not match *stokes_dim*.\n";
2480 os <<
"The WSV *instrument_pol* can not be empty.\n";
2485 if (error_found)
throw runtime_error(os.str());
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.";
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";
2504 if (error_found)
throw runtime_error(os.str());
2508 if (error_found)
throw runtime_error(os.str());
2512 if (iy_unit ==
"PlanckBT" || iy_unit ==
"RJBT") {
2518 Sparse Hpol(nfz * nnew, nin);
2519 Vector hrow(nin, 0);
2522 for (Index i = 0; i < nfz; i++) {
2523 Index col = i * npol;
2524 for (Index in = 0; in < nnew; in++) {
2532 hrow[Range(col, stokes_dim)], stokes_dim, instrument_pol[in],
w);
2534 Hpol.insert_row(row, hrow);
2544 Sparse Htmp = sensor_response;
2545 sensor_response.resize(Hpol.nrows(), Htmp.ncols());
2546 mult(sensor_response, Hpol, Htmp);
2549 sensor_response_pol_grid = instrument_pol;
2553 sensor_response_pol,
2554 sensor_response_dlos,
2555 sensor_response_f_grid,
2556 sensor_response_pol_grid,
2557 sensor_response_dlos_grid);
2564 const Vector& sensor_response_f_grid,
2566 const Matrix& sensor_response_dlos_grid,
2567 const Index& stokes_dim,
2568 const Vector& stokes_rotation,
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;
2582 bool error_found =
false;
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";
2593 if (stokes_dim < 3) {
2594 os <<
"To perform a rotation of the Stokes coordinate system,\n"
2595 <<
"*stokes_dim* must be >= 3.\n";
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";
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";
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";
2621 if (error_found)
throw runtime_error(os.str());
2626 Sparse H(sensor_response.nrows(), sensor_response.ncols());
2628 Sparse Hrot(npol, npol);
2629 Vector row(H.ncols(), 0);
2632 for (Index ilos = 0; ilos < nlos; ilos++) {
2636 for (Index ifr = 0; ifr < nf; ifr++) {
2637 for (Index ip = 0; ip < npol; ip++) {
2640 for (Index is = 0; is < npol; is++) {
2641 row[irow + is] = Hrot.ro(ip, is);
2643 H.insert_row(irow + ip, row);
2645 for (Index is = 0; is < npol; is++) {
2658 Sparse Htmp = sensor_response;
2659 sensor_response.resize(Htmp.nrows(), Htmp.ncols());
2660 mult(sensor_response, H, Htmp);
2669 Matrix& mblock_dlos,
2670 Sparse& sensor_response,
2671 Vector& sensor_response_f,
2673 Matrix& sensor_response_dlos,
2674 Vector& sensor_response_f_grid,
2676 Matrix& sensor_response_dlos_grid,
2679 const Index& atmosphere_dim,
2680 const Index& stokes_dim,
2681 const Matrix& sensor_description_amsu,
2683 const Numeric& spacing,
2686 const Index n = sensor_description_amsu.nrows();
2687 const Index m = sensor_description_amsu.ncols();
2697 if (5 > sensor_description_amsu.ncols()) {
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());
2704 ConstVectorView lo_multi = sensor_description_amsu(Range(joker), 0);
2705 ConstMatrixView offset = sensor_description_amsu(
2706 Range(joker), Range(1, m - 2));
2707 ConstVectorView verbosityVectIn =
2708 sensor_description_amsu(Range(joker), m - 1);
2709 ConstVectorView width = sensor_description_amsu(Range(joker), m - 2);
2713 const Numeric minRatioVerbosityVsFdiff =
2716 Vector verbosityVect(n);
2718 for (Index idx = 0; idx < n; ++idx) {
2719 if ((verbosityVectIn[idx] == 0) || (verbosityVectIn[idx] > width[idx])) {
2720 verbosityVect[idx] = ((Numeric)width[idx]) / 3;
2722 verbosityVect[idx] = verbosityVectIn[idx];
2730 for (Index i = 0; i < n; ++i) {
2732 for (Index j = 0; j < (m - 2); ++j) {
2734 if (offset(i, j) > 0) {
2739 numPB[i] = 1 << (int)numPB[i];
2740 if (numPB[i] == 1) {
2743 numPBpseudo[i] = numPB[i];
2750 os <<
"This function does currently not support more than 4 passbands per channel"
2752 throw runtime_error(os.str());
2759 ArrayOfVector f_backend_multi(n);
2760 for (Index i = 0; i < n; ++i) {
2762 Vector& f = f_backend_multi[i];
2764 f[0] = lo_multi[i] + 0.0 * width[i];
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;
2778 for (Index pbOffsetIdx = 0; pbOffsetIdx < numPBpseudo[i];
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;
2786 b_resp.
data[pbOffsetIdx * numVal + 0] = 0.0 / (Numeric)numPB[i];
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];
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;
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;
2805 if (pbOffsetIdx == 1) {
2806 pbOffset = 0.0 * width[i];
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;
2820 if (pbOffsetIdx == 0) {
2821 pbOffset = -offset(i, 0);
2823 if (pbOffsetIdx == 1) {
2824 pbOffset = offset(i, 0);
2827 if (numPB[i] == 4) {
2828 if (pbOffsetIdx == 0) {
2829 pbOffset = -offset(i, 0) - offset(i, 1);
2831 if (pbOffsetIdx == 1) {
2832 pbOffset = -offset(i, 0) + offset(i, 1);
2834 if (pbOffsetIdx == 2) {
2835 pbOffset = offset(i, 0) - offset(i, 1);
2837 if (pbOffsetIdx == 3) {
2838 pbOffset = offset(i, 0) + offset(i, 1);
2841 for (Index iii = 0; iii < numVal; ++iii) {
2842 f_range[pbOffsetIdx * numVal + iii] += 1 * pbOffset;
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]))
2851 if ((f_range[ii + 2] - f_range[ii - 1]) >
2854 f_range[ii + 1] = f_range[ii + 2] - verbosityVectIn[i] / 2;
2855 f_range[ii] = f_range[ii - 1] + verbosityVectIn[i] / 2;
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;
2862 f_range[ii - 1] + verbosityVectIn[i] / minRatioVerbosityVsFdiff;
2864 f_range[ii + 1] + verbosityVectIn[i] / minRatioVerbosityVsFdiff;
2874 for (Index i = 0; i < n; ++i) {
2876 r.
set_name(
"Sideband response function");
2877 r.
resize(numPBpseudo[i]);
2878 Vector f(numPBpseudo[i]);
2879 if (numPB[i] == 1) {
2881 f[0] = -.0 * width[i];
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];
2896 f[1] = -offset(i, 0) + offset(i, 1) - 0.5 * width[i];
2898 f[2] = +offset(i, 0) - offset(i, 1) + 0.5 * width[i];
2900 f[3] = +offset(i, 0) + offset(i, 1) + 0.5 * width[i];
2913 backend_channel_response_multi,
2919 AntennaOff(antenna_dim, mblock_dlos, verbosity);
2925 sensor_response_pol,
2926 sensor_response_dlos,
2927 sensor_response_f_grid,
2928 sensor_response_pol_grid,
2929 sensor_response_dlos_grid,
2939 Index numLO = lo_multi.nelem();
2943 ArrayOfVector srfgrid;
2945 const Index nlos = sensor_response_dlos_grid.nrows();
2948 for (Index idxLO = 0; idxLO < numLO; idxLO++) {
2949 Sparse sr1 = sensor_response;
2950 Vector srf1 = sensor_response_f;
2952 Matrix srdlos1 = sensor_response_dlos;
2953 Vector srfgrid1 = sensor_response_f_grid;
2962 sensor_response_pol_grid,
2963 sensor_response_dlos_grid,
2964 f_backend_multi[idxLO],
2965 backend_channel_response_multi[idxLO],
2971 srfgrid.push_back(srfgrid1);
2973 cumsumf[idxLO + 1] = cumsumf[idxLO] + srfgrid1.
nelem();
2978 const Index nfnew = cumsumf[numLO];
2979 sensor_response_f_grid.resize(nfnew);
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];
2987 const Index ncols = sr[0].ncols();
2988 const Index npolnew = sensor_response_pol_grid.
nelem();
2989 const Index nfpolnew = nfnew * npolnew;
2991 sensor_response.resize(nlos * nfpolnew, ncols);
2993 Vector dummy(ncols, 0.0);
2995 for (Index ilo = 0; ilo < numLO; ilo++) {
2996 const Index nfpolthis = (cumsumf[ilo + 1] - cumsumf[ilo]) * npolnew;
3001 for (Index ilos = 0; ilos < nlos; ilos++) {
3002 for (Index i = 0; i < nfpolthis; i++) {
3004 for (Index ic = 0; ic < ncols; ic++) {
3005 dummy[ic] = sr[ilo](ilos * nfpolthis + i, ic);
3008 sensor_response.insert_row(ilos * nfpolnew + cumsumf[ilo] * npolnew + i,
3015 sensor_response_pol,
3016 sensor_response_dlos,
3017 sensor_response_f_grid,
3018 sensor_response_pol_grid,
3019 sensor_response_dlos_grid);
3028 Matrix& mblock_dlos,
3029 Sparse& sensor_response,
3030 Vector& sensor_response_f,
3032 Matrix& sensor_response_dlos,
3033 Vector& sensor_response_f_grid,
3035 Matrix& sensor_response_dlos_grid,
3038 const Index& atmosphere_dim,
3039 const Index& stokes_dim,
3040 const Matrix& sensor_description_amsu,
3042 const Numeric& spacing,
3045 if (3 != sensor_description_amsu.ncols()) {
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());
3053 const Index n = sensor_description_amsu.nrows();
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);
3063 ArrayOfVector f_backend_multi(n);
3064 for (Index i = 0; i < n; ++i) {
3065 Vector& f = f_backend_multi[i];
3067 f[0] = lo_multi[i] + offset[i];
3072 for (Index i = 0; i < n; ++i) {
3073 backend_channel_response_multi[i].resize(1);
3075 r.
set_name(
"Backend channel response function");
3080 f[0] = -0.5 * width[i];
3081 f[1] = +0.5 * width[i];
3092 for (Index i = 0; i < n; ++i) {
3094 r.
set_name(
"Sideband response function");
3099 f[0] = -(offset[i] + 0.5 * width[i]);
3100 f[1] = +(offset[i] + 0.5 * width[i]);
3121 backend_channel_response_multi,
3125 AntennaOff(antenna_dim, mblock_dlos, verbosity);
3129 sensor_response_pol,
3130 sensor_response_dlos,
3131 sensor_response_f_grid,
3132 sensor_response_pol_grid,
3133 sensor_response_dlos_grid,
3144 sensor_response_pol,
3145 sensor_response_dlos,
3146 sensor_response_f_grid,
3147 sensor_response_pol_grid,
3148 sensor_response_dlos_grid,
3150 sideband_response_multi,
3151 sideband_mode_multi,
3153 backend_channel_response_multi,
3162 Sparse& wmrf_weights,
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());
3193 if (
min(wmrf_channels) < 0) {
3194 os <<
"Min(wmrf_channels) must be >= 0, but it is " <<
min(wmrf_channels)
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) <<
".";
3203 if (wmrf_channels.
nelem() == f_backend.nelem()) {
3205 out2 <<
" Retaining all channels.\n";
3207 out2 <<
" Reducing number of channels from " << f_backend.nelem() <<
" to "
3208 << wmrf_channels.
nelem() <<
".\n";
3214 Select(f_backend, f_backend, wmrf_channels, verbosity);
3219 Select(wmrf_weights, wmrf_weights, wmrf_channels, verbosity);
3230 selection.reserve(f_grid.nelem());
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) {
3238 for (i = 0; i < wmrf_weights.nrows(); ++i) {
3239 if (wmrf_weights(i, fi) != 0) {
3240 selection.push_back(fi);
3244 if (i == wmrf_weights.nrows()) {
3245 out3 <<
" The frequency with index " << fi
3246 <<
" is not used by any channel.\n";
3250 if (selection.
nelem() == f_grid.nelem()) {
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");
3256 out2 <<
" Reducing number of frequency grid points from " << f_grid.nelem()
3257 <<
" to " << selection.
nelem() <<
".\n";
3261 Select(f_grid, f_grid, selection, verbosity);
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);
3277 Sparse& sensor_response,
3278 Vector& sensor_response_f,
3280 Matrix& sensor_response_dlos,
3281 Vector& sensor_response_f_grid,
3284 const Matrix& sensor_response_dlos_grid,
3285 const Sparse& wmrf_weights,
3286 const Vector& f_backend,
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;
3298 bool error_found =
false;
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";
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";
3314 os <<
"One of f_grid, pol_grid, dlos_grid are empty. Sizes are: (" << nf
3315 <<
", " << npol <<
", " << nlos <<
")"
3322 const Index nrw = wmrf_weights.nrows();
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";
3334 const Index ncw = wmrf_weights.ncols();
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()
3347 if (error_found)
throw runtime_error(os.str());
3354 Sparse htmp = sensor_response;
3355 sensor_response.resize(wmrf_weights.nrows(), htmp.ncols());
3356 mult(sensor_response, wmrf_weights, htmp);
3359 out3 <<
" Size of *sensor_response*: " << sensor_response.nrows() <<
"x"
3360 << sensor_response.ncols() <<
"\n";
3363 sensor_response_f_grid = f_backend;
3367 sensor_response_pol,
3368 sensor_response_dlos,
3369 sensor_response_f_grid,
3370 sensor_response_pol_grid,
3371 sensor_response_dlos_grid);
3380 const Index& stokes_dim,
3381 const Vector& f_grid,
3385 const Index sensor_norm = 1, atmosphere_dim = 1;
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;
3395 AntennaOff(antenna_dim, mblock_dlos, verbosity);
3399 sensor_response_pol,
3400 sensor_response_dlos,
3401 sensor_response_f_grid,
3402 sensor_response_pol_grid,
3403 sensor_response_dlos_grid,
3414 linspace(f_backend, f_grid[0] + df / 2,
last(f_grid), df);
3423 sensor_response_pol,
3424 sensor_response_dlos,
3425 sensor_response_f_grid,
3426 sensor_response_pol_grid,
3427 sensor_response_dlos_grid,
3434 const Index nf = f_grid.nelem();
3435 const Index n = sensor_response.nrows();
3439 Vector iyb(nf * stokes_dim);
3441 for (Index is = 0; is < stokes_dim; is++) {
3442 iyb[Range(is, nf, stokes_dim)] = iy(joker, is);
3447 y_f = sensor_response_f;
3449 mult(y, sensor_response, iyb);
3460 ArrayOfVector& y_aux,
3463 const Index& stokes_dim,
3464 const Index& jacobian_do,
3465 const Matrix& sensor_pos,
3466 const Matrix& sensor_pol,
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;
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.");
3487 if (jacobian.nrows() != n1)
3488 throw runtime_error(
"Sizes of *y* and *jacobian* are inconsistent.");
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.");
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*.");
3515 const Vector y1 = y, y_f1 = y_f;
3516 const Matrix y_pos1 = y_pos, y_los1 = y_los, y_geo1 = y_geo;
3518 const ArrayOfVector y_aux1 = y_aux;
3519 Matrix jacobian1(0, 0);
3521 jacobian1 = jacobian;
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);
3533 jacobian.resize(n2, jacobian1.ncols());
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;
3541 const Numeric wq = cos(2 *
DEG2RAD * sensor_pol(r,
c));
3542 const Numeric wu = sin(2 *
DEG2RAD * sensor_pol(r,
c));
3545 y[iout] = y1[iin] + wq * y1[iin + 1] + wu * y1[iin + 2];
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);
3555 y_pol[iout] = (Index)sensor_pol(r,
c);
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];
base max(const Array< base > &x)
Max function.
base min(const Array< base > &x)
Min function.
The global header file for ARTS.
Constants of physical expressions as constexpr.
Index nelem() const ARTS_NOEXCEPT
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,...)
#define ARTS_USER_ERROR_IF(condition,...)
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 &)
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.
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.
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.
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.
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.
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.
constexpr Numeric SPEED_OF_LIGHT
constexpr Numeric DEG2RAD
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.
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.
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.
constexpr Numeric NAT_LOG_2
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.
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.
constexpr Numeric RAD2DEG
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.
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.
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.
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.
void backend_channel_responseFlat(ArrayOfGriddedField1 &r, const Numeric &resolution, const Verbosity &)
WORKSPACE METHOD: backend_channel_responseFlat.
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.
void AntennaOff(Index &antenna_dim, Matrix &mblock_dlos, const Verbosity &verbosity)
WORKSPACE METHOD: AntennaOff.
void WMRFSelectChannels(Vector &f_grid, Sparse &wmrf_weights, Vector &f_backend, const ArrayOfIndex &wmrf_channels, const Verbosity &verbosity)
WORKSPACE METHOD: WMRFSelectChannels.
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.
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.
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.
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.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
void flat(VectorView x, ConstMatrixView X)
flat
Numeric last(ConstVectorView x)
last
Declarations having to do with the four output streams.
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
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.
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
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
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
Sensor modelling functions.
Contains sorting routines.
Header file for special_interp.cc.
This file contains basic functions to handle XML data files.