00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00061 #include <cmath>
00062 #include <iostream>
00063 #include <stdexcept>
00064 #include "check_input.h"
00065 #include "interpolation.h"
00066 #include "math_funcs.h"
00067 #include "messages.h"
00068 #include "special_interp.h"
00069
00070
00071
00072
00073
00074
00075
00076
00078
00108 void fix_gridpos_at_boundary(
00109 ArrayOfGridPos& gp,
00110 const Index grid_size
00111 )
00112 {
00113
00114 const Numeric epsilon = 1e-6;
00115
00116 for( Index i = 0; i< gp.nelem(); i++)
00117 {
00118 if( gp[i].idx == -1 )
00119 {
00120 if (abs(gp[i].fd[0]-1)>epsilon)
00121 {
00122 out1 << " --- WARNING ---, fix_gridpos_at_boundary encountered a strange "
00123 << "Gridpos: idx = " << gp[i].idx << ", fd[0] = " << gp[i].fd[0]
00124 << ", fd[1] = " << gp[i].fd[1];
00125 }
00126 gp[i].idx += 1;
00127 gp[i].fd[0] = 0.;
00128 gp[i].fd[1] = 1.;
00129 }
00130 else if (gp[i].idx == grid_size-1)
00131 {
00132 if (abs(gp[i].fd[0])>epsilon)
00133 {
00134 out1 << " --- WARNING ---, fix_gridpos_at_boundary encountered a strange "
00135 << "Gridpos: idx = " << gp[i].idx << ", fd[0] = " << gp[i].fd[0]
00136 << ", fd[1] = " << gp[i].fd[1];
00137 }
00138 gp[i].idx -= 1;
00139 gp[i].fd[0] = 1.;
00140 gp[i].fd[1] = 0.;
00141 }
00142 if ((gp[i].idx < 0)||(gp[i].idx > grid_size-1))
00143 {
00144 ostringstream os;
00145 os << "Invalid GridPos: idx = " << gp[i].idx << ", fd[0] = " << gp[i].fd[0] << ", fd[1] = " << gp[i].fd[1];
00146 throw runtime_error(os.str());
00147 }
00148 assert(gp[i].idx > -1);
00149 }
00150 }
00151
00152
00154
00178 void interp_atmfield_gp2itw(
00179 Matrix& itw,
00180 const Index& atmosphere_dim,
00181 ConstVectorView p_grid,
00182 ConstVectorView lat_grid,
00183 ConstVectorView lon_grid,
00184 const ArrayOfGridPos& gp_p,
00185 const ArrayOfGridPos& gp_lat,
00186 const ArrayOfGridPos& gp_lon )
00187 {
00188 chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
00189
00190 const Index n = gp_p.nelem();
00191
00192 if( atmosphere_dim == 1 )
00193 {
00194 itw.resize(n,2);
00195 interpweights( itw, gp_p );
00196 }
00197
00198 else if( atmosphere_dim == 2 )
00199 {
00200 assert( gp_lat.nelem() == n );
00201 itw.resize(n,4);
00202 interpweights( itw, gp_p, gp_lat );
00203 }
00204
00205 else if( atmosphere_dim == 3 )
00206 {
00207 assert( gp_lat.nelem() == n );
00208 assert( gp_lon.nelem() == n );
00209 itw.resize(n,8);
00210 interpweights( itw, gp_p, gp_lat, gp_lon );
00211 }
00212 }
00213
00214
00215
00217
00246 void interp_atmfield_by_itw(
00247 VectorView x,
00248 const Index& atmosphere_dim,
00249 ConstVectorView p_grid,
00250 ConstVectorView lat_grid,
00251 ConstVectorView lon_grid,
00252 ConstTensor3View x_field,
00253 const ArrayOfGridPos& gp_p,
00254 const ArrayOfGridPos& gp_lat,
00255 const ArrayOfGridPos& gp_lon,
00256 ConstMatrixView itw )
00257 {
00258
00259 chk_atm_field( "", x_field, atmosphere_dim, p_grid, lat_grid,
00260 lon_grid );
00261 assert( x.nelem() == gp_p.nelem() );
00262
00263 if( atmosphere_dim == 1 )
00264 {
00265 assert( itw.ncols() == 2 );
00266 interp( x, itw, x_field(Range(joker),0,0), gp_p );
00267 }
00268
00269 else if( atmosphere_dim == 2 )
00270 {
00271 assert( itw.ncols() == 4 );
00272 interp( x, itw, x_field(Range(joker),Range(joker),0), gp_p, gp_lat );
00273 }
00274
00275 else if( atmosphere_dim == 3 )
00276 {
00277 assert( itw.ncols() == 8 );
00278 interp( x, itw, x_field, gp_p, gp_lat, gp_lon );
00279 }
00280 }
00281
00282
00283
00285
00311 void interp_atmfield_by_gp(
00312 VectorView x,
00313 const Index& atmosphere_dim,
00314 ConstVectorView p_grid,
00315 ConstVectorView lat_grid,
00316 ConstVectorView lon_grid,
00317 ConstTensor3View x_field,
00318 const ArrayOfGridPos& gp_p,
00319 const ArrayOfGridPos& gp_lat,
00320 const ArrayOfGridPos& gp_lon )
00321 {
00322 Matrix itw;
00323
00324 interp_atmfield_gp2itw( itw, atmosphere_dim, p_grid, lat_grid,
00325 lon_grid, gp_p, gp_lat, gp_lon );
00326
00327 interp_atmfield_by_itw( x, atmosphere_dim, p_grid, lat_grid, lon_grid,
00328 x_field, gp_p, gp_lat, gp_lon, itw );
00329 }
00330
00331
00332
00334
00340 Numeric interp_atmfield_by_gp(
00341 const Index& atmosphere_dim,
00342 ConstVectorView p_grid,
00343 ConstVectorView lat_grid,
00344 ConstVectorView lon_grid,
00345 ConstTensor3View x_field,
00346 const GridPos& gp_p,
00347 const GridPos& gp_lat,
00348 const GridPos& gp_lon )
00349 {
00350 ArrayOfGridPos agp_p(1), agp_lat(0), agp_lon(0);
00351
00352 gridpos_copy( agp_p[0], gp_p );
00353
00354 if( atmosphere_dim > 1 )
00355 {
00356 agp_lat.resize(1);
00357 gridpos_copy( agp_lat[0], gp_lat );
00358 }
00359
00360 if( atmosphere_dim > 2 )
00361 {
00362 agp_lon.resize(1);
00363 gridpos_copy( agp_lon[0], gp_lon );
00364 }
00365
00366 Vector x(1);
00367
00368 interp_atmfield_by_gp( x, atmosphere_dim, p_grid, lat_grid, lon_grid,
00369 x_field, agp_p, agp_lat, agp_lon );
00370
00371 return x[0];
00372 }
00373
00374
00375
00377
00397 void interp_atmsurface_gp2itw(
00398 Matrix& itw,
00399 const Index& atmosphere_dim,
00400 const ArrayOfGridPos& gp_lat,
00401 const ArrayOfGridPos& gp_lon )
00402 {
00403 if( atmosphere_dim == 1 )
00404 {
00405 itw.resize(1,1);
00406 itw = 1;
00407 }
00408
00409 else if( atmosphere_dim == 2 )
00410 {
00411 const Index n = gp_lat.nelem();
00412 itw.resize(n,2);
00413 interpweights( itw, gp_lat );
00414 }
00415
00416 else if( atmosphere_dim == 3 )
00417 {
00418 const Index n = gp_lat.nelem();
00419 assert( n == gp_lon.nelem() );
00420 itw.resize(n,4);
00421 interpweights( itw, gp_lat, gp_lon );
00422 }
00423 }
00424
00425
00426
00428
00455 void interp_atmsurface_by_itw(
00456 VectorView x,
00457 const Index& atmosphere_dim,
00458 ConstVectorView lat_grid,
00459 ConstVectorView lon_grid,
00460 ConstMatrixView x_surface,
00461 const ArrayOfGridPos& gp_lat,
00462 const ArrayOfGridPos& gp_lon,
00463 ConstMatrixView itw )
00464 {
00465
00466 chk_atm_surface( "", x_surface, atmosphere_dim, lat_grid,
00467 lon_grid );
00468
00469 if( atmosphere_dim == 1 )
00470 {
00471 assert( itw.ncols() == 1 );
00472 x = x_surface(0,0);
00473 }
00474
00475 else if( atmosphere_dim == 2 )
00476 {
00477 assert( x.nelem() == gp_lat.nelem() );
00478 assert( itw.ncols() == 2 );
00479 interp( x, itw, x_surface(Range(joker),0), gp_lat );
00480 }
00481
00482 else if( atmosphere_dim == 3 )
00483 {
00484 assert( x.nelem() == gp_lat.nelem() );
00485 assert( itw.ncols() == 4 );
00486 interp( x, itw, x_surface, gp_lat, gp_lon );
00487 }
00488 }
00489
00490
00491
00493
00517 void interp_atmsurface_by_gp(
00518 VectorView x,
00519 const Index& atmosphere_dim,
00520 ConstVectorView lat_grid,
00521 ConstVectorView lon_grid,
00522 ConstMatrixView x_surface,
00523 const ArrayOfGridPos& gp_lat,
00524 const ArrayOfGridPos& gp_lon )
00525 {
00526 Matrix itw;
00527
00528 interp_atmsurface_gp2itw( itw, atmosphere_dim, gp_lat, gp_lon );
00529
00530 interp_atmsurface_by_itw( x, atmosphere_dim, lat_grid, lon_grid,
00531 x_surface, gp_lat, gp_lon, itw );
00532 }
00533
00534
00535
00537
00543 Numeric interp_atmsurface_by_gp(
00544 const Index& atmosphere_dim,
00545 ConstVectorView lat_grid,
00546 ConstVectorView lon_grid,
00547 ConstMatrixView x_surface,
00548 const GridPos& gp_lat,
00549 const GridPos& gp_lon )
00550 {
00551 ArrayOfGridPos agp_lat(0), agp_lon(0);
00552
00553 if( atmosphere_dim > 1 )
00554 {
00555 agp_lat.resize(1);
00556 gridpos_copy( agp_lat[0], gp_lat );
00557 }
00558
00559 if( atmosphere_dim > 2 )
00560 {
00561 agp_lon.resize(1);
00562 gridpos_copy( agp_lon[0], gp_lon );
00563 }
00564
00565 Vector x(1);
00566
00567 interp_atmsurface_by_gp( x, atmosphere_dim, lat_grid, lon_grid, x_surface,
00568 agp_lat, agp_lon );
00569
00570 return x[0];
00571 }
00572
00573
00574
00576
00598 void interp_gfield3(
00599 Numeric& value,
00600 const GField3& gfield3,
00601 const Index& effective_dim,
00602 const Numeric& x,
00603 const Numeric& y,
00604 const Numeric& z,
00605 const String& dim0,
00606 const String& dim1,
00607 const String& dim2 )
00608 {
00609 chk_if_in_range( "effective_dim", effective_dim, 1, 3 );
00610
00611 ArrayOfGridPos gp0(1), gp1(1), gp2(1);
00612
00613
00614
00615 {
00616 const String gridname = gfield3.get_grid_name(0);
00617 const ConstVectorView grid = gfield3.get_numeric_grid(0);
00618
00619 if( dim0 != gridname )
00620 {
00621 ostringstream os;
00622 os << "Wrong quantity found for grid of first dimension:\n"
00623 << " expected quantity : " << dim0 << "\n"
00624 << " quantity in gfield : " << gridname << "\n";
00625 throw runtime_error( os.str() );
00626 }
00627
00628 chk_if_increasing( "first grid of gfield3", grid );
00629
00630 if( x < grid[0] || x > last(grid) )
00631 {
00632 ostringstream os;
00633 os << "Interpolation outside covered range for first dimension:\n"
00634 << " interpolation point : " << x << "\n"
00635 << " gfield grid range : [" << grid[0] << "," << last(grid)
00636 << "]\n";
00637 throw runtime_error( os.str() );
00638 }
00639
00640 gridpos( gp0, grid, Vector(1,x) );
00641 }
00642
00643
00644
00645 if( effective_dim >= 2 )
00646 {
00647 const String gridname = gfield3.get_grid_name(1);
00648 const ConstVectorView grid = gfield3.get_numeric_grid(1);
00649
00650 if( dim1 != gridname )
00651 {
00652 ostringstream os;
00653 os << "Wrong quantity found for grid of second dimension:\n"
00654 << " expected quantity : " << dim1 << "\n"
00655 << " quantity in gfield : " << gridname << "\n";
00656 throw runtime_error( os.str() );
00657 }
00658
00659 chk_if_increasing( "second grid of gfield3", grid );
00660
00661 if( y < grid[0] || y > last(grid) )
00662 {
00663 ostringstream os;
00664 os << "Interpolation outside covered range for second dimension:\n"
00665 << " interpolation point : " << y << "\n"
00666 << " gfield grid range : [" << grid[0] << "," << last(grid)
00667 << "]\n";
00668 throw runtime_error( os.str() );
00669 }
00670
00671 gridpos( gp1, grid, Vector(1,y) );
00672 }
00673
00674
00675
00676 if( effective_dim >= 3 )
00677 {
00678 const String gridname = gfield3.get_grid_name(2);
00679 const ConstVectorView grid = gfield3.get_numeric_grid(2);
00680
00681 if( dim2 != gridname )
00682 {
00683 ostringstream os;
00684 os << "Wrong quantity found for grid of third dimension:\n"
00685 << " expected quantity : " << dim2 << "\n"
00686 << " quantity in gfield : " << gridname << "\n";
00687 throw runtime_error( os.str() );
00688 }
00689
00690 chk_if_increasing( "third grid of gfield3", grid );
00691
00692 if( z < grid[0] || z > last(grid) )
00693 {
00694 ostringstream os;
00695 os << "Interpolation outside covered range for second dimension:\n"
00696 << " interpolation point : " << z << "\n"
00697 << " gfield grid range : [" << grid[0] << "," << last(grid)
00698 << "]\n";
00699 throw runtime_error( os.str() );
00700 }
00701
00702 gridpos( gp2, grid, Vector(1,z) );
00703 }
00704
00705
00706
00707 Vector result( 1 );
00708
00709 if( effective_dim == 1 )
00710 {
00711 if( gfield3.nrows() > 1 || gfield3.ncols() > 1 )
00712 {
00713 ostringstream os;
00714 os << "A 1D interpolation requested, but the provided gridded field "
00715 << "has an effective dimension of 2D or 3D.";
00716 throw runtime_error( os.str() );
00717 }
00718
00719 Matrix itw(1,2);
00720 interpweights( itw, gp0 );
00721 interp( result, itw, gfield3(joker,0,0), gp0 );
00722 }
00723
00724 else if( effective_dim == 2 )
00725 {
00726 if( gfield3.ncols() > 1 )
00727 {
00728 ostringstream os;
00729 os << "A 2D interpolation requested, but the provided gridded field "
00730 << "has an effective dimension of 3D.";
00731 throw runtime_error( os.str() );
00732 }
00733
00734 Matrix itw(1,4);
00735 interpweights( itw, gp0, gp1 );
00736 interp( result, itw, gfield3(joker,joker,0), gp0, gp1 );
00737 }
00738
00739 else if( effective_dim == 3 )
00740 {
00741 Matrix itw(1,8);
00742 interpweights( itw, gp0, gp1, gp2 );
00743 interp( result, itw, gfield3, gp0, gp1, gp2 );
00744 }
00745
00746 value = result[0];
00747 }
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00759
00780 void itw2p(
00781 VectorView p_values,
00782 ConstVectorView p_grid,
00783 const ArrayOfGridPos& gp,
00784 ConstMatrixView itw )
00785 {
00786 assert( itw.ncols() == 2 );
00787 assert( p_values.nelem() == itw.nrows() );
00788
00789
00790 Vector logpgrid( p_grid.nelem() );
00791
00792 transform( logpgrid, log, p_grid );
00793
00794 interp( p_values, itw, logpgrid, gp );
00795
00796 transform( p_values, exp, p_values );
00797 }
00798
00799
00800
00802
00827 void p2gridpos( ArrayOfGridPos& gp,
00828 ConstVectorView old_pgrid,
00829 ConstVectorView new_pgrid,
00830 const Numeric& extpolfac )
00831 {
00832
00833 Vector logold( old_pgrid.nelem() );
00834 Vector lognew( new_pgrid.nelem() );
00835
00836 transform( logold, log, old_pgrid );
00837 transform( lognew, log, new_pgrid );
00838
00839 gridpos( gp, logold, lognew, extpolfac );
00840 }
00841
00842
00843
00844
00845
00846
00847
00849
00866 void z_at_lat_2d(
00867 VectorView z,
00868 ConstVectorView p_grid,
00869
00870 #ifndef NDEBUG
00871 ConstVectorView lat_grid,
00872 #else
00873 ConstVectorView,
00874 #endif
00875 ConstMatrixView z_field,
00876 const GridPos& gp_lat )
00877 {
00878 const Index np = p_grid.nelem();
00879
00880 assert( z.nelem() == np );
00881 assert( z_field.nrows() == np );
00882 assert( z_field.ncols() == lat_grid.nelem() );
00883
00884 Matrix z_matrix(np,1);
00885 ArrayOfGridPos gp_z(np), agp_lat(1);
00886 Tensor3 itw(np,1,4);
00887
00888 gridpos_copy( agp_lat[0], gp_lat );
00889 gridpos( gp_z, p_grid, p_grid );
00890 interpweights( itw, gp_z, agp_lat );
00891 interp( z_matrix, itw, z_field, gp_z, agp_lat );
00892
00893 z = z_matrix(Range(joker),0);
00894 }
00895
00896
00897
00899
00918 void z_at_latlon(
00919 VectorView z,
00920 ConstVectorView p_grid,
00921
00922 #ifndef NDEBUG
00923 ConstVectorView lat_grid,
00924 ConstVectorView lon_grid,
00925 #else
00926 ConstVectorView,
00927 ConstVectorView,
00928 #endif
00929 ConstTensor3View z_field,
00930 const GridPos& gp_lat,
00931 const GridPos& gp_lon )
00932 {
00933 const Index np = p_grid.nelem();
00934
00935 assert( z.nelem() == np );
00936 assert( z_field.npages() == np );
00937 assert( z_field.nrows() == lat_grid.nelem() );
00938 assert( z_field.ncols() == lon_grid.nelem() );
00939
00940 Tensor3 z_tensor(np,1,1);
00941 ArrayOfGridPos agp_z(np), agp_lat(1), agp_lon(1);
00942 Tensor4 itw(np,1,1,8);
00943
00944 gridpos_copy( agp_lat[0], gp_lat );
00945 gridpos_copy( agp_lon[0], gp_lon );
00946 gridpos( agp_z, p_grid, p_grid );
00947 interpweights( itw, agp_z, agp_lat, agp_lon );
00948
00949 interp( z_tensor, itw, z_field, agp_z, agp_lat, agp_lon );
00950
00951 z = z_tensor(Range(joker),0,0);
00952 }
