ARTS 2.5.11 (git: 6827797f)
interpolation.cc
Go to the documentation of this file.
1
23#include "interpolation.h"
24#include <cmath>
25#include <iostream>
26#include "array.h"
27#include "check_input.h"
28#include "logic.h"
29
30// File-global constants:
31
33
44DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6;)
45
47
56const Numeric FD_TOL = 1.5e-3;
57
59
68#define LOOPIT(x) for (const Numeric* x = &t##x.fd[1]; x >= &t##x.fd[0]; --x)
69
71
79ostream& operator<<(ostream& os, const GridPos& gp) {
80 os << gp.idx << " " << gp.fd[0] << " " << gp.fd[1] << "\n";
81 return os;
82}
83
85
140 ConstVectorView old_grid,
141 ConstVectorView new_grid,
142 const Numeric& extpolfac) {
143 const Index n_old = old_grid.nelem();
144 const Index n_new = new_grid.nelem();
145
146 // Assert that gp has the right size:
147 ARTS_ASSERT(is_size(gp, n_new));
148
149 // Assert, that the old grid has more than one element
150 ARTS_ASSERT(1 < n_old);
151
152 // This function hast two parts, depending on whether old_grid is
153 // sorted in ascending or descending order. Maybe that's not too
154 // elegant, but it's the most efficient, because in this way there
155 // is no additional runtime overhead to handle both cases.
156
157 // We use only the first two elements to decide how the grid is
158 // sorted. (The rest of the grid is checked later by an ARTS_ASSERT.)
159 // If both values are the same, we still assume the grid is
160 // ascending. However, this will lead to an assertion fail later on,
161 // because the grid has to be strictly sorted.
162 bool ascending = (old_grid[0] <= old_grid[1]);
163
164 if (ascending) {
165 // So, old_grid should always be sorted in strictly ascending order.
166 // (No duplicate values.)
167 // Assert that this is so. This may depend on user input, however,
168 // inside this elementary function is not the place to check for
169 // that. There must be runtime checks on higher levels to insure
170 // that all grids are sorted. The assertion here is just as a last
171 // safety check.
172 ARTS_ASSERT(is_increasing(old_grid));
173
174 // Limits of extrapolation.
175 const Numeric og_min =
176 old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
177 const Numeric og_max =
178 old_grid[n_old - 1] +
179 extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
180
181 // We will make no firm assumptions about the new grid. But the case
182 // that we have in mind is that it is either also sorted, or at
183 // least partially sorted, for example like this:
184 // 5 4 3 2 2.5 3 4
185 // This kind of sequence should be typical if we interpolate
186 // atmospheric fields along a limb propagation path.
187
188 // Let's get some idea where the first point in the new grid is,
189 // relative to the old grid. We use linear interpolation between the
190 // maximum and the minimum of the old grid for this. (Taking
191 // into account the small allowed extrapolation.)
192 Numeric frac = (new_grid[0] - og_min) / (og_max - og_min);
193
194 // We are not checking if the value of frac is reasonable,
195 // because there is another assertion below to catch the
196 // consequences.
197
198 // Initialize current_position:
199 Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
200
201 // The above statement should satisfy
202 // 0 <= current_position <= n_old-2
203 // Assert that this is indeed the case.
204 // cout << "frac = " << frac << "\n";
205 // cout << "current_position = " << current_position << "\n";
206 ARTS_ASSERT(0 <= current_position);
207 ARTS_ASSERT(current_position <= n_old - 2);
208
209 // The variables lower and upper are used to remember the value of
210 // the old grid at current_position and one above current_position:
211 Numeric lower = old_grid[current_position];
212 Numeric upper = old_grid[current_position + 1];
213
214 // Loop over all points in the new grid:
215 for (Index i_new = 0; i_new < n_new; ++i_new) {
216 // Get a reference to the current element of gp:
217 GridPos& tgp = gp[i_new];
218 // And on the current value of the new grid:
219 const Numeric tng = new_grid[i_new];
220
221 // Verify that the new grid is within the limits of
222 // extrapolation that we have defined by og_min and og_max:
223 ARTS_ASSERT(og_min <= tng);
224 ARTS_ASSERT(tng <= og_max);
225
226 // cout << "lower / tng / upper = "
227 // << lower << " / "
228 // << tng << " / "
229 // << upper << "\n";
230
231 // Is current_position too high?
232 // (The current_position>0 condition is there so that the position
233 // stays 0 for extrapolation.)
234 if (tng < lower && current_position > 0) {
235 do {
236 --current_position;
237 lower = old_grid[current_position];
238 } while (tng < lower && current_position > 0);
239
240 upper = old_grid[current_position + 1];
241
242 tgp.idx = current_position;
243 tgp.fd[0] = (tng - lower) / (upper - lower);
244 tgp.fd[1] = 1.0 - tgp.fd[0];
245 } else {
246 // Is it too low?
247 // (The current_position<n_old condition is there so
248 // that uppers stays n_old-1 for extrapolation.)
249 if (tng >= upper && current_position < n_old - 2) {
250 do {
251 ++current_position;
252 upper = old_grid[current_position + 1];
253 } while (tng >= upper && current_position < n_old - 2);
254
255 lower = old_grid[current_position];
256
257 tgp.idx = current_position;
258 tgp.fd[0] = (tng - lower) / (upper - lower);
259 tgp.fd[1] = 1.0 - tgp.fd[0];
260 } else {
261 // None of the other two conditions were true. That means:
262 // lower <= tng < upper. The current_position is still
263 // valid.
264
265 // SAB 2010-04-28: Note that if a new grid point is exactly on
266 // top of an old grid point, then you are now guaranteed to get
267 // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
268
269 tgp.idx = current_position;
270 tgp.fd[0] = (tng - lower) / (upper - lower);
271 tgp.fd[1] = 1.0 - tgp.fd[0];
272 }
273 }
274
275 // cout << tgp.idx << " " << tgp.fd[0] << " " << tgp.fd[1] << endl;
276
277 // Safety check to ensure the above promise:
278 ARTS_ASSERT(old_grid[tgp.idx] <= tng || tgp.idx == 0);
279 }
280 } else // if (ascending)
281 {
282 // Now we are in the "descending old grid" part. We do exactly
283 // the same as in the other part, just accounting for the
284 // different order of things. Comments here refer only to
285 // interesting differences from the ascending case. See that
286 // case for more general comments.
287
288 // This time ensure strictly descending order:
289 ARTS_ASSERT(is_decreasing(old_grid));
290
291 // The max is now the first point, the min the last point!
292 // I think the sign is right here...
293 const Numeric og_max =
294 old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
295 const Numeric og_min =
296 old_grid[n_old - 1] +
297 extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
298
299 // We have to take 1- here, because we are starting from the
300 // high end.
301 Numeric frac = 1 - (new_grid[0] - og_min) / (og_max - og_min);
302
303 // We are not checking if the value of frac is reasonable,
304 // because there is another assertion below to catch the
305 // consequences.
306
307 Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
308
309 // The above statement should satisfy
310 // 0 <= current_position <= n_old-2
311 // Assert that this is indeed the case.
312 // cout << "frac = " << frac << "\n";
313 // cout << "current_position = " << current_position << "\n";
314 ARTS_ASSERT(0 <= current_position);
315 ARTS_ASSERT(current_position <= n_old - 2);
316
317 // Note, that old_grid[lower] has a higher numerical value than
318 // old_grid[upper]!
319 Numeric lower = old_grid[current_position];
320 Numeric upper = old_grid[current_position + 1];
321
322 for (Index i_new = 0; i_new < n_new; ++i_new) {
323 GridPos& tgp = gp[i_new];
324 const Numeric tng = new_grid[i_new];
325
326 // Verify that the new grid is within the limits of
327 // extrapolation that we have defined by og_min and og_max:
328 ARTS_ASSERT(og_min <= tng);
329 ARTS_ASSERT(tng <= og_max);
330
331 // cout << "lower / tng / upper = "
332 // << lower << " / "
333 // << tng << " / "
334 // << upper << "\n";
335
336 // Is current_position too high? (Sign of comparison changed
337 // compared to ascending case!)
338 if (tng > lower && current_position > 0) {
339 do {
340 --current_position;
341 lower = old_grid[current_position];
342 } while (tng > lower && current_position > 0);
343
344 upper = old_grid[current_position + 1];
345
346 tgp.idx = current_position;
347 tgp.fd[0] = (tng - lower) / (upper - lower);
348 tgp.fd[1] = 1.0 - tgp.fd[0];
349 } else {
350 // Is it too low? (Sign of comparison changed
351 // compared to ascending case!)
352 if (tng <= upper && current_position < n_old - 2) {
353 do {
354 ++current_position;
355 upper = old_grid[current_position + 1];
356 } while (tng <= upper && current_position < n_old - 2);
357
358 lower = old_grid[current_position];
359
360 tgp.idx = current_position;
361 tgp.fd[0] = (tng - lower) / (upper - lower);
362 tgp.fd[1] = 1.0 - tgp.fd[0];
363 } else {
364 // None of the other two conditions were true. That means:
365 // lower >= tng > upper. The current_position is still
366 // valid. (Note that upper and lower have switched
367 // place compared to the ascending case.)
368
369 // SAB 2010-04-28: Note that if a new grid point is exactly on
370 // top of an old grid point, then you are now guaranteed to get
371 // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
372
373 tgp.idx = current_position;
374 tgp.fd[0] = (tng - lower) / (upper - lower);
375 tgp.fd[1] = 1.0 - tgp.fd[0];
376 }
377 }
378
379 // Safety check to ensure the above promise:
380 ARTS_ASSERT(old_grid[tgp.idx] >= tng || tgp.idx == 0);
381 }
382 }
383}
384
386
408 ConstVectorView old_grid,
409 const Numeric& new_grid,
410 const Numeric& extpolfac) {
411 ArrayOfGridPos agp(1);
412 Vector v(1, new_grid);
413 gridpos(agp, old_grid, v, extpolfac);
414 gridpos_copy(gp, agp[0]);
415}
416
418
432void gridpos_1to1(ArrayOfGridPos& gp, ConstVectorView grid) {
433 const Index n = grid.nelem();
434 gp.resize(n);
435
436 for (Index i = 0; i < n - 1; i++) {
437 gp[i].idx = i;
438 gp[i].fd[0] = 0;
439 gp[i].fd[1] = 1;
440 }
441 gp[n - 1].idx = n - 2;
442 gp[n - 1].fd[0] = 1;
443 gp[n - 1].fd[1] = 0;
444}
445
447
456void gridpos_copy(GridPos& gp_new, const GridPos& gp_old) {
457 gp_new.idx = gp_old.idx;
458 gp_new.fd[0] = gp_old.fd[0];
459 gp_new.fd[1] = gp_old.fd[1];
460}
461
463
475Numeric fractional_gp(const GridPos& gp) {
476 return (Numeric(gp.idx) + gp.fd[0]);
477}
478
480
493 // Catch values that "must" be wrong
494 ARTS_ASSERT(gp.fd[0] > -FD_TOL);
495 ARTS_ASSERT(gp.fd[0] < 1.0 + FD_TOL);
496 ARTS_ASSERT(gp.fd[1] > -FD_TOL);
497 ARTS_ASSERT(gp.fd[1] < 1.0 + FD_TOL);
498
499 if (gp.fd[0] < 0.0) {
500 gp.fd[0] = 0.0;
501 gp.fd[1] = 1.0;
502 } else if (gp.fd[0] > 1.0) {
503 gp.fd[0] = 1.0;
504 gp.fd[1] = 0.0;
505 }
506
507 if (gp.fd[1] < 0.0) {
508 gp.fd[1] = 0.0;
509 gp.fd[0] = 1.0;
510 } else if (gp.fd[1] > 1.0) {
511 gp.fd[1] = 1.0;
512 gp.fd[0] = 0.0;
513 }
514}
515
517
538void gridpos_force_end_fd(GridPos& gp, const Index& n) {
539 ARTS_ASSERT(gp.idx >= 0);
540
541 // If fd=1, shift to grid index above
542 if (gp.fd[0] > 0.5) {
543 gp.idx += 1;
544 }
545 gp.fd[0] = 0;
546 gp.fd[1] = 1;
547
548 ARTS_ASSERT(gp.idx < n);
549
550 // End of complete grid range must be handled separately
551 if (gp.idx == n - 1) {
552 gp.idx -= 1;
553 gp.fd[0] = 1;
554 gp.fd[1] = 0;
555 }
556}
557
559
572void gridpos_upperend_check(GridPos& gp, const Index& ie) {
573 if (gp.idx == ie) {
574 ARTS_ASSERT(gp.fd[0] < 0.005); // To capture obviously bad cases
575 gp.idx -= 1;
576 gp.fd[0] = 1.0;
577 gp.fd[1] = 0.0;
578 }
579}
580
582
595void gridpos_upperend_check(ArrayOfGridPos& gp, const Index& ie) {
596 for (Index i = 0; i < gp.nelem(); i++) {
597 if (gp[i].idx == ie) {
598 ARTS_ASSERT(gp[i].fd[0] < 0.005); // To capture obviously bad cases
599 gp[i].idx -= 1;
600 gp[i].fd[0] = 1.0;
601 gp[i].fd[1] = 0.0;
602 }
603 }
604}
605
607
620 for (Index i = 0; i < gp.nelem(); i++) {
621 gp[i].idx = 0;
622 gp[i].fd[0] = 0;
623 gp[i].fd[1] = 1;
624 }
625}
626
628
641 const Index& i,
642 const bool& strict) {
643 if (strict) {
644 // Assume that gridpos_force_end_fd has been used. The expression 0==0
645 // should be safer than 1==1.
646 if (gp.idx == i && gp.fd[0] == 0) {
647 return true;
648 } else if (gp.idx == i - 1 && gp.fd[1] == 0) {
649 return true;
650 }
651 } else {
652 if (gp.idx == i && gp.fd[0] < FD_TOL) {
653 return true;
654 } else if (gp.idx == i - 1 && gp.fd[1] < FD_TOL) {
655 return true;
656 }
657 }
658 return false;
659}
660
662
680Index gridpos2gridrange(const GridPos& gp, const bool& upwards) {
681 ARTS_ASSERT(gp.fd[0] >= 0);
682 ARTS_ASSERT(gp.fd[1] >= 0);
683
684 // Not at a grid point
685 if (gp.fd[0] > 0 && gp.fd[1] > 0) {
686 return gp.idx;
687 }
688
689 // Fractional distance 0
690 else if (gp.fd[0] == 0) {
691 if (upwards)
692 return gp.idx;
693 else
694 return gp.idx - 1;
695 }
696
697 // Fractional distance 1
698 else {
699 if (upwards)
700 return gp.idx + 1;
701 else
702 return gp.idx;
703 }
704}
705
707// Red Interpolation
709
711
724void interpweights(VectorView itw, const GridPos& tc) {
725 ARTS_ASSERT(is_size(itw, 2)); // We must store 2 interpolation
726 // weights.
727
728 // Interpolation weights are stored in this order (l=lower
729 // u=upper, c=column):
730 // 1. l-c
731 // 2. u-c
732
733 Index iti = 0;
734
735 // We could use a straight-forward for loop here:
736 //
737 // for ( Index c=1; c>=0; --c )
738 // {
739 // ti[iti] = tc.fd[c];
740 // ++iti;
741 // }
742 //
743 // However, it is slightly faster to use pointers (I tried it,
744 // the speed gain is about 20%). So we should write something
745 // like:
746 //
747 // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
748 // {
749 // ti[iti] = *c;
750 // ++iti;
751 // }
752 //
753 // For higher dimensions we have to nest these loops. To avoid
754 // typos and safe typing, I use the LOOPIT macro, which expands
755 // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
756 // COMMAND!
757
758 LOOPIT(c) {
759 itw[iti] = *c;
760 ++iti;
761 }
762}
763
765
779void interpweights(VectorView itw, const GridPos& tr, const GridPos& tc) {
780 ARTS_ASSERT(is_size(itw, 4)); // We must store 4 interpolation
781 // weights.
782 Index iti = 0;
783
784 LOOPIT(r)
785 LOOPIT(c) {
786 itw[iti] = (*r) * (*c);
787 ++iti;
788 }
789}
790
792
807void interpweights(VectorView itw,
808 const GridPos& tp,
809 const GridPos& tr,
810 const GridPos& tc) {
811 ARTS_ASSERT(is_size(itw, 8)); // We must store 8 interpolation
812 // weights.
813 Index iti = 0;
814
815 LOOPIT(p)
816 LOOPIT(r)
817 LOOPIT(c) {
818 itw[iti] = (*p) * (*r) * (*c);
819 ++iti;
820 }
821}
822
824
840void interpweights(VectorView itw,
841 const GridPos& tb,
842 const GridPos& tp,
843 const GridPos& tr,
844 const GridPos& tc) {
845 ARTS_ASSERT(is_size(itw, 16)); // We must store 16 interpolation
846 // weights.
847 Index iti = 0;
848
849 LOOPIT(b)
850 LOOPIT(p)
851 LOOPIT(r)
852 LOOPIT(c) {
853 itw[iti] = (*b) * (*p) * (*r) * (*c);
854 ++iti;
855 }
856}
857
859
876void interpweights(VectorView itw,
877 const GridPos& ts,
878 const GridPos& tb,
879 const GridPos& tp,
880 const GridPos& tr,
881 const GridPos& tc) {
882 ARTS_ASSERT(is_size(itw, 32)); // We must store 32 interpolation
883 // weights.
884 Index iti = 0;
885
886 LOOPIT(s)
887 LOOPIT(b)
888 LOOPIT(p)
889 LOOPIT(r)
890 LOOPIT(c) {
891 itw[iti] = (*s) * (*b) * (*p) * (*r) * (*c);
892 ++iti;
893 }
894}
895
897
915void interpweights(VectorView itw,
916 const GridPos& tv,
917 const GridPos& ts,
918 const GridPos& tb,
919 const GridPos& tp,
920 const GridPos& tr,
921 const GridPos& tc) {
922 ARTS_ASSERT(is_size(itw, 64)); // We must store 64 interpolation
923 // weights.
924 Index iti = 0;
925
926 LOOPIT(v)
927 LOOPIT(s)
928 LOOPIT(b)
929 LOOPIT(p)
930 LOOPIT(r)
931 LOOPIT(c) {
932 itw[iti] = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
933 ++iti;
934 }
935}
936
938
953Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos& tc) {
954 ARTS_ASSERT(is_size(itw, 2)); // We need 2 interpolation
955 // weights.
956
957 // Check that interpolation weights are valid. The sum of all
958 // weights (last dimension) must always be approximately one.
959 ARTS_ASSERT(is_same_within_epsilon(sum(itw), 1, sum_check_epsilon));
960
961 // To store interpolated value:
962 Numeric tia = 0;
963
964 Index iti = 0;
965 for (Index c = 0; c < 2; ++c) {
966 tia += a[tc.idx + c] * itw[iti];
967 ++iti;
968 }
969
970 return tia;
971}
972
974
990Numeric interp(ConstVectorView itw,
991 ConstMatrixView a,
992 const GridPos& tr,
993 const GridPos& tc) {
994 ARTS_ASSERT(is_size(itw, 4)); // We need 4 interpolation
995 // weights.
996
997 // Check that interpolation weights are valid. The sum of all
998 // weights (last dimension) must always be approximately one.
999 ARTS_ASSERT(is_same_within_epsilon(sum(itw), 1, sum_check_epsilon));
1000
1001 // To store interpolated value:
1002 Numeric tia = 0;
1003
1004 Index iti = 0;
1005 for (Index r = 0; r < 2; ++r)
1006 for (Index c = 0; c < 2; ++c) {
1007 tia += a(tr.idx + r, tc.idx + c) * itw[iti];
1008 ++iti;
1009 }
1010
1011 return tia;
1012}
1013
1015
1032Numeric interp(ConstVectorView itw,
1033 ConstTensor3View a,
1034 const GridPos& tp,
1035 const GridPos& tr,
1036 const GridPos& tc) {
1037 ARTS_ASSERT(is_size(itw, 8)); // We need 8 interpolation
1038 // weights.
1039
1040 // Check that interpolation weights are valid. The sum of all
1041 // weights (last dimension) must always be approximately one.
1042 ARTS_ASSERT(is_same_within_epsilon(sum(itw), 1, sum_check_epsilon));
1043
1044 // To store interpolated value:
1045 Numeric tia = 0;
1046
1047 Index iti = 0;
1048 for (Index p = 0; p < 2; ++p)
1049 for (Index r = 0; r < 2; ++r)
1050 for (Index c = 0; c < 2; ++c) {
1051 tia += a(tp.idx + p, tr.idx + r, tc.idx + c) * itw[iti];
1052 ++iti;
1053 }
1054
1055 return tia;
1056}
1057
1059
1077Numeric interp(ConstVectorView itw,
1078 ConstTensor4View a,
1079 const GridPos& tb,
1080 const GridPos& tp,
1081 const GridPos& tr,
1082 const GridPos& tc) {
1083 ARTS_ASSERT(is_size(itw, 16)); // We need 16 interpolation
1084 // weights.
1085
1086 // Check that interpolation weights are valid. The sum of all
1087 // weights (last dimension) must always be approximately one.
1088 ARTS_ASSERT(is_same_within_epsilon(sum(itw), 1, sum_check_epsilon));
1089
1090 // To store interpolated value:
1091 Numeric tia = 0;
1092
1093 Index iti = 0;
1094 for (Index b = 0; b < 2; ++b)
1095 for (Index p = 0; p < 2; ++p)
1096 for (Index r = 0; r < 2; ++r)
1097 for (Index c = 0; c < 2; ++c) {
1098 tia += a(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
1099 itw[iti];
1100 ++iti;
1101 }
1102
1103 return tia;
1104}
1105
1107
1126Numeric interp(ConstVectorView itw,
1127 ConstTensor5View a,
1128 const GridPos& ts,
1129 const GridPos& tb,
1130 const GridPos& tp,
1131 const GridPos& tr,
1132 const GridPos& tc) {
1133 ARTS_ASSERT(is_size(itw, 32)); // We need 32 interpolation
1134 // weights.
1135
1136 // Check that interpolation weights are valid. The sum of all
1137 // weights (last dimension) must always be approximately one.
1138 ARTS_ASSERT(is_same_within_epsilon(sum(itw), 1, sum_check_epsilon));
1139
1140 // To store interpolated value:
1141 Numeric tia = 0;
1142
1143 Index iti = 0;
1144 for (Index s = 0; s < 2; ++s)
1145 for (Index b = 0; b < 2; ++b)
1146 for (Index p = 0; p < 2; ++p)
1147 for (Index r = 0; r < 2; ++r)
1148 for (Index c = 0; c < 2; ++c) {
1149 tia += a(ts.idx + s,
1150 tb.idx + b,
1151 tp.idx + p,
1152 tr.idx + r,
1153 tc.idx + c) *
1154 itw[iti];
1155 ++iti;
1156 }
1157
1158 return tia;
1159}
1160
1162
1182Numeric interp(ConstVectorView itw,
1183 ConstTensor6View a,
1184 const GridPos& tv,
1185 const GridPos& ts,
1186 const GridPos& tb,
1187 const GridPos& tp,
1188 const GridPos& tr,
1189 const GridPos& tc) {
1190 ARTS_ASSERT(is_size(itw, 64)); // We need 64 interpolation
1191 // weights.
1192
1193 // Check that interpolation weights are valid. The sum of all
1194 // weights (last dimension) must always be approximately one.
1195 ARTS_ASSERT(is_same_within_epsilon(sum(itw), 1, sum_check_epsilon));
1196
1197 // To store interpolated value:
1198 Numeric tia = 0;
1199
1200 Index iti = 0;
1201 for (Index v = 0; v < 2; ++v)
1202 for (Index s = 0; s < 2; ++s)
1203 for (Index b = 0; b < 2; ++b)
1204 for (Index p = 0; p < 2; ++p)
1205 for (Index r = 0; r < 2; ++r)
1206 for (Index c = 0; c < 2; ++c) {
1207 tia += a(tv.idx + v,
1208 ts.idx + s,
1209 tb.idx + b,
1210 tp.idx + p,
1211 tr.idx + r,
1212 tc.idx + c) *
1213 itw[iti];
1214 ++iti;
1215 }
1216
1217 return tia;
1218}
1219
1221// Blue interpolation
1223
1225
1240void interpweights(MatrixView itw, const ArrayOfGridPos& cgp) {
1241 Index n = cgp.nelem();
1242 ARTS_ASSERT(is_size(itw, n, 2)); // We must store 2 interpolation
1243 // weights for each position.
1244
1245 // We have to loop all the points in the sequence:
1246 for (Index i = 0; i < n; ++i) {
1247 // Current grid positions:
1248 const GridPos& tc = cgp[i];
1249
1250 // Interpolation weights are stored in this order (l=lower
1251 // u=upper, c=column):
1252 // 1. l-c
1253 // 2. u-c
1254
1255 Index iti = 0;
1256
1257 // We could use a straight-forward for loop here:
1258 //
1259 // for ( Index c=1; c>=0; --c )
1260 // {
1261 // ti[iti] = tc.fd[c];
1262 // ++iti;
1263 // }
1264 //
1265 // However, it is slightly faster to use pointers (I tried it,
1266 // the speed gain is about 20%). So we should write something
1267 // like:
1268 //
1269 // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
1270 // {
1271 // ti[iti] = *c;
1272 // ++iti;
1273 // }
1274 //
1275 // For higher dimensions we have to nest these loops. To avoid
1276 // typos and safe typing, I use the LOOPIT macro, which expands
1277 // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
1278 // COMMAND!
1279
1280 LOOPIT(c) {
1281 itw(i, iti) = *c;
1282 ++iti;
1283 }
1284 }
1285}
1286
1288
1309void interpweights(MatrixView itw,
1310 const ArrayOfGridPos& rgp,
1311 const ArrayOfGridPos& cgp) {
1312 Index n = cgp.nelem();
1313 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1314 ARTS_ASSERT(is_size(itw, n, 4)); // We must store 4 interpolation
1315 // weights for each position.
1316
1317 // We have to loop all the points in the sequence:
1318 for (Index i = 0; i < n; ++i) {
1319 // Current grid positions:
1320 const GridPos& tr = rgp[i];
1321 const GridPos& tc = cgp[i];
1322
1323 // Interpolation weights are stored in this order (l=lower
1324 // u=upper, r=row, c=column):
1325 // 1. l-r l-c
1326 // 2. l-r u-c
1327 // 3. u-r l-c
1328 // 4. u-r u-c
1329
1330 Index iti = 0;
1331
1332 LOOPIT(r)
1333 LOOPIT(c) {
1334 itw(i, iti) = (*r) * (*c);
1335 ++iti;
1336 }
1337 }
1338}
1339
1341
1363void interpweights(MatrixView itw,
1364 const ArrayOfGridPos& pgp,
1365 const ArrayOfGridPos& rgp,
1366 const ArrayOfGridPos& cgp) {
1367 Index n = cgp.nelem();
1368 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1369 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1370 ARTS_ASSERT(is_size(itw, n, 8)); // We must store 8 interpolation
1371 // weights for each position.
1372
1373 // We have to loop all the points in the sequence:
1374 for (Index i = 0; i < n; ++i) {
1375 // Current grid positions:
1376 const GridPos& tp = pgp[i];
1377 const GridPos& tr = rgp[i];
1378 const GridPos& tc = cgp[i];
1379
1380 Index iti = 0;
1381 LOOPIT(p)
1382 LOOPIT(r)
1383 LOOPIT(c) {
1384 itw(i, iti) = (*p) * (*r) * (*c);
1385 ++iti;
1386 }
1387 }
1388}
1389
1391
1414void interpweights(MatrixView itw,
1415 const ArrayOfGridPos& bgp,
1416 const ArrayOfGridPos& pgp,
1417 const ArrayOfGridPos& rgp,
1418 const ArrayOfGridPos& cgp) {
1419 Index n = cgp.nelem();
1420 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1421 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1422 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1423 ARTS_ASSERT(is_size(itw, n, 16)); // We must store 16 interpolation
1424 // weights for each position.
1425
1426 // We have to loop all the points in the sequence:
1427 for (Index i = 0; i < n; ++i) {
1428 // Current grid positions:
1429 const GridPos& tb = bgp[i];
1430 const GridPos& tp = pgp[i];
1431 const GridPos& tr = rgp[i];
1432 const GridPos& tc = cgp[i];
1433
1434 Index iti = 0;
1435 LOOPIT(b)
1436 LOOPIT(p)
1437 LOOPIT(r)
1438 LOOPIT(c) {
1439 itw(i, iti) = (*b) * (*p) * (*r) * (*c);
1440 ++iti;
1441 }
1442 }
1443}
1444
1446
1470void interpweights(MatrixView itw,
1471 const ArrayOfGridPos& sgp,
1472 const ArrayOfGridPos& bgp,
1473 const ArrayOfGridPos& pgp,
1474 const ArrayOfGridPos& rgp,
1475 const ArrayOfGridPos& cgp) {
1476 Index n = cgp.nelem();
1477 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1478 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1479 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1480 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1481 ARTS_ASSERT(is_size(itw, n, 32)); // We must store 32 interpolation
1482 // weights for each position.
1483
1484 // We have to loop all the points in the sequence:
1485 for (Index i = 0; i < n; ++i) {
1486 // Current grid positions:
1487 const GridPos& ts = sgp[i];
1488 const GridPos& tb = bgp[i];
1489 const GridPos& tp = pgp[i];
1490 const GridPos& tr = rgp[i];
1491 const GridPos& tc = cgp[i];
1492
1493 Index iti = 0;
1494 LOOPIT(s)
1495 LOOPIT(b)
1496 LOOPIT(p)
1497 LOOPIT(r)
1498 LOOPIT(c) {
1499 itw(i, iti) = (*s) * (*b) * (*p) * (*r) * (*c);
1500 ++iti;
1501 }
1502 }
1503}
1504
1506
1531void interpweights(MatrixView itw,
1532 const ArrayOfGridPos& vgp,
1533 const ArrayOfGridPos& sgp,
1534 const ArrayOfGridPos& bgp,
1535 const ArrayOfGridPos& pgp,
1536 const ArrayOfGridPos& rgp,
1537 const ArrayOfGridPos& cgp) {
1538 Index n = cgp.nelem();
1539 ARTS_ASSERT(is_size(vgp, n)); // vgp must have same size as cgp.
1540 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1541 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1542 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1543 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1544 ARTS_ASSERT(is_size(itw, n, 64)); // We must store 64 interpolation
1545 // weights for each position.
1546
1547 // We have to loop all the points in the sequence:
1548 for (Index i = 0; i < n; ++i) {
1549 // Current grid positions:
1550 const GridPos& tv = vgp[i];
1551 const GridPos& ts = sgp[i];
1552 const GridPos& tb = bgp[i];
1553 const GridPos& tp = pgp[i];
1554 const GridPos& tr = rgp[i];
1555 const GridPos& tc = cgp[i];
1556
1557 Index iti = 0;
1558 LOOPIT(v)
1559 LOOPIT(s)
1560 LOOPIT(b)
1561 LOOPIT(p)
1562 LOOPIT(r)
1563 LOOPIT(c) {
1564 itw(i, iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
1565 ++iti;
1566 }
1567 }
1568}
1569
1571
1587void interp(VectorView ia,
1588 ConstMatrixView itw,
1589 ConstVectorView a,
1590 const ArrayOfGridPos& cgp) {
1591 Index n = cgp.nelem();
1592 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1593 ARTS_ASSERT(is_size(itw, n, 2)); // We need 2 interpolation
1594 // weights for each position.
1595
1596 // Check that interpolation weights are valid. The sum of all
1597 // weights (last dimension) must always be approximately one. We
1598 // only check the first element.
1600 is_same_within_epsilon(sum(itw(0, Range(joker))), 1, sum_check_epsilon));
1601
1602 // We have to loop all the points in the sequence:
1603 for (Index i = 0; i < n; ++i) {
1604 // Current grid positions:
1605 const GridPos& tc = cgp[i];
1606
1607 // Get handle to current element of output vector and initialize
1608 // it to zero:
1609 Numeric& tia = ia[i];
1610 tia = 0;
1611
1612 Index iti = 0;
1613 for (Index c = 0; c < 2; ++c) {
1614 ARTS_ASSERT(tc.idx + c < a.nelem()); // Temporary !?
1615 tia += a(tc.idx + c) * itw(i, iti);
1616 ++iti;
1617 }
1618 }
1619}
1620
1622
1643void interp(VectorView ia,
1644 ConstMatrixView itw,
1645 ConstMatrixView a,
1646 const ArrayOfGridPos& rgp,
1647 const ArrayOfGridPos& cgp) {
1648 Index n = cgp.nelem();
1649 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1650 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1651 ARTS_ASSERT(is_size(itw, n, 4)); // We need 4 interpolation
1652 // weights for each position.
1653
1654 // Check that interpolation weights are valid. The sum of all
1655 // weights (last dimension) must always be approximately one. We
1656 // only check the first element.
1658 is_same_within_epsilon(sum(itw(0, Range(joker))), 1, sum_check_epsilon));
1659
1660 // We have to loop all the points in the sequence:
1661 for (Index i = 0; i < n; ++i) {
1662 // Current grid positions:
1663 const GridPos& tr = rgp[i];
1664 const GridPos& tc = cgp[i];
1665
1666 // Get handle to current element of output vector and initialize
1667 // it to zero:
1668 Numeric& tia = ia[i];
1669 tia = 0;
1670
1671 Index iti = 0;
1672 for (Index r = 0; r < 2; ++r)
1673 for (Index c = 0; c < 2; ++c) {
1674 tia += a(tr.idx + r, tc.idx + c) * itw(i, iti);
1675 ++iti;
1676 }
1677 }
1678}
1679
1681
1703void interp(VectorView ia,
1704 ConstMatrixView itw,
1705 ConstTensor3View a,
1706 const ArrayOfGridPos& pgp,
1707 const ArrayOfGridPos& rgp,
1708 const ArrayOfGridPos& cgp) {
1709 Index n = cgp.nelem();
1710 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1711 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1712 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1713 ARTS_ASSERT(is_size(itw, n, 8)); // We need 8 interpolation
1714 // weights for each position.
1715
1716 // Check that interpolation weights are valid. The sum of all
1717 // weights (last dimension) must always be approximately one. We
1718 // only check the first element.
1720 is_same_within_epsilon(sum(itw(0, Range(joker))), 1, sum_check_epsilon));
1721
1722 // We have to loop all the points in the sequence:
1723 for (Index i = 0; i < n; ++i) {
1724 // Current grid positions:
1725 const GridPos& tp = pgp[i];
1726 const GridPos& tr = rgp[i];
1727 const GridPos& tc = cgp[i];
1728
1729 // Get handle to current element of output vector and initialize
1730 // it to zero:
1731 Numeric& tia = ia[i];
1732 tia = 0;
1733
1734 Index iti = 0;
1735 for (Index p = 0; p < 2; ++p)
1736 for (Index r = 0; r < 2; ++r)
1737 for (Index c = 0; c < 2; ++c) {
1738 tia += a(tp.idx + p, tr.idx + r, tc.idx + c) * itw(i, iti);
1739 ++iti;
1740 }
1741 }
1742}
1743
1745
1768void interp(VectorView ia,
1769 ConstMatrixView itw,
1770 ConstTensor4View a,
1771 const ArrayOfGridPos& bgp,
1772 const ArrayOfGridPos& pgp,
1773 const ArrayOfGridPos& rgp,
1774 const ArrayOfGridPos& cgp) {
1775 Index n = cgp.nelem();
1776 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1777 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1778 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1779 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1780 ARTS_ASSERT(is_size(itw, n, 16)); // We need 16 interpolation
1781 // weights for each position.
1782
1783 // Check that interpolation weights are valid. The sum of all
1784 // weights (last dimension) must always be approximately one. We
1785 // only check the first element.
1787 is_same_within_epsilon(sum(itw(0, Range(joker))), 1, sum_check_epsilon));
1788
1789 // We have to loop all the points in the sequence:
1790 for (Index i = 0; i < n; ++i) {
1791 // Current grid positions:
1792 const GridPos& tb = bgp[i];
1793 const GridPos& tp = pgp[i];
1794 const GridPos& tr = rgp[i];
1795 const GridPos& tc = cgp[i];
1796
1797 // Get handle to current element of output vector and initialize
1798 // it to zero:
1799 Numeric& tia = ia[i];
1800 tia = 0;
1801
1802 Index iti = 0;
1803 for (Index b = 0; b < 2; ++b)
1804 for (Index p = 0; p < 2; ++p)
1805 for (Index r = 0; r < 2; ++r)
1806 for (Index c = 0; c < 2; ++c) {
1807 tia += a(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
1808 itw(i, iti);
1809 ++iti;
1810 }
1811 }
1812}
1813
1815
1839void interp(VectorView ia,
1840 ConstMatrixView itw,
1841 ConstTensor5View a,
1842 const ArrayOfGridPos& sgp,
1843 const ArrayOfGridPos& bgp,
1844 const ArrayOfGridPos& pgp,
1845 const ArrayOfGridPos& rgp,
1846 const ArrayOfGridPos& cgp) {
1847 Index n = cgp.nelem();
1848 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1849 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1850 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1851 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1852 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1853 ARTS_ASSERT(is_size(itw, n, 32)); // We need 32 interpolation
1854 // weights for each position.
1855
1856 // Check that interpolation weights are valid. The sum of all
1857 // weights (last dimension) must always be approximately one. We
1858 // only check the first element.
1860 is_same_within_epsilon(sum(itw(0, Range(joker))), 1, sum_check_epsilon));
1861
1862 // We have to loop all the points in the sequence:
1863 for (Index i = 0; i < n; ++i) {
1864 // Current grid positions:
1865 const GridPos& ts = sgp[i];
1866 const GridPos& tb = bgp[i];
1867 const GridPos& tp = pgp[i];
1868 const GridPos& tr = rgp[i];
1869 const GridPos& tc = cgp[i];
1870
1871 // Get handle to current element of output vector and initialize
1872 // it to zero:
1873 Numeric& tia = ia[i];
1874 tia = 0;
1875
1876 Index iti = 0;
1877 for (Index s = 0; s < 2; ++s)
1878 for (Index b = 0; b < 2; ++b)
1879 for (Index p = 0; p < 2; ++p)
1880 for (Index r = 0; r < 2; ++r)
1881 for (Index c = 0; c < 2; ++c) {
1882 tia += a(ts.idx + s,
1883 tb.idx + b,
1884 tp.idx + p,
1885 tr.idx + r,
1886 tc.idx + c) *
1887 itw(i, iti);
1888 ++iti;
1889 }
1890 }
1891}
1892
1894
1919void interp(VectorView ia,
1920 ConstMatrixView itw,
1921 ConstTensor6View a,
1922 const ArrayOfGridPos& vgp,
1923 const ArrayOfGridPos& sgp,
1924 const ArrayOfGridPos& bgp,
1925 const ArrayOfGridPos& pgp,
1926 const ArrayOfGridPos& rgp,
1927 const ArrayOfGridPos& cgp) {
1928 Index n = cgp.nelem();
1929 ARTS_ASSERT(is_size(ia, n)); // ia must have same size as cgp.
1930 ARTS_ASSERT(is_size(vgp, n)); // vgp must have same size as cgp.
1931 ARTS_ASSERT(is_size(sgp, n)); // sgp must have same size as cgp.
1932 ARTS_ASSERT(is_size(bgp, n)); // bgp must have same size as cgp.
1933 ARTS_ASSERT(is_size(pgp, n)); // pgp must have same size as cgp.
1934 ARTS_ASSERT(is_size(rgp, n)); // rgp must have same size as cgp.
1935 ARTS_ASSERT(is_size(itw, n, 64)); // We need 64 interpolation
1936 // weights for each position.
1937
1938 // Check that interpolation weights are valid. The sum of all
1939 // weights (last dimension) must always be approximately one. We
1940 // only check the first element.
1942 is_same_within_epsilon(sum(itw(0, Range(joker))), 1, sum_check_epsilon));
1943
1944 // We have to loop all the points in the sequence:
1945 for (Index i = 0; i < n; ++i) {
1946 // Current grid positions:
1947 const GridPos& tv = vgp[i];
1948 const GridPos& ts = sgp[i];
1949 const GridPos& tb = bgp[i];
1950 const GridPos& tp = pgp[i];
1951 const GridPos& tr = rgp[i];
1952 const GridPos& tc = cgp[i];
1953
1954 // Get handle to current element of output vector and initialize
1955 // it to zero:
1956 Numeric& tia = ia[i];
1957 tia = 0;
1958
1959 Index iti = 0;
1960 for (Index v = 0; v < 2; ++v)
1961 for (Index s = 0; s < 2; ++s)
1962 for (Index b = 0; b < 2; ++b)
1963 for (Index p = 0; p < 2; ++p)
1964 for (Index r = 0; r < 2; ++r)
1965 for (Index c = 0; c < 2; ++c) {
1966 tia += a(tv.idx + v,
1967 ts.idx + s,
1968 tb.idx + b,
1969 tp.idx + p,
1970 tr.idx + r,
1971 tc.idx + c) *
1972 itw(i, iti);
1973 ++iti;
1974 }
1975 }
1976}
1977
1979// Green interpolation
1981
1983
2004void interpweights(Tensor3View itw,
2005 const ArrayOfGridPos& rgp,
2006 const ArrayOfGridPos& cgp) {
2007 Index nr = rgp.nelem();
2008 Index nc = cgp.nelem();
2009 ARTS_ASSERT(is_size(itw, nr, nc, 4)); // We must store 4 interpolation
2010 // weights for each position.
2011
2012 // We have to loop all the points in the new grid:
2013 for (Index ir = 0; ir < nr; ++ir) {
2014 // Current grid position:
2015 const GridPos& tr = rgp[ir];
2016
2017 for (Index ic = 0; ic < nc; ++ic) {
2018 // Current grid position:
2019 const GridPos& tc = cgp[ic];
2020
2021 // Interpolation weights are stored in this order (l=lower
2022 // u=upper, r=row, c=column):
2023 // 1. l-r l-c
2024 // 2. l-r u-c
2025 // 3. u-r l-c
2026 // 4. u-r u-c
2027
2028 Index iti = 0;
2029
2030 LOOPIT(r)
2031 LOOPIT(c) {
2032 itw(ir, ic, iti) = (*r) * (*c);
2033 ++iti;
2034 }
2035 }
2036 }
2037}
2038
2040
2062void interpweights(Tensor4View itw,
2063 const ArrayOfGridPos& pgp,
2064 const ArrayOfGridPos& rgp,
2065 const ArrayOfGridPos& cgp) {
2066 Index np = pgp.nelem();
2067 Index nr = rgp.nelem();
2068 Index nc = cgp.nelem();
2069 // We must store 8 interpolation weights for each position:
2070 ARTS_ASSERT(is_size(itw, np, nr, nc, 8));
2071
2072 // We have to loop all the points in the new grid:
2073 for (Index ip = 0; ip < np; ++ip) {
2074 const GridPos& tp = pgp[ip];
2075 for (Index ir = 0; ir < nr; ++ir) {
2076 const GridPos& tr = rgp[ir];
2077 for (Index ic = 0; ic < nc; ++ic) {
2078 const GridPos& tc = cgp[ic];
2079
2080 Index iti = 0;
2081
2082 LOOPIT(p)
2083 LOOPIT(r)
2084 LOOPIT(c) {
2085 itw(ip, ir, ic, iti) = (*p) * (*r) * (*c);
2086 ++iti;
2087 }
2088 }
2089 }
2090 }
2091}
2092
2094
2117void interpweights(Tensor5View itw,
2118 const ArrayOfGridPos& bgp,
2119 const ArrayOfGridPos& pgp,
2120 const ArrayOfGridPos& rgp,
2121 const ArrayOfGridPos& cgp) {
2122 Index nb = bgp.nelem();
2123 Index np = pgp.nelem();
2124 Index nr = rgp.nelem();
2125 Index nc = cgp.nelem();
2126 // We must store 16 interpolation weights for each position:
2127 ARTS_ASSERT(is_size(itw, nb, np, nr, nc, 16));
2128
2129 // We have to loop all the points in the new grid:
2130 for (Index ib = 0; ib < nb; ++ib) {
2131 const GridPos& tb = bgp[ib];
2132 for (Index ip = 0; ip < np; ++ip) {
2133 const GridPos& tp = pgp[ip];
2134 for (Index ir = 0; ir < nr; ++ir) {
2135 const GridPos& tr = rgp[ir];
2136 for (Index ic = 0; ic < nc; ++ic) {
2137 const GridPos& tc = cgp[ic];
2138
2139 Index iti = 0;
2140
2141 LOOPIT(b)
2142 LOOPIT(p)
2143 LOOPIT(r)
2144 LOOPIT(c) {
2145 itw(ib, ip, ir, ic, iti) = (*b) * (*p) * (*r) * (*c);
2146 ++iti;
2147 }
2148 }
2149 }
2150 }
2151 }
2152}
2153
2155
2179void interpweights(Tensor6View itw,
2180 const ArrayOfGridPos& sgp,
2181 const ArrayOfGridPos& bgp,
2182 const ArrayOfGridPos& pgp,
2183 const ArrayOfGridPos& rgp,
2184 const ArrayOfGridPos& cgp) {
2185 Index ns = sgp.nelem();
2186 Index nb = bgp.nelem();
2187 Index np = pgp.nelem();
2188 Index nr = rgp.nelem();
2189 Index nc = cgp.nelem();
2190 // We must store 32 interpolation weights for each position:
2191 ARTS_ASSERT(is_size(itw, ns, nb, np, nr, nc, 32));
2192
2193 // We have to loop all the points in the new grid:
2194 for (Index is = 0; is < ns; ++is) {
2195 const GridPos& ts = sgp[is];
2196 for (Index ib = 0; ib < nb; ++ib) {
2197 const GridPos& tb = bgp[ib];
2198 for (Index ip = 0; ip < np; ++ip) {
2199 const GridPos& tp = pgp[ip];
2200 for (Index ir = 0; ir < nr; ++ir) {
2201 const GridPos& tr = rgp[ir];
2202 for (Index ic = 0; ic < nc; ++ic) {
2203 const GridPos& tc = cgp[ic];
2204
2205 Index iti = 0;
2206
2207 LOOPIT(s)
2208 LOOPIT(b)
2209 LOOPIT(p)
2210 LOOPIT(r)
2211 LOOPIT(c) {
2212 itw(is, ib, ip, ir, ic, iti) =
2213 (*s) * (*b) * (*p) * (*r) * (*c);
2214 ++iti;
2215 }
2216 }
2217 }
2218 }
2219 }
2220 }
2221}
2222
2224
2249void interpweights(Tensor7View itw,
2250 const ArrayOfGridPos& vgp,
2251 const ArrayOfGridPos& sgp,
2252 const ArrayOfGridPos& bgp,
2253 const ArrayOfGridPos& pgp,
2254 const ArrayOfGridPos& rgp,
2255 const ArrayOfGridPos& cgp) {
2256 Index nv = vgp.nelem();
2257 Index ns = sgp.nelem();
2258 Index nb = bgp.nelem();
2259 Index np = pgp.nelem();
2260 Index nr = rgp.nelem();
2261 Index nc = cgp.nelem();
2262 // We must store 64 interpolation weights for each position:
2263 ARTS_ASSERT(is_size(itw, nv, ns, nb, np, nr, nc, 64));
2264
2265 // We have to loop all the points in the new grid:
2266 for (Index iv = 0; iv < nv; ++iv) {
2267 const GridPos& tv = vgp[iv];
2268 for (Index is = 0; is < ns; ++is) {
2269 const GridPos& ts = sgp[is];
2270 for (Index ib = 0; ib < nb; ++ib) {
2271 const GridPos& tb = bgp[ib];
2272 for (Index ip = 0; ip < np; ++ip) {
2273 const GridPos& tp = pgp[ip];
2274 for (Index ir = 0; ir < nr; ++ir) {
2275 const GridPos& tr = rgp[ir];
2276 for (Index ic = 0; ic < nc; ++ic) {
2277 const GridPos& tc = cgp[ic];
2278
2279 Index iti = 0;
2280
2281 LOOPIT(v)
2282 LOOPIT(s)
2283 LOOPIT(b)
2284 LOOPIT(p)
2285 LOOPIT(r)
2286 LOOPIT(c) {
2287 itw(iv, is, ib, ip, ir, ic, iti) =
2288 (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
2289 ++iti;
2290 }
2291 }
2292 }
2293 }
2294 }
2295 }
2296 }
2297}
2298
2300
2321void interp(MatrixView ia,
2322 ConstTensor3View itw,
2323 ConstMatrixView a,
2324 const ArrayOfGridPos& rgp,
2325 const ArrayOfGridPos& cgp) {
2326 Index nr = rgp.nelem();
2327 Index nc = cgp.nelem();
2328 ARTS_ASSERT(is_size(ia, nr, nc));
2329 ARTS_ASSERT(is_size(itw, nr, nc, 4)); // We need 4 interpolation
2330 // weights for each position.
2331
2332 // Check that interpolation weights are valid. The sum of all
2333 // weights (last dimension) must always be approximately one. We
2334 // only check the first element.
2335 ARTS_ASSERT(is_same_within_epsilon(
2336 sum(itw(0, 0, Range(joker))), 1, sum_check_epsilon));
2337
2338 // We have to loop all the points in the new grid:
2339 for (Index ir = 0; ir < nr; ++ir) {
2340 // Current grid position:
2341 const GridPos& tr = rgp[ir];
2342
2343 for (Index ic = 0; ic < nc; ++ic) {
2344 // Current grid position:
2345 const GridPos& tc = cgp[ic];
2346
2347 // Get handle to current element of output tensor and initialize
2348 // it to zero:
2349 Numeric& tia = ia(ir, ic);
2350 tia = 0;
2351
2352 Index iti = 0;
2353 for (Index r = 0; r < 2; ++r)
2354 for (Index c = 0; c < 2; ++c) {
2355 ARTS_ASSERT(tr.idx + r < a.nrows()); // Temporary !?
2356 ARTS_ASSERT(tc.idx + c < a.ncols()); // Temporary !?
2357 tia += a(tr.idx + r, tc.idx + c) * itw(ir, ic, iti);
2358 ++iti;
2359 }
2360 }
2361 }
2362}
2363
2365
2387void interp(Tensor3View ia,
2388 ConstTensor4View itw,
2389 ConstTensor3View a,
2390 const ArrayOfGridPos& pgp,
2391 const ArrayOfGridPos& rgp,
2392 const ArrayOfGridPos& cgp) {
2393 Index np = pgp.nelem();
2394 Index nr = rgp.nelem();
2395 Index nc = cgp.nelem();
2396 ARTS_ASSERT(is_size(ia, np, nr, nc));
2397 ARTS_ASSERT(is_size(itw, np, nr, nc, 8));
2398
2399 // Check that interpolation weights are valid. The sum of all
2400 // weights (last dimension) must always be approximately one. We
2401 // only check the first element.
2402 ARTS_ASSERT(is_same_within_epsilon(
2403 sum(itw(0, 0, 0, Range(joker))), 1, sum_check_epsilon));
2404
2405 // We have to loop all the points in the new grid:
2406 for (Index ip = 0; ip < np; ++ip) {
2407 const GridPos& tp = pgp[ip];
2408 for (Index ir = 0; ir < nr; ++ir) {
2409 const GridPos& tr = rgp[ir];
2410 for (Index ic = 0; ic < nc; ++ic) {
2411 // Current grid position:
2412 const GridPos& tc = cgp[ic];
2413
2414 // Get handle to current element of output tensor and
2415 // initialize it to zero:
2416 Numeric& tia = ia(ip, ir, ic);
2417 tia = 0;
2418
2419 Index iti = 0;
2420 for (Index p = 0; p < 2; ++p)
2421 for (Index r = 0; r < 2; ++r)
2422 for (Index c = 0; c < 2; ++c) {
2423 ARTS_ASSERT(tp.idx + p < a.npages()); // Temporary !?
2424 ARTS_ASSERT(tr.idx + r < a.nrows()); // Temporary !?
2425 ARTS_ASSERT(tc.idx + c < a.ncols()); // Temporary !?
2426 tia += a(tp.idx + p, tr.idx + r, tc.idx + c) *
2427 itw(ip, ir, ic, iti);
2428 ++iti;
2429 }
2430 }
2431 }
2432 }
2433}
2434
2436
2459void interp(Tensor4View ia,
2460 ConstTensor5View itw,
2461 ConstTensor4View a,
2462 const ArrayOfGridPos& bgp,
2463 const ArrayOfGridPos& pgp,
2464 const ArrayOfGridPos& rgp,
2465 const ArrayOfGridPos& cgp) {
2466 Index nb = bgp.nelem();
2467 Index np = pgp.nelem();
2468 Index nr = rgp.nelem();
2469 Index nc = cgp.nelem();
2470 ARTS_ASSERT(is_size(ia, nb, np, nr, nc));
2471 ARTS_ASSERT(is_size(itw, nb, np, nr, nc, 16));
2472
2473 // Check that interpolation weights are valid. The sum of all
2474 // weights (last dimension) must always be approximately one. We
2475 // only check the first element.
2476 ARTS_ASSERT(is_same_within_epsilon(
2477 sum(itw(0, 0, 0, 0, Range(joker))), 1, sum_check_epsilon));
2478
2479 // We have to loop all the points in the new grid:
2480 for (Index ib = 0; ib < nb; ++ib) {
2481 const GridPos& tb = bgp[ib];
2482 for (Index ip = 0; ip < np; ++ip) {
2483 const GridPos& tp = pgp[ip];
2484 for (Index ir = 0; ir < nr; ++ir) {
2485 const GridPos& tr = rgp[ir];
2486 for (Index ic = 0; ic < nc; ++ic) {
2487 // Current grid position:
2488 const GridPos& tc = cgp[ic];
2489
2490 // Get handle to current element of output tensor and
2491 // initialize it to zero:
2492 Numeric& tia = ia(ib, ip, ir, ic);
2493 tia = 0;
2494
2495 Index iti = 0;
2496 for (Index b = 0; b < 2; ++b)
2497 for (Index p = 0; p < 2; ++p)
2498 for (Index r = 0; r < 2; ++r)
2499 for (Index c = 0; c < 2; ++c) {
2500 tia += a(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
2501 itw(ib, ip, ir, ic, iti);
2502 ++iti;
2503 }
2504 }
2505 }
2506 }
2507 }
2508}
2509
2511
2535void interp(Tensor5View ia,
2536 ConstTensor6View itw,
2537 ConstTensor5View a,
2538 const ArrayOfGridPos& sgp,
2539 const ArrayOfGridPos& bgp,
2540 const ArrayOfGridPos& pgp,
2541 const ArrayOfGridPos& rgp,
2542 const ArrayOfGridPos& cgp) {
2543 Index ns = sgp.nelem();
2544 Index nb = bgp.nelem();
2545 Index np = pgp.nelem();
2546 Index nr = rgp.nelem();
2547 Index nc = cgp.nelem();
2548 ARTS_ASSERT(is_size(ia, ns, nb, np, nr, nc));
2549 ARTS_ASSERT(is_size(itw, ns, nb, np, nr, nc, 32));
2550
2551 // Check that interpolation weights are valid. The sum of all
2552 // weights (last dimension) must always be approximately one. We
2553 // only check the first element.
2554 ARTS_ASSERT(is_same_within_epsilon(
2555 sum(itw(0, 0, 0, 0, 0, Range(joker))), 1, sum_check_epsilon));
2556
2557 // We have to loop all the points in the new grid:
2558 for (Index is = 0; is < ns; ++is) {
2559 const GridPos& ts = sgp[is];
2560 for (Index ib = 0; ib < nb; ++ib) {
2561 const GridPos& tb = bgp[ib];
2562 for (Index ip = 0; ip < np; ++ip) {
2563 const GridPos& tp = pgp[ip];
2564 for (Index ir = 0; ir < nr; ++ir) {
2565 const GridPos& tr = rgp[ir];
2566 for (Index ic = 0; ic < nc; ++ic) {
2567 // Current grid position:
2568 const GridPos& tc = cgp[ic];
2569
2570 // Get handle to current element of output tensor and
2571 // initialize it to zero:
2572 Numeric& tia = ia(is, ib, ip, ir, ic);
2573 tia = 0;
2574
2575 Index iti = 0;
2576 for (Index s = 0; s < 2; ++s)
2577 for (Index b = 0; b < 2; ++b)
2578 for (Index p = 0; p < 2; ++p)
2579 for (Index r = 0; r < 2; ++r)
2580 for (Index c = 0; c < 2; ++c) {
2581 tia += a(ts.idx + s,
2582 tb.idx + b,
2583 tp.idx + p,
2584 tr.idx + r,
2585 tc.idx + c) *
2586 itw(is, ib, ip, ir, ic, iti);
2587 ++iti;
2588 }
2589 }
2590 }
2591 }
2592 }
2593 }
2594}
2595
2597
2622void interp(Tensor6View ia,
2623 ConstTensor7View itw,
2624 ConstTensor6View a,
2625 const ArrayOfGridPos& vgp,
2626 const ArrayOfGridPos& sgp,
2627 const ArrayOfGridPos& bgp,
2628 const ArrayOfGridPos& pgp,
2629 const ArrayOfGridPos& rgp,
2630 const ArrayOfGridPos& cgp) {
2631 Index nv = vgp.nelem();
2632 Index ns = sgp.nelem();
2633 Index nb = bgp.nelem();
2634 Index np = pgp.nelem();
2635 Index nr = rgp.nelem();
2636 Index nc = cgp.nelem();
2637 ARTS_ASSERT(is_size(ia, nv, ns, nb, np, nr, nc));
2638 ARTS_ASSERT(is_size(itw, nv, ns, nb, np, nr, nc, 64));
2639
2640 // Check that interpolation weights are valid. The sum of all
2641 // weights (last dimension) must always be approximately one. We
2642 // only check the first element.
2643 ARTS_ASSERT(is_same_within_epsilon(
2644 sum(itw(0, 0, 0, 0, 0, 0, Range(joker))), 1, sum_check_epsilon));
2645
2646 // We have to loop all the points in the new grid:
2647 for (Index iv = 0; iv < nv; ++iv) {
2648 const GridPos& tv = vgp[iv];
2649 for (Index is = 0; is < ns; ++is) {
2650 const GridPos& ts = sgp[is];
2651 for (Index ib = 0; ib < nb; ++ib) {
2652 const GridPos& tb = bgp[ib];
2653 for (Index ip = 0; ip < np; ++ip) {
2654 const GridPos& tp = pgp[ip];
2655 for (Index ir = 0; ir < nr; ++ir) {
2656 const GridPos& tr = rgp[ir];
2657 for (Index ic = 0; ic < nc; ++ic) {
2658 // Current grid position:
2659 const GridPos& tc = cgp[ic];
2660
2661 // Get handle to current element of output tensor and
2662 // initialize it to zero:
2663 Numeric& tia = ia(iv, is, ib, ip, ir, ic);
2664 tia = 0;
2665
2666 Index iti = 0;
2667 for (Index v = 0; v < 2; ++v)
2668 for (Index s = 0; s < 2; ++s)
2669 for (Index b = 0; b < 2; ++b)
2670 for (Index p = 0; p < 2; ++p)
2671 for (Index r = 0; r < 2; ++r)
2672 for (Index c = 0; c < 2; ++c) {
2673 tia += a(tv.idx + v,
2674 ts.idx + s,
2675 tb.idx + b,
2676 tp.idx + p,
2677 tr.idx + r,
2678 tc.idx + c) *
2679 itw(iv, is, ib, ip, ir, ic, iti);
2680 ++iti;
2681 }
2682 }
2683 }
2684 }
2685 }
2686 }
2687 }
2688}
2689
2690
2692
2711void polint(Numeric& y_int,
2712 Numeric& dy_int,
2713 ConstVectorView xa,
2714 ConstVectorView ya,
2715 const Index& n,
2716 const Numeric& x) {
2717 Index ns = 1;
2718 Numeric den, dif, dift, ho, hp, w;
2719
2720 dif = abs(x - xa[0]);
2721
2722 Vector c(n);
2723 Vector d(n);
2724
2725 // Find the index of the closest table entry
2726 for (Index i = 0; i < n; i++) {
2727 if ((dift = abs(x - xa[i])) < dif) {
2728 ns = i;
2729 dif = dift;
2730 }
2731 // Initialize c and d
2732 c[i] = ya[i];
2733 d[i] = ya[i];
2734 }
2735 // Initial approximation to y
2736 y_int = ya[ns--];
2737
2738 for (Index m = 1; m < n; m++) {
2739 for (Index i = 0; i < n - m; i++) {
2740 ho = xa[i] - x;
2741 hp = xa[i + m] - x;
2742 w = c[i + 1] - d[i];
2743 den = ho - hp;
2744 // This error occurs when two input xa's are identical.
2745 ARTS_ASSERT(den != 0.);
2746 den = w / den;
2747 d[i] = hp * den;
2748 c[i] = ho * den;
2749 }
2750 y_int += (dy_int = (2 * (ns + 1) < (n - m) ? c[ns + 1] : d[ns--]));
2751 }
2752}
2753
2755
2770Numeric interp_poly(ConstVectorView x,
2771 ConstVectorView y,
2772 const Numeric& x_i,
2773 const GridPos& gp) {
2774 Index N_x = x.nelem();
2775
2776 ARTS_ASSERT(N_x == y.nelem());
2777 ARTS_ASSERT(N_x > 2);
2778
2779 Vector xa(4), ya(4);
2780 Numeric y_int;
2781 Numeric dy_int;
2782 y_int = 0.;
2783
2784 // 1 - polynomial interpolation (3 points) with grid position search
2785 // 2 - polynomial interpolation (3 points) without grid position search
2786 // 3 - polynomial interpolation (4 points)
2787
2788 Index interp_method = 1;
2789
2790 if (interp_method == 1) {
2791 // Pick out three points for interpolation
2792 if ((gp.fd[0] <= 0.5 && gp.idx > 0) || gp.idx == N_x - 2) {
2793 xa[0] = x[gp.idx - 1];
2794 xa[1] = x[gp.idx];
2795 xa[2] = x[gp.idx + 1];
2796
2797 ya[0] = y[gp.idx - 1];
2798 ya[1] = y[gp.idx];
2799 ya[2] = y[gp.idx + 1];
2800 }
2801
2802 else if ((gp.fd[0] > 0.5 && gp.idx < N_x - 2) || gp.idx == 0) {
2803 xa[0] = x[gp.idx];
2804 xa[1] = x[gp.idx + 1];
2805 xa[2] = x[gp.idx + 2];
2806
2807 ya[0] = y[gp.idx];
2808 ya[1] = y[gp.idx + 1];
2809 ya[2] = y[gp.idx + 2];
2810 }
2811
2812 else if (gp.idx == N_x - 1) {
2813 xa[0] = x[N_x - 2];
2814 xa[1] = x[N_x - 1];
2815 xa[2] = x[N_x];
2816
2817 ya[0] = y[N_x - 2];
2818 ya[1] = y[N_x - 1];
2819 ya[2] = y[N_x];
2820 } else {
2821 ARTS_ASSERT(false);
2822 arts_exit();
2823 }
2824
2825 polint(y_int, dy_int, xa, ya, 3, x_i);
2826
2827 }
2828
2829 else if (interp_method == 2) {
2830 if (gp.idx == 0) {
2831 xa[0] = x[gp.idx];
2832 xa[1] = x[gp.idx + 1];
2833 xa[2] = x[gp.idx + 2];
2834
2835 ya[0] = y[gp.idx];
2836 ya[1] = y[gp.idx + 1];
2837 ya[2] = y[gp.idx + 2];
2838 } else if (gp.idx == N_x - 1) {
2839 xa[0] = x[gp.idx - 2];
2840 xa[1] = x[gp.idx - 1];
2841 xa[2] = x[gp.idx];
2842
2843 ya[0] = y[gp.idx - 2];
2844 ya[1] = y[gp.idx - 1];
2845 ya[2] = y[gp.idx];
2846 } else {
2847 xa[0] = x[gp.idx - 1];
2848 xa[1] = x[gp.idx];
2849 xa[2] = x[gp.idx + 1];
2850
2851 ya[0] = y[gp.idx - 1];
2852 ya[1] = y[gp.idx];
2853 ya[2] = y[gp.idx + 1];
2854 }
2855
2856 // Polynominal interpolation, n = 3
2857 polint(y_int, dy_int, xa, ya, 3, x_i);
2858 }
2859
2860 else if (interp_method == 3) {
2861 // Take 4 points
2862 if (gp.idx == 0) {
2863 xa[0] = -x[gp.idx + 1];
2864 xa[1] = x[gp.idx + 0];
2865 xa[2] = x[gp.idx + 1];
2866 xa[3] = x[gp.idx + 2];
2867
2868 ya[0] = y[gp.idx + 1];
2869 ya[1] = y[gp.idx + 0];
2870 ya[2] = y[gp.idx + 1];
2871 ya[3] = y[gp.idx + 2];
2872 } else if (gp.idx == N_x - 1) {
2873 xa[0] = x[gp.idx - 1];
2874 xa[1] = x[gp.idx - 0];
2875 xa[2] = 2 * x[gp.idx] - x[gp.idx - 1];
2876 xa[3] = 2 * x[gp.idx] - x[gp.idx - 2];
2877
2878 ya[0] = y[gp.idx - 1];
2879 ya[1] = y[gp.idx - 0];
2880 ya[2] = y[gp.idx - 1];
2881 ya[3] = y[gp.idx - 2];
2882 } else if (gp.idx == N_x - 2) {
2883 xa[0] = x[gp.idx - 2];
2884 xa[1] = x[gp.idx - 1];
2885 xa[2] = x[gp.idx];
2886 xa[3] = x[gp.idx + 1];
2887
2888 ya[0] = y[gp.idx - 2];
2889 ya[1] = y[gp.idx - 1];
2890 ya[2] = y[gp.idx];
2891 ya[3] = y[gp.idx + 1];
2892 } else {
2893 xa[0] = x[gp.idx - 1];
2894 xa[1] = x[gp.idx];
2895 xa[2] = x[gp.idx + 1];
2896 xa[3] = x[gp.idx + 2];
2897
2898 ya[0] = y[gp.idx - 1];
2899 ya[1] = y[gp.idx];
2900 ya[2] = y[gp.idx + 1];
2901 ya[3] = y[gp.idx + 2];
2902 }
2903 // Polinominal interpolation, n = 4
2904 polint(y_int, dy_int, xa, ya, 4, x_i);
2905 }
2906
2907 return y_int;
2908}
This file contains the definition of Array.
void arts_exit(int status)
This is the exit function of ARTS.
Definition: arts.cc:25
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
#define DEBUG_ONLY(...)
Definition: debug.h:53
#define ARTS_ASSERT(condition,...)
Definition: debug.h:84
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.
Structure to store a grid position.
Definition: interpolation.h:56
std::array< Numeric, 2 > fd
Definition: interpolation.h:58
Index idx
Definition: interpolation.h:57
#define d
#define v
#define w
#define a
#define c
#define b