redistrib_jd

 This function redistributes the function f defined on the cells limited
 by the grid x onto the cells limited bythe grid y

 INPUTS : 
   - x : grid delimiting the original cells [nx,1];
   - f : function evaluated on the original cells [nx-1,1];
   - y : grid delimiting the output cells [ny,1];
   - testsum : (1) test that integral is conserved (0) do not test [1,1]; 

 OUTPUTS : 
   - g : function evaluated on the output cells [ny-1,1];

 by J. Decker, E. Nilsson and Y. Peysson

0001 function g = redistrib_jd(x,f,y,testsum)
0002 %
0003 % This function redistributes the function f defined on the cells limited
0004 % by the grid x onto the cells limited bythe grid y
0005 %
0006 % INPUTS :
0007 %   - x : grid delimiting the original cells [nx,1];
0008 %   - f : function evaluated on the original cells [nx-1,1];
0009 %   - y : grid delimiting the output cells [ny,1];
0010 %   - testsum : (1) test that integral is conserved (0) do not test [1,1];
0011 %
0012 % OUTPUTS :
0013 %   - g : function evaluated on the output cells [ny-1,1];
0014 %
0015 % by J. Decker, E. Nilsson and Y. Peysson
0016 %
0017 if nargin < 3,
0018     error('Not enough input arguments')
0019 end
0020 %
0021 x = x(:);
0022 f = f(:);
0023 %
0024 nx = length(x);
0025 %
0026 if nx < 2,
0027     error('The grid x should delimit at least one cell');
0028 end
0029 if length(f) ~= nx - 1
0030     error('The length of f must be nx - 1');
0031 end
0032 %
0033 if nargin < 3,% assumes y=x
0034     g = f;
0035     return
0036 end
0037 %
0038 y = y(:);
0039 ny = length(y);
0040 %
0041 dx = diff(x);
0042 dy = diff(y);
0043 %
0044 if any(dx <= 0),
0045     error('The grid x should be monotonically increasing.');
0046 end
0047 if any(dy <= 0),
0048     error('The grid y should be monotonically increasing.');
0049 end
0050 %
0051 if y(1) > x(1) || y(end) < x(end),
0052     disp('WARNING : some contributions of f will be lost as the y grid does not completely overlap the x grid')
0053 end
0054 %
0055 X = repmat(x.',[ny-1,1]);
0056 Y = repmat(y,[1,nx-1]);
0057 %
0058 maskYminp = Y(2:end,:) < X(:,2:end);
0059 maskYminm = Y(1:end-1,:) < X(:,1:end-1);
0060 %
0061 D = (Y(2:end,:).*maskYminp + X(:,2:end).*~maskYminp) - (Y(1:end-1,:).*~maskYminm + X(:,1:end-1).*maskYminm);% min(Xp,Yp) - max(Xm,Ym)
0062 %
0063 D(D < 0) = 0;
0064 %
0065 g = D*f./dy;
0066 %
0067 if nargin == 4 && testsum,
0068     %
0069     sx = sum(f.*dx);
0070     sy = sum(g.*dy);
0071     %
0072     disp(['The relative difference on the integral is : ',num2str(abs((sy - sx)/sx))])
0073     %
0074 end
0075 
0076 
0077 
0078