octave_packages/optiminterp-0.3.3/example_optiminterp.m

   1 %% Copyright (C) 2008 Alexander Barth
   2 %%
   3 %% This program is free software; you can redistribute it and/or modify
   4 %% it under the terms of the GNU General Public License as published by
   5 %% the Free Software Foundation; either version 2 of the License, or
   6 %% (at your option) any 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, see <http://www.gnu.org/licenses/>.
  15
  16 % Example program of the optimal interpolation toolbox
  17
  18
  19 % the grid onto which the observations are interpolated
  20
  21 [xi,yi] = ndgrid(linspace(0,1,100));
  22
  23 % background estimate or first guess
  24 xb = 10 + xi;
  25
  26 % number of observations to interpolate
  27
  28 on = 200;
  29
  30 % create randomly located observations within
  31 % the square [0 1] x [0 1]
  32
  33 x = rand(1,on);
  34 y = rand(1,on);
  35
  36 % the underlying function to interpolate
  37
  38 yo = 10 + x + sin(6*x) .* cos(6*y);
  39
  40 % the error variance of the observations divided by the error
  41 % variance of the background field
  42
  43 var = 0.1 * ones(on,1);
  44
  45 % the correlation length in x and y direction
  46
  47 lenx = 0.1;
  48 leny = 0.1;
  49
  50 % number of influential observations
  51
  52 m = 30;
  53
  54 % subtract the first guess from the observations
  55 % (DON'T FORGET THIS - THIS IS VERY IMPORTANT)
  56
  57 Hxb = interp2(xi(:,1),yi(1,:),xb',x,y);
  58 f = yo - Hxb;
  59
  60 % run the optimal interpolation
  61 % fi is the interpolated field and vari is its error variance
  62
  63 [fi,vari] = optiminterp2(x,y,f,var,lenx,leny,m,xi,yi);
  64
  65 % Add the first guess back
  66
  67 xa = fi + xb;