26#include "matpack_data.h"
27#include "matpack_math.h"
49 const Index& antenna_dim,
51 const Index& antenna_dim
_U_,
53 ConstVectorView antenna_dza,
55 ConstVectorView za_grid,
56 ConstVectorView f_grid,
58 const Index do_norm) {
60 const Index n_za = za_grid.nelem();
61 const Index n_f = f_grid.nelem();
64 const Index n_ant = antenna_dza.nelem();
73 const Index n_ar_pol =
75 ConstVectorView aresponse_f_grid =
77 ConstVectorView aresponse_za_grid =
83 const Index n_ar_f = aresponse_f_grid.nelem();
84 const Index n_ar_za = aresponse_za_grid.nelem();
85 const Index pol_step = n_ar_pol > 1;
97 const Index nfpol = n_f * n_pol;
100 H.resize(n_ant * nfpol, n_za * nfpol);
103 Vector hrow(H.ncols(), 0.0);
104 Vector hza(n_za, 0.0);
107 Vector aresponse(n_ar_za, 0.0);
110 for (Index ia = 0; ia < n_ant; ia++) {
111 Vector shifted_aresponse_za_grid{aresponse_za_grid};
112 shifted_aresponse_za_grid += antenna_dza[ia];
118 for (Index f = 0; f < n_f; f++) {
120 for (Index ip = 0; ip < n_pol; ip++) {
125 Index new_antenna = 1;
131 aresponse = antenna_response.
data(ip, 0, joker, 0);
132 }
else if (f == 0 && ip == 0)
134 aresponse = antenna_response.
data(0, 0, joker, 0);
140 if (ip == 0 || pol_step)
144 gridpos(gp_f, aresponse_f_grid, Vector(1, f_grid[f]));
145 gridpos(gp_za, aresponse_za_grid, aresponse_za_grid);
146 Tensor3 itw(1, n_ar_za, 4);
148 Matrix aresponse_matrix(1, n_ar_za);
151 antenna_response.
data(ip, joker, joker, 0),
154 aresponse = aresponse_matrix(0, joker);
164 hza, aresponse, shifted_aresponse_za_grid, za_grid);
174 const Index ii = f * n_pol + ip;
176 hrow[Range(ii, n_za, nfpol)] = hza;
178 H.insert_row(ia * nfpol + ii, hrow);
190 const Index& antenna_dim,
192 const Index& antenna_dim
_U_,
194 ConstMatrixView antenna_dlos,
196 ConstMatrixView mblock_dlos,
197 ConstVectorView f_grid,
201 const Index n_dlos = mblock_dlos.nrows();
202 const Index n_f = f_grid.nelem();
206 "For the gridded_dlos option, *mblock_dlos* "
207 "must have two columns.");
211 for(Index i=0; i<n_dlos-1 && mblock_dlos(i+1,0) > mblock_dlos(i,0); i++ ) {
215 "For the gridded_dlos option, the number of za angles "
216 "(among dlos directions) must be >= 2.");
218 "For the gridded_dlos option, the number of dlos angles "
219 "must be a product of two integers.");
220 const Index naa = n_dlos / nza;
221 const Vector za_grid{mblock_dlos(Range(0,nza),0)};
222 const Vector aa_grid{mblock_dlos(Range(0,naa,nza),1)};
223 for(Index i=0; i<n_dlos; i++ ) {
225 "Zenith angle of dlos ", i,
" (0-based) differs to zenith "
226 "angle of dlos ", i-nza,
", while they need to be equal "
227 "to form rectangular grid.")
229 "Azimuth angle of dlos ", i,
" (0-based) differs to azimuth "
230 "angle ", (i/nza)*nza ,
", while they need to be equal "
231 "to form rectangular grid.")
235 const Index n_ant = antenna_dlos.nrows();
243 const Index n_ar_pol =
245 ConstVectorView aresponse_f_grid =
247 ConstVectorView aresponse_za_grid =
249 ConstVectorView aresponse_aa_grid =
252 const Index n_ar_f = aresponse_f_grid.nelem();
253 const Index n_ar_za = aresponse_za_grid.nelem();
254 const Index n_ar_aa = aresponse_aa_grid.nelem();
255 const Index pol_step = n_ar_pol > 1;
269 Tensor4 aresponse_with_cos(n_ar_pol, n_ar_f, n_ar_za, n_ar_aa);
270 for(Index i3=0; i3<n_ar_za; i3++) {
271 const Numeric t = cos(
DEG2RAD * aresponse_za_grid[i3]);
272 for(Index i4=0; i4<n_ar_aa; i4++) {
273 for(Index i2=0; i2<n_ar_f; i2++) {
274 for(Index i1=0; i1<n_ar_pol; i1++) {
275 aresponse_with_cos(i1,i2,i3,i4) = t * antenna_response.
data(i1,i2,i3,i4);
282 const Index nfpol = n_f * n_pol;
285 H.resize(n_ant * nfpol, n_dlos * nfpol);
288 Vector hrow(H.ncols(), 0.0);
289 Vector hza(n_dlos, 0.0);
292 Matrix aresponse(n_ar_za, n_ar_aa, 0.0);
295 for (Index ia = 0; ia < n_ant; ia++) {
300 for (Index f = 0; f < n_f; f++) {
302 for (Index ip = 0; ip < n_pol; ip++) {
307 Index new_antenna = 1;
313 aresponse = aresponse_with_cos(ip, 0, joker, joker);
314 }
else if (f == 0 && ip == 0)
316 aresponse = aresponse_with_cos(0, 0, joker, joker);
322 if (ip == 0 || pol_step) {
325 gridpos(gp_f, aresponse_f_grid, Vector(1, f_grid[f]));
326 gridpos(gp_za, aresponse_za_grid, aresponse_za_grid);
327 gridpos(gp_aa, aresponse_aa_grid, aresponse_aa_grid);
328 Tensor4 itw(1, n_ar_za, n_ar_aa, 8);
330 Tensor3 aresponse_matrix(1, n_ar_za, n_ar_aa);
333 aresponse_with_cos(ip, joker, joker, joker),
337 aresponse = aresponse_matrix(0, joker, joker);
348 Vector zas{aresponse_za_grid};
349 zas += antenna_dlos(ia, 0);
351 "The zenith angle grid in *mblock_dlos* is too narrow. "
352 "It must be extended downwards with at least ",
353 za_grid[0]-zas[0],
" deg.")
355 "The zenith angle grid in *mblock_dlos* is too narrow. "
356 "It must be extended upwards with at least ",
357 zas[n_ar_za-1] - za_grid[nza-1],
" deg.")
363 Vector aas{aresponse_aa_grid};
364 if (antenna_dlos.ncols() > 1) { aas += antenna_dlos(ia, 1); }
366 "The azimuth angle grid in *mblock_dlos* is too narrow. "
367 "It must be extended downwards with at least ",
368 aa_grid[0]-aas[0],
" deg.")
370 "The azimuth angle grid in *mblock_dlos* is too narrow. "
371 "It must be extended upwards with at least ",
372 aas[n_ar_aa-1] - aa_grid[naa-1],
" deg.")
379 Tensor3 itw(n_ar_za, n_ar_za, 4);
386 for (Index iaa = 0; iaa < n_ar_aa; iaa++) {
387 const Index
a = gp_aa[iaa].idx;
388 const Index
b =
a + 1;
390 for (Index iza = 0; iza < n_ar_za; iza++) {
392 const Index z = gp_za[iza].idx;
393 const Index x = z + 1;
395 if( itw(iza,iaa,0) > 1e-9 ) {
396 hza[
a*nza+z] += aresponse(iza,iaa) * itw(iza,iaa,0);
398 if( itw(iza,iaa,1) > 1e-9 ) {
399 hza[
b*nza+z] += aresponse(iza,iaa) * itw(iza,iaa,1);
401 if( itw(iza,iaa,2) > 1e-9 ) {
402 hza[
a*nza+x] += aresponse(iza,iaa) * itw(iza,iaa,2);
404 if( itw(iza,iaa,3) > 1e-9 ) {
405 hza[
b*nza+x] += aresponse(iza,iaa) * itw(iza,iaa,3);
416 const Index ii = f * n_pol + ip;
418 hrow[Range(ii, n_dlos, nfpol)] = hza;
420 H.insert_row(ia * nfpol + ii, hrow);
432 const Index& antenna_dim,
434 const Index& antenna_dim
_U_,
436 ConstMatrixView antenna_dlos,
438 ConstMatrixView mblock_dlos,
439 ConstVectorView solid_angles,
440 ConstVectorView f_grid,
444 const Index n_dlos = mblock_dlos.nrows();
445 const Index n_f = f_grid.nelem();
448 const Index n_ant = antenna_dlos.nrows();
457 const Index n_ar_pol =
459 ConstVectorView aresponse_f_grid =
461 ConstVectorView aresponse_za_grid =
463 ConstVectorView aresponse_aa_grid =
466 const Index n_ar_f = aresponse_f_grid.nelem();
467 const Index n_ar_za = aresponse_za_grid.nelem();
468 const Index n_ar_aa = aresponse_aa_grid.nelem();
469 const Index pol_step = n_ar_pol > 1;
482 Tensor4 aresponse_with_cos(n_ar_pol, n_ar_f, n_ar_za, n_ar_aa);
483 for(Index i3=0; i3<n_ar_za; i3++) {
484 const Numeric t = cos(
DEG2RAD * aresponse_za_grid[i3]);
485 for(Index i4=0; i4<n_ar_aa; i4++) {
486 for(Index i2=0; i2<n_ar_f; i2++) {
487 for(Index i1=0; i1<n_ar_pol; i1++) {
488 aresponse_with_cos(i1,i2,i3,i4) = t * antenna_response.
data(i1,i2,i3,i4);
495 const Index nfpol = n_f * n_pol;
498 H.resize(n_ant * nfpol, n_dlos * nfpol);
501 Vector hrow(H.ncols(), 0.0);
502 Vector hdlos(n_dlos, 0.0);
505 Matrix aresponse(n_ar_za, n_ar_aa, 0.0);
511 for (Index ia = 0; ia < n_ant; ia++) {
516 for (Index f = 0; f < n_f; f++) {
518 for (Index ip = 0; ip < n_pol; ip++) {
523 Index new_antenna = 1;
529 aresponse = antenna_response.
data(ip, 0, joker, joker);
530 }
else if (f == 0 && ip == 0)
532 aresponse = antenna_response.
data(0, 0, joker, joker);
538 if (ip == 0 || pol_step) {
541 gridpos(gp_f, aresponse_f_grid, Vector(1, f_grid[f]));
542 gridpos(gp_za, aresponse_za_grid, aresponse_za_grid);
543 gridpos(gp_aa, aresponse_aa_grid, aresponse_aa_grid);
544 Tensor4 itw(1, n_ar_za, n_ar_aa, 8);
546 Tensor3 aresponse_matrix(1, n_ar_za, n_ar_aa);
549 antenna_response.
data(ip, joker, joker, joker),
553 aresponse = aresponse_matrix(0, joker, joker);
562 for (Index l = 0; l < n_dlos; l++) {
563 const Numeric za = mblock_dlos(l, 0) - antenna_dlos(ia, 0);
565 if (mblock_dlos.ncols() > 1) {
566 aa += mblock_dlos(l, 1);
568 if (antenna_dlos.ncols() > 1) {
569 aa -= antenna_dlos(ia, 1);
574 if (za < aresponse_za_grid[0] ||
575 za > aresponse_za_grid[n_ar_za - 1] ||
576 aa < aresponse_aa_grid[0] ||
577 aa > aresponse_aa_grid[n_ar_aa - 1]) {
583 gridpos(gp_za, aresponse_za_grid, Vector(1, za));
584 gridpos(gp_aa, aresponse_aa_grid, Vector(1, aa));
588 interp(value, itw, aresponse, gp_za, gp_aa);
589 hdlos[l] = solid_angles[l] * value[0];
599 const Index ii = f * n_pol + ip;
601 hrow[Range(ii, n_dlos, nfpol)] = hdlos;
603 H.insert_row(ia * nfpol + ii, hrow);
617 ConstVectorView f_grid,
620 const Index& do_norm) {
657 const Numeric lim_low = 0;
658 const Numeric lim_high = -filter_grid[0];
662 list<Numeric> l_mixer;
663 for (Index i = 0; i < f_grid.nelem(); i++) {
664 if (fabs(f_grid[i] - lo) >= lim_low && fabs(f_grid[i] - lo) <= lim_high) {
665 l_mixer.push_back(fabs(f_grid[i] - lo));
668 l_mixer.push_back(lim_high);
671 f_mixer.resize((Index)l_mixer.size());
673 for (list<Numeric>::iterator li = l_mixer.begin(); li != l_mixer.end();
680 H.resize(f_mixer.nelem() * n_pol * n_sp, f_grid.nelem() * n_pol * n_sp);
685 Vector row_temp(f_grid.nelem());
686 Vector row_final(f_grid.nelem() * n_pol * n_sp);
688 Vector if_grid{f_grid};
691 for (Index i = 0; i < f_mixer.nelem(); i++) {
693 row_temp, filter.
data, filter_grid, if_grid, f_mixer[i], -f_mixer[i]);
696 if (do_norm) row_temp /= sum(row_temp);
699 for (Index p = 0; p < n_pol; p++) {
701 for (Index
a = 0;
a < n_sp;
a++) {
705 a * f_grid.nelem() * n_pol + p, f_grid.nelem(), n_pol)] = row_temp;
706 H.insert_row(
a * f_mixer.nelem() * n_pol + p + i * n_pol, row_final);
717 const Index stokes_dim,
723 if (iy_unit ==
"PlanckBT" || iy_unit ==
"RJBT") {
729 const Index nch = mm_pol.
nelem();
732 for (Index i = 0; i < nch; i++) {
733 if (mm_pol[i] ==
"AMSU-H") {
736 }
else if (mm_pol[i] ==
"AMSU-V") {
739 }
else if (mm_pol[i] ==
"ISMAR-H") {
742 }
else if (mm_pol[i] ==
"ISMAR-V") {
745 }
else if (mm_pol[i] ==
"MARSS-H") {
748 }
else if (mm_pol[i] ==
"MARSS-V") {
751 }
else if (mm_pol[i] ==
"H" || mm_pol[i] ==
"V" || mm_pol[i] ==
"LHC" ||
752 mm_pol[i] ==
"RHC") {
767 H = Sparse(nch, nch * stokes_dim);
769 for (Index i = 0; i < nch; i++) {
788 }
else if (pol[i] ==
"H") {
790 }
else if (pol[i] ==
"LHC")
793 }
else if (pol[i] ==
"RHC")
812 if (rot[i] ==
"none") {
814 Vector hrow(nch * stokes_dim, 0.0);
818 stokes2pol(hrow[Range(i * stokes_dim, stokes_dim)], stokes_dim, ipol,
w);
819 H.insert_row(i, hrow);
825 Sparse Hrot(stokes_dim, stokes_dim);
826 if (rot[i] ==
"AMSU") {
830 }
else if (rot[i] ==
"ISMAR") {
834 }
else if (rot[i] ==
"MARSS") {
848 Sparse Hpol(1, stokes_dim);
853 Vector hrow(stokes_dim);
855 Hpol.insert_row(0, hrow);
859 Sparse Hc(1, stokes_dim);
860 mult(Hc, Hpol, Hrot);
863 Vector hrow(nch * stokes_dim, 0.0);
864 const Index i0 = i * stokes_dim;
865 for (Index s = 0; s < stokes_dim; s++) {
866 hrow[i0 + s] = Hc(0, s);
868 H.insert_row(i, hrow);
877 Matrix& sensor_response_dlos,
878 ConstVectorView sensor_response_f_grid,
880 ConstMatrixView sensor_response_dlos_grid) {
882 const Index nf = sensor_response_f_grid.nelem();
883 const Index npol = sensor_response_pol_grid.
nelem();
884 const Index nlos = sensor_response_dlos_grid.nrows();
885 const Index n = nf * npol * nlos;
888 sensor_response_f.resize(n);
889 sensor_response_pol.resize(n);
890 sensor_response_dlos.resize(n, sensor_response_dlos_grid.ncols());
893 for (Index ilos = 0; ilos < nlos; ilos++) {
894 const Index i2 = ilos * nf * npol;
896 for (Index ifr = 0; ifr < nf; ifr++) {
897 const Index i3 = i2 + ifr * npol;
899 for (Index ip = 0; ip < npol; ip++) {
900 const Index i = i3 + ip;
902 sensor_response_f[i] = sensor_response_f_grid[ifr];
903 sensor_response_pol[i] = sensor_response_pol_grid[ip];
904 sensor_response_dlos(i, joker) = sensor_response_dlos_grid(ilos, joker);
913 ConstVectorView ch_f,
915 ConstVectorView sensor_f,
918 const Index& do_norm) {
924 Index freq_full = ch_response.
nelem() > 1;
931 const Index nin_f = sensor_f.nelem();
932 const Index nout_f = ch_f.nelem();
933 const Index nin = n_sp * nin_f * n_pol;
934 const Index nout = n_sp * nout_f * n_pol;
941 Vector ch_response_f;
942 Vector weights(nin_f);
943 Vector weights_long(nin, 0.0);
945 for (Index ifr = 0; ifr < nout_f; ifr++) {
946 const Index irp = ifr * freq_full;
949 ch_response_f = ch_response[irp].get_numeric_grid(GFIELD1_F_GRID);
950 ch_response_f += ch_f[ifr];
954 weights, ch_response[irp].data, ch_response_f, sensor_f);
957 if (do_norm) weights /= sum(weights);
961 for (Index sp = 0; sp < n_sp; sp++) {
962 for (Index pol = 0; pol < n_pol; pol++) {
964 weights_long[Range(sp * nin_f * n_pol + pol, nin_f, n_pol)] = weights;
967 H.insert_row(sp * nout_f * n_pol + ifr * n_pol + pol, weights_long);
979 const Index& stokes_dim,
980 const Index& ipol_1based,
985 "Valid polarization indices are 1 to 10 (1-based).");
987 ArrayOfVector s2p(10);
992 s2p[3] = {0, 0, 0, 1};
995 s2p[6] = {nv, 0, nv};
996 s2p[7] = {nv, 0, -nv};
997 s2p[8] = {nv, 0, 0, nv};
998 s2p[9] = {nv, 0, 0, -nv};
1000 const Index l = s2p[ipol_1based - 1].nelem();
1002 "You have selected polarization with 1-based index: ", ipol_1based,
1004 "but this polarization demands stokes_dim >= ", l,
"\n",
1005 "while the actual values of stokes_dim is ", stokes_dim)
1007 w[Range(0, l)] = s2p[ipol_1based - 1];
1008 if (l < stokes_dim) {
1009 w[Range(l, stokes_dim - l)] = 0;
1045 const Index nf = fmin.nelem();
1058 if (fmax[i] >= fmin[j]) {
1063 if (fmax[j] > fmax[i]) fmax[i] = fmax[j];
1068 Index n_behind = nf - 1 - j;
1070 Vector dummy = fmin;
1071 fmin.resize(nf - 1);
1072 fmin[Range(0, j)] = dummy[Range(0, j)];
1073 if (n_behind > 0) fmin[Range(j, n_behind)] = dummy[Range(j + 1, n_behind)];
1076 fmax.resize(nf - 1);
1077 fmax[Range(0, j)] = dummy[Range(0, j)];
1078 if (n_behind > 0) fmax[Range(j, n_behind)] = dummy[Range(j + 1, n_behind)];
1090 const Vector& f_backend,
1092 const Numeric& delta,
1097 const Index n_chan = f_backend.nelem();
1103 "There must be at least one channel.\n"
1104 "(The vector *f_backend* must have at least one element.)")
1108 "Variables *f_backend_multi* and *backend_channel_response_multi*\n"
1109 "must have same number of bands for each LO.")
1112 for (Index i = 0; i < n_chan; ++i) {
1114 const Vector& backend_f_grid =
1115 backend_channel_response[i].get_numeric_grid(0);
1118 "The frequency grid for the backend channel response of\n"
1119 "channel ", i,
" is not strictly increasing.\n"
1120 "It is: ", backend_f_grid,
"\n")
1125 out2 <<
" Original channel characteristics:\n"
1126 <<
" min nominal max (all in Hz):\n";
1133 for (Index idx = 0; idx < n_chan; ++idx) {
1134 const Vector& backend_filter = backend_channel_response[idx].data;
1135 if (backend_filter.nelem() >
1137 for (Index idy = 1; idy < backend_filter.nelem(); ++idy) {
1138 if ((backend_filter[idy] > 0) && (backend_filter[idy - 1] == 0)) {
1148 "No passbands found.\n"
1149 "*backend_channel_response* must be zero around the passbands.\n"
1150 "backend_channel_response.data = [0, >0, >0, 0]\n"
1151 "Borders between passbands are identified as [...0,0...]");
1153 Vector fmin_pb(numPB);
1154 Vector fmax_pb(numPB);
1157 for (Index idx = 0; idx < n_chan; ++idx) {
1190 const Vector& backend_f_grid =
1191 backend_channel_response[idx].get_numeric_grid(0);
1192 const Vector& backend_filter = backend_channel_response[idx].data;
1193 if (backend_filter.nelem() >=
1196 for (Index idy = 1; idy < backend_filter.nelem(); ++idy) {
1198 fmin_pb[pbIdx] = f_backend[idx] + backend_f_grid[0];
1199 }
else if ((backend_filter[idy] > 0) &&
1200 (backend_filter[idy - 1] == 0)) {
1201 fmin_pb[pbIdx] = f_backend[idx] + backend_f_grid[idy - 1] - delta;
1203 if ((backend_filter[idy] == 0) && (backend_filter[idy - 1] > 0)) {
1204 fmax_pb[pbIdx] = f_backend[idx] + backend_f_grid[idy] + delta;
1206 <<
"fmin_pb " << fmin_pb[pbIdx] <<
" "
1207 <<
"f_backend" << f_backend[idx] <<
" "
1208 <<
"fmax_pb " << fmax_pb[pbIdx] <<
" "
1209 <<
"diff " << fmax_pb[pbIdx] - fmin_pb[pbIdx] <<
"\n";
1213 fmax_pb[pbIdx - 1] = f_backend[idx] +
last(backend_f_grid);
1216 fmin_pb[pbIdx] = f_backend[idx] + backend_f_grid[0] - delta;
1217 fmax_pb[pbIdx] = f_backend[idx] +
1218 backend_f_grid[backend_f_grid.nelem() - 1] +
1221 <<
"fmin_pb " << fmin_pb[pbIdx] <<
" "
1222 <<
"f_backend" << f_backend[pbIdx] <<
" "
1223 <<
"fmax_pb " << fmax_pb[pbIdx] <<
"\n";
1247 out2 <<
" resize numPb " << numPB <<
"\n";
1248 for (Index idx = 0; idx < numPB; ++idx) {
1249 fmin[idx] = fmin_pb[isorted[idx]];
1250 fmax[idx] = fmax_pb[isorted[idx]];
1259 for (Index i = 0; i < fmin.nelem() - 1; ++i) {
1260 bool continue_checking =
true;
1264 while (continue_checking && i < fmin.nelem() - 1) {
1272 out2 <<
" New channel characteristics:\n"
1273 <<
" min max (all in Hz):\n";
1274 for (Index i = 0; i < fmin.nelem(); ++i)
1275 out2 <<
" " << fmin[i] <<
" " << fmax[i] <<
"\n";
1286 ConstVectorView x_f_in,
1287 ConstVectorView x_g_in) {
1289 const Index nf = x_f_in.nelem();
1290 const Index ng = x_g_in.nelem();
1296 ARTS_ASSERT(is_increasing(x_g_in) || is_decreasing(x_g_in));
1300 Numeric xfmin = x_f_in[0];
1301 Numeric xfmax = x_f_in[nf - 1];
1304 const bool xg_reversed = is_decreasing(x_g_in);
1305 Vector x_g{xg_reversed ? reverse(x_g_in) : Vector{x_g_in}};
1312 const Numeric df = xfmax - xfmin;
1315 for (Index i = 0; i < nf; i++) {
1316 x_f[i] = (x_f_in[i] - xfmin) / df;
1318 for (Index i = 0; i < ng; i++) {
1319 x_g[i] = (x_g[i] - xfmin) / df;
1331 for (Index it = 0; it < nf; it++) {
1332 l_x.push_back(x_f[it]);
1334 for (Index it = 0; it < ng; it++) {
1335 if (x_g[it] > xfmin && x_g[it] < xfmax) {
1336 l_x.push_back(x_g[it]);
1342 Vector x_ref(l_x.size());
1344 for (list<Numeric>::iterator li = l_x.begin(); li != l_x.end(); li++) {
1354 Numeric dx, a0, b0, c0, a1, b1, c1, x3, x2, x1;
1356 for (Index i = 0; i < x_ref.nelem() - 1; i++) {
1359 while (x_g[i_g + 1] <= x_ref[i]) {
1362 while (x_f[i_f + 1] <= x_ref[i]) {
1368 if (x_ref[i] >= xfmin && x_ref[i] < xfmax) {
1370 dx = (x_f[i_f + 1] - x_f[i_f]) * (x_g[i_g + 1] - x_g[i_g]);
1373 a0 = (f[i_f] - f[i_f + 1]) / 3.0;
1374 b0 = (-f[i_f] * (x_g[i_g + 1] + x_f[i_f + 1]) +
1375 f[i_f + 1] * (x_g[i_g + 1] + x_f[i_f])) /
1377 c0 = x_g[i_g + 1] * (f[i_f] * x_f[i_f + 1] - f[i_f + 1] * x_f[i_f]);
1380 b1 = (f[i_f] * (x_g[i_g] + x_f[i_f + 1]) -
1381 f[i_f + 1] * (x_g[i_g] + x_f[i_f])) /
1383 c1 = x_g[i_g] * (-f[i_f] * x_f[i_f + 1] + f[i_f + 1] * x_f[i_f]);
1385 x1 = x_ref[i + 1] - x_ref[i];
1390 x2 = x1 * (2 * x_ref[i] + x1);
1391 x3 = x1 * (3 * x_ref[i] * (x_ref[i] + x1) + x1 * x1);
1395 h[i_g] += df * (a0 * x3 + b0 * x2 + c0 * x1) / dx;
1396 h[i_g + 1] += df * (a1 * x3 + b1 * x2 + c1 * x1) / dx;
1402 Vector tmp = reverse(h);
1410 "Significant negative response value obtained, "
1411 "despite sensor reponse is strictly positive. This "
1412 "indicates numerical problems. Is there any very "
1413 "small spacing of the sensor response grid?"
1414 "Please, send a report to Patrick if you see this!");
1420 ConstVectorView x_g_in,
1421 const Numeric& limit1,
1422 const Numeric& limit2) {
1424 const Index ng = x_g_in.nelem();
1432 const bool xg_reversed = is_decreasing(x_g_in);
1433 const Vector x_g{xg_reversed ? reverse(x_g_in) : Vector{x_g_in}};
1441 if (limit1 == limit2) {
1446 else if (limit1 == x_g[0] && limit2 == x_g[ng - 1]) {
1447 h[0] = (x_g[1] - x_g[0]) / 2.0;
1448 for (Index i = 1; i < ng - 1; i++) {
1449 h[i] = (x_g[i + 1] - x_g[i - 1]) / 2.0;
1451 h[ng - 1] = (x_g[ng - 1] - x_g[ng - 2]) / 2.0;
1459 for (Index i = 0; i < ng; i++) {
1460 bool inside =
false;
1463 if (limit1 < x_g[1]) {
1466 x2 =
min(limit2, x_g[1]);
1468 }
else if (i == ng - 1) {
1469 if (limit2 > x_g[ng - 2]) {
1471 x1 =
max(limit1, x_g[ng - 2]);
1475 if ((limit1 < x_g[i + 1] && limit2 > x_g[i - 1]) ||
1476 (limit2 > x_g[i - 1] && limit1 < x_g[i + 1])) {
1478 x1 =
max(limit1, x_g[i - 1]);
1479 x2 =
min(limit2, x_g[i + 1]);
1489 const Numeric r = 1.0 / (x_g[i] - x_g[i - 1]);
1490 const Numeric y1 = r * (x1 - x_g[i - 1]);
1491 const Numeric dx =
min(x2, x_g[i]) - x1;
1492 const Numeric y2 = y1 + r * dx;
1493 h[i] = 0.5 * dx * (y1 + y2);
1500 const Numeric r = 1.0 / (x_g[i + 1] - x_g[i]);
1501 const Numeric y2 = r * (x_g[i + 1] - x2);
1502 const Numeric dx = x2 -
max(x1, x_g[i]);
1503 const Numeric y1 = y2 + r * dx;
1504 h[i] += 0.5 * dx * (y1 + y2);
1511 Vector tmp = reverse(h);
1519 ConstVectorView x_f,
1520 ConstVectorView x_g,
1536 gridpos(gp1g, x_g, ExhaustiveConstVectorView{x1});
1537 gridpos(gp1f, x_f, ExhaustiveConstVectorView{x1});
1541 interp(ExhaustiveVectorView{f1}, itw1, f, gp1f);
1545 gridpos(gp2g, x_g, ExhaustiveConstVectorView{x2});
1546 gridpos(gp2f, x_f, ExhaustiveConstVectorView{x2});
1550 interp(ExhaustiveVectorView{f2}, itw2, f, gp2f);
1554 h[gp1g[0].idx] += f1 * gp1g[0].fd[1];
1555 h[gp1g[0].idx + 1] += f1 * gp1g[0].fd[0];
1556 h[gp2g[0].idx] += f2 * gp2g[0].fd[1];
1557 h[gp2g[0].idx + 1] += f2 * gp2g[0].fd[0];
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
const ArrayOfString & get_string_grid(Index i) const
Get a string 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(...)
#define ARTS_USER_ERROR_IF(condition,...)
Implementation of gridded fields.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
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 auto deg2rad(auto x) noexcept
Converts degrees to radians.
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.
void muellersparse_rotation(Sparse &H, const Index &stokes_dim, const Numeric &rotangle)
muellersparse_rotation
Declaration of functions in rte.cc.
void antenna2d_gridded_dlos(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_dlos, const GriddedField4 &antenna_response, ConstMatrixView mblock_dlos, ConstVectorView f_grid, const Index n_pol)
antenna2d_interp_gridded_dlos
void find_effective_channel_boundaries(Vector &fmin, Vector &fmax, const Vector &f_backend, const ArrayOfGriddedField1 &backend_channel_response, const Numeric &delta, const Verbosity &verbosity)
Calculate channel boundaries from instrument response functions.
void integration_bin_by_vecmult(VectorView h, ConstVectorView x_g_in, const Numeric &limit1, const Numeric &limit2)
integration_bin_by_vecmult
constexpr Numeric DEG2RAD
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.
constexpr Numeric NAT_LOG_2
void summation_by_vecmult(VectorView h, ConstVectorView f, ConstVectorView x_f, ConstVectorView x_g, const Numeric x1, const Numeric x2)
summation_by_vecmult
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
bool test_and_merge_two_channels(Vector &fmin, Vector &fmax, Index i, Index j)
Test if two instrument channels overlap, and if so, merge them.
void spectrometer_matrix(Sparse &H, ConstVectorView ch_f, const ArrayOfGriddedField1 &ch_response, ConstVectorView sensor_f, const Index &n_pol, const Index &n_sp, const Index &do_norm)
spectrometer_matrix
void antenna2d_interp_response(Sparse &H, const Index &antenna_dim, ConstMatrixView antenna_dlos, const GriddedField4 &antenna_response, ConstMatrixView mblock_dlos, ConstVectorView solid_angles, ConstVectorView f_grid, const Index n_pol)
antenna2d_interp_response
void integration_func_by_vecmult(VectorView h, ConstVectorView f, ConstVectorView x_f_in, ConstVectorView x_g_in)
integration_func_by_vecmult
void antenna1d_matrix(Sparse &H, const Index &antenna_dim, ConstVectorView antenna_dza, const GriddedField4 &antenna_response, ConstVectorView za_grid, ConstVectorView f_grid, const Index n_pol, const Index do_norm)
antenna1d_matrix
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.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes