Home > atmlab > gridcreation > p_merge.m

p_merge

PURPOSE ^

+ P_MERGE: Merges two vertical profiles weighted by the uncertainties of the

SYNOPSIS ^

function [vprofile,si] = p_merge(p1,p2,g1,g2,ngrid,p1u,p2u,verbose,debug)

DESCRIPTION ^

+ P_MERGE: Merges two vertical profiles weighted by the uncertainties of the
  individual profiles

  Program interpolates the two profiles and uncertainties unto a common
  grid and merges them using the formula for weigthed mean:
  profile = (np1./(np1u.*np1u) + np2./(np2u.*np2u)) ./ ...
            (1./(np1u.*np1u) + 1./(np2u.*np2u));

  The uncertainties are also merged: 1/(u)^2 = 1/(u1)^2 + 1/(u2)^2

  OUT       new vertical profile and the error/uncertainties per level
  IN      p1 = profile 1, p2 = profile 2, g1 = old grid of profile 1, g2 =
  old grid of profile 2, ngrid = new common grid, p1u = uncertainties of
  profile 1, uncertainties of profile 2.

  Requirements: max(g2) must equal max(ngrid).
                min(g1) must equal min(ngrid).
                All grids must monotonic decreasing with max value in
                index 1.
                The grids to be merged must overlap!
                The uncertainties must be in absolute terms!

  Options: verbose: set to true or false if you want to see the data
                    inside the function
           debug:   set to true or false if you want to stop inside the
                    function using "keyboard"
+

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

p_merge.m

SOURCE CODE ^

