ARTS  2.2.66
interpolation.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
40 #include <iostream>
41 #include <cmath>
42 #include "array.h"
43 #include "check_input.h"
44 #include "interpolation.h"
45 #include "logic.h"
46 
47 
48 
49 // File-global constants:
50 
52 
64 
65 
67 
76 const Numeric FD_TOL = 1.5e-3;
77 
78 
79 
80 
82 
91 #define LOOPIT(x) for ( const Numeric* x=&t##x.fd[1]; x>=&t##x.fd[0]; --x )
92 
93 
94 
95 
97 
105 ostream& operator<<(ostream& os, const GridPos& gp)
106 {
107  os << gp.idx << " " << gp.fd[0] << " " << gp.fd[1] << "\n";
108  return os;
109 }
110 
111 
113 
168  ConstVectorView old_grid,
169  ConstVectorView new_grid,
170  const Numeric& extpolfac )
171 {
172  const Index n_old = old_grid.nelem();
173  const Index n_new = new_grid.nelem();
174 
175  // Assert that gp has the right size:
176  assert( is_size(gp,n_new) );
177 
178  // Assert, that the old grid has more than one element
179  assert( 1 < n_old );
180 
181  // This function hast two parts, depending on whether old_grid is
182  // sorted in ascending or descending order. Maybe that's not too
183  // elegant, but it's the most efficient, because in this way there
184  // is no additional runtime overhead to handle both cases.
185 
186  // We use only the first two elements to decide how the grid is
187  // sorted. (The rest of the grid is checked later by an assert.)
188  // If both values are the same, we still assume the grid is
189  // ascending. However, this will lead to an assertion fail later on,
190  // because the grid has to be strictly sorted.
191  bool ascending = ( old_grid[0] <= old_grid[1] );
192 
193  if (ascending)
194  {
195  // So, old_grid should always be sorted in strictly ascending order.
196  // (No duplicate values.)
197  // Assert that this is so. This may depend on user input, however,
198  // inside this elementary function is not the place to check for
199  // that. There must be runtime checks on higher levels to insure
200  // that all grids are sorted. The assertion here is just as a last
201  // safety check.
202  assert( is_increasing(old_grid) );
203 
204  // Limits of extrapolation.
205  const Numeric og_min = old_grid[0] -
206  extpolfac * ( old_grid[1] - old_grid[0] );
207  const Numeric og_max = old_grid[n_old-1] +
208  extpolfac * ( old_grid[n_old-1] - old_grid[n_old-2] );
209 
210  // We will make no firm assumptions about the new grid. But the case
211  // that we have in mind is that it is either also sorted, or at
212  // least partially sorted, for example like this:
213  // 5 4 3 2 2.5 3 4
214  // This kind of sequence should be typical if we interpolate
215  // atmospheric fields along a limb propagation path.
216 
217  // Let's get some idea where the first point in the new grid is,
218  // relative to the old grid. We use linear interpolation between the
219  // maximum and the minimum of the old grid for this. (Taking
220  // into account the small allowed extrapolation.)
221  Numeric frac = (new_grid[0]-og_min)/(og_max-og_min);
222 
223  // We are not checking if the value of frac is reasonable,
224  // because there is another assertion below to catch the
225  // consequences.
226 
227  // Initialize current_position:
228  Index current_position = (Index) rint(frac*(Numeric)(n_old-2));
229 
230  // The above statement should satisfy
231  // 0 <= current_position <= n_old-2
232  // Assert that this is indeed the case.
233 // cout << "frac = " << frac << "\n";
234 // cout << "current_position = " << current_position << "\n";
235  assert( 0<= current_position );
236  assert( current_position <= n_old-2 );
237 
238  // The variables lower and upper are used to remember the value of
239  // the old grid at current_position and one above current_position:
240  Numeric lower = old_grid[current_position];
241  Numeric upper = old_grid[current_position+1];
242 
243  // Loop over all points in the new grid:
244  for ( Index i_new=0; i_new<n_new; ++i_new )
245  {
246  // Get a reference to the current element of gp:
247  GridPos& tgp = gp[i_new];
248  // And on the current value of the new grid:
249  const Numeric tng = new_grid[i_new];
250 
251  // Verify that the new grid is within the limits of
252  // extrapolation that we have defined by og_min and og_max:
253  assert( og_min <= tng );
254  assert( tng <= og_max );
255 
256 // cout << "lower / tng / upper = "
257 // << lower << " / "
258 // << tng << " / "
259 // << upper << "\n";
260 
261  // Is current_position too high?
262  // (The current_position>0 condition is there so that the position
263  // stays 0 for extrapolation.)
264  if ( tng < lower && current_position > 0 )
265  {
266  do
267  {
268  --current_position;
269  lower = old_grid[current_position];
270  }
271  while ( tng < lower && current_position > 0 );
272 
273  upper = old_grid[current_position+1];
274 
275  tgp.idx = current_position;
276  tgp.fd[0] = (tng-lower)/(upper-lower);
277  tgp.fd[1] = 1.0 - tgp.fd[0];
278  }
279  else
280  {
281  // Is it too low?
282  // (The current_position<n_old condition is there so
283  // that uppers stays n_old-1 for extrapolation.)
284  if ( tng >= upper && current_position < n_old-2 )
285  {
286  do
287  {
288  ++current_position;
289  upper = old_grid[current_position+1];
290  }
291  while ( tng >= upper && current_position < n_old-2 );
292 
293  lower = old_grid[current_position];
294 
295  tgp.idx = current_position;
296  tgp.fd[0] = (tng-lower)/(upper-lower);
297  tgp.fd[1] = 1.0 - tgp.fd[0];
298  }
299  else
300  {
301  // None of the other two conditions were true. That means:
302  // lower <= tng < upper. The current_position is still
303  // valid.
304 
305  // SAB 2010-04-28: Note that if a new grid point is exactly on
306  // top of an old grid point, then you are now guaranteed to get
307  // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
308 
309  tgp.idx = current_position;
310  tgp.fd[0] = (tng-lower)/(upper-lower);
311  tgp.fd[1] = 1.0 - tgp.fd[0];
312  }
313  }
314 
315 // cout << tgp.idx << " " << tgp.fd[0] << " " << tgp.fd[1] << endl;
316 
317  // Safety check to ensure the above promise:
318  assert(old_grid[tgp.idx]<=tng || tgp.idx==0);
319  }
320  }
321  else // if (ascending)
322  {
323  // Now we are in the "descending old grid" part. We do exactly
324  // the same as in the other part, just accounting for the
325  // different order of things. Comments here refer only to
326  // interesting differences from the ascending case. See that
327  // case for more general comments.
328 
329  // This time ensure strictly descending order:
330  assert( is_decreasing(old_grid) );
331 
332  // The max is now the first point, the min the last point!
333  // I think the sign is right here...
334  const Numeric og_max = old_grid[0] -
335  extpolfac * ( old_grid[1] - old_grid[0] );
336  const Numeric og_min = old_grid[n_old-1] +
337  extpolfac * ( old_grid[n_old-1] - old_grid[n_old-2] );
338 
339  // We have to take 1- here, because we are starting from the
340  // high end.
341  Numeric frac = 1 - (new_grid[0]-og_min)/(og_max-og_min);
342 
343  // We are not checking if the value of frac is reasonable,
344  // because there is another assertion below to catch the
345  // consequences.
346 
347  Index current_position = (Index) rint(frac*(Numeric)(n_old-2));
348 
349  // The above statement should satisfy
350  // 0 <= current_position <= n_old-2
351  // Assert that this is indeed the case.
352 // cout << "frac = " << frac << "\n";
353 // cout << "current_position = " << current_position << "\n";
354  assert( 0<= current_position );
355  assert( current_position <= n_old-2 );
356 
357  // Note, that old_grid[lower] has a higher numerical value than
358  // old_grid[upper]!
359  Numeric lower = old_grid[current_position];
360  Numeric upper = old_grid[current_position+1];
361 
362  for ( Index i_new=0; i_new<n_new; ++i_new )
363  {
364  GridPos& tgp = gp[i_new];
365  const Numeric tng = new_grid[i_new];
366 
367  // Verify that the new grid is within the limits of
368  // extrapolation that we have defined by og_min and og_max:
369  assert( og_min <= tng );
370  assert( tng <= og_max );
371 
372 // cout << "lower / tng / upper = "
373 // << lower << " / "
374 // << tng << " / "
375 // << upper << "\n";
376 
377  // Is current_position too high? (Sign of comparison changed
378  // compared to ascending case!)
379  if ( tng > lower && current_position > 0 )
380  {
381  do
382  {
383  --current_position;
384  lower = old_grid[current_position];
385  }
386  while ( tng > lower && current_position > 0 );
387 
388  upper = old_grid[current_position+1];
389 
390  tgp.idx = current_position;
391  tgp.fd[0] = (tng-lower)/(upper-lower);
392  tgp.fd[1] = 1.0 - tgp.fd[0];
393  }
394  else
395  {
396  // Is it too low? (Sign of comparison changed
397  // compared to ascending case!)
398  if ( tng <= upper && current_position < n_old-2 )
399  {
400  do
401  {
402  ++current_position;
403  upper = old_grid[current_position+1];
404  }
405  while ( tng <= upper && current_position < n_old-2 );
406 
407  lower = old_grid[current_position];
408 
409  tgp.idx = current_position;
410  tgp.fd[0] = (tng-lower)/(upper-lower);
411  tgp.fd[1] = 1.0 - tgp.fd[0];
412  }
413  else
414  {
415  // None of the other two conditions were true. That means:
416  // lower >= tng > upper. The current_position is still
417  // valid. (Note that upper and lower have switched
418  // place compared to the ascending case.)
419 
420  // SAB 2010-04-28: Note that if a new grid point is exactly on
421  // top of an old grid point, then you are now guaranteed to get
422  // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
423 
424  tgp.idx = current_position;
425  tgp.fd[0] = (tng-lower)/(upper-lower);
426  tgp.fd[1] = 1.0 - tgp.fd[0];
427  }
428  }
429 
430  // Safety check to ensure the above promise:
431  assert(old_grid[tgp.idx]>=tng || tgp.idx==0);
432  }
433  }
434 }
435 
436 
437 
439 
460 void gridpos( GridPos& gp,
461  ConstVectorView old_grid,
462  const Numeric& new_grid,
463  const Numeric& extpolfac )
464 {
465  ArrayOfGridPos agp(1);
466  Vector v( 1, new_grid );
467  gridpos( agp, old_grid, v, extpolfac );
468  gridpos_copy( gp, agp[0] );
469 }
470 
471 
472 
474 
489  ArrayOfGridPos& gp,
490  ConstVectorView grid )
491 {
492  const Index n = grid.nelem();
493  gp.resize( n );
494 
495  for( Index i=0; i<n-1; i++ )
496  {
497  gp[i].idx = i;
498  gp[i].fd[0] = 0;
499  gp[i].fd[1] = 1;
500  }
501  gp[n-1].idx = n-2;
502  gp[n-1].fd[0] = 1;
503  gp[n-1].fd[1] = 0;
504 }
505 
506 
508 
517 void gridpos_copy( GridPos& gp_new, const GridPos& gp_old )
518 {
519  gp_new.idx = gp_old.idx;
520  gp_new.fd[0] = gp_old.fd[0];
521  gp_new.fd[1] = gp_old.fd[1];
522 }
523 
524 
525 
527 
540 {
541  return ( Numeric(gp.idx) + gp.fd[0] );
542 }
543 
544 
545 
547 
560 {
561  // Catch values that "must" be wrong
562  assert( gp.fd[0] > -FD_TOL );
563  assert( gp.fd[0] < 1.0 + FD_TOL );
564  assert( gp.fd[1] > -FD_TOL );
565  assert( gp.fd[1] < 1.0 + FD_TOL );
566 
567  if( gp.fd[0] < 0.0 )
568  { gp.fd[0] = 0.0; gp.fd[1] = 1.0; }
569  else if( gp.fd[0] > 1.0 )
570  { gp.fd[0] = 1.0; gp.fd[1] = 0.0; }
571 
572  if( gp.fd[1] < 0.0 )
573  { gp.fd[1] = 0.0; gp.fd[0] = 1.0; }
574  else if( gp.fd[1] > 1.0 )
575  { gp.fd[1] = 1.0; gp.fd[0] = 0.0; }
576 }
577 
578 
579 
581 
603  GridPos& gp,
604  const Index& n )
605 {
606  assert( gp.idx >= 0 );
607 
608  // If fd=1, shift to grid index above
609  if( gp.fd[0] > 0.5 )
610  {
611  gp.idx += 1;
612  }
613  gp.fd[0] = 0;
614  gp.fd[1] = 1;
615 
616  assert( gp.idx < n );
617 
618  // End of complete grid range must be handled separately
619  if( gp.idx == n-1 )
620  {
621  gp.idx -= 1;
622  gp.fd[0] = 1;
623  gp.fd[1] = 0;
624  }
625 }
626 
627 
628 
630 
644  GridPos& gp,
645  const Index& ie )
646 {
647  if( gp.idx == ie )
648  {
649  assert( gp.fd[0] < 0.005 ); // To capture obviously bad cases
650  gp.idx -= 1;
651  gp.fd[0] = 1.0;
652  gp.fd[1] = 0.0;
653  }
654 }
655 
656 
658 
672  ArrayOfGridPos& gp,
673  const Index& ie )
674 {
675  for (Index i = 0; i < gp.nelem(); i++ )
676  {
677  if( gp[i].idx == ie )
678  {
679  assert( gp[i].fd[0] < 0.005 ); // To capture obviously bad cases
680  gp[i].idx -= 1;
681  gp[i].fd[0] = 1.0;
682  gp[i].fd[1] = 0.0;
683  }
684  }
685 }
686 
687 
688 
690 
703  const GridPos& gp,
704  const Index& i,
705  const bool& strict )
706 {
707  if( strict )
708  {
709  // Assume that gridpos_force_end_fd has been used. The expression 0==0
710  // should be safer than 1==1.
711  if( gp.idx == i && gp.fd[0] == 0 )
712  { return true; }
713  else if( gp.idx == i-1 && gp.fd[1] == 0 )
714  { return true; }
715  }
716  else
717  {
718  if( gp.idx == i && gp.fd[0] < FD_TOL )
719  { return true; }
720  else if( gp.idx == i-1 && gp.fd[1] < FD_TOL )
721  { return true; }
722  }
723  return false;
724 }
725 
726 
727 
729 
748  const GridPos& gp,
749  const bool& upwards )
750 {
751  assert( gp.fd[0] >= 0 );
752  assert( gp.fd[1] >= 0 );
753 
754  // Not at a grid point
755  if( gp.fd[0] > 0 && gp.fd[1] > 0 )
756  {
757  return gp.idx;
758  }
759 
760  // Fractional distance 0
761  else if( gp.fd[0] == 0 )
762  {
763  if( upwards )
764  return gp.idx;
765  else
766  return gp.idx - 1;
767  }
768 
769  // Fractional distance 1
770  else
771  {
772  if( upwards )
773  return gp.idx + 1;
774  else
775  return gp.idx;
776  }
777 }
778 
779 
780 
781 
782 
783 
785 // Red Interpolation
787 
789 
803  const GridPos& tc )
804 {
805  assert(is_size(itw,2)); // We must store 2 interpolation
806  // weights.
807 
808  // Interpolation weights are stored in this order (l=lower
809  // u=upper, c=column):
810  // 1. l-c
811  // 2. u-c
812 
813  Index iti = 0;
814 
815  // We could use a straight-forward for loop here:
816  //
817  // for ( Index c=1; c>=0; --c )
818  // {
819  // ti[iti] = tc.fd[c];
820  // ++iti;
821  // }
822  //
823  // However, it is slightly faster to use pointers (I tried it,
824  // the speed gain is about 20%). So we should write something
825  // like:
826  //
827  // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
828  // {
829  // ti[iti] = *c;
830  // ++iti;
831  // }
832  //
833  // For higher dimensions we have to nest these loops. To avoid
834  // typos and safe typing, I use the LOOPIT macro, which expands
835  // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
836  // COMMAND!
837 
838  LOOPIT(c)
839  {
840  itw.get(iti) = *c;
841  ++iti;
842  }
843 }
844 
846 
861  const GridPos& tr,
862  const GridPos& tc )
863 {
864  assert(is_size(itw,4)); // We must store 4 interpolation
865  // weights.
866  Index iti = 0;
867 
868  LOOPIT(r)
869  LOOPIT(c)
870  {
871  itw.get(iti) = (*r) * (*c);
872  ++iti;
873  }
874 }
875 
877 
893  const GridPos& tp,
894  const GridPos& tr,
895  const GridPos& tc )
896 {
897  assert(is_size(itw,8)); // We must store 8 interpolation
898  // weights.
899  Index iti = 0;
900 
901  LOOPIT(p)
902  LOOPIT(r)
903  LOOPIT(c)
904  {
905  itw.get(iti) = (*p) * (*r) * (*c);
906  ++iti;
907  }
908 }
909 
911 
928  const GridPos& tb,
929  const GridPos& tp,
930  const GridPos& tr,
931  const GridPos& tc )
932 {
933  assert(is_size(itw,16)); // We must store 16 interpolation
934  // weights.
935  Index iti = 0;
936 
937  LOOPIT(b)
938  LOOPIT(p)
939  LOOPIT(r)
940  LOOPIT(c)
941  {
942  itw.get(iti) = (*b) * (*p) * (*r) * (*c);
943  ++iti;
944  }
945 }
946 
948 
966  const GridPos& ts,
967  const GridPos& tb,
968  const GridPos& tp,
969  const GridPos& tr,
970  const GridPos& tc )
971 {
972  assert(is_size(itw,32)); // We must store 32 interpolation
973  // weights.
974  Index iti = 0;
975 
976  LOOPIT(s)
977  LOOPIT(b)
978  LOOPIT(p)
979  LOOPIT(r)
980  LOOPIT(c)
981  {
982  itw.get(iti) = (*s) * (*b) * (*p) * (*r) * (*c);
983  ++iti;
984  }
985 }
986 
988 
1007  const GridPos& tv,
1008  const GridPos& ts,
1009  const GridPos& tb,
1010  const GridPos& tp,
1011  const GridPos& tr,
1012  const GridPos& tc )
1013 {
1014  assert(is_size(itw,64)); // We must store 64 interpolation
1015  // weights.
1016  Index iti = 0;
1017 
1018  LOOPIT(v)
1019  LOOPIT(s)
1020  LOOPIT(b)
1021  LOOPIT(p)
1022  LOOPIT(r)
1023  LOOPIT(c)
1024  {
1025  itw.get(iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
1026  ++iti;
1027  }
1028 }
1029 
1031 
1047  ConstVectorView a,
1048  const GridPos& tc )
1049 {
1050  assert(is_size(itw,2)); // We need 2 interpolation
1051  // weights.
1052 
1053  // Check that interpolation weights are valid. The sum of all
1054  // weights (last dimension) must always be approximately one.
1055  assert( is_same_within_epsilon( itw.sum(),
1056  1,
1057  sum_check_epsilon ) );
1058 
1059  // To store interpolated value:
1060  Numeric tia = 0;
1061 
1062  Index iti = 0;
1063  for ( Index c=0; c<2; ++c )
1064  {
1065  tia += a.get(tc.idx+c) * itw.get(iti);
1066  ++iti;
1067  }
1068 
1069  return tia;
1070 }
1071 
1073 
1090  ConstMatrixView a,
1091  const GridPos& tr,
1092  const GridPos& tc )
1093 {
1094  assert(is_size(itw,4)); // We need 4 interpolation
1095  // weights.
1096 
1097  // Check that interpolation weights are valid. The sum of all
1098  // weights (last dimension) must always be approximately one.
1099  assert( is_same_within_epsilon( itw.sum(),
1100  1,
1101  sum_check_epsilon ) );
1102 
1103  // To store interpolated value:
1104  Numeric tia = 0;
1105 
1106  Index iti = 0;
1107  for ( Index r=0; r<2; ++r )
1108  for ( Index c=0; c<2; ++c )
1109  {
1110  tia += a.get(tr.idx+r,
1111  tc.idx+c) * itw.get(iti);
1112  ++iti;
1113  }
1114 
1115  return tia;
1116 }
1117 
1119 
1137  ConstTensor3View a,
1138  const GridPos& tp,
1139  const GridPos& tr,
1140  const GridPos& tc )
1141 {
1142  assert(is_size(itw,8)); // We need 8 interpolation
1143  // weights.
1144 
1145  // Check that interpolation weights are valid. The sum of all
1146  // weights (last dimension) must always be approximately one.
1147  assert( is_same_within_epsilon( itw.sum(),
1148  1,
1149  sum_check_epsilon ) );
1150 
1151  // To store interpolated value:
1152  Numeric tia = 0;
1153 
1154  Index iti = 0;
1155  for ( Index p=0; p<2; ++p )
1156  for ( Index r=0; r<2; ++r )
1157  for ( Index c=0; c<2; ++c )
1158  {
1159  tia += a.get(tp.idx+p,
1160  tr.idx+r,
1161  tc.idx+c) * itw.get(iti);
1162  ++iti;
1163  }
1164 
1165  return tia;
1166 }
1167 
1169 
1188  ConstTensor4View a,
1189  const GridPos& tb,
1190  const GridPos& tp,
1191  const GridPos& tr,
1192  const GridPos& tc )
1193 {
1194  assert(is_size(itw,16)); // We need 16 interpolation
1195  // weights.
1196 
1197  // Check that interpolation weights are valid. The sum of all
1198  // weights (last dimension) must always be approximately one.
1199  assert( is_same_within_epsilon( itw.sum(),
1200  1,
1201  sum_check_epsilon ) );
1202 
1203  // To store interpolated value:
1204  Numeric tia = 0;
1205 
1206  Index iti = 0;
1207  for ( Index b=0; b<2; ++b )
1208  for ( Index p=0; p<2; ++p )
1209  for ( Index r=0; r<2; ++r )
1210  for ( Index c=0; c<2; ++c )
1211  {
1212  tia += a.get(tb.idx+b,
1213  tp.idx+p,
1214  tr.idx+r,
1215  tc.idx+c) * itw.get(iti);
1216  ++iti;
1217  }
1218 
1219  return tia;
1220 }
1221 
1223 
1243  ConstTensor5View a,
1244  const GridPos& ts,
1245  const GridPos& tb,
1246  const GridPos& tp,
1247  const GridPos& tr,
1248  const GridPos& tc )
1249 {
1250  assert(is_size(itw,32)); // We need 32 interpolation
1251  // weights.
1252 
1253  // Check that interpolation weights are valid. The sum of all
1254  // weights (last dimension) must always be approximately one.
1255  assert( is_same_within_epsilon( itw.sum(),
1256  1,
1257  sum_check_epsilon ) );
1258 
1259  // To store interpolated value:
1260  Numeric tia = 0;
1261 
1262  Index iti = 0;
1263  for ( Index s=0; s<2; ++s )
1264  for ( Index b=0; b<2; ++b )
1265  for ( Index p=0; p<2; ++p )
1266  for ( Index r=0; r<2; ++r )
1267  for ( Index c=0; c<2; ++c )
1268  {
1269  tia += a.get(ts.idx+s,
1270  tb.idx+b,
1271  tp.idx+p,
1272  tr.idx+r,
1273  tc.idx+c) * itw.get(iti);
1274  ++iti;
1275  }
1276 
1277  return tia;
1278 }
1279 
1281 
1302  ConstTensor6View a,
1303  const GridPos& tv,
1304  const GridPos& ts,
1305  const GridPos& tb,
1306  const GridPos& tp,
1307  const GridPos& tr,
1308  const GridPos& tc )
1309 {
1310  assert(is_size(itw,64)); // We need 64 interpolation
1311  // weights.
1312 
1313  // Check that interpolation weights are valid. The sum of all
1314  // weights (last dimension) must always be approximately one.
1315  assert( is_same_within_epsilon( itw.sum(),
1316  1,
1317  sum_check_epsilon ) );
1318 
1319  // To store interpolated value:
1320  Numeric tia = 0;
1321 
1322  Index iti = 0;
1323  for ( Index v=0; v<2; ++v )
1324  for ( Index s=0; s<2; ++s )
1325  for ( Index b=0; b<2; ++b )
1326  for ( Index p=0; p<2; ++p )
1327  for ( Index r=0; r<2; ++r )
1328  for ( Index c=0; c<2; ++c )
1329  {
1330  tia += a.get(tv.idx+v,
1331  ts.idx+s,
1332  tb.idx+b,
1333  tp.idx+p,
1334  tr.idx+r,
1335  tc.idx+c) * itw.get(iti);
1336  ++iti;
1337  }
1338 
1339  return tia;
1340 }
1341 
1342 
1344 // Blue interpolation
1346 
1348 
1364  const ArrayOfGridPos& cgp )
1365 {
1366  Index n = cgp.nelem();
1367  assert(is_size(itw,n,2)); // We must store 2 interpolation
1368  // weights for each position.
1369 
1370  // We have to loop all the points in the sequence:
1371  for ( Index i=0; i<n; ++i )
1372  {
1373  // Current grid positions:
1374  const GridPos& tc = cgp[i];
1375 
1376  // Interpolation weights are stored in this order (l=lower
1377  // u=upper, c=column):
1378  // 1. l-c
1379  // 2. u-c
1380 
1381  Index iti = 0;
1382 
1383  // We could use a straight-forward for loop here:
1384  //
1385  // for ( Index c=1; c>=0; --c )
1386  // {
1387  // ti[iti] = tc.fd[c];
1388  // ++iti;
1389  // }
1390  //
1391  // However, it is slightly faster to use pointers (I tried it,
1392  // the speed gain is about 20%). So we should write something
1393  // like:
1394  //
1395  // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
1396  // {
1397  // ti[iti] = *c;
1398  // ++iti;
1399  // }
1400  //
1401  // For higher dimensions we have to nest these loops. To avoid
1402  // typos and safe typing, I use the LOOPIT macro, which expands
1403  // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
1404  // COMMAND!
1405 
1406  LOOPIT(c)
1407  {
1408  itw.get(i,iti) = *c;
1409  ++iti;
1410  }
1411  }
1412 }
1413 
1415 
1437  const ArrayOfGridPos& rgp,
1438  const ArrayOfGridPos& cgp )
1439 {
1440  Index n = cgp.nelem();
1441  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1442  assert(is_size(itw,n,4)); // We must store 4 interpolation
1443  // weights for each position.
1444 
1445  // We have to loop all the points in the sequence:
1446  for ( Index i=0; i<n; ++i )
1447  {
1448  // Current grid positions:
1449  const GridPos& tr = rgp[i];
1450  const GridPos& tc = cgp[i];
1451 
1452  // Interpolation weights are stored in this order (l=lower
1453  // u=upper, r=row, c=column):
1454  // 1. l-r l-c
1455  // 2. l-r u-c
1456  // 3. u-r l-c
1457  // 4. u-r u-c
1458 
1459  Index iti = 0;
1460 
1461  LOOPIT(r)
1462  LOOPIT(c)
1463  {
1464  itw.get(i,iti) = (*r) * (*c);
1465  ++iti;
1466  }
1467  }
1468 }
1469 
1471 
1494  const ArrayOfGridPos& pgp,
1495  const ArrayOfGridPos& rgp,
1496  const ArrayOfGridPos& cgp )
1497 {
1498  Index n = cgp.nelem();
1499  assert(is_size(pgp,n)); // pgp must have same size as cgp.
1500  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1501  assert(is_size(itw,n,8)); // We must store 8 interpolation
1502  // weights for each position.
1503 
1504  // We have to loop all the points in the sequence:
1505  for ( Index i=0; i<n; ++i )
1506  {
1507  // Current grid positions:
1508  const GridPos& tp = pgp[i];
1509  const GridPos& tr = rgp[i];
1510  const GridPos& tc = cgp[i];
1511 
1512  Index iti = 0;
1513  LOOPIT(p)
1514  LOOPIT(r)
1515  LOOPIT(c)
1516  {
1517  itw.get(i,iti) = (*p) * (*r) * (*c);
1518  ++iti;
1519  }
1520  }
1521 }
1522 
1524 
1548  const ArrayOfGridPos& bgp,
1549  const ArrayOfGridPos& pgp,
1550  const ArrayOfGridPos& rgp,
1551  const ArrayOfGridPos& cgp )
1552 {
1553  Index n = cgp.nelem();
1554  assert(is_size(bgp,n)); // bgp must have same size as cgp.
1555  assert(is_size(pgp,n)); // pgp must have same size as cgp.
1556  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1557  assert(is_size(itw,n,16)); // We must store 16 interpolation
1558  // weights for each position.
1559 
1560  // We have to loop all the points in the sequence:
1561  for ( Index i=0; i<n; ++i )
1562  {
1563  // Current grid positions:
1564  const GridPos& tb = bgp[i];
1565  const GridPos& tp = pgp[i];
1566  const GridPos& tr = rgp[i];
1567  const GridPos& tc = cgp[i];
1568 
1569  Index iti = 0;
1570  LOOPIT(b)
1571  LOOPIT(p)
1572  LOOPIT(r)
1573  LOOPIT(c)
1574  {
1575  itw.get(i,iti) = (*b) * (*p) * (*r) * (*c);
1576  ++iti;
1577  }
1578  }
1579 }
1580 
1582 
1607  const ArrayOfGridPos& sgp,
1608  const ArrayOfGridPos& bgp,
1609  const ArrayOfGridPos& pgp,
1610  const ArrayOfGridPos& rgp,
1611  const ArrayOfGridPos& cgp )
1612 {
1613  Index n = cgp.nelem();
1614  assert(is_size(sgp,n)); // sgp must have same size as cgp.
1615  assert(is_size(bgp,n)); // bgp must have same size as cgp.
1616  assert(is_size(pgp,n)); // pgp must have same size as cgp.
1617  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1618  assert(is_size(itw,n,32)); // We must store 32 interpolation
1619  // weights for each position.
1620 
1621  // We have to loop all the points in the sequence:
1622  for ( Index i=0; i<n; ++i )
1623  {
1624  // Current grid positions:
1625  const GridPos& ts = sgp[i];
1626  const GridPos& tb = bgp[i];
1627  const GridPos& tp = pgp[i];
1628  const GridPos& tr = rgp[i];
1629  const GridPos& tc = cgp[i];
1630 
1631  Index iti = 0;
1632  LOOPIT(s)
1633  LOOPIT(b)
1634  LOOPIT(p)
1635  LOOPIT(r)
1636  LOOPIT(c)
1637  {
1638  itw.get(i,iti) = (*s) * (*b) * (*p) * (*r) * (*c);
1639  ++iti;
1640  }
1641  }
1642 }
1643 
1645 
1671  const ArrayOfGridPos& vgp,
1672  const ArrayOfGridPos& sgp,
1673  const ArrayOfGridPos& bgp,
1674  const ArrayOfGridPos& pgp,
1675  const ArrayOfGridPos& rgp,
1676  const ArrayOfGridPos& cgp )
1677 {
1678  Index n = cgp.nelem();
1679  assert(is_size(vgp,n)); // vgp must have same size as cgp.
1680  assert(is_size(sgp,n)); // sgp must have same size as cgp.
1681  assert(is_size(bgp,n)); // bgp must have same size as cgp.
1682  assert(is_size(pgp,n)); // pgp must have same size as cgp.
1683  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1684  assert(is_size(itw,n,64)); // We must store 64 interpolation
1685  // weights for each position.
1686 
1687  // We have to loop all the points in the sequence:
1688  for ( Index i=0; i<n; ++i )
1689  {
1690  // Current grid positions:
1691  const GridPos& tv = vgp[i];
1692  const GridPos& ts = sgp[i];
1693  const GridPos& tb = bgp[i];
1694  const GridPos& tp = pgp[i];
1695  const GridPos& tr = rgp[i];
1696  const GridPos& tc = cgp[i];
1697 
1698  Index iti = 0;
1699  LOOPIT(v)
1700  LOOPIT(s)
1701  LOOPIT(b)
1702  LOOPIT(p)
1703  LOOPIT(r)
1704  LOOPIT(c)
1705  {
1706  itw.get(i,iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
1707  ++iti;
1708  }
1709  }
1710 }
1711 
1713 
1730  ConstMatrixView itw,
1731  ConstVectorView a,
1732  const ArrayOfGridPos& cgp)
1733 {
1734  Index n = cgp.nelem();
1735  assert(is_size(ia,n)); // ia must have same size as cgp.
1736  assert(is_size(itw,n,2)); // We need 2 interpolation
1737  // weights for each position.
1738 
1739  // Check that interpolation weights are valid. The sum of all
1740  // weights (last dimension) must always be approximately one. We
1741  // only check the first element.
1742  assert( is_same_within_epsilon( itw(0,Range(joker)).sum(),
1743  1,
1744  sum_check_epsilon ) );
1745 
1746  // We have to loop all the points in the sequence:
1747  for ( Index i=0; i<n; ++i )
1748  {
1749  // Current grid positions:
1750  const GridPos& tc = cgp[i];
1751 
1752  // Get handle to current element of output vector and initialize
1753  // it to zero:
1754  Numeric& tia = ia[i];
1755  tia = 0;
1756 
1757  Index iti = 0;
1758  for ( Index c=0; c<2; ++c )
1759  {
1760  tia += a.get(tc.idx+c) * itw.get(i,iti);
1761  ++iti;
1762  }
1763  }
1764 }
1765 
1767 
1789  ConstMatrixView itw,
1790  ConstMatrixView a,
1791  const ArrayOfGridPos& rgp,
1792  const ArrayOfGridPos& cgp)
1793 {
1794  Index n = cgp.nelem();
1795  assert(is_size(ia,n)); // ia must have same size as cgp.
1796  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1797  assert(is_size(itw,n,4)); // We need 4 interpolation
1798  // weights for each position.
1799 
1800  // Check that interpolation weights are valid. The sum of all
1801  // weights (last dimension) must always be approximately one. We
1802  // only check the first element.
1803  assert( is_same_within_epsilon( itw(0,Range(joker)).sum(),
1804  1,
1805  sum_check_epsilon ) );
1806 
1807  // We have to loop all the points in the sequence:
1808  for ( Index i=0; i<n; ++i )
1809  {
1810  // Current grid positions:
1811  const GridPos& tr = rgp[i];
1812  const GridPos& tc = cgp[i];
1813 
1814  // Get handle to current element of output vector and initialize
1815  // it to zero:
1816  Numeric& tia = ia[i];
1817  tia = 0;
1818 
1819  Index iti = 0;
1820  for ( Index r=0; r<2; ++r )
1821  for ( Index c=0; c<2; ++c )
1822  {
1823  tia += a.get(tr.idx+r,
1824  tc.idx+c) * itw.get(i,iti);
1825  ++iti;
1826  }
1827  }
1828 }
1829 
1831 
1854  ConstMatrixView itw,
1855  ConstTensor3View a,
1856  const ArrayOfGridPos& pgp,
1857  const ArrayOfGridPos& rgp,
1858  const ArrayOfGridPos& cgp)
1859 {
1860  Index n = cgp.nelem();
1861  assert(is_size(ia,n)); // ia must have same size as cgp.
1862  assert(is_size(pgp,n)); // pgp must have same size as cgp.
1863  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1864  assert(is_size(itw,n,8)); // We need 8 interpolation
1865  // weights for each position.
1866 
1867  // Check that interpolation weights are valid. The sum of all
1868  // weights (last dimension) must always be approximately one. We
1869  // only check the first element.
1870  assert( is_same_within_epsilon( itw(0,Range(joker)).sum(),
1871  1,
1872  sum_check_epsilon ) );
1873 
1874  // We have to loop all the points in the sequence:
1875  for ( Index i=0; i<n; ++i )
1876  {
1877  // Current grid positions:
1878  const GridPos& tp = pgp[i];
1879  const GridPos& tr = rgp[i];
1880  const GridPos& tc = cgp[i];
1881 
1882  // Get handle to current element of output vector and initialize
1883  // it to zero:
1884  Numeric& tia = ia[i];
1885  tia = 0;
1886 
1887  Index iti = 0;
1888  for ( Index p=0; p<2; ++p )
1889  for ( Index r=0; r<2; ++r )
1890  for ( Index c=0; c<2; ++c )
1891  {
1892  tia += a.get(tp.idx+p,
1893  tr.idx+r,
1894  tc.idx+c) * itw.get(i,iti);
1895  ++iti;
1896  }
1897  }
1898 }
1899 
1901 
1925  ConstMatrixView itw,
1926  ConstTensor4View a,
1927  const ArrayOfGridPos& bgp,
1928  const ArrayOfGridPos& pgp,
1929  const ArrayOfGridPos& rgp,
1930  const ArrayOfGridPos& cgp)
1931 {
1932  Index n = cgp.nelem();
1933  assert(is_size(ia,n)); // ia must have same size as cgp.
1934  assert(is_size(bgp,n)); // bgp must have same size as cgp.
1935  assert(is_size(pgp,n)); // pgp must have same size as cgp.
1936  assert(is_size(rgp,n)); // rgp must have same size as cgp.
1937  assert(is_size(itw,n,16)); // We need 16 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.
1943  assert( is_same_within_epsilon( itw(0,Range(joker)).sum(),
1944  1,
1945  sum_check_epsilon ) );
1946 
1947  // We have to loop all the points in the sequence:
1948  for ( Index i=0; i<n; ++i )
1949  {
1950  // Current grid positions:
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 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  {
1967  tia += a.get(tb.idx+b,
1968  tp.idx+p,
1969  tr.idx+r,
1970  tc.idx+c) * itw.get(i,iti);
1971  ++iti;
1972  }
1973  }
1974 }
1975 
1977 
2002  ConstMatrixView itw,
2003  ConstTensor5View a,
2004  const ArrayOfGridPos& sgp,
2005  const ArrayOfGridPos& bgp,
2006  const ArrayOfGridPos& pgp,
2007  const ArrayOfGridPos& rgp,
2008  const ArrayOfGridPos& cgp)
2009 {
2010  Index n = cgp.nelem();
2011  assert(is_size(ia,n)); // ia must have same size as cgp.
2012  assert(is_size(sgp,n)); // sgp must have same size as cgp.
2013  assert(is_size(bgp,n)); // bgp must have same size as cgp.
2014  assert(is_size(pgp,n)); // pgp must have same size as cgp.
2015  assert(is_size(rgp,n)); // rgp must have same size as cgp.
2016  assert(is_size(itw,n,32)); // We need 32 interpolation
2017  // weights for each position.
2018 
2019  // Check that interpolation weights are valid. The sum of all
2020  // weights (last dimension) must always be approximately one. We
2021  // only check the first element.
2022  assert( is_same_within_epsilon( itw(0,Range(joker)).sum(),
2023  1,
2024  sum_check_epsilon ) );
2025 
2026  // We have to loop all the points in the sequence:
2027  for ( Index i=0; i<n; ++i )
2028  {
2029  // Current grid positions:
2030  const GridPos& ts = sgp[i];
2031  const GridPos& tb = bgp[i];
2032  const GridPos& tp = pgp[i];
2033  const GridPos& tr = rgp[i];
2034  const GridPos& tc = cgp[i];
2035 
2036  // Get handle to current element of output vector and initialize
2037  // it to zero:
2038  Numeric& tia = ia[i];
2039  tia = 0;
2040 
2041  Index iti = 0;
2042  for ( Index s=0; s<2; ++s )
2043  for ( Index b=0; b<2; ++b )
2044  for ( Index p=0; p<2; ++p )
2045  for ( Index r=0; r<2; ++r )
2046  for ( Index c=0; c<2; ++c )
2047  {
2048  tia += a.get(ts.idx+s,
2049  tb.idx+b,
2050  tp.idx+p,
2051  tr.idx+r,
2052  tc.idx+c) * itw.get(i,iti);
2053  ++iti;
2054  }
2055  }
2056 }
2057 
2059 
2085  ConstMatrixView itw,
2086  ConstTensor6View a,
2087  const ArrayOfGridPos& vgp,
2088  const ArrayOfGridPos& sgp,
2089  const ArrayOfGridPos& bgp,
2090  const ArrayOfGridPos& pgp,
2091  const ArrayOfGridPos& rgp,
2092  const ArrayOfGridPos& cgp)
2093 {
2094  Index n = cgp.nelem();
2095  assert(is_size(ia,n)); // ia must have same size as cgp.
2096  assert(is_size(vgp,n)); // vgp must have same size as cgp.
2097  assert(is_size(sgp,n)); // sgp must have same size as cgp.
2098  assert(is_size(bgp,n)); // bgp must have same size as cgp.
2099  assert(is_size(pgp,n)); // pgp must have same size as cgp.
2100  assert(is_size(rgp,n)); // rgp must have same size as cgp.
2101  assert(is_size(itw,n,64)); // We need 64 interpolation
2102  // weights for each position.
2103 
2104  // Check that interpolation weights are valid. The sum of all
2105  // weights (last dimension) must always be approximately one. We
2106  // only check the first element.
2107  assert( is_same_within_epsilon( itw(0,Range(joker)).sum(),
2108  1,
2109  sum_check_epsilon ) );
2110 
2111  // We have to loop all the points in the sequence:
2112  for ( Index i=0; i<n; ++i )
2113  {
2114  // Current grid positions:
2115  const GridPos& tv = vgp[i];
2116  const GridPos& ts = sgp[i];
2117  const GridPos& tb = bgp[i];
2118  const GridPos& tp = pgp[i];
2119  const GridPos& tr = rgp[i];
2120  const GridPos& tc = cgp[i];
2121 
2122  // Get handle to current element of output vector and initialize
2123  // it to zero:
2124  Numeric& tia = ia[i];
2125  tia = 0;
2126 
2127  Index iti = 0;
2128  for ( Index v=0; v<2; ++v )
2129  for ( Index s=0; s<2; ++s )
2130  for ( Index b=0; b<2; ++b )
2131  for ( Index p=0; p<2; ++p )
2132  for ( Index r=0; r<2; ++r )
2133  for ( Index c=0; c<2; ++c )
2134  {
2135  tia += a.get(tv.idx+v,
2136  ts.idx+s,
2137  tb.idx+b,
2138  tp.idx+p,
2139  tr.idx+r,
2140  tc.idx+c) * itw.get(i,iti);
2141  ++iti;
2142  }
2143  }
2144 }
2145 
2147 // Green interpolation
2149 
2151 
2173  const ArrayOfGridPos& rgp,
2174  const ArrayOfGridPos& cgp )
2175 {
2176  Index nr = rgp.nelem();
2177  Index nc = cgp.nelem();
2178  assert(is_size(itw,nr,nc,4)); // We must store 4 interpolation
2179  // weights for each position.
2180 
2181  // We have to loop all the points in the new grid:
2182  for ( Index ir=0; ir<nr; ++ir )
2183  {
2184  // Current grid position:
2185  const GridPos& tr = rgp[ir];
2186 
2187  for ( Index ic=0; ic<nc; ++ic )
2188  {
2189  // Current grid position:
2190  const GridPos& tc = cgp[ic];
2191 
2192  // Interpolation weights are stored in this order (l=lower
2193  // u=upper, r=row, c=column):
2194  // 1. l-r l-c
2195  // 2. l-r u-c
2196  // 3. u-r l-c
2197  // 4. u-r u-c
2198 
2199  Index iti = 0;
2200 
2201  LOOPIT(r)
2202  LOOPIT(c)
2203  {
2204  itw.get(ir,ic,iti) = (*r) * (*c);
2205  ++iti;
2206  }
2207  }
2208  }
2209 }
2210 
2212 
2235  const ArrayOfGridPos& pgp,
2236  const ArrayOfGridPos& rgp,
2237  const ArrayOfGridPos& cgp )
2238 {
2239  Index np = pgp.nelem();
2240  Index nr = rgp.nelem();
2241  Index nc = cgp.nelem();
2242  // We must store 8 interpolation weights for each position:
2243  assert(is_size(itw,np,nr,nc,8));
2244 
2245  // We have to loop all the points in the new grid:
2246  for ( Index ip=0; ip<np; ++ip )
2247  {
2248  const GridPos& tp = pgp[ip];
2249  for ( Index ir=0; ir<nr; ++ir )
2250  {
2251  const GridPos& tr = rgp[ir];
2252  for ( Index ic=0; ic<nc; ++ic )
2253  {
2254  const GridPos& tc = cgp[ic];
2255 
2256  Index iti = 0;
2257 
2258  LOOPIT(p)
2259  LOOPIT(r)
2260  LOOPIT(c)
2261  {
2262  itw.get(ip,ir,ic,iti) =
2263  (*p) * (*r) * (*c);
2264  ++iti;
2265  }
2266  }
2267  }
2268  }
2269 }
2270 
2272 
2296  const ArrayOfGridPos& bgp,
2297  const ArrayOfGridPos& pgp,
2298  const ArrayOfGridPos& rgp,
2299  const ArrayOfGridPos& cgp )
2300 {
2301  Index nb = bgp.nelem();
2302  Index np = pgp.nelem();
2303  Index nr = rgp.nelem();
2304  Index nc = cgp.nelem();
2305  // We must store 16 interpolation weights for each position:
2306  assert(is_size(itw,nb,np,nr,nc,16));
2307 
2308  // We have to loop all the points in the new grid:
2309  for ( Index ib=0; ib<nb; ++ib )
2310  {
2311  const GridPos& tb = bgp[ib];
2312  for ( Index ip=0; ip<np; ++ip )
2313  {
2314  const GridPos& tp = pgp[ip];
2315  for ( Index ir=0; ir<nr; ++ir )
2316  {
2317  const GridPos& tr = rgp[ir];
2318  for ( Index ic=0; ic<nc; ++ic )
2319  {
2320  const GridPos& tc = cgp[ic];
2321 
2322  Index iti = 0;
2323 
2324  LOOPIT(b)
2325  LOOPIT(p)
2326  LOOPIT(r)
2327  LOOPIT(c)
2328  {
2329  itw.get(ib,ip,ir,ic,iti) =
2330  (*b) * (*p) * (*r) * (*c);
2331  ++iti;
2332  }
2333  }
2334  }
2335  }
2336  }
2337 }
2338 
2340 
2365  const ArrayOfGridPos& sgp,
2366  const ArrayOfGridPos& bgp,
2367  const ArrayOfGridPos& pgp,
2368  const ArrayOfGridPos& rgp,
2369  const ArrayOfGridPos& cgp )
2370 {
2371  Index ns = sgp.nelem();
2372  Index nb = bgp.nelem();
2373  Index np = pgp.nelem();
2374  Index nr = rgp.nelem();
2375  Index nc = cgp.nelem();
2376  // We must store 32 interpolation weights for each position:
2377  assert(is_size(itw,ns,nb,np,nr,nc,32));
2378 
2379  // We have to loop all the points in the new grid:
2380  for ( Index is=0; is<ns; ++is )
2381  {
2382  const GridPos& ts = sgp[is];
2383  for ( Index ib=0; ib<nb; ++ib )
2384  {
2385  const GridPos& tb = bgp[ib];
2386  for ( Index ip=0; ip<np; ++ip )
2387  {
2388  const GridPos& tp = pgp[ip];
2389  for ( Index ir=0; ir<nr; ++ir )
2390  {
2391  const GridPos& tr = rgp[ir];
2392  for ( Index ic=0; ic<nc; ++ic )
2393  {
2394  const GridPos& tc = cgp[ic];
2395 
2396  Index iti = 0;
2397 
2398  LOOPIT(s)
2399  LOOPIT(b)
2400  LOOPIT(p)
2401  LOOPIT(r)
2402  LOOPIT(c)
2403  {
2404  itw.get(is,ib,ip,ir,ic,iti) =
2405  (*s) * (*b) * (*p) * (*r) * (*c);
2406  ++iti;
2407  }
2408  }
2409  }
2410  }
2411  }
2412  }
2413 }
2414 
2416 
2442  const ArrayOfGridPos& vgp,
2443  const ArrayOfGridPos& sgp,
2444  const ArrayOfGridPos& bgp,
2445  const ArrayOfGridPos& pgp,
2446  const ArrayOfGridPos& rgp,
2447  const ArrayOfGridPos& cgp )
2448 {
2449  Index nv = vgp.nelem();
2450  Index ns = sgp.nelem();
2451  Index nb = bgp.nelem();
2452  Index np = pgp.nelem();
2453  Index nr = rgp.nelem();
2454  Index nc = cgp.nelem();
2455  // We must store 64 interpolation weights for each position:
2456  assert(is_size(itw,nv,ns,nb,np,nr,nc,64));
2457 
2458  // We have to loop all the points in the new grid:
2459  for ( Index iv=0; iv<nv; ++iv )
2460  {
2461  const GridPos& tv = vgp[iv];
2462  for ( Index is=0; is<ns; ++is )
2463  {
2464  const GridPos& ts = sgp[is];
2465  for ( Index ib=0; ib<nb; ++ib )
2466  {
2467  const GridPos& tb = bgp[ib];
2468  for ( Index ip=0; ip<np; ++ip )
2469  {
2470  const GridPos& tp = pgp[ip];
2471  for ( Index ir=0; ir<nr; ++ir )
2472  {
2473  const GridPos& tr = rgp[ir];
2474  for ( Index ic=0; ic<nc; ++ic )
2475  {
2476  const GridPos& tc = cgp[ic];
2477 
2478  Index iti = 0;
2479 
2480  LOOPIT(v)
2481  LOOPIT(s)
2482  LOOPIT(b)
2483  LOOPIT(p)
2484  LOOPIT(r)
2485  LOOPIT(c)
2486  {
2487  itw.get(iv,is,ib,ip,ir,ic,iti) =
2488  (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
2489  ++iti;
2490  }
2491  }
2492  }
2493  }
2494  }
2495  }
2496  }
2497 }
2498 
2500 
2522  ConstTensor3View itw,
2523  ConstMatrixView a,
2524  const ArrayOfGridPos& rgp,
2525  const ArrayOfGridPos& cgp)
2526 {
2527  Index nr = rgp.nelem();
2528  Index nc = cgp.nelem();
2529  assert(is_size(ia,nr,nc));
2530  assert(is_size(itw,nr,nc,4)); // We need 4 interpolation
2531  // weights for each position.
2532 
2533  // Check that interpolation weights are valid. The sum of all
2534  // weights (last dimension) must always be approximately one. We
2535  // only check the first element.
2536  assert( is_same_within_epsilon( itw(0,0,Range(joker)).sum(),
2537  1,
2538  sum_check_epsilon ) );
2539 
2540  // We have to loop all the points in the new grid:
2541  for ( Index ir=0; ir<nr; ++ir )
2542  {
2543  // Current grid position:
2544  const GridPos& tr = rgp[ir];
2545 
2546  for ( Index ic=0; ic<nc; ++ic )
2547  {
2548  // Current grid position:
2549  const GridPos& tc = cgp[ic];
2550 
2551  // Get handle to current element of output tensor and initialize
2552  // it to zero:
2553  Numeric& tia = ia(ir,ic);
2554  tia = 0;
2555 
2556  Index iti = 0;
2557  for ( Index r=0; r<2; ++r )
2558  for ( Index c=0; c<2; ++c )
2559  {
2560  tia += a.get(tr.idx+r,
2561  tc.idx+c) * itw.get(ir,ic,iti);
2562  ++iti;
2563  }
2564  }
2565  }
2566 }
2567 
2569 
2592  ConstTensor4View itw,
2593  ConstTensor3View a,
2594  const ArrayOfGridPos& pgp,
2595  const ArrayOfGridPos& rgp,
2596  const ArrayOfGridPos& cgp)
2597 {
2598  Index np = pgp.nelem();
2599  Index nr = rgp.nelem();
2600  Index nc = cgp.nelem();
2601  assert(is_size(ia,
2602  np,nr,nc));
2603  assert(is_size(itw,
2604  np,nr,nc,
2605  8));
2606 
2607  // Check that interpolation weights are valid. The sum of all
2608  // weights (last dimension) must always be approximately one. We
2609  // only check the first element.
2610  assert( is_same_within_epsilon( itw(0,0,0,Range(joker)).sum(),
2611  1,
2612  sum_check_epsilon ) );
2613 
2614  // We have to loop all the points in the new grid:
2615  for ( Index ip=0; ip<np; ++ip )
2616  {
2617  const GridPos& tp = pgp[ip];
2618  for ( Index ir=0; ir<nr; ++ir )
2619  {
2620  const GridPos& tr = rgp[ir];
2621  for ( Index ic=0; ic<nc; ++ic )
2622  {
2623  // Current grid position:
2624  const GridPos& tc = cgp[ic];
2625 
2626  // Get handle to current element of output tensor and
2627  // initialize it to zero:
2628  Numeric& tia = ia(ip,ir,ic);
2629  tia = 0;
2630 
2631  Index iti = 0;
2632  for ( Index p=0; p<2; ++p )
2633  for ( Index r=0; r<2; ++r )
2634  for ( Index c=0; c<2; ++c )
2635  {
2636  tia += a.get(tp.idx+p,
2637  tr.idx+r,
2638  tc.idx+c) * itw.get(ip,ir,ic,
2639  iti);
2640  ++iti;
2641  }
2642  }
2643  }
2644  }
2645 }
2646 
2648 
2672  ConstTensor5View itw,
2673  ConstTensor4View a,
2674  const ArrayOfGridPos& bgp,
2675  const ArrayOfGridPos& pgp,
2676  const ArrayOfGridPos& rgp,
2677  const ArrayOfGridPos& cgp)
2678 {
2679  Index nb = bgp.nelem();
2680  Index np = pgp.nelem();
2681  Index nr = rgp.nelem();
2682  Index nc = cgp.nelem();
2683  assert(is_size(ia,
2684  nb,np,nr,nc));
2685  assert(is_size(itw,
2686  nb,np,nr,nc,
2687  16));
2688 
2689  // Check that interpolation weights are valid. The sum of all
2690  // weights (last dimension) must always be approximately one. We
2691  // only check the first element.
2692  assert( is_same_within_epsilon( itw(0,0,0,0,Range(joker)).sum(),
2693  1,
2694  sum_check_epsilon ) );
2695 
2696  // We have to loop all the points in the new grid:
2697  for ( Index ib=0; ib<nb; ++ib )
2698  {
2699  const GridPos& tb = bgp[ib];
2700  for ( Index ip=0; ip<np; ++ip )
2701  {
2702  const GridPos& tp = pgp[ip];
2703  for ( Index ir=0; ir<nr; ++ir )
2704  {
2705  const GridPos& tr = rgp[ir];
2706  for ( Index ic=0; ic<nc; ++ic )
2707  {
2708  // Current grid position:
2709  const GridPos& tc = cgp[ic];
2710 
2711  // Get handle to current element of output tensor and
2712  // initialize it to zero:
2713  Numeric& tia = ia(ib,ip,ir,ic);
2714  tia = 0;
2715 
2716  Index iti = 0;
2717  for ( Index b=0; b<2; ++b )
2718  for ( Index p=0; p<2; ++p )
2719  for ( Index r=0; r<2; ++r )
2720  for ( Index c=0; c<2; ++c )
2721  {
2722  tia += a.get(tb.idx+b,
2723  tp.idx+p,
2724  tr.idx+r,
2725  tc.idx+c) * itw.get(ib,ip,ir,ic,
2726  iti);
2727  ++iti;
2728  }
2729  }
2730  }
2731  }
2732  }
2733 }
2734 
2736 
2761  ConstTensor6View itw,
2762  ConstTensor5View a,
2763  const ArrayOfGridPos& sgp,
2764  const ArrayOfGridPos& bgp,
2765  const ArrayOfGridPos& pgp,
2766  const ArrayOfGridPos& rgp,
2767  const ArrayOfGridPos& cgp)
2768 {
2769  Index ns = sgp.nelem();
2770  Index nb = bgp.nelem();
2771  Index np = pgp.nelem();
2772  Index nr = rgp.nelem();
2773  Index nc = cgp.nelem();
2774  assert(is_size(ia,
2775  ns,nb,np,nr,nc));
2776  assert(is_size(itw,
2777  ns,nb,np,nr,nc,
2778  32));
2779 
2780  // Check that interpolation weights are valid. The sum of all
2781  // weights (last dimension) must always be approximately one. We
2782  // only check the first element.
2783  assert( is_same_within_epsilon( itw(0,0,0,0,0,Range(joker)).sum(),
2784  1,
2785  sum_check_epsilon ) );
2786 
2787  // We have to loop all the points in the new grid:
2788  for ( Index is=0; is<ns; ++is )
2789  {
2790  const GridPos& ts = sgp[is];
2791  for ( Index ib=0; ib<nb; ++ib )
2792  {
2793  const GridPos& tb = bgp[ib];
2794  for ( Index ip=0; ip<np; ++ip )
2795  {
2796  const GridPos& tp = pgp[ip];
2797  for ( Index ir=0; ir<nr; ++ir )
2798  {
2799  const GridPos& tr = rgp[ir];
2800  for ( Index ic=0; ic<nc; ++ic )
2801  {
2802  // Current grid position:
2803  const GridPos& tc = cgp[ic];
2804 
2805  // Get handle to current element of output tensor and
2806  // initialize it to zero:
2807  Numeric& tia = ia(is,ib,ip,ir,ic);
2808  tia = 0;
2809 
2810  Index iti = 0;
2811  for ( Index s=0; s<2; ++s )
2812  for ( Index b=0; b<2; ++b )
2813  for ( Index p=0; p<2; ++p )
2814  for ( Index r=0; r<2; ++r )
2815  for ( Index c=0; c<2; ++c )
2816  {
2817  tia += a.get(ts.idx+s,
2818  tb.idx+b,
2819  tp.idx+p,
2820  tr.idx+r,
2821  tc.idx+c) * itw.get(is,ib,ip,ir,ic,
2822  iti);
2823  ++iti;
2824  }
2825  }
2826  }
2827  }
2828  }
2829  }
2830 }
2831 
2833 
2859  ConstTensor7View itw,
2860  ConstTensor6View a,
2861  const ArrayOfGridPos& vgp,
2862  const ArrayOfGridPos& sgp,
2863  const ArrayOfGridPos& bgp,
2864  const ArrayOfGridPos& pgp,
2865  const ArrayOfGridPos& rgp,
2866  const ArrayOfGridPos& cgp)
2867 {
2868  Index nv = vgp.nelem();
2869  Index ns = sgp.nelem();
2870  Index nb = bgp.nelem();
2871  Index np = pgp.nelem();
2872  Index nr = rgp.nelem();
2873  Index nc = cgp.nelem();
2874  assert(is_size(ia,
2875  nv,ns,nb,np,nr,nc));
2876  assert(is_size(itw,
2877  nv,ns,nb,np,nr,nc,
2878  64));
2879 
2880  // Check that interpolation weights are valid. The sum of all
2881  // weights (last dimension) must always be approximately one. We
2882  // only check the first element.
2883  assert( is_same_within_epsilon( itw(0,0,0,0,0,0,Range(joker)).sum(),
2884  1,
2885  sum_check_epsilon ) );
2886 
2887  // We have to loop all the points in the new grid:
2888  for ( Index iv=0; iv<nv; ++iv )
2889  {
2890  const GridPos& tv = vgp[iv];
2891  for ( Index is=0; is<ns; ++is )
2892  {
2893  const GridPos& ts = sgp[is];
2894  for ( Index ib=0; ib<nb; ++ib )
2895  {
2896  const GridPos& tb = bgp[ib];
2897  for ( Index ip=0; ip<np; ++ip )
2898  {
2899  const GridPos& tp = pgp[ip];
2900  for ( Index ir=0; ir<nr; ++ir )
2901  {
2902  const GridPos& tr = rgp[ir];
2903  for ( Index ic=0; ic<nc; ++ic )
2904  {
2905  // Current grid position:
2906  const GridPos& tc = cgp[ic];
2907 
2908  // Get handle to current element of output tensor and
2909  // initialize it to zero:
2910  Numeric& tia = ia(iv,is,ib,ip,ir,ic);
2911  tia = 0;
2912 
2913  Index iti = 0;
2914  for ( Index v=0; v<2; ++v )
2915  for ( Index s=0; s<2; ++s )
2916  for ( Index b=0; b<2; ++b )
2917  for ( Index p=0; p<2; ++p )
2918  for ( Index r=0; r<2; ++r )
2919  for ( Index c=0; c<2; ++c )
2920  {
2921  tia += a.get(tv.idx+v,
2922  ts.idx+s,
2923  tb.idx+b,
2924  tp.idx+p,
2925  tr.idx+r,
2926  tc.idx+c) * itw.get(iv,is,ib,ip,ir,ic,
2927  iti);
2928  ++iti;
2929  }
2930  }
2931  }
2932  }
2933  }
2934  }
2935  }
2936 }
2937 
2938 
2940 
2956  ConstVectorView y,
2957  const Numeric& x_i,
2958  const GridPos& gp)
2959 {
2960  Index N_x = x.nelem();
2961 
2962  assert(N_x == y.nelem());
2963  assert(N_x > 2);
2964 
2965  Vector xa(4), ya(4);
2966  Numeric y_int;
2967  Numeric dy_int;
2968  y_int = 0.;
2969 
2970  // 1 - polynomial interpolation (3 points) with grid position search
2971  // 2 - polynomial interpolation (3 points) without grid position search
2972  // 3 - polynomial interpolation (4 points)
2973 
2974  Index interp_method = 1;
2975 
2976  if (interp_method == 1)
2977  {
2978 
2979  // Pick out three points for interpolation
2980  if((gp.fd[0] <= 0.5 && gp.idx > 0) || gp.idx == N_x-2 )
2981  {
2982  xa[0] = x[gp.idx - 1];
2983  xa[1] = x[gp.idx];
2984  xa[2] = x[gp.idx + 1];
2985 
2986  ya[0] = y[gp.idx - 1];
2987  ya[1] = y[gp.idx];
2988  ya[2] = y[gp.idx + 1];
2989  }
2990 
2991  else if((gp.fd[0] > 0.5 && gp.idx < N_x-2) || gp.idx == 0 )
2992  {
2993  xa[0] = x[gp.idx];
2994  xa[1] = x[gp.idx + 1];
2995  xa[2] = x[gp.idx + 2];
2996 
2997  ya[0] = y[gp.idx];
2998  ya[1] = y[gp.idx + 1];
2999  ya[2] = y[gp.idx + 2];
3000  }
3001 
3002  else if(gp.idx == N_x-1)
3003  {
3004  xa[0] = x[N_x - 2];
3005  xa[1] = x[N_x - 1];
3006  xa[2] = x[N_x];
3007 
3008  ya[0] = y[N_x - 2];
3009  ya[1] = y[N_x - 1];
3010  ya[2] = y[N_x];
3011  }
3012  else
3013  {
3014  assert(false);
3015  arts_exit();
3016  }
3017 
3018  polint(y_int, dy_int, xa, ya, 3, x_i);
3019 
3020  }
3021 
3022  else if (interp_method == 2)
3023  {
3024  if( gp.idx == 0 )
3025  {
3026  xa[0] = x[gp.idx];
3027  xa[1] = x[gp.idx + 1];
3028  xa[2] = x[gp.idx + 2];
3029 
3030  ya[0] = y[gp.idx];
3031  ya[1] = y[gp.idx + 1];
3032  ya[2] = y[gp.idx + 2];
3033  }
3034  else if(gp.idx == N_x-1)
3035  {
3036  xa[0] = x[gp.idx - 2];
3037  xa[1] = x[gp.idx - 1];
3038  xa[2] = x[gp.idx];
3039 
3040  ya[0] = y[gp.idx - 2];
3041  ya[1] = y[gp.idx - 1];
3042  ya[2] = y[gp.idx];
3043  }
3044  else
3045  {
3046  xa[0] = x[gp.idx - 1];
3047  xa[1] = x[gp.idx];
3048  xa[2] = x[gp.idx + 1];
3049 
3050  ya[0] = y[gp.idx - 1];
3051  ya[1] = y[gp.idx];
3052  ya[2] = y[gp.idx + 1];
3053  }
3054 
3055  // Polinominal interpolation, n = 3
3056  polint(y_int, dy_int, xa, ya, 3, x_i);
3057  }
3058 
3059  else if (interp_method == 3)
3060  {
3061  // Take 4 points
3062  if( gp.idx == 0 )
3063  {
3064  xa[0] = - x[gp.idx + 1];
3065  xa[1] = x[gp.idx + 0];
3066  xa[2] = x[gp.idx + 1];
3067  xa[3] = x[gp.idx + 2];
3068 
3069  ya[0] = y[gp.idx + 1];
3070  ya[1] = y[gp.idx + 0];
3071  ya[2] = y[gp.idx + 1];
3072  ya[3] = y[gp.idx + 2];
3073  }
3074  else if(gp.idx == N_x-1)
3075  {
3076  xa[0] = x[gp.idx - 1];
3077  xa[1] = x[gp.idx - 0];
3078  xa[2] = 2*x[gp.idx] - x[gp.idx-1];
3079  xa[3] = 2*x[gp.idx] - x[gp.idx-2];
3080 
3081  ya[0] = y[gp.idx - 1];
3082  ya[1] = y[gp.idx - 0];
3083  ya[2] = y[gp.idx - 1];
3084  ya[3] = y[gp.idx - 2];
3085  }
3086  else if(gp.idx == N_x-2)
3087  {
3088  xa[0] = x[gp.idx - 2];
3089  xa[1] = x[gp.idx - 1];
3090  xa[2] = x[gp.idx ];
3091  xa[3] = x[gp.idx + 1];
3092 
3093  ya[0] = y[gp.idx - 2];
3094  ya[1] = y[gp.idx - 1];
3095  ya[2] = y[gp.idx];
3096  ya[3] = y[gp.idx + 1];
3097  }
3098  else
3099  {
3100  xa[0] = x[gp.idx - 1];
3101  xa[1] = x[gp.idx];
3102  xa[2] = x[gp.idx + 1];
3103  xa[3] = x[gp.idx + 2];
3104 
3105  ya[0] = y[gp.idx - 1];
3106  ya[1] = y[gp.idx];
3107  ya[2] = y[gp.idx + 1];
3108  ya[3] = y[gp.idx + 2];
3109  }
3110  // Polinominal interpolation, n = 4
3111  polint(y_int, dy_int, xa, ya, 4, x_i);
3112  }
3113 
3114 
3115  return y_int;
3116 }
3117 
3118 
3119 
3120 
3122 
3141 void polint(Numeric& y_int,
3142  Numeric& dy_int,
3143  ConstVectorView xa,
3144  ConstVectorView ya,
3145  const Index& n,
3146  const Numeric& x)
3147 {
3148  Index ns = 1;
3149  Numeric den, dif, dift, ho, hp, w;
3150 
3151  dif = abs(x-xa[0]);
3152 
3153  Vector c(n);
3154  Vector d(n);
3155 
3156  // Find the index of the closest table entry
3157  for(Index i=0; i<n; i++)
3158  {
3159  if( (dift = abs(x-xa[i])) < dif)
3160  {
3161  ns = i;
3162  dif = dift;
3163  }
3164  // Initialize c and d
3165  c[i] = ya[i];
3166  d[i] = ya[i];
3167  }
3168  // Initial approximation to y
3169  y_int = ya[ns--];
3170 
3171  for(Index m=1; m<n; m++)
3172  {
3173  for(Index i=0; i < n-m; i++)
3174  {
3175  ho = xa[i] - x;
3176  hp = xa[i+m] - x;
3177  w = c[i+1] - d[i];
3178  den = ho-hp;
3179  // This error occurs when two input xa's are identical.
3180  assert(den != 0.);
3181  den = w/den;
3182  d[i] = hp * den;
3183  c[i] = ho * den;
3184  }
3185  y_int += (dy_int = (2*(ns+1) < (n-m) ? c[ns+1] : d[ns--] ));
3186  }
3187 }
3188 
3189 
3190 
3191 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
The VectorView class.
Definition: matpackI.h:372
The Tensor4View class.
Definition: matpackIV.h:243
Index nelem() const
Number of elements.
Definition: array.h:176
void arts_exit(int status)
This is the exit function of ARTS.
Definition: arts.cc:42
A constant view of a Tensor7.
Definition: matpackVII.h:162
The Tensor7View class.
Definition: matpackVII.h:780
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
Numeric get(Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVI.h:386
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
A constant view of a Tensor6.
Definition: matpackVI.h:159
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:275
The range class.
Definition: matpackI.h:148
Numeric get(Index r, Index c) const
Get element implementation without assertions.
Definition: matpackI.h:618
Numeric fractional_gp(const GridPos &gp)
fractional_gp
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
Numeric get(Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVI.h:662
The Tensor6View class.
Definition: matpackVI.h:449
Numeric get(Index r, Index c) const
Get element implementation without assertions.
Definition: matpackI.h:691
Header file for interpolation.cc.
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Numeric get(Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackV.h:322
Structure to store a grid position.
Definition: interpolation.h:74
A constant view of a Tensor4.
Definition: matpackIV.h:141
void polint(Numeric &y_int, Numeric &dy_int, ConstVectorView xa, ConstVectorView ya, const Index &n, const Numeric &x)
Polynomial interpolation.
This file contains the definition of Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
ostream & operator<<(ostream &os, const GridPos &gp)
Output operator for GridPos.
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:446
Numeric fd[2]
Definition: interpolation.h:76
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:186
The Tensor3View class.
Definition: matpackIII.h:232
Numeric get(Index n) const
Get element implementation without assertions.
Definition: matpackI.h:311
Numeric get(Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIII.h:252
const Numeric sum_check_epsilon
The maximum difference from 1 that we allow for a sum check.
#define LOOPIT(x)
Macro for interpolation weight loops.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
A constant view of a Tensor5.
Definition: matpackV.h:152
#define abs(x)
Definition: continua.cc:20458
const Joker joker
Index idx
Definition: interpolation.h:75
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Numeric get(Index l, Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVII.h:1315
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Numeric get(Index n) const
Get element implementation without assertions.
Definition: matpackI.h:387
The Tensor5View class.
Definition: matpackV.h:276
Header file for logic.cc.
#define ns
Definition: continua.cc:21931
Numeric get(Index l, Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVII.h:713
A constant view of a Tensor3.
Definition: matpackIII.h:139
A constant view of a Vector.
Definition: matpackI.h:292
A constant view of a Matrix.
Definition: matpackI.h:596
Numeric get(Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIV.h:272
Numeric get(Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackV.h:216
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:324
void gridpos_1to1(ArrayOfGridPos &gp, ConstVectorView grid)
gridpos_1to1
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:91
const Numeric FD_TOL
Allowed tolerance for fractional distance values.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric get(Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIII.h:176
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
Numeric get(Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIV.h:185