function [Bl,dels,r,lon,lat,lmcosi]=geoboxcap(L,dom,N,degres,act)
% [Bl,dels,r,lon,lat,lmcosi]=GEOBOXCAP(L,dom,N,degres,act)
%
% Returns the spherical harmonic power spectrum of an all-or-nothing
% boxcar over a particular geographical region.
%
% INPUT:
% 
% L          Bandwidth of the output spectrum
% dom        'africa',  'eurasia',  'namerica', 'australia', 'greenland',
%            'samerica', 'amazon', 'orinoco', 'gpsnamerica', 'antarctica'
%            OR: [lon lat] an ordered list defining a closed curve [degrees]
% N          The splining smoothness of the boundary [default: 10] 
% degres     The resolution of the underlying spatial grid [default: Nyquist]
% act        1 Actually perform the calculations [default]
%            0 Don't, but simply return the mask function
%
% OUTPUT:
%
% Bl         The power spectrum 
% dels       The spherical harmonic degrees
% r          The "mask" with the "continent function"
% lon,lat    The grid on which  this is defined
% lmcosi     The complete set of spherical harmonic expansion coefficients
%
% EXAMPLE:
%
% [Bl1,dels1]=geoboxcap(18,'australia',[],[]);
% [Bl2,dels2]=geoboxcap(18,'australia',[],1);
%
% SEE ALSO:
%
% BPBOXCAP, KERNELC
% 
% Last modified by fjsimons-at-alum.mit.edu, 10/10/2009

% Default inputs
defval('L',18)
defval('dom','australia')
defval('N',10)
defval('act',1)

if isstr(dom)
  % Run the named function to return the coordinates
  XY=eval(sprintf('%s(%i)',dom,N));
else
  % Get the coordinates as define from the input in degrees
  XY=dom; 
end

% Make sure the coordinates make sense
XY(:,1)=XY(:,1)-360*[XY(:,1)>360];

% Make a grid of ones and zeroes depending on the desired resolution
degN=180/sqrt(L*(L+1));
defval('degres',degN);

% Default grid is all of the planet
defval('c11cmn',[0 90 360 -90]);

% The number of longitude and latitude grid points that will be computed
nlon=ceil([c11cmn(3)-c11cmn(1)]/degres+1);
nlat=ceil([c11cmn(2)-c11cmn(4)]/degres+1);

% Longitude grid vector in degrees
lon=linspace(c11cmn(1),c11cmn(3),nlon);
% Latitudex grid vector in degrees
lat=linspace(c11cmn(2),c11cmn(4),nlat);
% Make the input grid
r=repmat(0,nlat,nlon);
[LON,LAT]=meshgrid(lon,lat);

% Now decide is we're inside or outside of the region
r(inpolygon(LON,LAT,XY(:,1),XY(:,2)))=1;

if act==1
  % And now do the spherical harmonic transform but only to L
  lmcosi=xyz2plm(r,L);
  
  % Calculate the power spectral density
  [Bl,dels]=plm2spec(lmcosi,2);
else
  [Bl,dels,lmcosi]=deal(NaN);
end

% If these ever used to build loops etc they must be a row
dels=dels(:)';
% Provide output as desired
varns={Bl,dels,r,lon,lat,lmcosi};
varargout=varns(1:nargout);