Home > atmlab > covmat > covmat3d.m

covmat3d

PURPOSE ^

COVMAT3D Creation of 1D-3D covariance matrices

SYNOPSIS ^

function S = covmat3d(dim,D,grid1,grid2,grid3,dtype)

DESCRIPTION ^

 COVMAT3D   Creation of 1D-3D covariance matrices

   Converts specification data to a covariance matrix. The covariance is
   described by the structure *D*. Correlation up to 3 dimensions are
   handled. The dimensions are here denoted as 1, 2 and 3, in order to be
   general. The physical meaning of these dimensions differ. For e.g. gas
   species they correspond to pressure/altitude, latitude and longitude.  For
   atmospheric data *dtype* shall be set to 'atm', which result in that
   log10(p) is used as vertical coordinate (assumed to be dim 1). All pressure
   grids are converted (including e.g. CL1_GRID1). A consequence is that
   correlation lengths shall be given in pressure decades. For the Earth's
   atmosphere one pressure decade is roughly 16 km.

   The data are assumed to be stored with dimension 1 as innermost "loop",
   dimension 2 next etc. For 3D case, the data order is as follows:
 
       [(x1,y1,z1) (x2,y1,z2) ...(xn,y1,z1) (x1,y2,z) ... (x1,y1,z2) ...]'

   The main field of D is FORMAT, specifying the format used (only one
   option implemented so far):

      D.FORMAT = 'param' 
      ---
      The covariance matrix is here described following the format of 
      *covmat3d_from_cfun*. The function *covmat3d_from_cfun* is also
      used for handling most general 2D and 3D cases. For 1D, 
      *covmat1d_from_cfun* is used. The D-format of *covmat3d_from_cfun*
      is used strictly for these two cases.
      It is further possible to add the field SEPERABLE to D. If this
      field is true, the function *covmat_partstat_corr* is used.
      Correlation lengths (D.CL1 etc.) are here throughout described as 
      for 1D. Otherwise as above, such as scalar SI and CL1 are allowed.
      A simple example:
         D.FORMAT    = 'param'; 
         D.SEPERABLE = 1;
         D.CCO       = 0.01; 
         D.CFUN1     = 'exp';
         D.CL1       = 0.2;
         D.CFUN2     = 'lin';
         D.CL2       = 1;          
   
 FORMAT   S = covmat3d(dim,D,grid1,grid2,grid3,dtype)
              or
          S = covmat3d(dim,D,G,dtype)

 OUT   S       Covariance matrix.
 IN    dim     Dimensionality.
       D       Specification data.
       grid1   Retrieval grid in dimension 1.
       G       In this argument pattern, this argument shall be an
               ArrayOfVectors, where first vector is used as grid1 etc.
               Can hold more grids than needed considering *dim*.
 OPT   grid2   Retrieval grid in dimension 2. Not needed if *dim* < 2.
       grid3   Retrieval grid in dimension 3. Not needed if *dim* < 3.
       dtype   Type of data. Only specified options are [] and 'atm'. See
               above for details. Default is [].

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

covmat3d.m

SOURCE CODE ^

