toy_models/ice_thick_distrib/runmean.m

function [fldOut]=runmean(fldIn,halfWindow,dim,varargin);
%object:    compute running mean window ('rmw') over a dimension
%input:     fldIn is the field to which the rmw will be applied
%           halfWindow is the half width of the rmw
%           dim is the dimension over which the rmw will be applied
%optional:  doCycle states whether the boundary condition is cyclic (1) or
%               not (0; default). If doCycle==0, the no. of averaged points
%               decreases from 1+2*halfWindow to halfWindow at the edges,
%               and we mask all of those edge points with NaNs.
%output:    fldOut is the resulting field
%
%notes:     - NaNs are discarded in the rmw, implying that an average is
%           computed if the rmw contains at least 1 valid point.
%           - setting halfWindow to 0 implies fldOut=fldIn.

%determine and check a couple things
fld_isa_gcmfaces=isa(fldIn,'gcmfaces');
if fld_isa_gcmfaces;
    if dim<3;
        error('for gcmfaces objects runmean excludes dim=1 or 2');
    end;
end;

if nargin==3; doCycle=0; else; doCycle=varargin{1}; end;

%switch to array format if needed
if fld_isa_gcmfaces; fldIn=convert2array(fldIn); end;

%switch dim to first dimension
sizeIn=size(fldIn); 

perm1to2=[1:length(sizeIn)]; 
perm1to2=[dim perm1to2(find(perm1to2~=dim))];
perm2to1=[[1:dim-1]+1 1 [dim+1:length(sizeIn)]]; 
sizeIn2=sizeIn(perm1to2);

fldIn=permute(fldIn,perm1to2); 

sizeCur=size(fldIn);
if ~doCycle;
    %add NaNs at both edges
    sizeCur(1)=sizeCur(1)+2*halfWindow;
    fldIn2=NaN*ones(sizeCur);
    fldIn2(halfWindow+1:end-halfWindow,:,:,:)=fldIn;
    fldIn=fldIn2; clear fldIn2;
end;

%create mask and remove NaNs:
fldMsk=~isnan(fldIn);
fldIn(isnan(fldIn))=0;
fldCnt=0*fldIn;

%apply the running mean
fldOut=zeros(sizeCur);
for tcur=-halfWindow:halfWindow
    %To have halfWindow*2 coeffs rather than halfWindow*2+1, centered to the current 
    %point, it is convenient to reduce the weight of the farthest points to 1/2.
    %This became necessary to get proper annual means, from monthly data, with halfWindow=6.
    if abs(tcur)==halfWindow; fac=1/2; else; fac=1; end;
    tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
    fldOut=fldOut+fac*tmp1;
    tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
    fldCnt=fldCnt+fac*tmp1;
end

% fldCnt(fldCnt<2*halfWindow)=NaN;
fldCnt(fldCnt<halfWindow/2)=NaN; fprintf('attention: modified runmean\n');
fldOut=fldOut./fldCnt;

if ~doCycle;
    fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:);
    %consistent with old version bug (one point offset)
    %     fldOut=fldOut(halfWindow:end-halfWindow,:,:,:);
end;

%switch dimensions order back to original order
fldOut=permute(fldOut,perm2to1);

%switch back to gcmfaces format if needed
if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;


1	gforget	1.1	function [fldOut]=runmean(fldIn,halfWindow,dim,varargin);
2			%object: compute running mean window ('rmw') over a dimension
3			%input: fldIn is the field to which the rmw will be applied
4			% halfWindow is the half width of the rmw
5			% dim is the dimension over which the rmw will be applied
6			%optional: doCycle states whether the boundary condition is cyclic (1) or
7			% not (0; default). If doCycle==0, the no. of averaged points
8			% decreases from 1+2*halfWindow to halfWindow at the edges,
9			% and we mask all of those edge points with NaNs.
10			%output: fldOut is the resulting field
11			%
12			%notes: - NaNs are discarded in the rmw, implying that an average is
13			% computed if the rmw contains at least 1 valid point.
14			% - setting halfWindow to 0 implies fldOut=fldIn.
15
16			%determine and check a couple things
17			fld_isa_gcmfaces=isa(fldIn,'gcmfaces');
18			if fld_isa_gcmfaces;
19			if dim<3;
20			error('for gcmfaces objects runmean excludes dim=1 or 2');
21			end;
22			end;
23
24			if nargin==3; doCycle=0; else; doCycle=varargin{1}; end;
25
26			%switch to array format if needed
27			if fld_isa_gcmfaces; fldIn=convert2array(fldIn); end;
28
29			%switch dim to first dimension
30			sizeIn=size(fldIn);
31
32			perm1to2=[1:length(sizeIn)];
33			perm1to2=[dim perm1to2(find(perm1to2~=dim))];
34			perm2to1=[[1:dim-1]+1 1 [dim+1:length(sizeIn)]];
35			sizeIn2=sizeIn(perm1to2);
36
37			fldIn=permute(fldIn,perm1to2);
38
39			sizeCur=size(fldIn);
40			if ~doCycle;
41			%add NaNs at both edges
42			sizeCur(1)=sizeCur(1)+2*halfWindow;
43			fldIn2=NaN*ones(sizeCur);
44			fldIn2(halfWindow+1:end-halfWindow,:,:,:)=fldIn;
45			fldIn=fldIn2; clear fldIn2;
46			end;
47
48			%create mask and remove NaNs:
49			fldMsk=~isnan(fldIn);
50			fldIn(isnan(fldIn))=0;
51			fldCnt=0*fldIn;
52
53			%apply the running mean
54			fldOut=zeros(sizeCur);
55			for tcur=-halfWindow:halfWindow
56			%To have halfWindow2 coeffs rather than halfWindow2+1, centered to the current
57			%point, it is convenient to reduce the weight of the farthest points to 1/2.
58			%This became necessary to get proper annual means, from monthly data, with halfWindow=6.
59			if abs(tcur)==halfWindow; fac=1/2; else; fac=1; end;
60			tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
61			fldOut=fldOut+fac*tmp1;
62			tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
63			fldCnt=fldCnt+fac*tmp1;
64			end
65
66			% fldCnt(fldCnt<2*halfWindow)=NaN;
67			fldCnt(fldCnt<halfWindow/2)=NaN; fprintf('attention: modified runmean\n');
68			fldOut=fldOut./fldCnt;
69
70			if ~doCycle;
71			fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:);
72			%consistent with old version bug (one point offset)
73			% fldOut=fldOut(halfWindow:end-halfWindow,:,:,:);
74			end;
75
76			%switch dimensions order back to original order
77			fldOut=permute(fldOut,perm2to1);
78
79			%switch back to gcmfaces format if needed
80			if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;
81
82
83