%filter method from Swenson & Wahr(2006) "Post-processing removal of correlated errors in GRACE data, Geophys. Res. Lett., 33, L08402, doi:10.1029/2005GL025285.".

%Gaussian filter
r=500000000;      %%%%%%%Radius
R=6378140000;
Degree=60;

W=zeros(1,Degree+1); 
b=log(2)/(1-cos(r/R));
W(1)=1;
W(2)=((1+exp(-2*b))/(1-exp(-2*b))-1/b);
for I=2:Degree
   W(I+1)=-(2*I-1)*W(I)/b+W(I-1);
end


function mm=filter(CC)
%format long 
double mm;
mm=CC;
[b,c]=size(CC);

for deg=8:b
    for l=9:50
L=zeros(3,3);      
  for i=0:2
       for j=0:2
            for n=l-8:2:l+8
               L(i+1,j+1)=L(i+1,j+1)+(n^i)*(n^j);
           end
       end 
  end 
M=zeros(3,1);
  m=0;
  for i=0:2
       for n=l-8:2:l+8
         M(m+1)=(n^i)*mm(n,deg)+M(m+1);
       end
       m=m+1;
  end
A=inv(L);
mm(l,deg)=[1 l l^2]*A*M;
    end 
end


end