0001 % COVMAT3D   Creation of 1D-3D covariance matrices
0002 %
0003 %   Converts specification data to a covariance matrix. The covariance is
0004 %   described by the structure *D*. Correlation up to 3 dimensions are
0005 %   handled. The dimensions are here denoted as 1, 2 and 3, in order to be
0006 %   general. The physical meaning of these dimensions differ. For e.g. gas
0007 %   species they correspond to pressure/altitude, latitude and longitude.  For
0008 %   atmospheric data *dtype* shall be set to 'atm', which result in that
0009 %   log10(p) is used as vertical coordinate (assumed to be dim 1). All pressure
0010 %   grids are converted (including e.g. CL1_GRID1). A consequence is that
0011 %   correlation lengths shall be given in pressure decades. For the Earth's
0012 %   atmosphere one pressure decade is roughly 16 km.
0013 %
0014 %   The data are assumed to be stored with dimension 1 as innermost "loop",
0015 %   dimension 2 next etc. For 3D case, the data order is as follows:
0016 %
0017 %       [(x1,y1,z1) (x2,y1,z2) ...(xn,y1,z1) (x1,y2,z) ... (x1,y1,z2) ...]'
0018 %
0019 %   The main field of D is FORMAT, specifying the format used (only one
0020 %   option implemented so far):
0021 %
0022 %      D.FORMAT = 'param'
0023 %      ---
0024 %      The covariance matrix is here described following the format of
0025 %      *covmat3d_from_cfun*. The function *covmat3d_from_cfun* is also
0026 %      used for handling most general 2D and 3D cases. For 1D,
0027 %      *covmat1d_from_cfun* is used. The D-format of *covmat3d_from_cfun*
0028 %      is used strictly for these two cases.
0029 %      It is further possible to add the field SEPERABLE to D. If this
0030 %      field is true, the function *covmat_partstat_corr* is used.
0031 %      Correlation lengths (D.CL1 etc.) are here throughout described as
0032 %      for 1D. Otherwise as above, such as scalar SI and CL1 are allowed.
0033 %      A simple example:
0034 %         D.FORMAT    = 'param';
0035 %         D.SEPERABLE = 1;
0036 %         D.CCO       = 0.01;
0037 %         D.CFUN1     = 'exp';
0038 %         D.CL1       = 0.2;
0039 %         D.CFUN2     = 'lin';
0040 %         D.CL2       = 1;
0041 %
0042 % FORMAT   S = covmat3d(dim,D,grid1,grid2,grid3,dtype)
0043 %              or
0044 %          S = covmat3d(dim,D,G,dtype)
0045 %
0046 % OUT   S       Covariance matrix.
0047 % IN    dim     Dimensionality.
0048 %       D       Specification data.
0049 %       grid1   Retrieval grid in dimension 1.
0050 %       G       In this argument pattern, this argument shall be an
0051 %               ArrayOfVectors, where first vector is used as grid1 etc.
0052 %               Can hold more grids than needed considering *dim*.
0053 % OPT   grid2   Retrieval grid in dimension 2. Not needed if *dim* < 2.
0054 %       grid3   Retrieval grid in dimension 3. Not needed if *dim* < 3.
0055 %       dtype   Type of data. Only specified options are [] and 'atm'. See
0056 %               above for details. Default is [].
0057 
0058 % 2006-08-21   Created by Patrick Eriksson.
0059 
0060 function S = covmat3d(dim,D,grid1,grid2,grid3,dtype)
0061                                                                           %&%
0062                                                                           %&%
0063 %--- Simple checks of input                                               %&%
0064 %                                                                         %&%
0065 rqre_nargin( 3, nargin );                                                 %&%
0066 rqre_alltypes( dim, {@istensor0,@iswhole} );                              %&%
0067 rqre_in_range( dim, 1, 3 );                                               %&%
0068 rqre_datatype( D, @isstruct );                                            %&%
0069 %
0070 if iscell( grid1 )    %--- Argument pattern 2 --------------------
0071   if nargin > 3
0072     dtype = grid2;
0073     if nargin > 4                                                         %&%
0074       error('To many input arguments for argument pattern 2.');           %&%
0075     end                                                                   %&%
0076   else
0077     dtype = [];
0078   end
0079   G = grid1;
0080 elseif isvector( grid1 )   %--- Argument pattern 1 ---------------
0081   rqre_nargin( 2+dim, nargin );                                           %&%
0082   rqre_datatype( grid1, @istensor1 );                                     %&%
0083   G{1} = grid1;
0084   if dim >= 2
0085     rqre_datatype( grid2, @istensor1 );                                   %&%
0086     G{2} = grid2;
0087     if dim == 3
0088       rqre_datatype( grid3, @istensor1 );                                 %&%
0089       G{3} = grid3;
0090     end
0091   end
0092   if nargin < 6
0093     dtype = [];
0094   end
0095 else                                                                      %&%
0096   error('Unknown type of input argument 3.');                             %&%
0097 end
0098 
0099 
0100 %--- Unit conversions
0101 %
0102 if isempty(dtype)
0103   %
0104 elseif strcmp(dtype,'atm')
0105   G{1} = log10(G{1});
0106   if isfield(D,'SI_GRID1'),  D.SI_GRID1  = log10(D.SI_GRID1);   end
0107   if isfield(D,'CL1_GRID1'), D.CL1_GRID1 = log10(D.CL1_GRID1);  end
0108   if isfield(D,'CL2_GRID1'), D.CL2_GRID1 = log10(D.CL2_GRID1);  end
0109   if isfield(D,'CL3_GRID1'), D.CL3_GRID1 = log10(D.CL3_GRID1);  end
0110 else
0111   error(sprintf('Unknown option for *dtype* (%s).',dtype));
0112 end
0113 
0114 
0115 
0116 %--- Parameterised version ---------------------------------------------------
0117 %
0118 % Some extra flexibility is offered for 1D and SEPERABLE option. Both row and
0119 % column vectors are allowed.
0120 %
0121 if strcmp( D.FORMAT, 'param' )
0122 
0123   %- Check consistency inside D
0124   %
0125   for id = 1 : length(G)
0126     if ~isscalar(D.SI) & size(D.SI,id)~=length(D.(sprintf('SI_GRID%d',id)))
0127       error( sprintf( 'Inconsistency between D.SI and D.SI_GRID%d.', id ) );
0128     end
0129     if ~isfield( D, 'SEPERABLE' )  |  ~D.SEPERABLE
0130     if ~isscalar(D.CL1) & size(D.CL1,id)~=length(D.(sprintf('CL1_GRID%d',id)))
0131       error( sprintf( 'Inconsistency between D.CL1 and D.CL1_GRID%d.', id ) );
0132     end
0133     if length(G)>=2  &  ...
0134        ~isscalar(D.CL2) & size(D.CL2,id)~=length(D.(sprintf('CL2_GRID%d',id)))
0135       error( sprintf( 'Inconsistency between D.CL2 and D.CL2_GRID%d.', id ) );
0136     end          
0137     if length(G)>=3  &  ...
0138        ~isscalar(D.CL3) & size(D.CL3,id)~=length(D.(sprintf('CL3_GRID%d',id)))
0139       error( sprintf( 'Inconsistency between D.CL3 and D.CL3_GRID%d.', id ) );
0140     end          
0141     end
0142   end
0143   
0144 
0145   if dim == 1               % 1D is handled separately.
0146     %
0147     if isscalar( D.SI )  
0148       Std = D.SI;
0149     elseif isvector( D.SI )
0150       Std = [ vec2col(D.SI_GRID1) vec2col(D.SI) ];
0151     else
0152       error( 'Not allowed format for D.SI.' );
0153     end
0154     %
0155     if isscalar( D.CL1 )
0156       Cl = D.CL1;
0157     elseif isvector( D.CL1 )
0158       Cl = [ vec2col(D.CL1_GRID1) vec2col(D.CL1) ];
0159     else
0160       error( 'Not allowed format of D.CL1 for 1D.' );
0161     end
0162     %
0163     S = covmat1d_from_cfun( G{1}, Std, D.CFUN1, Cl, D.CCO );
0164     %
0165   else    % 2D and 3D
0166     %
0167     % Most general option:
0168     if ~( isfield( D, 'SEPERABLE' )  &  D.SEPERABLE )
0169       S = covmat3d_from_cfun( dim, D, G{1:dim} );
0170       %
0171     else % Seperable option:
0172       %
0173       if isscalar( D.CL1 )    % Correlation for dim 1
0174         Cl = D.CL1;
0175       elseif isvector( D.CL1 )
0176         Cl = [ vec2col(D.CL1_GRID1) vec2col(D.CL1) ];
0177       else
0178         error( 'Not allowed format of D.CL1 for 2D/3D and D.SEPERABLE.' );
0179       end
0180       C{1} = covmat1d_from_cfun( G{1}, [], D.CFUN1, Cl, D.CCO );
0181       %
0182       if isscalar( D.CL2 )    % Correlation for dim 2
0183         Cl = D.CL2;
0184       elseif isvector( D.CL2 )
0185         Cl = [ vec2col(D.CL2_GRID1) vec2col(D.CL2) ];
0186       else
0187         error( 'Not allowed format of D.CL2 for 2D/3D and D.SEPERABLE.' );
0188       end
0189       C{2} = covmat1d_from_cfun( G{2}, [], D.CFUN2, Cl, D.CCO );
0190       %
0191       if dim == 3             % Correlation for dim 3
0192         if isscalar( D.CL3 )
0193           Cl = D.CL3;
0194         elseif isvector( D.CL3 )
0195           Cl = [ vec2col(D.CL3_GRID1) vec2col(D.CL3) ];
0196         else
0197           error( 'Not allowed format of D.CL3 for 3D and D.SEPERABLE.' );
0198         end
0199         C{3} = covmat1d_from_cfun( G{3}, [], D.CFUN3, Cl, D.CCO );
0200       end
0201       %
0202       if isscalar( D.SI )     % Standard devaition
0203         si = D.SI;
0204       else
0205         x = repmat( vec2col(grid1), length(grid2), 1 );          
0206         y = repmat( vec2row(grid2), length(grid1), 1 );
0207         y = y(:);
0208         if dim == 2
0209           si = interpd( [x,y], D.SI, D.SI_GRID1, D.SI_GRID2, 'linear' );
0210         else
0211           if dim == 3
0212             x = repmat( x, length(grid3), 1 );
0213             y = repmat( y, length(grid3), 1 );
0214             z = repmat( vec2row(grid3), length(grid1)*length(grid2), 1 );
0215             z = z(:);
0216             si = interpd( [x,y,z], D.SI, D.SI_GRID1, D.SI_GRID2, ...
0217                                                   D.SI_GRID3, 'linear' );
0218           end
0219         end
0220       end
0221       %
0222       S  = covmat_partstat_corr( si, C{1:dim} );
0223     end
0224   end
0225   
0226 
0227 %--- Unknown version
0228 else
0229   error( sprintf( 'Unknown choice for D.FORMAT (%s).', D.FORMAT ) );
0230 end

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