function lmcosif=plmfilt(lmcosi,L,wlen)
% lmcosif=plmfilt(lmcosi,L,wlen)
%
% A very simple filtering routine for spherical harmonics
%
% INPUT:
%
% lmcosi    Standard [degree order cosine sine] coefficients, not
%           necessarily starting at zero degree
% L         Single degree with lowpass degree of halfway point
%           or two degrees for a bandpass filter, NaN for nothing 
% wlen      The odd length of the transition band [default: 5]
%
% OUTPUT:
%
% lmcosid   Filtered equivalent of the input array
%
% Last modified by fjsimons-at-alum.mit.edu, 06/22/2011

% Specify defaults
defval('wlen',5)
defval('L',180)

if isnan(L)
  lmcosif=lmcosi;
  return
end

% Make sure it's sorted descending; only lowpass or highpass just yet
L=sort(L,'descend');

if ~mod(wlen,2)
  error('Total filter length must be odd')
end

% What is the half length of the filter
whalf=(wlen-1)/2;

% Generate the window as filter coefficients
[w,wl,wr]=fhanning(2*wlen);

% Initialize
lmcosif=lmcosi;

% Where does it all start, that determines the offset
lofs=addmup(lmcosi(1)-1);

% Repeat them the right number of times for the center degree
for index=1:length(L)
  if L(index)+whalf>lmcosi(end,1) | L(index)-whalf<lmcosi(1,1)
    error('Filter too long or bandlimit too high/low for input')
  end
  degfilt=L(index)+[-whalf:whalf];
  % Cut to the right
  if index==1 
    wrr=gamini(wr,degfilt+1)';
  end
  % Cut to the left if there is a LOWER degree coming
  if index==2
    wrr=flipud(gamini(wr,degfilt+1)');
  end
  
  disp(sprintf('Filtering degrees %i - (%i) - %i',...
	       degfilt(1),L(index),degfilt(end)))
  % Perform the filtering 
  degind=[addmup(degfilt(1)-1)+1-lofs:...
	  addmup(degfilt(end))-lofs]';
  lmcosif(degind,3)=lmcosif(degind,3).*wrr;
  lmcosif(degind,4)=lmcosif(degind,4).*wrr;
end

% Truncate to the new maximal resolution
%lmcosif=lmcosif(1-lofs:addmup(L(1)+whalf),:);
lmcosif=lmcosif(1:addmup(L(1)+whalf)-lofs,:);

% Truncate from the new minimal resolution
if length(L)==2
  %lmcosif=lmcosif(addmup(L(2)-whalf-1)+1:end,:);
  lmcosif=lmcosif(addmup(L(2)-whalf-1)-lofs+1:end,:);
end