Home > atmlab > gridcreation > gridselect2D.m

gridselect2D

PURPOSE ^

GRIDSELECT2D Simple selection of 2D grids

SYNOPSIS ^

function [xc,yc,Ac] = gridselect2D( xf, yf, Af, tol, repr, abs_or_rel )

DESCRIPTION ^

 GRIDSELECT2D   Simple selection of 2D grids

    The function needs a set of realisations of the quantity for which the
    grids shall be selected. These realisations are given along dim. 3 of the
    3D variable *Af*, are provided on a sufficient fine grid that 
    representation errors can be neglected and shall span the possible cases 
    that shall be represented. 

    The dimension of Af shall be [ xf, yf, number of realisation ]

    The grid is selected by first including the end points of *xf* and *yf*. 
    The coarse representation is transformed to the fine one, and a new grid 
    point is included where the maximum deviation is found (where all
    realisations are considered in parallel). The point is included in
    either *xc* or *yc*, depending on which that most effectively reduces
    the deviation. This procedure is repeated until the coarse and 
    representations deviates less then *tol*.

    See also GRIDSELECT1D and GRIDSELECT3D.

 FORMAT   [xc,yc,Ac] = gridselect2D( xf, yf, Af, tol, repr, abs_or_rel )

 OUT   xc          Obtained coarse grid for x dimension.
       yc          Obtained coarse grid for y dimension.
       Ac          Af at grid points xc and yc.
 IN    xf          Reference fine grid for x dimension.
       yf          Reference fine grid for y dimension.
       Af          Reference data (a matrix) on fine grid.
       tol         Acceptable tolarance.
       repr        String describing representation used. Possible choices
                   are equal to allowed 'methods' for *interp1*.
 OPT   abs_or_rel  String telling if tolarance is in absolute ('abs') or
                   relative ('rel') units. Default is 'abs'.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

gridselect2D.m

SOURCE CODE ^

