ARTS 2.5.9 (git: 825fa5f2)
special_interp.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
62#include <cmath>
63#include <iostream>
64#include <stdexcept>
65#include "auto_md.h"
66#include "check_input.h"
67#include "math_funcs.h"
68#include "messages.h"
69#include "special_interp.h"
70
71/*===========================================================================
72 === Point interpolation functions for atmospheric grids and fields
73 ===========================================================================*/
74
76 const Index& atmosphere_dim,
77 const ArrayOfGridPos& gp_p,
78 const ArrayOfGridPos& gp_lat,
79 const ArrayOfGridPos& gp_lon) {
80 const Index n = gp_p.nelem();
81
82 if (atmosphere_dim == 1) {
83 itw.resize(n, 2);
84 interpweights(itw, gp_p);
85 }
86
87 else if (atmosphere_dim == 2) {
88 ARTS_ASSERT(gp_lat.nelem() == n);
89 itw.resize(n, 4);
90 interpweights(itw, gp_p, gp_lat);
91 }
92
93 else if (atmosphere_dim == 3) {
94 ARTS_ASSERT(gp_lat.nelem() == n);
95 ARTS_ASSERT(gp_lon.nelem() == n);
96 itw.resize(n, 8);
97 interpweights(itw, gp_p, gp_lat, gp_lon);
98 }
99}
100
102 const Index& atmosphere_dim,
103 ConstTensor3View x_field,
104 const ArrayOfGridPos& gp_p,
105 const ArrayOfGridPos& gp_lat,
106 const ArrayOfGridPos& gp_lon,
107 ConstMatrixView itw) {
108 ARTS_ASSERT(x.nelem() == gp_p.nelem());
109
110 if (atmosphere_dim == 1) {
111 ARTS_ASSERT(itw.ncols() == 2);
112 interp(x, itw, x_field(Range(joker), 0, 0), gp_p);
113 }
114
115 else if (atmosphere_dim == 2) {
116 ARTS_ASSERT(itw.ncols() == 4);
117 interp(x, itw, x_field(Range(joker), Range(joker), 0), gp_p, gp_lat);
118 }
119
120 else if (atmosphere_dim == 3) {
121 ARTS_ASSERT(itw.ncols() == 8);
122 interp(x, itw, x_field, gp_p, gp_lat, gp_lon);
123 }
124}
125
127 const Index& atmosphere_dim,
128 ConstTensor3View x_field,
129 const ArrayOfGridPos& gp_p,
130 const ArrayOfGridPos& gp_lat,
131 const ArrayOfGridPos& gp_lon) {
132 Matrix itw;
133
134 interp_atmfield_gp2itw(itw, atmosphere_dim, gp_p, gp_lat, gp_lon);
135
136 interp_atmfield_by_itw(x, atmosphere_dim, x_field, gp_p, gp_lat, gp_lon, itw);
137}
138
139Numeric interp_atmfield_by_gp(const Index& atmosphere_dim,
140 ConstTensor3View x_field,
141 const GridPos& gp_p,
142 const GridPos& gp_lat,
143 const GridPos& gp_lon) {
144 ArrayOfGridPos agp_p(1), agp_lat(0), agp_lon(0);
145
146 gridpos_copy(agp_p[0], gp_p);
147
148 if (atmosphere_dim > 1) {
149 agp_lat.resize(1);
150 gridpos_copy(agp_lat[0], gp_lat);
151 }
152
153 if (atmosphere_dim > 2) {
154 agp_lon.resize(1);
155 gridpos_copy(agp_lon[0], gp_lon);
156 }
157
158 Vector x(1);
159
160 interp_atmfield_by_gp(x, atmosphere_dim, x_field, agp_p, agp_lat, agp_lon);
161
162 return x[0];
163}
164
166 GridPos& gp_p_out,
167 GridPos& gp_lat_out,
168 GridPos& gp_lon_out,
169 const GridPos& gp_p_in,
170 const GridPos& gp_lat_in,
171 const GridPos& gp_lon_in,
172 const Index& atmosphere_dim,
173 const ArrayOfIndex& cloudbox_limits) {
174 // Shift grid positions to cloud box grids
175 if (atmosphere_dim == 1) {
176 gridpos_copy(gp_p_out, gp_p_in);
177 gp_p_out.idx -= cloudbox_limits[0];
178 gridpos_upperend_check(gp_p_out, cloudbox_limits[1] - cloudbox_limits[0]);
179 ARTS_ASSERT(itw.nelem() == 2);
180 interpweights(itw, gp_p_out);
181 } else if (atmosphere_dim == 2) {
182 gridpos_copy(gp_p_out, gp_p_in);
183 gridpos_copy(gp_lat_out, gp_lat_in);
184 gp_p_out.idx -= cloudbox_limits[0];
185 gp_lat_out.idx -= cloudbox_limits[2];
186 gridpos_upperend_check(gp_p_out, cloudbox_limits[1] - cloudbox_limits[0]);
187 gridpos_upperend_check(gp_lat_out, cloudbox_limits[3] - cloudbox_limits[2]);
188 ARTS_ASSERT(itw.nelem() == 4);
189 interpweights(itw, gp_p_out, gp_lat_out);
190 } else {
191 gridpos_copy(gp_p_out, gp_p_in);
192 gridpos_copy(gp_lat_out, gp_lat_in);
193 gridpos_copy(gp_lon_out, gp_lon_in);
194 gp_p_out.idx -= cloudbox_limits[0];
195 gp_lat_out.idx -= cloudbox_limits[2];
196 gp_lon_out.idx -= cloudbox_limits[4];
197 gridpos_upperend_check(gp_p_out, cloudbox_limits[1] - cloudbox_limits[0]);
198 gridpos_upperend_check(gp_lat_out, cloudbox_limits[3] - cloudbox_limits[2]);
199 gridpos_upperend_check(gp_lon_out, cloudbox_limits[5] - cloudbox_limits[4]);
200 ARTS_ASSERT(itw.nelem() == 8);
201 interpweights(itw, gp_p_out, gp_lat_out, gp_lon_out);
202 }
203}
204
225 const Index& atmosphere_dim,
226 const ArrayOfGridPos& gp_lat,
227 const ArrayOfGridPos& gp_lon) {
228 if (atmosphere_dim == 1) {
229 itw.resize(1, 1);
230 itw = 1;
231 }
232
233 else if (atmosphere_dim == 2) {
234 const Index n = gp_lat.nelem();
235 itw.resize(n, 2);
236 interpweights(itw, gp_lat);
237 }
238
239 else if (atmosphere_dim == 3) {
240 const Index n = gp_lat.nelem();
241 ARTS_ASSERT(n == gp_lon.nelem());
242 itw.resize(n, 4);
243 interpweights(itw, gp_lat, gp_lon);
244 }
245}
246
248 const Index& atmosphere_dim,
249 ConstMatrixView x_surface,
250 const ArrayOfGridPos& gp_lat,
251 const ArrayOfGridPos& gp_lon,
252 ConstMatrixView itw) {
253 if (atmosphere_dim == 1) {
254 ARTS_ASSERT(itw.ncols() == 1);
255 x = x_surface(0, 0);
256 }
257
258 else if (atmosphere_dim == 2) {
259 ARTS_ASSERT(x.nelem() == gp_lat.nelem());
260 ARTS_ASSERT(itw.ncols() == 2);
261 interp(x, itw, x_surface(Range(joker), 0), gp_lat);
262 }
263
264 else if (atmosphere_dim == 3) {
265 ARTS_ASSERT(x.nelem() == gp_lat.nelem());
266 ARTS_ASSERT(itw.ncols() == 4);
267 interp(x, itw, x_surface, gp_lat, gp_lon);
268 }
269}
270
272 const Index& atmosphere_dim,
273 ConstMatrixView x_surface,
274 const ArrayOfGridPos& gp_lat,
275 const ArrayOfGridPos& gp_lon) {
276 Matrix itw;
277
278 interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
279
280 interp_atmsurface_by_itw(x, atmosphere_dim, x_surface, gp_lat, gp_lon, itw);
281}
282
284 ConstMatrixView x_surface,
285 const GridPos& gp_lat,
286 const GridPos& gp_lon) {
287 ArrayOfGridPos agp_lat(0), agp_lon(0);
288
289 if (atmosphere_dim > 1) {
290 agp_lat.resize(1);
291 gridpos_copy(agp_lat[0], gp_lat);
292 }
293
294 if (atmosphere_dim > 2) {
295 agp_lon.resize(1);
296 gridpos_copy(agp_lon[0], gp_lon);
297 }
298
299 Vector x(1);
300
301 interp_atmsurface_by_gp(x, atmosphere_dim, x_surface, agp_lat, agp_lon);
302
303 return x[0];
304}
305
306/*===========================================================================
307 === Regridding
308 ===========================================================================*/
309
311 const Index& atmosphere_dim,
312 ConstTensor3View field_old,
313 const ArrayOfGridPos& gp_p,
314 const ArrayOfGridPos& gp_lat,
315 const ArrayOfGridPos& gp_lon) {
316 const Index n1 = gp_p.nelem();
317
318 if (atmosphere_dim == 1) {
319 field_new.resize(n1, 1, 1);
320 Matrix itw(n1, 2);
321 interpweights(itw, gp_p);
322 interp(field_new(joker, 0, 0), itw, field_old(joker, 0, 0), gp_p);
323 } else if (atmosphere_dim == 2) {
324 const Index n2 = gp_lat.nelem();
325 field_new.resize(n1, n2, 1);
326 Tensor3 itw(n1, n2, 4);
327 interpweights(itw, gp_p, gp_lat);
328 interp(field_new(joker, joker, 0),
329 itw,
330 field_old(joker, joker, 0),
331 gp_p,
332 gp_lat);
333 } else if (atmosphere_dim == 3) {
334 const Index n2 = gp_lat.nelem();
335 const Index n3 = gp_lon.nelem();
336 field_new.resize(n1, n2, n3);
337 Tensor4 itw(n1, n2, n3, 8);
338 interpweights(itw, gp_p, gp_lat, gp_lon);
339 interp(field_new, itw, field_old, gp_p, gp_lat, gp_lon);
340 }
341}
342
344 const Index& atmosphere_dim,
345 ConstMatrixView field_old,
346 const ArrayOfGridPos& gp_lat,
347 const ArrayOfGridPos& gp_lon) {
348 if (atmosphere_dim == 1) {
349 field_new = field_old;
350 } else if (atmosphere_dim == 2) {
351 const Index n1 = gp_lat.nelem();
352 field_new.resize(n1, 1);
353 Matrix itw(n1, 2);
354 interpweights(itw, gp_lat);
355 interp(field_new(joker, 0), itw, field_old(joker, 0), gp_lat);
356 } else if (atmosphere_dim == 3) {
357 const Index n1 = gp_lat.nelem();
358 const Index n2 = gp_lon.nelem();
359 field_new.resize(n1, n2);
360 Tensor3 itw(n1, n2, 4);
361 interpweights(itw, gp_lat, gp_lon);
362 interp(field_new, itw, field_old, gp_lat, gp_lon);
363 }
364}
365
367 ArrayOfGridPos& gp_lat,
368 ArrayOfGridPos& gp_lon,
369 const RetrievalQuantity& rq,
370 const Index& atmosphere_dim,
371 const Vector& p_grid,
372 const Vector& lat_grid,
373 const Vector& lon_grid) {
374 gp_p.resize(rq.Grids()[0].nelem());
375 p2gridpos(gp_p, p_grid, rq.Grids()[0], 0);
376 //
377 if (atmosphere_dim >= 2) {
378 gp_lat.resize(rq.Grids()[1].nelem());
379 gridpos(gp_lat, lat_grid, rq.Grids()[1], 0);
380 } else {
381 gp_lat.resize(0);
382 }
383 //
384 if (atmosphere_dim >= 3) {
385 gp_lon.resize(rq.Grids()[2].nelem());
386 gridpos(gp_lon, lon_grid, rq.Grids()[2], 0);
387 } else {
388 gp_lon.resize(0);
389 }
390}
391
393 ArrayOfGridPos& gp_lon,
394 const RetrievalQuantity& rq,
395 const Index& atmosphere_dim,
396 const Vector& lat_grid,
397 const Vector& lon_grid) {
398 if (atmosphere_dim >= 2) {
399 gp_lat.resize(rq.Grids()[0].nelem());
400 gridpos(gp_lat, lat_grid, rq.Grids()[0], 0);
401 } else {
402 gp_lat.resize(0);
403 }
404 //
405 if (atmosphere_dim >= 3) {
406 gp_lon.resize(rq.Grids()[1].nelem());
407 gridpos(gp_lon, lon_grid, rq.Grids()[1], 0);
408 } else {
409 gp_lon.resize(0);
410 }
411}
412
414 ArrayOfGridPos& gp_lat,
415 ArrayOfGridPos& gp_lon,
416 Index& n_p,
417 Index& n_lat,
418 Index& n_lon,
419 const ArrayOfVector& ret_grids,
420 const Index& atmosphere_dim,
421 const Vector& p_grid,
422 const Vector& lat_grid,
423 const Vector& lon_grid) {
424 // We want here an extrapolation to infinity ->
425 // extremly high extrapolation factor
426 const Numeric inf_proxy = 1.0e99;
427
428 gp_p.resize(p_grid.nelem());
429 n_p = ret_grids[0].nelem();
430 if (n_p > 1) {
431 p2gridpos(gp_p, ret_grids[0], p_grid, inf_proxy);
433 } else {
434 gp4length1grid(gp_p);
435 }
436
437 if (atmosphere_dim >= 2) {
438 gp_lat.resize(lat_grid.nelem());
439 n_lat = ret_grids[1].nelem();
440 if (n_lat > 1) {
441 gridpos(gp_lat, ret_grids[1], lat_grid, inf_proxy);
443 } else {
444 gp4length1grid(gp_lat);
445 }
446 } else {
447 gp_lat.resize(0);
448 n_lat = 1;
449 }
450 //
451 if (atmosphere_dim >= 3) {
452 gp_lon.resize(lon_grid.nelem());
453 n_lon = ret_grids[2].nelem();
454 if (n_lon > 1) {
455 gridpos(gp_lon, ret_grids[2], lon_grid, inf_proxy);
457 } else {
458 gp4length1grid(gp_lon);
459 }
460 } else {
461 gp_lon.resize(0);
462 n_lon = 1;
463 }
464}
465
467 ArrayOfGridPos& gp_lon,
468 Index& n_lat,
469 Index& n_lon,
470 const ArrayOfVector& ret_grids,
471 const Index& atmosphere_dim,
472 const Vector& lat_grid,
473 const Vector& lon_grid) {
474 // We want here an extrapolation to infinity ->
475 // extremly high extrapolation factor
476 const Numeric inf_proxy = 1.0e99;
477
478 if (atmosphere_dim >= 2) {
479 gp_lat.resize(lat_grid.nelem());
480 n_lat = ret_grids[0].nelem();
481 if (n_lat > 1) {
482 gridpos(gp_lat, ret_grids[0], lat_grid, inf_proxy);
484 } else {
485 gp4length1grid(gp_lat);
486 }
487 } else {
488 gp_lat.resize(0);
489 n_lat = 1;
490 }
491 //
492 if (atmosphere_dim >= 3) {
493 gp_lon.resize(lon_grid.nelem());
494 n_lon = ret_grids[1].nelem();
495 if (n_lon > 1) {
496 gridpos(gp_lon, ret_grids[1], lon_grid, inf_proxy);
498 } else {
499 gp4length1grid(gp_lon);
500 }
501 } else {
502 gp_lon.resize(0);
503 n_lon = 1;
504 }
505}
506
508 const Index& atmosphere_dim,
509 ConstTensor3View field_old,
510 const ArrayOfGridPos& gp_p,
511 const ArrayOfGridPos& gp_lat,
512 const ArrayOfGridPos& gp_lon) {
513 const Index n1 = gp_p.nelem();
514
515 const bool np_is1 = field_old.npages() == 1 ? true : false;
516 const bool nlat_is1 =
517 atmosphere_dim > 1 && field_old.nrows() == 1 ? true : false;
518 const bool nlon_is1 =
519 atmosphere_dim > 2 && field_old.ncols() == 1 ? true : false;
520
521 // If no length 1, we can use standard function
522 if (!np_is1 && !nlat_is1 && !nlon_is1) {
524 field_new, atmosphere_dim, field_old, gp_p, gp_lat, gp_lon);
525 } else {
526 //--- 1D (1 possibilities left) -------------------------------------------
527 if (atmosphere_dim == 1) { // 1: No interpolation at all
528 field_new.resize(n1, 1, 1);
529 field_new(joker, 0, 0) = field_old(0, 0, 0);
530 }
531
532 //--- 2D (3 possibilities left) -------------------------------------------
533 else if (atmosphere_dim == 2) {
534 const Index n2 = gp_lat.nelem();
535 field_new.resize(n1, n2, 1);
536 //
537 if (np_is1 && nlat_is1) // 1: No interpolation at all
538 {
539 // Here we need no interpolation at all
540 field_new(joker, joker, 0) = field_old(0, 0, 0);
541 } else if (np_is1) // 2: Latitude interpolation
542 {
543 Matrix itw(n2, 2);
544 interpweights(itw, gp_lat);
545 Vector tmp(n2);
546 interp(tmp, itw, field_old(0, joker, 0), gp_lat);
547 for (Index p = 0; p < n1; p++) {
548 ARTS_ASSERT(gp_p[p].fd[0] < 1e-6);
549 field_new(p, joker, 0) = tmp;
550 }
551 } else // 3: Pressure interpolation
552 {
553 Matrix itw(n1, 2);
554 interpweights(itw, gp_p);
555 Vector tmp(n1);
556 interp(tmp, itw, field_old(joker, 0, 0), gp_p);
557 for (Index lat = 0; lat < n2; lat++) {
558 ARTS_ASSERT(gp_lat[lat].fd[0] < 1e-6);
559 field_new(joker, lat, 0) = tmp;
560 }
561 }
562 }
563
564 //--- 3D (7 possibilities left) -------------------------------------------
565 else if (atmosphere_dim == 3) {
566 const Index n2 = gp_lat.nelem();
567 const Index n3 = gp_lon.nelem();
568 field_new.resize(n1, n2, n3);
569 //
570 if (np_is1 && nlat_is1 && nlon_is1) // 1: No interpolation at all
571 {
572 field_new(joker, joker, joker) = field_old(0, 0, 0);
573 }
574
575 else if (np_is1) // No pressure interpolation --------------
576 {
577 if (nlat_is1) // 2: Just longitude interpolation
578 {
579 Matrix itw(n3, 2);
580 interpweights(itw, gp_lon);
581 Vector tmp(n3);
582 interp(tmp, itw, field_old(0, 0, joker), gp_lon);
583 for (Index p = 0; p < n1; p++) {
584 ARTS_ASSERT(gp_p[p].fd[0] < 1e-6);
585 for (Index lat = 0; lat < n2; lat++) {
586 ARTS_ASSERT(gp_lat[lat].fd[0] < 1e-6);
587 field_new(p, lat, joker) = tmp;
588 }
589 }
590 } else if (nlon_is1) // 3: Just latitude interpolation
591 {
592 Matrix itw(n2, 2);
593 interpweights(itw, gp_lat);
594 Vector tmp(n2);
595 interp(tmp, itw, field_old(0, joker, 0), gp_lat);
596 for (Index p = 0; p < n1; p++) {
597 ARTS_ASSERT(gp_p[p].fd[0] < 1e-6);
598 for (Index lon = 0; lon < n3; lon++) {
599 ARTS_ASSERT(gp_lon[lon].fd[0] < 1e-6);
600 field_new(p, joker, lon) = tmp;
601 }
602 }
603 } else // 4: Both lat and lon interpolation
604 {
605 Tensor3 itw(n2, n3, 4);
606 interpweights(itw, gp_lat, gp_lon);
607 Matrix tmp(n2, n3);
608 interp(tmp, itw, field_old(0, joker, joker), gp_lat, gp_lon);
609 for (Index p = 0; p < n1; p++) {
610 ARTS_ASSERT(gp_p[p].fd[0] < 1e-6);
611 field_new(p, joker, joker) = tmp;
612 }
613 }
614 }
615
616 else // Pressure interpolation --------------
617 {
618 if (nlat_is1 && nlon_is1) // 5: Just pressure interpolatiom
619 {
620 Matrix itw(n1, 2);
621 interpweights(itw, gp_p);
622 Vector tmp(n1);
623 interp(tmp, itw, field_old(joker, 0, 0), gp_p);
624 for (Index lat = 0; lat < n2; lat++) {
625 ARTS_ASSERT(gp_lat[lat].fd[0] < 1e-6);
626 for (Index lon = 0; lon < n3; lon++) {
627 ARTS_ASSERT(gp_lon[lon].fd[0] < 1e-6);
628 field_new(joker, lat, lon) = tmp;
629 }
630 }
631 } else if (nlat_is1) // 6: Both p and lon interpolation
632 {
633 Tensor3 itw(n1, n3, 4);
634 interpweights(itw, gp_p, gp_lon);
635 Matrix tmp(n1, n3);
636 interp(tmp, itw, field_old(joker, 0, joker), gp_p, gp_lon);
637 for (Index lat = 0; lat < n2; lat++) {
638 ARTS_ASSERT(gp_lat[lat].fd[0] < 1e-6);
639 field_new(joker, lat, joker) = tmp;
640 }
641 } else // 7: Both p and lat interpolation
642 {
643 Tensor3 itw(n1, n2, 4);
644 interpweights(itw, gp_p, gp_lat);
645 Matrix tmp(n1, n2);
646 interp(tmp, itw, field_old(joker, joker, 0), gp_p, gp_lat);
647 for (Index lon = 0; lon < n3; lon++) {
648 ARTS_ASSERT(gp_lon[lon].fd[0] < 1e-6);
649 field_new(joker, joker, lon) = tmp;
650 }
651 }
652 }
653 }
654 }
655}
656
657/* So far just a temporary test */
659 const Index& atmosphere_dim,
660 ConstMatrixView field_old,
661 const ArrayOfGridPos& gp_lat,
662 const ArrayOfGridPos& gp_lon) {
663 // As 1D is so simple, let's do it here and not go to standard function
664 if (atmosphere_dim == 1) {
665 field_new = field_old;
666 } else {
667 const bool nlat_is1 = field_old.nrows() == 1 ? true : false;
668 const bool nlon_is1 =
669 atmosphere_dim > 2 && field_old.ncols() == 1 ? true : false;
670
671 // If no length 1, we can use standard function
672 if (!nlat_is1 && !nlon_is1) {
674 field_new, atmosphere_dim, field_old, gp_lat, gp_lon);
675 } else {
676 if (atmosphere_dim == 2) { // 1: No interpolation at all
677 const Index n1 = gp_lat.nelem();
678 field_new.resize(n1, 1);
679 field_new(joker, 0) = field_old(0, 0);
680 } else {
681 const Index n1 = gp_lat.nelem();
682 const Index n2 = gp_lon.nelem();
683 field_new.resize(n1, n2);
684 //
685 if (nlat_is1 && nlon_is1) // 1: No interpolation at all
686 {
687 field_new(joker, joker) = field_old(0, 0);
688 } else if (nlon_is1) // 2: Just latitude interpolation
689 {
690 Matrix itw(n1, 2);
691 interpweights(itw, gp_lat);
692 Vector tmp(n1);
693 interp(tmp, itw, field_old(joker, 0), gp_lat);
694 for (Index lon = 0; lon < n2; lon++) {
695 ARTS_ASSERT(gp_lon[lon].fd[0] < 1e-6);
696 field_new(joker, lon) = tmp;
697 }
698 } else // 2: Just longitude interpolation
699 {
700 Matrix itw(n2, 2);
701 interpweights(itw, gp_lon);
702 Vector tmp(n2);
703 interp(tmp, itw, field_old(0, joker), gp_lon);
704 for (Index lat = 0; lat < n1; lat++) {
705 ARTS_ASSERT(gp_lat[lat].fd[0] < 1e-6);
706 field_new(lat, joker) = tmp;
707 }
708 }
709 }
710 }
711 }
712}
713
714/*===========================================================================
715 === Conversion altitudes / pressure
716 ===========================================================================*/
717
718void itw2p(VectorView p_values,
719 ConstVectorView p_grid,
720 const ArrayOfGridPos& gp,
721 ConstMatrixView itw) {
722 ARTS_ASSERT(itw.ncols() == 2);
723 ARTS_ASSERT(p_values.nelem() == itw.nrows());
724
725 // Local variable to store log of the pressure grid:
726 Vector logpgrid(p_grid.nelem());
727
728 transform(logpgrid, log, p_grid);
729
730 interp(p_values, itw, logpgrid, gp);
731
732 transform(p_values, exp, p_values);
733}
734
760 ConstVectorView old_pgrid,
761 ConstVectorView new_pgrid,
762 const Numeric& extpolfac) {
763 // Local variable to store log of the pressure grids
764 Vector logold(old_pgrid.nelem());
765 Vector lognew(new_pgrid.nelem());
766
767 transform(logold, log, old_pgrid);
768 transform(lognew, log, new_pgrid);
769
770 gridpos(gp, logold, lognew, extpolfac);
771}
772
774 GridPos& gp_lat,
775 GridPos& gp_lon,
776 const Index& atmosphere_dim,
777 ConstVectorView p_grid,
778 ConstVectorView lat_grid,
779 ConstVectorView lon_grid,
780 ConstTensor3View z_field,
781 ConstVectorView rte_pos) {
782 chk_rte_pos(atmosphere_dim, rte_pos);
783
784 if (atmosphere_dim == 1) {
786 "Altitude interpolation", z_field(joker, 0, 0), rte_pos[0]);
787 gridpos(gp_p, z_field(joker, 0, 0), rte_pos[0]);
788 } else {
789 // Determine z at lat/lon (z_grid) by blue interpolation
790 const Index np = p_grid.nelem();
791 Vector z_grid(np);
792 ArrayOfGridPos agp_z, agp_lat(np);
793 //
794 gridpos_1to1(agp_z, p_grid);
795 //
796 chk_interpolation_grids("Latitude interpolation", lat_grid, rte_pos[1]);
797 gridpos(gp_lat, lat_grid, rte_pos[1]);
798
799 if (atmosphere_dim == 2) {
800 for (Index i = 0; i < np; i++) {
801 agp_lat[i] = gp_lat;
802 }
803 Matrix itw(np, 4);
804 interpweights(itw, agp_z, agp_lat);
805 interp(z_grid, itw, z_field(joker, joker, 0), agp_z, agp_lat);
806 } else {
807 chk_interpolation_grids("Longitude interpolation", lon_grid, rte_pos[2]);
808 gridpos(gp_lon, lon_grid, rte_pos[2]);
809 ArrayOfGridPos agp_lon(np);
810 for (Index i = 0; i < np; i++) {
811 agp_lat[i] = gp_lat;
812 agp_lon[i] = gp_lon;
813 }
814 Matrix itw(np, 8);
815 interpweights(itw, agp_z, agp_lat, agp_lon);
816 interp(z_grid, itw, z_field, agp_z, agp_lat, agp_lon);
817 }
818
819 // And use z_grid to get gp_p (gp_al and gp_lon determined above)
820 chk_interpolation_grids("Altitude interpolation", z_grid, rte_pos[0]);
821 gridpos(gp_p, z_grid, rte_pos[0]);
822 }
823}
824
826 GridPos& gp_lon,
827 const Index& atmosphere_dim,
828 ConstVectorView lat_grid,
829 ConstVectorView lon_grid,
830 ConstVectorView rte_pos) {
831 chk_rte_pos(atmosphere_dim, rte_pos);
832
833 if (atmosphere_dim == 1) {
834 } else {
835 chk_interpolation_grids("Latitude interpolation", lat_grid, rte_pos[1]);
836 gridpos(gp_lat, lat_grid, rte_pos[1]);
837
838 if (atmosphere_dim == 3) {
839 chk_interpolation_grids("Longitude interpolation", lon_grid, rte_pos[2]);
840 gridpos(gp_lon, lon_grid, rte_pos[2]);
841 }
842 }
843}
844
846 ConstVectorView p_grid,
847// FIXME only used in assertion
848#ifndef NDEBUG
849 ConstVectorView lat_grid,
850#else
852#endif
853 ConstMatrixView z_field,
854 const GridPos& gp_lat) {
855 const Index np = p_grid.nelem();
856
857 ARTS_ASSERT(z.nelem() == np);
858 ARTS_ASSERT(z_field.nrows() == np);
859 ARTS_ASSERT(z_field.ncols() == lat_grid.nelem());
860
861 Matrix z_matrix(np, 1);
862 ArrayOfGridPos gp_z(np), agp_lat(1);
863 Tensor3 itw(np, 1, 4);
864
865 gridpos_copy(agp_lat[0], gp_lat);
866 gridpos(gp_z, p_grid, p_grid);
867 interpweights(itw, gp_z, agp_lat);
868 interp(z_matrix, itw, z_field, gp_z, agp_lat);
869
870 z = z_matrix(Range(joker), 0);
871}
872
874 ConstVectorView p_grid,
875//FIXME only used in assertion
876#ifndef NDEBUG
877 ConstVectorView lat_grid,
878 ConstVectorView lon_grid,
879#else
882#endif
883 ConstTensor3View z_field,
884 const GridPos& gp_lat,
885 const GridPos& gp_lon) {
886 const Index np = p_grid.nelem();
887
888 ARTS_ASSERT(z.nelem() == np);
889 ARTS_ASSERT(z_field.npages() == np);
890 ARTS_ASSERT(z_field.nrows() == lat_grid.nelem());
891 ARTS_ASSERT(z_field.ncols() == lon_grid.nelem());
892
893 Tensor3 z_tensor(np, 1, 1);
894 ArrayOfGridPos agp_z(np), agp_lat(1), agp_lon(1);
895 Tensor4 itw(np, 1, 1, 8);
896
897 gridpos_copy(agp_lat[0], gp_lat);
898 gridpos_copy(agp_lon[0], gp_lon);
899 gridpos(agp_z, p_grid, p_grid);
900 interpweights(itw, agp_z, agp_lat, agp_lon);
901
902 interp(z_tensor, itw, z_field, agp_z, agp_lat, agp_lon);
903
904 z = z_tensor(Range(joker), 0, 0);
905}
906
907/*===========================================================================
908 === Interpolation of GriddedFields
909 ===========================================================================*/
910
912 MatrixView n_imag,
913 const GriddedField3& complex_n,
914 const String& varname,
915 ConstVectorView f_grid,
916 ConstVectorView t_grid) {
917 // Set expected order of grids
918 Index gfield_fID = 0;
919 Index gfield_tID = 1;
920 Index gfield_compID = 2;
921
922 // Check of complex_n
923 //
924 complex_n.checksize_strict();
925 //
926 chk_griddedfield_gridname(complex_n, gfield_fID, "Frequency");
927 chk_griddedfield_gridname(complex_n, gfield_tID, "Temperature");
928 chk_griddedfield_gridname(complex_n, gfield_compID, "Complex");
929 //
930 if (complex_n.data.ncols() != 2) {
931 ostringstream os;
932 os << "The data in *" << varname
933 << "* must have exactly two pages. One page "
934 << "each\nfor the real and imaginary part of the complex refractive index.";
935 }
936
937 // Frequency and temperature grid sizes
938 const Index nf_in = complex_n.data.npages();
939 const Index nt_in = complex_n.data.nrows();
940 const Index nf_out = f_grid.nelem();
941 const Index nt_out = t_grid.nelem();
942
943 //Assert size of output variables
944 ARTS_ASSERT(n_real.nrows() == nf_out && n_real.ncols() == nt_out);
945 ARTS_ASSERT(n_imag.nrows() == nf_out && n_imag.ncols() == nt_out);
946
947 const Vector& f_grid_in = complex_n.get_numeric_grid(gfield_fID);
948 const Vector& t_grid_in = complex_n.get_numeric_grid(gfield_tID);
949
950 // Expand/interpolate in frequency dimension
951 //
952 Matrix nrf(nf_out, nt_in), nif(nf_out, nt_in);
953 //
954 if (nf_in == 1) {
955 for (Index i = 0; i < nf_out; i++) {
956 nrf(i, joker) = complex_n.data(0, joker, 0);
957 nif(i, joker) = complex_n.data(0, joker, 1);
958 }
959 } else {
960 chk_interpolation_grids("Frequency interpolation", f_grid_in, f_grid);
961 //
962 ArrayOfGridPos gp(nf_out);
963 Matrix itw(nf_out, 2);
964 gridpos(gp, f_grid_in, f_grid);
965 interpweights(itw, gp);
966 for (Index i = 0; i < nt_in; i++) {
967 interp(nrf(joker, i), itw, complex_n.data(joker, i, 0), gp);
968 interp(nif(joker, i), itw, complex_n.data(joker, i, 1), gp);
969 }
970 }
971
972 // Expand/interpolate in temperature dimension
973 //
974 if (nt_in == 1) {
975 for (Index i = 0; i < nt_out; i++) {
976 n_real(joker, i) = nrf(joker, 0);
977 n_imag(joker, i) = nif(joker, 0);
978 }
979 } else {
980 chk_interpolation_grids("Temperature interpolation", t_grid_in, t_grid);
981 //
982 ArrayOfGridPos gp(nt_out);
983 Matrix itw(nt_out, 2);
984 gridpos(gp, t_grid_in, t_grid);
985 interpweights(itw, gp);
986 for (Index i = 0; i < nf_out; i++) {
987 interp(n_real(i, joker), itw, nrf(i, joker), gp);
988 interp(n_imag(i, joker), itw, nif(i, joker), gp);
989 }
990 }
991}
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:906
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
void chk_griddedfield_gridname(const GriddedField &gf, const Index gridindex, const String &gridname)
Check name of grid in GriddedField.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Tensor3.
Definition: matpackIII.h:133
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
void checksize_strict() const final
Strict consistency check.
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
The range class.
Definition: matpackI.h:160
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:325
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:408
The Tensor3 class.
Definition: matpackIII.h:352
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:435
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
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
void jacobian_type_extrapol(ArrayOfGridPos &gp)
Adopts grid positions to extrapolation used for jacobians.
Definition: jacobian.cc:836
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1433
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
void complex_n_interp(MatrixView n_real, MatrixView n_imag, const GriddedField3 &complex_n, const String &varname, ConstVectorView f_grid, ConstVectorView t_grid)
General function for interpolating data of complex n type.
void regrid_atmfield_by_gp(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric field, for precalculated grid positions.
void interp_atmsurface_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of a surface-type variable.
void regrid_atmfield_by_gp_oem(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of atmospheric field OEM-type.
void get_gp_atmsurf_to_rq(ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric surfaces to retrieval grids.
void regrid_atmsurf_by_gp(Matrix &field_new, const Index &atmosphere_dim, ConstMatrixView field_old, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric surface, for precalculated grid positions.
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void z_at_latlon(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, const GridPos &gp_lat, const GridPos &gp_lon)
Returns the geomtrical altitudes of p_grid for one latitude and one longitude.
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
Converts atmospheric a grid position to weights for interpolation of a field defined ONLY inside the ...
void get_gp_rq_to_atmgrids(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, Index &n_p, Index &n_lat, Index &n_lon, const ArrayOfVector &ret_grids, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
void regrid_atmsurf_by_gp_oem(Matrix &field_new, const Index &atmosphere_dim, ConstMatrixView field_old, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of surface field OEM-type.
void get_gp_atmgrids_to_rq(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
void interp_atmfield_by_gp(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates an atmospheric field given the grid positions.
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
Converts a geographical position (rte_pos) to grid positions for p, lat and lon.
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
void z_at_lat_2d(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstMatrixView z_field, const GridPos &gp_lat)
Returns the geomtrical altitudes of p_grid for one latitude.
void interp_atmsurface_by_gp(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates a surface-type variable given the grid positions.
void interp_atmsurface_by_itw(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates a surface-type variable with pre-calculated weights by interp_atmsurface_gp2itw.
Header file for special_interp.cc.
Structure to store a grid position.
Definition: interpolation.h:73
Index idx
Definition: interpolation.h:74