function [Bl,dels]=bpboxcap(TH,Lmax,meth,evens,sord,xver)
% [Bl,dels]=BPBOXCAP(TH,Lmax,meth,evens,sord,xver)
%
% Returns the spherical harmonic power spectrum of cylindrical 
% boxcar polar caps: 1 in the region and 0 elsewhere.
%
% INPUT:
%
% TH       Colatitudinal radius of the cap, in degrees, may be vector
%          Colatitudinal halfwidth of the belt, degrees, may be vector
% Lmax     Maximum degree of this non-bandlimited function
% meth     1 Analytically, using unnormalized Legendre polynomials [default]
%          2 Numerically, integrating normalized Legendre polynomials 
% evens    0 For all degrees [default]
%          1 Only return the even degrees
% sord     1 Single cap of diameter 2TH [default]
%          2 Double cap left by subtracting belt of width 2TH
%          3 Equatorial belt of width 2TH
% xver     1 Provides extra verification
%          0 Doesn't [default]
%
% OUTPUT:
%
% Bl       The power spectrum 
% dels     The spherical harmonic degrees
%
% EXAMPLE:
%
% bpboxcap('demo') % should return no error messages and pass all checks
%
% Last modified by fjsimons-at-alum.mit.edu, 04/25/2007

% For all the others, sum((2*dels+1).*Bl)=4*pi does NOT hold since they
% are not normalized that way.

if ~isstr(TH)
  defval('TH',100)
  defval('Lmax',100)
  defval('meth',1)
  defval('sord',1)
  defval('xver',0)
  
  t=cputime;
  
  % If double boxcar cap, always symmetric, always faster, reinsert at the end
  if sord==2 || sord==3
    waseven=evens;
    defval('evens',1);
  else
    defval('evens',0)
  end
  
  % Returns the spherical harmonic degrees
  dels=0:evens+1:Lmax;
  dels1=0:Lmax+1;
  
  switch meth 
   case 1
    if sord==1
      % Get the unnormalized spherical harmonics
      [pel,mu,norms]=plm(dels1,0,cos(TH*pi/180),xver);
    else
      [pel,mu,norms]=plm(dels1,0,sin(TH*pi/180),xver);
    end
    % Get the power spectrum by application of the analytical formula
    Bl=pi*([repmat(1,1,length(TH)) ; pel(1+evens:evens+1:end-2,:)]...
	   -[pel(2:evens+1:end,:)]).^2 ...
       ./repmat((2*dels'+1).^2,1,length(TH));
    if sord==2 || sord==3
      Bl=Bl.*[1+(-1).^repmat(dels',1,length(TH))].^2;
    end  
    if sord==3 % Need extra term for the zeroth order
	       % Here you cannot take out the prefactors
	       Bl(1,:)=pi*(2-2*[repmat(1,1,length(TH))-pel(2,:)]).^2;
    end
   case 2
    switch sord
     case 1
      intv=[cos(TH(:)*pi/180) repmat(1,length(TH),1)];
     case 2
      intv=[cos(pi/2-TH(:)*pi/180) repmat(1,length(TH),1)];
     case 3
      intv=[cos(pi/2+TH(:)*pi/180) cos(pi/2-TH(:)*pi/180)];
    end
    % Overdo the GL weights for all but the last one
    [wGL,xGL]=gausslegendrecof(Lmax,[],intv)
    % Initialize this matrix as a potentially 3D object
    xel=repmat(NaN,[length(dels) size(xGL)]);
    % Do Gauss-Legendre on a whole matrix at the same time
    % Writing (:,:,) explicitly helps when there is only one degree
    xel(:,:,:)=xlm(dels,0,acos(xGL),xver);
    % Initialize the results matrix
    Bl=repmat(NaN,length(dels),length(TH));
    for index=1:length(TH)
      % Get the power spectrum by numerical integration
      Bl(:,index)=...
	  [(1+[sord==2])*2*pi*wGL(:,index)'*xel(:,:,index)'].^2./(2*dels+1);
    end
  end
  
  if (sord==2 || sord==3) 
    if evens==1 && waseven==0
      Bll=zeros(length(dels1)-1,length(TH));
      dels=dels1(1:end-1);
      Bll(1:2:end,:)=Bl;
      Bl=Bll;
    elseif evens==0 && waseven==0
      % Strictly speaking, for method 2 this is only necessary for sord=2
      Bl(2:2:end,:)=0;
    end
  end
  
  goods='BPBOXCAP area check passed';
  % Check the B0 term which should equal the area^2 divided by 4pi
  A=4*pi*spharea(TH,sord);
  difer(Bl(1,:)-A.^2/4/pi,[],[],goods)
  
  %disp(sprintf('Elapsed CPU time %8.3f',cputime-t))

elseif strcmp(TH,'demo')
  [b1,l1]=bpboxcap([20 30 40],50,1,1,1);
  [b2,l2]=bpboxcap([20 30 40],50,2,1,1); difer(b1-b2)
  [b1,l1]=bpboxcap([20 30 40],50,1,0,1);
  [b2,l2]=bpboxcap([20 30 40],50,2,0,1); difer(b1-b2)
  [b1,l1]=bpboxcap([20 30 40],50,1,1,2);
  [b2,l2]=bpboxcap([20 30 40],50,2,1,2); difer(b1-b2)
  [b1,l1]=bpboxcap([20 30 40],50,1,0,2);
  [b2,l2]=bpboxcap([20 30 40],50,2,0,2); difer(b1-b2)
  [b1,l1]=bpboxcap([20 30 40],50,1,0,3);
  [b2,l2]=bpboxcap([20 30 40],50,2,0,3); difer(b1-b2)
  [b1,l1]=bpboxcap([20 30 40],50,1,1,3);
  [b2,l2]=bpboxcap([20 30 40],50,2,1,3); difer(b1-b2)
  b=bpboxcap(180,ceil(rand*200),1);
  difer(2*sum(b)-b(1)-4*pi)
end