0001 % GRIDSELECT2D   Simple selection of 2D grids
0002 %
0003 %    The function needs a set of realisations of the quantity for which the
0004 %    grids shall be selected. These realisations are given along dim. 3 of the
0005 %    3D variable *Af*, are provided on a sufficient fine grid that
0006 %    representation errors can be neglected and shall span the possible cases
0007 %    that shall be represented.
0008 %
0009 %    The dimension of Af shall be [ xf, yf, number of realisation ]
0010 %
0011 %    The grid is selected by first including the end points of *xf* and *yf*.
0012 %    The coarse representation is transformed to the fine one, and a new grid
0013 %    point is included where the maximum deviation is found (where all
0014 %    realisations are considered in parallel). The point is included in
0015 %    either *xc* or *yc*, depending on which that most effectively reduces
0016 %    the deviation. This procedure is repeated until the coarse and
0017 %    representations deviates less then *tol*.
0018 %
0019 %    See also GRIDSELECT1D and GRIDSELECT3D.
0020 %
0021 % FORMAT   [xc,yc,Ac] = gridselect2D( xf, yf, Af, tol, repr, abs_or_rel )
0022 %
0023 % OUT   xc          Obtained coarse grid for x dimension.
0024 %       yc          Obtained coarse grid for y dimension.
0025 %       Ac          Af at grid points xc and yc.
0026 % IN    xf          Reference fine grid for x dimension.
0027 %       yf          Reference fine grid for y dimension.
0028 %       Af          Reference data (a matrix) on fine grid.
0029 %       tol         Acceptable tolarance.
0030 %       repr        String describing representation used. Possible choices
0031 %                   are equal to allowed 'methods' for *interp1*.
0032 % OPT   abs_or_rel  String telling if tolarance is in absolute ('abs') or
0033 %                   relative ('rel') units. Default is 'abs'.
0034 
0035 % 2004-09-23   Created by Patrick Eriksson.
0036 
0037 
0038 function [xc,yc,Ac] = gridselect2D( xf, yf, Af, tol, repr, abs_or_rel )
0039 
0040 if nargin < 6
0041   abs_or_rel = 'abs';
0042 end
0043 
0044 nx = length( xf );
0045 ny = length( yf );
0046 
0047 
0048 %=== Check input
0049 %
0050 if ~isvector( xf )
0051   error( 'Input argument *xf* must be a vector.' );
0052 end
0053 %
0054 if ~isvector( yf )
0055   error( 'Input argument *yf* must be a vector.' );
0056 end
0057 %
0058 if ~( ndims( Af ) == 2  |  ndims( Af ) == 3 )  |  ~isnumeric( Af )
0059   error( 'Input argument *Af* must be a 3D matrix.' );
0060 end
0061 %
0062 if size( Af, 1 ) ~= nx  &&  size( Af, 2 ) ~= ny
0063   error( 'Sizes of *xf*/*yf* and *Af* do not match.' );
0064 end
0065 %
0066 if strcmp( abs_or_rel, 'rel' )  &  any( Af == 0 )
0067    error( 'With *abs_or_rel* = ''rel'', *Af can not contain zeros.' );
0068 end
0069 
0070 xind = [1 nx];
0071 yind = [1 ny];
0072 
0073 xf = vec2row( xf );
0074 yf = vec2row( yf );
0075 
0076 Ac   = zeros( size( Af ) );
0077 
0078 Ac = do_interp( xf(xind), yf(yind), Af(xind,yind,:), xf, yf, repr );
0079 [e,xip,yip] = get_max_dev( Af, Ac, abs_or_rel );
0080 
0081 while e > tol
0082 
0083   %= Test to put in xip
0084   if ~any( xip == xind )
0085     xindtmp = sort( [ xind xip ] );
0086     Ac = do_interp( xf(xindtmp), yf(yind), Af(xindtmp,yind,:), xf, yf, repr );
0087     [e1,xip1,yip1] = get_max_dev( Af, Ac, abs_or_rel );
0088   else
0089     e1 = Inf;
0090   end
0091 
0092   %= Test to put in yip
0093   if ~any( yip == yind )
0094     yindtmp = sort( [ yind yip ] );
0095     Ac = do_interp( xf(xind), yf(yindtmp), Af(xind,yindtmp,:), xf, yf, repr );
0096     [e2,xip2,yip2] = get_max_dev( Af, Ac, abs_or_rel );
0097   else
0098     e2 = Inf;
0099   end
0100   
0101   if e1 <= e2
0102     xind = xindtmp;
0103     xip  = xip1;
0104     yip  = yip1;
0105     e    = e1;
0106   else
0107     yind = yindtmp;
0108     xip  = xip2;
0109     yip  = yip2;
0110     e    = e2;
0111   end
0112 end
0113 
0114 xc = xf(xind);
0115 yc = yf(yind);
0116 
0117 if nargout > 2
0118   Ac = Af(xind,yind,:);
0119 end
0120 
0121 
0122 return
0123 
0124 
0125 function Ac = do_interp( xc, yc, A, xf, yf, repr )
0126   %
0127   for i = 1 : size(A,3)
0128     Ac(:,:,i) = interp2( yc, xc, A(:,:,i), yf, xf', repr );
0129   end
0130   %
0131 return
0132 
0133 
0134 function [e,xip,yip] = get_max_dev( Af, Ac, abs_or_rel );
0135   %
0136   if strcmp( abs_or_rel, 'abs' )
0137     [e,xip] = max( max( max( abs( Ac - Af ), [], 3 ), [], 2 ) );
0138     [e,yip] = max( max( max( abs( Ac - Af ), [], 3 ), [], 1 ) );
0139   elseif strcmp( abs_or_rel, 'rel' )
0140     [e,xip] = max( max( max( abs( (Ac - Af)./Af ), [], 3 ), [], 2 ) );
0141     [e,yip] = max( max( max( abs( (Ac - Af)./Af ), [], 3 ), [], 1 ) );
0142   else
0143     error( ...
0144       sprintf('Unknown error type (%s). Can only be ''abs'' or ''rel''.', ...
0145                                                                 abs_or_rel ) );
0146   end
0147   %
0148 return

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005