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