0001 %+ P_MERGE: Merges two vertical profiles weighted by the uncertainties of the
0002 %  individual profiles
0003 %
0004 %  Program interpolates the two profiles and uncertainties unto a common
0005 %  grid and merges them using the formula for weigthed mean:
0006 %  profile = (np1./(np1u.*np1u) + np2./(np2u.*np2u)) ./ ...
0007 %            (1./(np1u.*np1u) + 1./(np2u.*np2u));
0008 %
0009 %  The uncertainties are also merged: 1/(u)^2 = 1/(u1)^2 + 1/(u2)^2
0010 %
0011 %  OUT       new vertical profile and the error/uncertainties per level
0012 %  IN      p1 = profile 1, p2 = profile 2, g1 = old grid of profile 1, g2 =
0013 %  old grid of profile 2, ngrid = new common grid, p1u = uncertainties of
0014 %  profile 1, uncertainties of profile 2.
0015 %
0016 %  Requirements: max(g2) must equal max(ngrid).
0017 %                min(g1) must equal min(ngrid).
0018 %                All grids must monotonic decreasing with max value in
0019 %                index 1.
0020 %                The grids to be merged must overlap!
0021 %                The uncertainties must be in absolute terms!
0022 %
0023 %  Options: verbose: set to true or false if you want to see the data
0024 %                    inside the function
0025 %           debug:   set to true or false if you want to stop inside the
0026 %                    function using "keyboard"
0027 %+
0028 
0029 % 2010-10-22 Created by Marston Johnston
0030 
0031 function [vprofile,si] = p_merge(p1,p2,g1,g2,ngrid,p1u,p2u,verbose,debug)
0032 
0033 if verbose || debug
0034     disp('Input data...');
0035     for p = 1:numel(p1), fprintf('g1(%d) = %f -> %f -> %f\n',p,g1(p),p1(p),p1u(p)); end
0036     for p = 1:numel(p2), fprintf('p2(%d) = %f %f -> %f\n',p,g2(p),p2(p),p2u(p)); end
0037 end
0038 vprofile = NaN(size(ngrid));
0039 si       = NaN(size(ngrid));
0040 
0041 % Find the degree of overlapp between the vectors. The grids must be
0042 % monotonic decreasing
0043 
0044 if all(diff(g1)<0) == true
0045     if verbose || debug, disp('Grid 1 is good'); end
0046 else
0047     error('Grid 1 must be monotonic decreasing');
0048 end
0049 
0050 if all(diff(g2)<0) == true
0051     if verbose || debug, disp('Grid 2 is good'); end
0052 else
0053     error('Grid 2 must be monotonic decreasing');
0054 end
0055 
0056 if all(diff(ngrid)<0) == true
0057     if verbose || debug, disp('ngrid is good'); end
0058 else
0059     error('Ngrid must be monotonic decreasing');
0060 end
0061 
0062 % Find where each layer is in relation to the other
0063 if max(round(g2)) == max(round(ngrid))
0064     layer1 = vec2col(g2);
0065     values1 = vec2col(p2);
0066     uncertainty1 = vec2col(p2u);
0067     layer2 = vec2col(g1);
0068     values2 = vec2col(g1);
0069     uncertainty2 = vec2col(p1u);
0070 elseif max(round(g1)) == max(round(ngrid))
0071     layer1 = vec2col(g1);
0072     values1 = vec2col(p1);
0073     uncertainty1 = vec2col(p1u);
0074     layer2 = vec2col(g2);
0075     values2 = vec2col(p2);
0076     uncertainty2 = vec2col(p2u);
0077 end
0078 % Find the overlapping areas
0079 if min(layer1) < max(layer2)
0080     if verbose, disp('There is an overlap!'); end
0081     i = min(layer1) <= ngrid & max(layer2) >= ngrid;
0082     iu1 = layer1 > max(layer2);
0083     iu2 = layer2 < min(layer1);
0084     % Find must be used here to get the indices
0085     ic1 = find(layer1 <= max(ngrid(i)));
0086     ic2 = find(layer2 >= min(ngrid(i)));
0087     % Add to the end points to ensure full data coverage. Making sure it
0088     % does not go past the end points!
0089     if ic1(1) > 1 && ic1(end) < numel(layer1)
0090         ic1 = [ic1(1)-1; ic1; ic1(end)+1];
0091     elseif ic1(1) == 1 && ic1(end) < numel(layer1)
0092         ic1 = [ic1; ic1(end)+1];
0093     elseif ic1(1) > 1 && ic1(end) == numel(layer1)
0094         ic1 = [ic1(1)-1; ic1];
0095     end
0096     if ic2(1) > 1 && ic2(end) < numel(layer2)
0097         ic2 = [ic2(1)-1; ic2; ic2(end)+1];
0098     elseif ic2(1) == 1 && ic2(end) < numel(layer2)
0099         ic2 = [ic2; ic2(end)+1];
0100     elseif ic2(1) > 1 && ic2(end) == numel(layer2)
0101         ic2 = [ic2(1)-1; ic2];
0102     end
0103 else
0104     error('Warning: There are no overlapping areas!');
0105 end
0106 % Copy the unique parts
0107 v1 = ismember(ngrid,layer1(iu1));
0108 v2 = ismember(ngrid,layer2(iu2));
0109 vprofile(v1) = values1(iu1);
0110 vprofile(v2) = values2(iu2);
0111 si(v1) = uncertainty1(iu1);
0112 si(v2) = uncertainty2(iu2);
0113 if verbose || debug
0114     disp('After copying the unique parts...');
0115     for p = 1:numel(values1), fprintf('layer1(%d) = %f %f\n',p,layer1(p),values1(p)); end
0116     for p = 1:numel(values2), fprintf('layer2(%d) = %f %f\n',p,layer2(p),values2(p)); end
0117     for p = 1:numel(vprofile), fprintf('vprofile(%d) = %f\n',p,vprofile(p)); end
0118     for p = 1:numel(si), fprintf('si(%d) = %f\n',p,si(p)); end
0119 end
0120 % Interpolate the common overlapp
0121 if sum(i) > 0
0122     values_common1      = interp1( layer1(ic1), values1(ic1), ngrid(i) );
0123     uncertainty_common1 = interp1( layer1(ic1), uncertainty1(ic1), ngrid(i) );
0124     values_common2      = interp1( layer2(ic2), values2(ic2), ngrid(i) );
0125     uncertainty_common2 = interp1( layer2(ic2), uncertainty2(ic2), ngrid(i) );
0126     % Calculates the weigths
0127     w1 = 1 ./ (uncertainty_common1.^2);
0128     w2 = 1 ./ (uncertainty_common2.^2);
0129     ws = w1 + w2;
0130     % Calculate the the new profile values and uncertainty values
0131     vprofile(i) = ( w1.*values_common1 + w2.*values_common2 ) ./ ws;
0132     si(i)      = sqrt( 1 ./ ws );
0133 end
0134 if verbose || debug
0135     disp('After copying the common parts...');
0136     for p = 1:numel(values1), fprintf('layer1(%d) = %f %f\n',p,layer1(p),values1(p)); end
0137     for p = 1:numel(values2), fprintf('layer2(%d) = %f %f\n',p,layer2(p),values2(p)); end
0138     for p = 1:numel(vprofile), fprintf('vprofile(%d) = %f\n',p,vprofile(p)); end
0139     for p = 1:numel(si), fprintf('si(%d) = %f\n',p,si(p)); end
0140 end
0141 
0142 if debug, disp('Stopping in debug mode. Type return to continue!'), keyboard; end
0143 
0144 % Some error checks!
0145 if any(isnan(vprofile))
0146     for p = 1:numel(vprofile), fprintf('vprofile(%d) = %f\n',p,vprofile(p)); end
0147     error('Resulting profile includes nans');
0148 end
0149 if any(isnan(si))
0150     for p = 1:numel(si), fprintf('si(%d) = %f\n',p,si(p)); end
0151     error('Resulting uncertainties includes nans');
0152 end

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