ARTS 2.5.0 (git: 9ee3ac6c)
interpolation.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
40#include "interpolation.h"
41#include <cmath>
42#include <iostream>
43#include "array.h"
44#include "check_input.h"
45#include "logic.h"
46
47// File-global constants:
48
50
61DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6;)
62
64
73const Numeric FD_TOL = 1.5e-3;
74
76
85#define LOOPIT(x) for (const Numeric* x = &t##x.fd[1]; x >= &t##x.fd[0]; --x)
86
88
96ostream& operator<<(ostream& os, const GridPos& gp) {
97 os << gp.idx << " " << gp.fd[0] << " " << gp.fd[1] << "\n";
98 return os;
99}
100
102
157 ConstVectorView old_grid,
158 ConstVectorView new_grid,
159 const Numeric& extpolfac) {
160 const Index n_old = old_grid.nelem();
161 const Index n_new = new_grid.nelem();
162
163 // Assert that gp has the right size:
164 ARTS_ASSERT(is_size(gp, n_new));
165
166 // Assert, that the old grid has more than one element
167 ARTS_ASSERT(1 < n_old);
168
169 // This function hast two parts, depending on whether old_grid is
170 // sorted in ascending or descending order. Maybe that's not too
171 // elegant, but it's the most efficient, because in this way there
172 // is no additional runtime overhead to handle both cases.
173
174 // We use only the first two elements to decide how the grid is
175 // sorted. (The rest of the grid is checked later by an ARTS_ASSERT.)
176 // If both values are the same, we still assume the grid is
177 // ascending. However, this will lead to an assertion fail later on,
178 // because the grid has to be strictly sorted.
179 bool ascending = (old_grid[0] <= old_grid[1]);
180
181 if (ascending) {
182 // So, old_grid should always be sorted in strictly ascending order.
183 // (No duplicate values.)
184 // Assert that this is so. This may depend on user input, however,
185 // inside this elementary function is not the place to check for
186 // that. There must be runtime checks on higher levels to insure
187 // that all grids are sorted. The assertion here is just as a last
188 // safety check.
189 ARTS_ASSERT(is_increasing(old_grid));
190
191 // Limits of extrapolation.
192 const Numeric og_min =
193 old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
194 const Numeric og_max =
195 old_grid[n_old - 1] +
196 extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
197
198 // We will make no firm assumptions about the new grid. But the case
199 // that we have in mind is that it is either also sorted, or at
200 // least partially sorted, for example like this:
201 // 5 4 3 2 2.5 3 4
202 // This kind of sequence should be typical if we interpolate
203 // atmospheric fields along a limb propagation path.
204
205 // Let's get some idea where the first point in the new grid is,
206 // relative to the old grid. We use linear interpolation between the
207 // maximum and the minimum of the old grid for this. (Taking
208 // into account the small allowed extrapolation.)
209 Numeric frac = (new_grid[0] - og_min) / (og_max - og_min);
210
211 // We are not checking if the value of frac is reasonable,
212 // because there is another assertion below to catch the
213 // consequences.
214
215 // Initialize current_position:
216 Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
217
218 // The above statement should satisfy
219 // 0 <= current_position <= n_old-2
220 // Assert that this is indeed the case.
221 // cout << "frac = " << frac << "\n";
222 // cout << "current_position = " << current_position << "\n";
223 ARTS_ASSERT(0 <= current_position);
224 ARTS_ASSERT(current_position <= n_old - 2);
225
226 // The variables lower and upper are used to remember the value of
227 // the old grid at current_position and one above current_position:
228 Numeric lower = old_grid[current_position];
229 Numeric upper = old_grid[current_position + 1];
230
231 // Loop over all points in the new grid:
232 for (Index i_new = 0; i_new < n_new; ++i_new) {
233 // Get a reference to the current element of gp:
234 GridPos& tgp = gp[i_new];
235 // And on the current value of the new grid:
236 const Numeric tng = new_grid[i_new];
237
238 // Verify that the new grid is within the limits of
239 // extrapolation that we have defined by og_min and og_max:
240 ARTS_ASSERT(og_min <= tng);
241 ARTS_ASSERT(tng <= og_max);
242
243 // cout << "lower / tng / upper = "
244 // << lower << " / "
245 // << tng << " / "
246 // << upper << "\n";
247
248 // Is current_position too high?
249 // (The current_position>0 condition is there so that the position
250 // stays 0 for extrapolation.)
251 if (tng < lower && current_position > 0) {
252 do {
253 --current_position;
254 lower = old_grid[current_position];
255 } while (tng < lower && current_position > 0);
256
257 upper = old_grid[current_position + 1];
258
259 tgp.idx = current_position;
260 tgp.fd[0] = (tng - lower) / (upper - lower);
261 tgp.fd[1] = 1.0 - tgp.fd[0];
262 } else {
263 // Is it too low?
264 // (The current_position<n_old condition is there so
265 // that uppers stays n_old-1 for extrapolation.)
266 if (tng >= upper && current_position < n_old - 2) {
267 do {
268 ++current_position;
269 upper = old_grid[current_position + 1];
270 } while (tng >= upper && current_position < n_old - 2);
271
272 lower = old_grid[current_position];
273
274 tgp.idx = current_position;
275 tgp.fd[0] = (tng - lower) / (upper - lower);
276 tgp.fd[1] = 1.0 - tgp.fd[0];
277 } else {
278 // None of the other two conditions were true. That means:
279 // lower <= tng < upper. The current_position is still
280 // valid.
281
282 // SAB 2010-04-28: Note that if a new grid point is exactly on
283 // top of an old grid point, then you are now guaranteed to get
284 // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
285
286 tgp.idx = current_position;
287 tgp.fd[0] = (tng - lower) / (upper - lower);
288 tgp.fd[1] = 1.0 - tgp.fd[0];
289 }
290 }
291
292 // cout << tgp.idx << " " << tgp.fd[0] << " " << tgp.fd[1] << endl;
293
294 // Safety check to ensure the above promise:
295 ARTS_ASSERT(old_grid[tgp.idx] <= tng || tgp.idx == 0);
296 }
297 } else // if (ascending)
298 {
299 // Now we are in the "descending old grid" part. We do exactly
300 // the same as in the other part, just accounting for the
301 // different order of things. Comments here refer only to
302 // interesting differences from the ascending case. See that
303 // case for more general comments.
304
305 // This time ensure strictly descending order:
306 ARTS_ASSERT(is_decreasing(old_grid));
307
308 // The max is now the first point, the min the last point!
309 // I think the sign is right here...
310 const Numeric og_max =
311 old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
312 const Numeric og_min =
313 old_grid[n_old - 1] +
314 extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
315
316 // We have to take 1- here, because we are starting from the
317 // high end.
318 Numeric frac = 1 - (new_grid[0] - og_min) / (og_max - og_min);
319
320 // We are not checking if the value of frac is reasonable,
321 // because there is another assertion below to catch the
322 // consequences.
323
324 Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
325
326 // The above statement should satisfy
327 // 0 <= current_position <= n_old-2
328 // Assert that this is indeed the case.
329 // cout << "frac = " << frac << "\n";
330 // cout << "current_position = " << current_position << "\n";
331 ARTS_ASSERT(0 <= current_position);
332 ARTS_ASSERT(current_position <= n_old - 2);
333
334 // Note, that old_grid[lower] has a higher numerical value than
335 // old_grid[upper]!
336 Numeric lower = old_grid[current_position];
337 Numeric upper = old_grid[current_position + 1];
338
339 for (Index i_new = 0; i_new < n_new; ++i_new) {
340 GridPos& tgp = gp[i_new];
341 const Numeric tng = new_grid[i_new];
342
343 // Verify that the new grid is within the limits of
344 // extrapolation that we have defined by og_min and og_max:
345 ARTS_ASSERT(og_min <= tng);
346 ARTS_ASSERT(tng <= og_max);
347
348 // cout << "lower / tng / upper = "
349 // << lower << " / "
350 // << tng << " / "
351 // << upper << "\n";
352
353 // Is current_position too high? (Sign of comparison changed
354 // compared to ascending case!)
355 if (tng > lower && current_position > 0) {
356 do {
357 --current_position;
358 lower = old_grid[current_position];
359 } while (tng > lower && current_position > 0);
360
361 upper = old_grid[current_position + 1];
362
363 tgp.idx = current_position;
364 tgp.fd[0] = (tng - lower) / (upper - lower);
365 tgp.fd[1] = 1.0 - tgp.fd[0];
366 } else {
367 // Is it too low? (Sign of comparison changed
368 // compared to ascending case!)
369 if (tng <= upper && current_position < n_old - 2) {
370 do {
371 ++current_position;
372 upper = old_grid[current_position + 1];
373 } while (tng <= upper && current_position < n_old - 2);
374
375 lower = old_grid[current_position];
376
377 tgp.idx = current_position;
378 tgp.fd[0] = (tng - lower) / (upper - lower);
379 tgp.fd[1] = 1.0 - tgp.fd[0];
380 } else {
381 // None of the other two conditions were true. That means:
382 // lower >= tng > upper. The current_position is still
383 // valid. (Note that upper and lower have switched
384 // place compared to the ascending case.)
385
386 // SAB 2010-04-28: Note that if a new grid point is exactly on
387 // top of an old grid point, then you are now guaranteed to get
388 // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
389
390 tgp.idx = current_position;
391 tgp.fd[0] = (tng - lower) / (upper - lower);
392 tgp.fd[1] = 1.0 - tgp.fd[0];
393 }
394 }
395
396 // Safety check to ensure the above promise:
397 ARTS_ASSERT(old_grid[tgp.idx] >= tng || tgp.idx == 0);
398 }
399 }
400}
401
403
425 ConstVectorView old_grid,
426 const Numeric& new_grid,
427 const Numeric& extpolfac) {
428 ArrayOfGridPos agp(1);
429 Vector v(1, new_grid);
430 gridpos(agp, old_grid, v, extpolfac);
431 gridpos_copy(gp, agp[0]);
432}
433
435
450 const Index n = grid.nelem();
451 gp.resize(n);
452
453 for (Index i = 0; i < n - 1; i++) {
454 gp[i].idx = i;
455 gp[i].fd[0] = 0;
456 gp[i].fd[1] = 1;
457 }
458 gp[n - 1].idx = n - 2;
459 gp[n - 1].fd[0] = 1;
460 gp[n - 1].fd[1] = 0;
461}
462
464
473void gridpos_copy(GridPos& gp_new, const GridPos& gp_old) {
474 gp_new.idx = gp_old.idx;
475 gp_new.fd[0] = gp_old.fd[0];
476 gp_new.fd[1] = gp_old.fd[1];
477}
478
480
493 return (Numeric(gp.idx) + gp.fd[0]);
494}
495
497
510 // Catch values that "must" be wrong
511 ARTS_ASSERT(gp.fd[0] > -FD_TOL);
512 ARTS_ASSERT(gp.fd[0] < 1.0 + FD_TOL);
513 ARTS_ASSERT(gp.fd[1] > -FD_TOL);
514 ARTS_ASSERT(gp.fd[1] < 1.0 + FD_TOL);
515
516 if (gp.fd[0] < 0.0) {
517 gp.fd[0] = 0.0;
518 gp.fd[1] = 1.0;
519 } else if (gp.fd[0] > 1.0) {
520 gp.fd[0] = 1.0;
521 gp.fd[1] = 0.0;
522 }
523
524 if (gp.fd[1] < 0.0) {
525 gp.fd[1] = 0.0;
526 gp.fd[0] = 1.0;
527 } else if (gp.fd[1] > 1.0) {
528 gp.fd[1] = 1.0;
529 gp.fd[0] = 0.0;
530 }
531}
532
534
555void gridpos_force_end_fd(GridPos& gp, const Index& n) {
556 ARTS_ASSERT(gp.idx >= 0);
557
558 // If fd=1, shift to grid index above
559 if (gp.fd[0] > 0.5) {
560 gp.idx += 1;
561 }
562 gp.fd[0] = 0;
563 gp.fd[1] = 1;
564
565 ARTS_ASSERT(gp.idx < n);
566
567 // End of complete grid range must be handled separately
568 if (gp.idx == n - 1) {
569 gp.idx -= 1;
570 gp.fd[0] = 1;
571 gp.fd[1] = 0;
572 }
573}
574
576
590 if (gp.idx == ie) {
591 ARTS_ASSERT(gp.fd[0] < 0.005); // To capture obviously bad cases
592 gp.idx -= 1;
593 gp.fd[0] = 1.0;
594 gp.fd[1] = 0.0;
595 }
596}
597
599
613 for (Index i = 0; i < gp.nelem(); i++) {
614 if (gp[i].idx == ie) {
615 ARTS_ASSERT(gp[i].fd[0] < 0.005); // To capture obviously bad cases
616 gp[i].idx -= 1;
617 gp[i].fd[0] = 1.0;
618 gp[i].fd[1] = 0.0;
619 }
620 }
621}
622
624
637 for (Index i = 0; i < gp.nelem(); i++) {
638 gp[i].idx = 0;
639 gp[i].fd[0] = 0;
640 gp[i].fd[1] = 1;
641 }
642}
643
645
658 const Index& i,
659 const bool& strict) {
660 if (strict) {
661 // Assume that gridpos_force_end_fd has been used. The expression 0==0
662 // should be safer than 1==1.
663 if (gp.idx == i && gp.fd[0] == 0) {
664 return true;
665 } else if (gp.idx == i - 1 && gp.fd[1] == 0) {
666 return true;
667 }
668 } else {
669 if (gp.idx == i && gp.fd[0] < FD_TOL) {
670 return true;
671 } else if (gp.idx == i - 1 && gp.fd[1] < FD_TOL) {
672 return true;
673 }
674 }
675 return false;
676}
677
679
697Index gridpos2gridrange(const GridPos& gp, const bool& upwards) {
698 ARTS_ASSERT(gp.fd[0] >= 0);
699 ARTS_ASSERT(gp.fd[1] >= 0);
700
701 // Not at a grid point
702 if (gp.fd[0] > 0 && gp.fd[1] > 0) {
703 return gp.idx;
704 }
705
706 // Fractional distance 0
707 else if (gp.fd[0] == 0) {
708 if (upwards)
709 return gp.idx;
710 else
711 return gp.idx - 1;
712 }
713
714 // Fractional distance 1
715 else {
716 if (upwards)
717 return gp.idx + 1;
718 else
719 return gp.idx;
720 }
721}
722
724// Red Interpolation
726
728
741void interpweights(VectorView itw, const GridPos& tc) {
742 ARTS_ASSERT(is_size(itw, 2)); // We must store 2 interpolation
743 // weights.
744
745 // Interpolation weights are stored in this order (l=lower
746 // u=upper, c=column):
747 // 1. l-c
748 // 2. u-c
749
750 Index iti = 0;
751
752 // We could use a straight-forward for loop here:
753 //
754 // for ( Index c=1; c>=0; --c )
755 // {
756 // ti[iti] = tc.fd[c];
757 // ++iti;
758 // }
759 //
760 // However, it is slightly faster to use pointers (I tried it,
761 // the speed gain is about 20%). So we should write something
762 // like:
763 //
764 // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
765 // {
766 // ti[iti] = *c;
767 // ++iti;
768 // }
769 //
770 // For higher dimensions we have to nest these loops. To avoid
771 // typos and safe typing, I use the LOOPIT macro, which expands
772 // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
773 // COMMAND!
774
775 LOOPIT(c) {
776 itw.get(iti) = *c;
777 ++iti;
778 }
779}
780
782
796void interpweights(VectorView itw, const GridPos& tr, const GridPos& tc) {
797 ARTS_ASSERT(is_size(itw, 4)); // We must store 4 interpolation
798 // weights.
799 Index iti = 0;
800
801 LOOPIT(r)
802 LOOPIT(c) {
803 itw.get(iti) = (*r) * (*c);
804 ++iti;
805 }
806}
807
809
825 const GridPos& tp,
826 const GridPos& tr,
827 const GridPos& tc) {
828 ARTS_ASSERT(is_size(itw, 8)); // We must store 8 interpolation
829 // weights.
830 Index iti = 0;
831
832 LOOPIT(p)
833 LOOPIT(r)
834 LOOPIT(c) {
835 itw.get(iti) = (*p) * (*r) * (*c);
836 ++iti;
837 }
838}
839
841
858 const GridPos& tb,
859 const GridPos& tp,
860 const GridPos& tr,
861 const GridPos& tc) {
862 ARTS_ASSERT(is_size(itw, 16)); // We must store 16 interpolation
863 // weights.
864 Index iti = 0;
865
866 LOOPIT(b)
867 LOOPIT(p)
868 LOOPIT(r)
869 LOOPIT(c) {
870 itw.get(iti) = (*b) * (*p) * (*r) * (*c);
871 ++iti;
872 }
873}
874
876
894 const GridPos& ts,
895 const GridPos& tb,
896 const GridPos& tp,
897 const GridPos& tr,
898 const GridPos& tc) {
899 ARTS_ASSERT(is_size(itw, 32)); // We must store 32 interpolation
900 // weights.
901 Index iti = 0;
902
903 LOOPIT(s)
904 LOOPIT(b)
905 LOOPIT(p)
906 LOOPIT(r)
907 LOOPIT(c) {
908 itw.get(iti) = (*s) * (*b) * (*p) * (*r) * (*c);
909 ++iti;
910 }
911}
912
914
933 const GridPos& tv,
934 const GridPos& ts,
935 const GridPos& tb,
936 const GridPos& tp,
937 const GridPos& tr,
938 const GridPos& tc) {
939 ARTS_ASSERT(is_size(itw, 64)); // We must store 64 interpolation
940 // weights.
941 Index iti = 0;
942
943 LOOPIT(v)
944 LOOPIT(s)
945 LOOPIT(b)
946 LOOPIT(p)
947 LOOPIT(r)
948 LOOPIT(c) {
949 itw.get(iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
950 ++iti;
951 }
952}
953
955
971 ARTS_ASSERT(is_size(itw, 2)); // We need 2 interpolation
972 // weights.
973
974 // Check that interpolation weights are valid. The sum of all
975 // weights (last dimension) must always be approximately one.
976 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
977
978 // To store interpolated value:
979 Numeric tia = 0;
980
981 Index iti = 0;
982 for (Index c = 0; c < 2; ++c) {
983 tia += a.get(tc.idx + c) * itw.get(iti);
984 ++iti;
985 }
986
987 return tia;
988}
989
991
1009 const GridPos& tr,
1010 const GridPos& tc) {
1011 ARTS_ASSERT(is_size(itw, 4)); // We need 4 interpolation
1012 // weights.
1013
1014 // Check that interpolation weights are valid. The sum of all
1015 // weights (last dimension) must always be approximately one.
1016 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1017
1018 // To store interpolated value:
1019 Numeric tia = 0;
1020
1021 Index iti = 0;
1022 for (Index r = 0; r < 2; ++r)
1023 for (Index c = 0; c < 2; ++c) {
1024 tia += a.get(tr.idx + r, tc.idx + c) * itw.get(iti);
1025 ++iti;
1026 }
1027
1028 return tia;
1029}
1030
1032
1051 const GridPos& tp,
1052 const GridPos& tr,
1053 const GridPos& tc) {
1054 ARTS_ASSERT(is_size(itw, 8)); // We need 8 interpolation
1055 // weights.
1056
1057 // Check that interpolation weights are valid. The sum of all
1058 // weights (last dimension) must always be approximately one.
1059 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1060
1061 // To store interpolated value:
1062 Numeric tia = 0;
1063
1064 Index iti = 0;
1065 for (Index p = 0; p < 2; ++p)
1066 for (Index r = 0; r < 2; ++r)
1067 for (Index c = 0; c < 2; ++c) {
1068 tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) * itw.get(iti);
1069 ++iti;
1070 }
1071
1072 return tia;
1073}
1074
1076
1096 const GridPos& tb,
1097 const GridPos& tp,
1098 const GridPos& tr,
1099 const GridPos& tc) {
1100 ARTS_ASSERT(is_size(itw, 16)); // We need 16 interpolation
1101 // weights.
1102
1103 // Check that interpolation weights are valid. The sum of all
1104 // weights (last dimension) must always be approximately one.
1105 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1106
1107 // To store interpolated value:
1108 Numeric tia = 0;
1109
1110 Index iti = 0;
1111 for (Index b = 0; b < 2; ++b)
1112 for (Index p = 0; p < 2; ++p)
1113 for (Index r = 0; r < 2; ++r)
1114 for (Index c = 0; c < 2; ++c) {
1115 tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
1116 itw.get(iti);
1117 ++iti;
1118 }
1119
1120 return tia;
1121}
1122
1124
1145 const GridPos& ts,
1146 const GridPos& tb,
1147 const GridPos& tp,
1148 const GridPos& tr,
1149 const GridPos& tc) {
1150 ARTS_ASSERT(is_size(itw, 32)); // We need 32 interpolation
1151 // weights.
1152
1153 // Check that interpolation weights are valid. The sum of all
1154 // weights (last dimension) must always be approximately one.
1155 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1156
1157 // To store interpolated value:
1158 Numeric tia = 0;
1159
1160 Index iti = 0;
1161 for (Index s = 0; s < 2; ++s)
1162 for (Index b = 0; b < 2; ++b)
1163 for (Index p = 0; p < 2; ++p)
1164 for (Index r = 0; r < 2; ++r)
1165 for (Index c = 0; c < 2; ++c) {
1166 tia += a.get(ts.idx + s,
1167 tb.idx + b,
1168 tp.idx + p,
1169 tr.idx + r,
1170 tc.idx + c) *
1171 itw.get(iti);
1172 ++iti;
1173 }
1174
1175 return tia;
1176}
1177
1179
1201 const GridPos& tv,
1202 const GridPos& ts,
1203 const GridPos& tb,
1204 const GridPos& tp,
1205 const GridPos& tr,
1206 const GridPos& tc) {
1207 ARTS_ASSERT(is_size(itw, 64)); // We need 64 interpolation
1208 // weights.
1209
1210 // Check that interpolation weights are valid. The sum of all
1211 // weights (last dimension) must always be approximately one.
1212 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1213
1214 // To store interpolated value:
1215 Numeric tia = 0;
1216
1217 Index iti = 0;
1218 for (Index v = 0; v < 2; ++v)
1219 for (Index s = 0; s < 2; ++s)
1220 for (Index b = 0; b < 2; ++b)
1221 for (Index p = 0; p < 2; ++p)
1222 for (Index r = 0; r < 2; ++r)
1223 for (Index c = 0; c < 2; ++c) {
1224 tia += a.get(tv.idx + v,
1225 ts.idx + s,
1226 tb.idx + b,
1227 tp.idx + p,
1228 tr.idx + r,
1229 tc.idx + c) *
1230 itw.get(iti);
1231 ++iti;
1232 }
1233
1234 return tia;
1235}
1236
1238// Blue interpolation
1240
1242
1258 Index n = cgp.nelem();
1259 ARTS_ASSERT(is_size(itw, n, 2)); // We must store 2 interpolation
1260 // weights for each position.
1261
1262 // We have to loop all the points in the sequence:
1263 for (Index i = 0; i < n; ++i) {
1264 // Current grid positions:
1265 const GridPos& tc = cgp[i];
1266
1267 // Interpolation weights are stored in this order (l=lower
1268 // u=upper, c=column):
1269 // 1. l-c
1270 // 2. u-c
1271
1272 Index iti = 0;
1273
1274 // We could use a straight-forward for loop here:
1275 //
1276 // for ( Index c=1; c>=0; --c )
1277 // {
1278 // ti[iti] = tc.fd[c];
1279 // ++iti;
1280 // }
1281 //
1282 // However, it is slightly faster to use pointers (I tried it,
1283 // the speed gain is about 20%). So we should write something
1284 // like:
1285 //
1286 // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
1287 // {
1288 // ti[iti] = *c;
1289 // ++iti;
1290 // }
1291 //
1292 // For higher dimensions we have to nest these loops. To avoid
1293 // typos and safe typing, I use the LOOPIT macro, which expands
1294 // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
1295 // COMMAND!
1296
1297 LOOPIT(c) {
1298 itw.get(i, iti) = *c;
1299 ++iti;
1300 }
1301 }
1302}
1303
1305
1327 const ArrayOfGridPos& rgp,
1328 const ArrayOfGridPos& cgp) {
1329 Index n = cgp.nelem();
1330 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1331 ARTS_ASSERT(is_size(itw, n, 4)); // We must store 4 interpolation
1332 // weights for each position.
1333
1334 // We have to loop all the points in the sequence:
1335 for (Index i = 0; i < n; ++i) {
1336 // Current grid positions:
1337 const GridPos& tr = rgp[i];
1338 const GridPos& tc = cgp[i];
1339
1340 // Interpolation weights are stored in this order (l=lower
1341 // u=upper, r=row, c=column):
1342 // 1. l-r l-c
1343 // 2. l-r u-c
1344 // 3. u-r l-c
1345 // 4. u-r u-c
1346
1347 Index iti = 0;
1348
1349 LOOPIT(r)
1350 LOOPIT(c) {
1351 itw.get(i, iti) = (*r) * (*c);
1352 ++iti;
1353 }
1354 }
1355}
1356
1358
1381 const ArrayOfGridPos& pgp,
1382 const ArrayOfGridPos& rgp,
1383 const ArrayOfGridPos& cgp) {
1384 Index n = cgp.nelem();
1385 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1386 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1387 ARTS_ASSERT(is_size(itw, n, 8)); // We must store 8 interpolation
1388 // weights for each position.
1389
1390 // We have to loop all the points in the sequence:
1391 for (Index i = 0; i < n; ++i) {
1392 // Current grid positions:
1393 const GridPos& tp = pgp[i];
1394 const GridPos& tr = rgp[i];
1395 const GridPos& tc = cgp[i];
1396
1397 Index iti = 0;
1398 LOOPIT(p)
1399 LOOPIT(r)
1400 LOOPIT(c) {
1401 itw.get(i, iti) = (*p) * (*r) * (*c);
1402 ++iti;
1403 }
1404 }
1405}
1406
1408
1432 const ArrayOfGridPos& bgp,
1433 const ArrayOfGridPos& pgp,
1434 const ArrayOfGridPos& rgp,
1435 const ArrayOfGridPos& cgp) {
1436 Index n = cgp.nelem();
1437 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1438 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1439 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1440 ARTS_ASSERT(is_size(itw, n, 16)); // We must store 16 interpolation
1441 // weights for each position.
1442
1443 // We have to loop all the points in the sequence:
1444 for (Index i = 0; i < n; ++i) {
1445 // Current grid positions:
1446 const GridPos& tb = bgp[i];
1447 const GridPos& tp = pgp[i];
1448 const GridPos& tr = rgp[i];
1449 const GridPos& tc = cgp[i];
1450
1451 Index iti = 0;
1452 LOOPIT(b)
1453 LOOPIT(p)
1454 LOOPIT(r)
1455 LOOPIT(c) {
1456 itw.get(i, iti) = (*b) * (*p) * (*r) * (*c);
1457 ++iti;
1458 }
1459 }
1460}
1461
1463
1488 const ArrayOfGridPos& sgp,
1489 const ArrayOfGridPos& bgp,
1490 const ArrayOfGridPos& pgp,
1491 const ArrayOfGridPos& rgp,
1492 const ArrayOfGridPos& cgp) {
1493 Index n = cgp.nelem();
1494 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1495 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1496 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1497 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1498 ARTS_ASSERT(is_size(itw, n, 32)); // We must store 32 interpolation
1499 // weights for each position.
1500
1501 // We have to loop all the points in the sequence:
1502 for (Index i = 0; i < n; ++i) {
1503 // Current grid positions:
1504 const GridPos& ts = sgp[i];
1505 const GridPos& tb = bgp[i];
1506 const GridPos& tp = pgp[i];
1507 const GridPos& tr = rgp[i];
1508 const GridPos& tc = cgp[i];
1509
1510 Index iti = 0;
1511 LOOPIT(s)
1512 LOOPIT(b)
1513 LOOPIT(p)
1514 LOOPIT(r)
1515 LOOPIT(c) {
1516 itw.get(i, iti) = (*s) * (*b) * (*p) * (*r) * (*c);
1517 ++iti;
1518 }
1519 }
1520}
1521
1523
1549 const ArrayOfGridPos& vgp,
1550 const ArrayOfGridPos& sgp,
1551 const ArrayOfGridPos& bgp,
1552 const ArrayOfGridPos& pgp,
1553 const ArrayOfGridPos& rgp,
1554 const ArrayOfGridPos& cgp) {
1555 Index n = cgp.nelem();
1556 ARTS_ASSERT(is_size(vgp, n)); // vgp must have same size as cgp.
1557 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1558 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1559 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1560 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1561 ARTS_ASSERT(is_size(itw, n, 64)); // We must store 64 interpolation
1562 // weights for each position.
1563
1564 // We have to loop all the points in the sequence:
1565 for (Index i = 0; i < n; ++i) {
1566 // Current grid positions:
1567 const GridPos& tv = vgp[i];
1568 const GridPos& ts = sgp[i];
1569 const GridPos& tb = bgp[i];
1570 const GridPos& tp = pgp[i];
1571 const GridPos& tr = rgp[i];
1572 const GridPos& tc = cgp[i];
1573
1574 Index iti = 0;
1575 LOOPIT(v)
1576 LOOPIT(s)
1577 LOOPIT(b)
1578 LOOPIT(p)
1579 LOOPIT(r)
1580 LOOPIT(c) {
1581 itw.get(i, iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
1582 ++iti;
1583 }
1584 }
1585}
1586
1588
1605 ConstMatrixView itw,
1607 const ArrayOfGridPos& cgp) {
1608 Index n = cgp.nelem();
1609 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1610 ARTS_ASSERT(is_size(itw, n, 2)); // We need 2 interpolation
1611 // weights for each position.
1612
1613 // Check that interpolation weights are valid. The sum of all
1614 // weights (last dimension) must always be approximately one. We
1615 // only check the first element.
1617 is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1618
1619 // We have to loop all the points in the sequence:
1620 for (Index i = 0; i < n; ++i) {
1621 // Current grid positions:
1622 const GridPos& tc = cgp[i];
1623
1624 // Get handle to current element of output vector and initialize
1625 // it to zero:
1626 Numeric& tia = ia[i];
1627 tia = 0;
1628
1629 Index iti = 0;
1630 for (Index c = 0; c < 2; ++c) {
1631 ARTS_ASSERT(tc.idx + c < a.nelem()); // Temporary !?
1632 tia += a.get(tc.idx + c) * itw.get(i, iti);
1633 ++iti;
1634 }
1635 }
1636}
1637
1639
1661 ConstMatrixView itw,
1663 const ArrayOfGridPos& rgp,
1664 const ArrayOfGridPos& cgp) {
1665 Index n = cgp.nelem();
1666 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1667 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1668 ARTS_ASSERT(is_size(itw, n, 4)); // We need 4 interpolation
1669 // weights for each position.
1670
1671 // Check that interpolation weights are valid. The sum of all
1672 // weights (last dimension) must always be approximately one. We
1673 // only check the first element.
1675 is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1676
1677 // We have to loop all the points in the sequence:
1678 for (Index i = 0; i < n; ++i) {
1679 // Current grid positions:
1680 const GridPos& tr = rgp[i];
1681 const GridPos& tc = cgp[i];
1682
1683 // Get handle to current element of output vector and initialize
1684 // it to zero:
1685 Numeric& tia = ia[i];
1686 tia = 0;
1687
1688 Index iti = 0;
1689 for (Index r = 0; r < 2; ++r)
1690 for (Index c = 0; c < 2; ++c) {
1691 tia += a.get(tr.idx + r, tc.idx + c) * itw.get(i, iti);
1692 ++iti;
1693 }
1694 }
1695}
1696
1698
1721 ConstMatrixView itw,
1723 const ArrayOfGridPos& pgp,
1724 const ArrayOfGridPos& rgp,
1725 const ArrayOfGridPos& cgp) {
1726 Index n = cgp.nelem();
1727 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1728 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1729 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1730 ARTS_ASSERT(is_size(itw, n, 8)); // We need 8 interpolation
1731 // weights for each position.
1732
1733 // Check that interpolation weights are valid. The sum of all
1734 // weights (last dimension) must always be approximately one. We
1735 // only check the first element.
1737 is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1738
1739 // We have to loop all the points in the sequence:
1740 for (Index i = 0; i < n; ++i) {
1741 // Current grid positions:
1742 const GridPos& tp = pgp[i];
1743 const GridPos& tr = rgp[i];
1744 const GridPos& tc = cgp[i];
1745
1746 // Get handle to current element of output vector and initialize
1747 // it to zero:
1748 Numeric& tia = ia[i];
1749 tia = 0;
1750
1751 Index iti = 0;
1752 for (Index p = 0; p < 2; ++p)
1753 for (Index r = 0; r < 2; ++r)
1754 for (Index c = 0; c < 2; ++c) {
1755 tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) * itw.get(i, iti);
1756 ++iti;
1757 }
1758 }
1759}
1760
1762
1786 ConstMatrixView itw,
1788 const ArrayOfGridPos& bgp,
1789 const ArrayOfGridPos& pgp,
1790 const ArrayOfGridPos& rgp,
1791 const ArrayOfGridPos& cgp) {
1792 Index n = cgp.nelem();
1793 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1794 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1795 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1796 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1797 ARTS_ASSERT(is_size(itw, n, 16)); // We need 16 interpolation
1798 // weights for each position.
1799
1800 // Check that interpolation weights are valid. The sum of all
1801 // weights (last dimension) must always be approximately one. We
1802 // only check the first element.
1804 is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1805
1806 // We have to loop all the points in the sequence:
1807 for (Index i = 0; i < n; ++i) {
1808 // Current grid positions:
1809 const GridPos& tb = bgp[i];
1810 const GridPos& tp = pgp[i];
1811 const GridPos& tr = rgp[i];
1812 const GridPos& tc = cgp[i];
1813
1814 // Get handle to current element of output vector and initialize
1815 // it to zero:
1816 Numeric& tia = ia[i];
1817 tia = 0;
1818
1819 Index iti = 0;
1820 for (Index b = 0; b < 2; ++b)
1821 for (Index p = 0; p < 2; ++p)
1822 for (Index r = 0; r < 2; ++r)
1823 for (Index c = 0; c < 2; ++c) {
1824 tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
1825 itw.get(i, iti);
1826 ++iti;
1827 }
1828 }
1829}
1830
1832
1857 ConstMatrixView itw,
1859 const ArrayOfGridPos& sgp,
1860 const ArrayOfGridPos& bgp,
1861 const ArrayOfGridPos& pgp,
1862 const ArrayOfGridPos& rgp,
1863 const ArrayOfGridPos& cgp) {
1864 Index n = cgp.nelem();
1865 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1866 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1867 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1868 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1869 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1870 ARTS_ASSERT(is_size(itw, n, 32)); // We need 32 interpolation
1871 // weights for each position.
1872
1873 // Check that interpolation weights are valid. The sum of all
1874 // weights (last dimension) must always be approximately one. We
1875 // only check the first element.
1877 is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1878
1879 // We have to loop all the points in the sequence:
1880 for (Index i = 0; i < n; ++i) {
1881 // Current grid positions:
1882 const GridPos& ts = sgp[i];
1883 const GridPos& tb = bgp[i];
1884 const GridPos& tp = pgp[i];
1885 const GridPos& tr = rgp[i];
1886 const GridPos& tc = cgp[i];
1887
1888 // Get handle to current element of output vector and initialize
1889 // it to zero:
1890 Numeric& tia = ia[i];
1891 tia = 0;
1892
1893 Index iti = 0;
1894 for (Index s = 0; s < 2; ++s)
1895 for (Index b = 0; b < 2; ++b)
1896 for (Index p = 0; p < 2; ++p)
1897 for (Index r = 0; r < 2; ++r)
1898 for (Index c = 0; c < 2; ++c) {
1899 tia += a.get(ts.idx + s,
1900 tb.idx + b,
1901 tp.idx + p,
1902 tr.idx + r,
1903 tc.idx + c) *
1904 itw.get(i, iti);
1905 ++iti;
1906 }
1907 }
1908}
1909
1911
1937 ConstMatrixView itw,
1939 const ArrayOfGridPos& vgp,
1940 const ArrayOfGridPos& sgp,
1941 const ArrayOfGridPos& bgp,
1942 const ArrayOfGridPos& pgp,
1943 const ArrayOfGridPos& rgp,
1944 const ArrayOfGridPos& cgp) {
1945 Index n = cgp.nelem();
1946 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1947 ARTS_ASSERT(is_size(vgp, n)); // vgp must have same size as cgp.
1948 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1949 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1950 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1951 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1952 ARTS_ASSERT(is_size(itw, n, 64)); // We need 64 interpolation
1953 // weights for each position.
1954
1955 // Check that interpolation weights are valid. The sum of all
1956 // weights (last dimension) must always be approximately one. We
1957 // only check the first element.
1959 is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1960
1961 // We have to loop all the points in the sequence:
1962 for (Index i = 0; i < n; ++i) {
1963 // Current grid positions:
1964 const GridPos& tv = vgp[i];
1965 const GridPos& ts = sgp[i];
1966 const GridPos& tb = bgp[i];
1967 const GridPos& tp = pgp[i];
1968 const GridPos& tr = rgp[i];
1969 const GridPos& tc = cgp[i];
1970
1971 // Get handle to current element of output vector and initialize
1972 // it to zero:
1973 Numeric& tia = ia[i];
1974 tia = 0;
1975
1976 Index iti = 0;
1977 for (Index v = 0; v < 2; ++v)
1978 for (Index s = 0; s < 2; ++s)
1979 for (Index b = 0; b < 2; ++b)
1980 for (Index p = 0; p < 2; ++p)
1981 for (Index r = 0; r < 2; ++r)
1982 for (Index c = 0; c < 2; ++c) {
1983 tia += a.get(tv.idx + v,
1984 ts.idx + s,
1985 tb.idx + b,
1986 tp.idx + p,
1987 tr.idx + r,
1988 tc.idx + c) *
1989 itw.get(i, iti);
1990 ++iti;
1991 }
1992 }
1993}
1994
1996// Green interpolation
1998
2000
2022 const ArrayOfGridPos& rgp,
2023 const ArrayOfGridPos& cgp) {
2024 Index nr = rgp.nelem();
2025 Index nc = cgp.nelem();
2026 ARTS_ASSERT(is_size(itw, nr, nc, 4)); // We must store 4 interpolation
2027 // weights for each position.
2028
2029 // We have to loop all the points in the new grid:
2030 for (Index ir = 0; ir < nr; ++ir) {
2031 // Current grid position:
2032 const GridPos& tr = rgp[ir];
2033
2034 for (Index ic = 0; ic < nc; ++ic) {
2035 // Current grid position:
2036 const GridPos& tc = cgp[ic];
2037
2038 // Interpolation weights are stored in this order (l=lower
2039 // u=upper, r=row, c=column):
2040 // 1. l-r l-c
2041 // 2. l-r u-c
2042 // 3. u-r l-c
2043 // 4. u-r u-c
2044
2045 Index iti = 0;
2046
2047 LOOPIT(r)
2048 LOOPIT(c) {
2049 itw.get(ir, ic, iti) = (*r) * (*c);
2050 ++iti;
2051 }
2052 }
2053 }
2054}
2055
2057
2080 const ArrayOfGridPos& pgp,
2081 const ArrayOfGridPos& rgp,
2082 const ArrayOfGridPos& cgp) {
2083 Index np = pgp.nelem();
2084 Index nr = rgp.nelem();
2085 Index nc = cgp.nelem();
2086 // We must store 8 interpolation weights for each position:
2087 ARTS_ASSERT(is_size(itw, np, nr, nc, 8));
2088
2089 // We have to loop all the points in the new grid:
2090 for (Index ip = 0; ip < np; ++ip) {
2091 const GridPos& tp = pgp[ip];
2092 for (Index ir = 0; ir < nr; ++ir) {
2093 const GridPos& tr = rgp[ir];
2094 for (Index ic = 0; ic < nc; ++ic) {
2095 const GridPos& tc = cgp[ic];
2096
2097 Index iti = 0;
2098
2099 LOOPIT(p)
2100 LOOPIT(r)
2101 LOOPIT(c) {
2102 itw.get(ip, ir, ic, iti) = (*p) * (*r) * (*c);
2103 ++iti;
2104 }
2105 }
2106 }
2107 }
2108}
2109
2111
2135 const ArrayOfGridPos& bgp,
2136 const ArrayOfGridPos& pgp,
2137 const ArrayOfGridPos& rgp,
2138 const ArrayOfGridPos& cgp) {
2139 Index nb = bgp.nelem();
2140 Index np = pgp.nelem();
2141 Index nr = rgp.nelem();
2142 Index nc = cgp.nelem();
2143 // We must store 16 interpolation weights for each position:
2144 ARTS_ASSERT(is_size(itw, nb, np, nr, nc, 16));
2145
2146 // We have to loop all the points in the new grid:
2147 for (Index ib = 0; ib < nb; ++ib) {
2148 const GridPos& tb = bgp[ib];
2149 for (Index ip = 0; ip < np; ++ip) {
2150 const GridPos& tp = pgp[ip];
2151 for (Index ir = 0; ir < nr; ++ir) {
2152 const GridPos& tr = rgp[ir];
2153 for (Index ic = 0; ic < nc; ++ic) {
2154 const GridPos& tc = cgp[ic];
2155
2156 Index iti = 0;
2157
2158 LOOPIT(b)
2159 LOOPIT(p)
2160 LOOPIT(r)
2161 LOOPIT(c) {
2162 itw.get(ib, ip, ir, ic, iti) = (*b) * (*p) * (*r) * (*c);
2163 ++iti;
2164 }
2165 }
2166 }
2167 }
2168 }
2169}
2170
2172
2197 const ArrayOfGridPos& sgp,
2198 const ArrayOfGridPos& bgp,
2199 const ArrayOfGridPos& pgp,
2200 const ArrayOfGridPos& rgp,
2201 const ArrayOfGridPos& cgp) {
2202 Index ns = sgp.nelem();
2203 Index nb = bgp.nelem();
2204 Index np = pgp.nelem();
2205 Index nr = rgp.nelem();
2206 Index nc = cgp.nelem();
2207 // We must store 32 interpolation weights for each position:
2208 ARTS_ASSERT(is_size(itw, ns, nb, np, nr, nc, 32));
2209
2210 // We have to loop all the points in the new grid:
2211 for (Index is = 0; is < ns; ++is) {
2212 const GridPos& ts = sgp[is];
2213 for (Index ib = 0; ib < nb; ++ib) {
2214 const GridPos& tb = bgp[ib];
2215 for (Index ip = 0; ip < np; ++ip) {
2216 const GridPos& tp = pgp[ip];
2217 for (Index ir = 0; ir < nr; ++ir) {
2218 const GridPos& tr = rgp[ir];
2219 for (Index ic = 0; ic < nc; ++ic) {
2220 const GridPos& tc = cgp[ic];
2221
2222 Index iti = 0;
2223
2224 LOOPIT(s)
2225 LOOPIT(b)
2226 LOOPIT(p)
2227 LOOPIT(r)
2228 LOOPIT(c) {
2229 itw.get(is, ib, ip, ir, ic, iti) =
2230 (*s) * (*b) * (*p) * (*r) * (*c);
2231 ++iti;
2232 }
2233 }
2234 }
2235 }
2236 }
2237 }
2238}
2239
2241
2267 const ArrayOfGridPos& vgp,
2268 const ArrayOfGridPos& sgp,
2269 const ArrayOfGridPos& bgp,
2270 const ArrayOfGridPos& pgp,
2271 const ArrayOfGridPos& rgp,
2272 const ArrayOfGridPos& cgp) {
2273 Index nv = vgp.nelem();
2274 Index ns = sgp.nelem();
2275 Index nb = bgp.nelem();
2276 Index np = pgp.nelem();
2277 Index nr = rgp.nelem();
2278 Index nc = cgp.nelem();
2279 // We must store 64 interpolation weights for each position:
2280 ARTS_ASSERT(is_size(itw, nv, ns, nb, np, nr, nc, 64));
2281
2282 // We have to loop all the points in the new grid:
2283 for (Index iv = 0; iv < nv; ++iv) {
2284 const GridPos& tv = vgp[iv];
2285 for (Index is = 0; is < ns; ++is) {
2286 const GridPos& ts = sgp[is];
2287 for (Index ib = 0; ib < nb; ++ib) {
2288 const GridPos& tb = bgp[ib];
2289 for (Index ip = 0; ip < np; ++ip) {
2290 const GridPos& tp = pgp[ip];
2291 for (Index ir = 0; ir < nr; ++ir) {
2292 const GridPos& tr = rgp[ir];
2293 for (Index ic = 0; ic < nc; ++ic) {
2294 const GridPos& tc = cgp[ic];
2295
2296 Index iti = 0;
2297
2298 LOOPIT(v)
2299 LOOPIT(s)
2300 LOOPIT(b)
2301 LOOPIT(p)
2302 LOOPIT(r)
2303 LOOPIT(c) {
2304 itw.get(iv, is, ib, ip, ir, ic, iti) =
2305 (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
2306 ++iti;
2307 }
2308 }
2309 }
2310 }
2311 }
2312 }
2313 }
2314}
2315
2317
2339 ConstTensor3View itw,
2341 const ArrayOfGridPos& rgp,
2342 const ArrayOfGridPos& cgp) {
2343 Index nr = rgp.nelem();
2344 Index nc = cgp.nelem();
2345 ARTS_ASSERT(is_size(ia, nr, nc));
2346 ARTS_ASSERT(is_size(itw, nr, nc, 4)); // We need 4 interpolation
2347 // weights for each position.
2348
2349 // Check that interpolation weights are valid. The sum of all
2350 // weights (last dimension) must always be approximately one. We
2351 // only check the first element.
2353 itw(0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2354
2355 // We have to loop all the points in the new grid:
2356 for (Index ir = 0; ir < nr; ++ir) {
2357 // Current grid position:
2358 const GridPos& tr = rgp[ir];
2359
2360 for (Index ic = 0; ic < nc; ++ic) {
2361 // Current grid position:
2362 const GridPos& tc = cgp[ic];
2363
2364 // Get handle to current element of output tensor and initialize
2365 // it to zero:
2366 Numeric& tia = ia(ir, ic);
2367 tia = 0;
2368
2369 Index iti = 0;
2370 for (Index r = 0; r < 2; ++r)
2371 for (Index c = 0; c < 2; ++c) {
2372 ARTS_ASSERT(tr.idx + r < a.nrows()); // Temporary !?
2373 ARTS_ASSERT(tc.idx + c < a.ncols()); // Temporary !?
2374 tia += a.get(tr.idx + r, tc.idx + c) * itw.get(ir, ic, iti);
2375 ++iti;
2376 }
2377 }
2378 }
2379}
2380
2382
2405 ConstTensor4View itw,
2407 const ArrayOfGridPos& pgp,
2408 const ArrayOfGridPos& rgp,
2409 const ArrayOfGridPos& cgp) {
2410 Index np = pgp.nelem();
2411 Index nr = rgp.nelem();
2412 Index nc = cgp.nelem();
2413 ARTS_ASSERT(is_size(ia, np, nr, nc));
2414 ARTS_ASSERT(is_size(itw, np, nr, nc, 8));
2415
2416 // Check that interpolation weights are valid. The sum of all
2417 // weights (last dimension) must always be approximately one. We
2418 // only check the first element.
2420 itw(0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2421
2422 // We have to loop all the points in the new grid:
2423 for (Index ip = 0; ip < np; ++ip) {
2424 const GridPos& tp = pgp[ip];
2425 for (Index ir = 0; ir < nr; ++ir) {
2426 const GridPos& tr = rgp[ir];
2427 for (Index ic = 0; ic < nc; ++ic) {
2428 // Current grid position:
2429 const GridPos& tc = cgp[ic];
2430
2431 // Get handle to current element of output tensor and
2432 // initialize it to zero:
2433 Numeric& tia = ia(ip, ir, ic);
2434 tia = 0;
2435
2436 Index iti = 0;
2437 for (Index p = 0; p < 2; ++p)
2438 for (Index r = 0; r < 2; ++r)
2439 for (Index c = 0; c < 2; ++c) {
2440 ARTS_ASSERT(tp.idx + p < a.npages()); // Temporary !?
2441 ARTS_ASSERT(tr.idx + r < a.nrows()); // Temporary !?
2442 ARTS_ASSERT(tc.idx + c < a.ncols()); // Temporary !?
2443 tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) *
2444 itw.get(ip, ir, ic, iti);
2445 ++iti;
2446 }
2447 }
2448 }
2449 }
2450}
2451
2453
2477 ConstTensor5View itw,
2479 const ArrayOfGridPos& bgp,
2480 const ArrayOfGridPos& pgp,
2481 const ArrayOfGridPos& rgp,
2482 const ArrayOfGridPos& cgp) {
2483 Index nb = bgp.nelem();
2484 Index np = pgp.nelem();
2485 Index nr = rgp.nelem();
2486 Index nc = cgp.nelem();
2487 ARTS_ASSERT(is_size(ia, nb, np, nr, nc));
2488 ARTS_ASSERT(is_size(itw, nb, np, nr, nc, 16));
2489
2490 // Check that interpolation weights are valid. The sum of all
2491 // weights (last dimension) must always be approximately one. We
2492 // only check the first element.
2494 itw(0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2495
2496 // We have to loop all the points in the new grid:
2497 for (Index ib = 0; ib < nb; ++ib) {
2498 const GridPos& tb = bgp[ib];
2499 for (Index ip = 0; ip < np; ++ip) {
2500 const GridPos& tp = pgp[ip];
2501 for (Index ir = 0; ir < nr; ++ir) {
2502 const GridPos& tr = rgp[ir];
2503 for (Index ic = 0; ic < nc; ++ic) {
2504 // Current grid position:
2505 const GridPos& tc = cgp[ic];
2506
2507 // Get handle to current element of output tensor and
2508 // initialize it to zero:
2509 Numeric& tia = ia(ib, ip, ir, ic);
2510 tia = 0;
2511
2512 Index iti = 0;
2513 for (Index b = 0; b < 2; ++b)
2514 for (Index p = 0; p < 2; ++p)
2515 for (Index r = 0; r < 2; ++r)
2516 for (Index c = 0; c < 2; ++c) {
2517 tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
2518 itw.get(ib, ip, ir, ic, iti);
2519 ++iti;
2520 }
2521 }
2522 }
2523 }
2524 }
2525}
2526
2528
2553 ConstTensor6View itw,
2555 const ArrayOfGridPos& sgp,
2556 const ArrayOfGridPos& bgp,
2557 const ArrayOfGridPos& pgp,
2558 const ArrayOfGridPos& rgp,
2559 const ArrayOfGridPos& cgp) {
2560 Index ns = sgp.nelem();
2561 Index nb = bgp.nelem();
2562 Index np = pgp.nelem();
2563 Index nr = rgp.nelem();
2564 Index nc = cgp.nelem();
2565 ARTS_ASSERT(is_size(ia, ns, nb, np, nr, nc));
2566 ARTS_ASSERT(is_size(itw, ns, nb, np, nr, nc, 32));
2567
2568 // Check that interpolation weights are valid. The sum of all
2569 // weights (last dimension) must always be approximately one. We
2570 // only check the first element.
2572 itw(0, 0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2573
2574 // We have to loop all the points in the new grid:
2575 for (Index is = 0; is < ns; ++is) {
2576 const GridPos& ts = sgp[is];
2577 for (Index ib = 0; ib < nb; ++ib) {
2578 const GridPos& tb = bgp[ib];
2579 for (Index ip = 0; ip < np; ++ip) {
2580 const GridPos& tp = pgp[ip];
2581 for (Index ir = 0; ir < nr; ++ir) {
2582 const GridPos& tr = rgp[ir];
2583 for (Index ic = 0; ic < nc; ++ic) {
2584 // Current grid position:
2585 const GridPos& tc = cgp[ic];
2586
2587 // Get handle to current element of output tensor and
2588 // initialize it to zero:
2589 Numeric& tia = ia(is, ib, ip, ir, ic);
2590 tia = 0;
2591
2592 Index iti = 0;
2593 for (Index s = 0; s < 2; ++s)
2594 for (Index b = 0; b < 2; ++b)
2595 for (Index p = 0; p < 2; ++p)
2596 for (Index r = 0; r < 2; ++r)
2597 for (Index c = 0; c < 2; ++c) {
2598 tia += a.get(ts.idx + s,
2599 tb.idx + b,
2600 tp.idx + p,
2601 tr.idx + r,
2602 tc.idx + c) *
2603 itw.get(is, ib, ip, ir, ic, iti);
2604 ++iti;
2605 }
2606 }
2607 }
2608 }
2609 }
2610 }
2611}
2612
2614
2640 ConstTensor7View itw,
2642 const ArrayOfGridPos& vgp,
2643 const ArrayOfGridPos& sgp,
2644 const ArrayOfGridPos& bgp,
2645 const ArrayOfGridPos& pgp,
2646 const ArrayOfGridPos& rgp,
2647 const ArrayOfGridPos& cgp) {
2648 Index nv = vgp.nelem();
2649 Index ns = sgp.nelem();
2650 Index nb = bgp.nelem();
2651 Index np = pgp.nelem();
2652 Index nr = rgp.nelem();
2653 Index nc = cgp.nelem();
2654 ARTS_ASSERT(is_size(ia, nv, ns, nb, np, nr, nc));
2655 ARTS_ASSERT(is_size(itw, nv, ns, nb, np, nr, nc, 64));
2656
2657 // Check that interpolation weights are valid. The sum of all
2658 // weights (last dimension) must always be approximately one. We
2659 // only check the first element.
2661 itw(0, 0, 0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2662
2663 // We have to loop all the points in the new grid:
2664 for (Index iv = 0; iv < nv; ++iv) {
2665 const GridPos& tv = vgp[iv];
2666 for (Index is = 0; is < ns; ++is) {
2667 const GridPos& ts = sgp[is];
2668 for (Index ib = 0; ib < nb; ++ib) {
2669 const GridPos& tb = bgp[ib];
2670 for (Index ip = 0; ip < np; ++ip) {
2671 const GridPos& tp = pgp[ip];
2672 for (Index ir = 0; ir < nr; ++ir) {
2673 const GridPos& tr = rgp[ir];
2674 for (Index ic = 0; ic < nc; ++ic) {
2675 // Current grid position:
2676 const GridPos& tc = cgp[ic];
2677
2678 // Get handle to current element of output tensor and
2679 // initialize it to zero:
2680 Numeric& tia = ia(iv, is, ib, ip, ir, ic);
2681 tia = 0;
2682
2683 Index iti = 0;
2684 for (Index v = 0; v < 2; ++v)
2685 for (Index s = 0; s < 2; ++s)
2686 for (Index b = 0; b < 2; ++b)
2687 for (Index p = 0; p < 2; ++p)
2688 for (Index r = 0; r < 2; ++r)
2689 for (Index c = 0; c < 2; ++c) {
2690 tia += a.get(tv.idx + v,
2691 ts.idx + s,
2692 tb.idx + b,
2693 tp.idx + p,
2694 tr.idx + r,
2695 tc.idx + c) *
2696 itw.get(iv, is, ib, ip, ir, ic, iti);
2697 ++iti;
2698 }
2699 }
2700 }
2701 }
2702 }
2703 }
2704 }
2705}
2706
2707
2709
2728void polint(Numeric& y_int,
2729 Numeric& dy_int,
2730 ConstVectorView xa,
2731 ConstVectorView ya,
2732 const Index& n,
2733 const Numeric& x) {
2734 Index ns = 1;
2735 Numeric den, dif, dift, ho, hp, w;
2736
2737 dif = abs(x - xa[0]);
2738
2739 Vector c(n);
2740 Vector d(n);
2741
2742 // Find the index of the closest table entry
2743 for (Index i = 0; i < n; i++) {
2744 if ((dift = abs(x - xa[i])) < dif) {
2745 ns = i;
2746 dif = dift;
2747 }
2748 // Initialize c and d
2749 c[i] = ya[i];
2750 d[i] = ya[i];
2751 }
2752 // Initial approximation to y
2753 y_int = ya[ns--];
2754
2755 for (Index m = 1; m < n; m++) {
2756 for (Index i = 0; i < n - m; i++) {
2757 ho = xa[i] - x;
2758 hp = xa[i + m] - x;
2759 w = c[i + 1] - d[i];
2760 den = ho - hp;
2761 // This error occurs when two input xa's are identical.
2762 ARTS_ASSERT(den != 0.);
2763 den = w / den;
2764 d[i] = hp * den;
2765 c[i] = ho * den;
2766 }
2767 y_int += (dy_int = (2 * (ns + 1) < (n - m) ? c[ns + 1] : d[ns--]));
2768 }
2769}
2770
2772
2789 const Numeric& x_i,
2790 const GridPos& gp) {
2791 Index N_x = x.nelem();
2792
2793 ARTS_ASSERT(N_x == y.nelem());
2794 ARTS_ASSERT(N_x > 2);
2795
2796 Vector xa(4), ya(4);
2797 Numeric y_int;
2798 Numeric dy_int;
2799 y_int = 0.;
2800
2801 // 1 - polynomial interpolation (3 points) with grid position search
2802 // 2 - polynomial interpolation (3 points) without grid position search
2803 // 3 - polynomial interpolation (4 points)
2804
2805 Index interp_method = 1;
2806
2807 if (interp_method == 1) {
2808 // Pick out three points for interpolation
2809 if ((gp.fd[0] <= 0.5 && gp.idx > 0) || gp.idx == N_x - 2) {
2810 xa[0] = x[gp.idx - 1];
2811 xa[1] = x[gp.idx];
2812 xa[2] = x[gp.idx + 1];
2813
2814 ya[0] = y[gp.idx - 1];
2815 ya[1] = y[gp.idx];
2816 ya[2] = y[gp.idx + 1];
2817 }
2818
2819 else if ((gp.fd[0] > 0.5 && gp.idx < N_x - 2) || gp.idx == 0) {
2820 xa[0] = x[gp.idx];
2821 xa[1] = x[gp.idx + 1];
2822 xa[2] = x[gp.idx + 2];
2823
2824 ya[0] = y[gp.idx];
2825 ya[1] = y[gp.idx + 1];
2826 ya[2] = y[gp.idx + 2];
2827 }
2828
2829 else if (gp.idx == N_x - 1) {
2830 xa[0] = x[N_x - 2];
2831 xa[1] = x[N_x - 1];
2832 xa[2] = x[N_x];
2833
2834 ya[0] = y[N_x - 2];
2835 ya[1] = y[N_x - 1];
2836 ya[2] = y[N_x];
2837 } else {
2838 ARTS_ASSERT(false);
2839 arts_exit();
2840 }
2841
2842 polint(y_int, dy_int, xa, ya, 3, x_i);
2843
2844 }
2845
2846 else if (interp_method == 2) {
2847 if (gp.idx == 0) {
2848 xa[0] = x[gp.idx];
2849 xa[1] = x[gp.idx + 1];
2850 xa[2] = x[gp.idx + 2];
2851
2852 ya[0] = y[gp.idx];
2853 ya[1] = y[gp.idx + 1];
2854 ya[2] = y[gp.idx + 2];
2855 } else if (gp.idx == N_x - 1) {
2856 xa[0] = x[gp.idx - 2];
2857 xa[1] = x[gp.idx - 1];
2858 xa[2] = x[gp.idx];
2859
2860 ya[0] = y[gp.idx - 2];
2861 ya[1] = y[gp.idx - 1];
2862 ya[2] = y[gp.idx];
2863 } else {
2864 xa[0] = x[gp.idx - 1];
2865 xa[1] = x[gp.idx];
2866 xa[2] = x[gp.idx + 1];
2867
2868 ya[0] = y[gp.idx - 1];
2869 ya[1] = y[gp.idx];
2870 ya[2] = y[gp.idx + 1];
2871 }
2872
2873 // Polynominal interpolation, n = 3
2874 polint(y_int, dy_int, xa, ya, 3, x_i);
2875 }
2876
2877 else if (interp_method == 3) {
2878 // Take 4 points
2879 if (gp.idx == 0) {
2880 xa[0] = -x[gp.idx + 1];
2881 xa[1] = x[gp.idx + 0];
2882 xa[2] = x[gp.idx + 1];
2883 xa[3] = x[gp.idx + 2];
2884
2885 ya[0] = y[gp.idx + 1];
2886 ya[1] = y[gp.idx + 0];
2887 ya[2] = y[gp.idx + 1];
2888 ya[3] = y[gp.idx + 2];
2889 } else if (gp.idx == N_x - 1) {
2890 xa[0] = x[gp.idx - 1];
2891 xa[1] = x[gp.idx - 0];
2892 xa[2] = 2 * x[gp.idx] - x[gp.idx - 1];
2893 xa[3] = 2 * x[gp.idx] - x[gp.idx - 2];
2894
2895 ya[0] = y[gp.idx - 1];
2896 ya[1] = y[gp.idx - 0];
2897 ya[2] = y[gp.idx - 1];
2898 ya[3] = y[gp.idx - 2];
2899 } else if (gp.idx == N_x - 2) {
2900 xa[0] = x[gp.idx - 2];
2901 xa[1] = x[gp.idx - 1];
2902 xa[2] = x[gp.idx];
2903 xa[3] = x[gp.idx + 1];
2904
2905 ya[0] = y[gp.idx - 2];
2906 ya[1] = y[gp.idx - 1];
2907 ya[2] = y[gp.idx];
2908 ya[3] = y[gp.idx + 1];
2909 } else {
2910 xa[0] = x[gp.idx - 1];
2911 xa[1] = x[gp.idx];
2912 xa[2] = x[gp.idx + 1];
2913 xa[3] = x[gp.idx + 2];
2914
2915 ya[0] = y[gp.idx - 1];
2916 ya[1] = y[gp.idx];
2917 ya[2] = y[gp.idx + 1];
2918 ya[3] = y[gp.idx + 2];
2919 }
2920 // Polinominal interpolation, n = 4
2921 polint(y_int, dy_int, xa, ya, 4, x_i);
2922 }
2923
2924 return y_int;
2925}
This file contains the definition of Array.
void arts_exit(int status)
This is the exit function of ARTS.
Definition: arts.cc:42
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
A constant view of a Matrix.
Definition: matpackI.h:1014
Numeric get(Index r, Index c) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1041
A constant view of a Tensor3.
Definition: matpackIII.h:132
Numeric get(Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIII.h:180
A constant view of a Tensor4.
Definition: matpackIV.h:133
Numeric get(Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIV.h:220
A constant view of a Tensor5.
Definition: matpackV.h:143
Numeric get(Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackV.h:265
A constant view of a Tensor6.
Definition: matpackVI.h:149
Numeric get(Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVI.h:550
A constant view of a Tensor7.
Definition: matpackVII.h:147
Numeric get(Index l, Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVII.h:1211
A constant view of a Vector.
Definition: matpackI.h:489
Numeric get(Index n) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:530
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:52
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The MatrixView class.
Definition: matpackI.h:1125
Numeric & get(Index r, Index c) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1150
The range class.
Definition: matpackI.h:165
The Tensor3View class.
Definition: matpackIII.h:239
Numeric & get(Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackIII.h:275
The Tensor4View class.
Definition: matpackIV.h:284
Numeric & get(Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackIV.h:348
The Tensor5View class.
Definition: matpackV.h:333
Numeric & get(Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackV.h:431
The Tensor6View class.
Definition: matpackVI.h:621
Numeric & get(Index v, Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackVI.h:1015
The Tensor7View class.
Definition: matpackVII.h:1286
Numeric & get(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackVII.h:2308
The VectorView class.
Definition: matpackI.h:626
Numeric & get(Index n) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:654
The Vector class.
Definition: matpackI.h:876
#define DEBUG_ONLY(...)
Definition: debug.h:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define LOOPIT(x)
Macro for interpolation weight loops.
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
const Numeric FD_TOL
The maximum difference from 1 that we allow for a sum check.
ostream & operator<<(ostream &os, const GridPos &gp)
Output operator for GridPos.
void polint(Numeric &y_int, Numeric &dy_int, ConstVectorView xa, ConstVectorView ya, const Index &n, const Numeric &x)
Polynomial interpolation.
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
void gridpos_1to1(ArrayOfGridPos &gp, ConstVectorView grid)
gridpos_1to1
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Header file for interpolation.cc.
#define abs(x)
#define ns
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:351
Header file for logic.cc.
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
void sum(PropagationMatrix &pm, const ComplexVectorView &abs, const PolarizationVector &polvec, const bool do_phase) ARTS_NOEXCEPT
Sums the Zeeman components into a propagation matrix.
Definition: zeemandata.cc:375
#define d
#define v
#define w
#define a
#define c
#define b
Structure to store a grid position.
Definition: interpolation.h:73
Numeric fd[2]
Definition: interpolation.h:75
Index idx
Definition: interpolation.h:74