function [Klmlmp,XY]=kernelc(Lmax,dom,pars,ngl)
% [Klmlmp,XY]=KERNELC(Lmax,dom,pars,ngl)
%
% Calculation of the localization matrix for some domain on the sphere.
% NOT FOR POLAR PATCHES! AND NOT GOOD FOR NEAR-POLAR PATCHES!
% NOT FOR REGIONS CONTAINING THE NORTH POLE OR THE SOUTH POLE!
%
% INPUT:
%
% Lmax       Maximum angular degree (bandwidth)
% dom        'patch'   Spherical patch [default], with specs in 'pars'
%                      Note: better use GRUNBAUM / PLM2ROT in this case
%            'sqpatch' Square patch with [thN thS phW phE] in 'pars'
%            'england', 'eurasia',  'namerica', 'australia', 'greenland'
%            'africa', 'samerica', 'amazon', 'orinoco', with specs in 'pars'
% pars       [th0,ph0,thR] for 'patch'
%                 th0  Colatitude of the cap center, in radians
%                 ph0  Longitude of the cap center, in radians
%                 thR  Radius of the cap, in radians
%            N  splining smoothness for geographical regions [default: 10]
% ngl        The degree of the Gauss-Legendre integration [default: 200]
%
% OUTPUT:
%
% Klmlmp     The localization kernel whose eigenfunctions we want,
%            indexed as: degree  [0  1  1  1  2  2  2  2  2]
%                        order   [0  0 -1  1  0 -1  1 -2  2]
%            The function LOCALIZATION later reindexes this in LMCOSI
%            fashion. Note you can use ADDMOUT and ADDMON to modify.
% XY         The outlines of the region into which you are localizing
% 
%
% EXAMPLE:
%
% L=18;
% [Klmlmp,XY]=kernelc(L,'australia');
%
% See also LOCALIZATION, SDWREGIONS, GLMALPHA, DLMLMP
%
% Last modified by fjsimons-at-alum.mit.edu, 08/20/2007

t0=clock;

defval('Lmax',18); 
defval('dom','patch')
defval('ngl',200)

if ~strcmp(dom,'sqpatch')
  fnpl=sprintf('%s/WREG-%s-%i.mat',...
	       fullfile(getenv('IFILES'),'KERNELC'),dom,Lmax);
else
  fnpl=sprintf('%s/%s-%i-%i-%i-%i-%i.mat',...
	       fullfile(getenv('IFILES'),'KERNELC'),dom,Lmax,...
	       round(pars(1)*180/pi),round(pars(2)*180/pi),...
	       round(pars(3)*180/pi),round(pars(4)*180/pi));
end

if exist(fnpl,'file')==2 && ~strcmp(dom,'patch')
  eval(sprintf('load %s',fnpl))
  disp(sprintf('%s loaded by KERNELC',fnpl))
else
  if strcmp(dom,'patch')
    defval('pars',[90 75 30]*pi/180);
    % For future reference 
    th0=pars(1); ph0=pars(2); thR=pars(3);
    if th0==0
       error('Not for polar caps! Use GRUNBAUM or SDWCAP instead')
    end
    if thR>th0
       error('Not for near-polar caps! Use GRUNBAUM, SDWCAP, then rotate')
    end
    % Convert all angles to degrees for CAPLOC only
    [lon,lat]=caploc([ph0 pi/2-th0]/pi*180,thR/pi*180,100,1);
    % Northern and Southern points, in radians
    thN=(th0-thR);
    thS=(th0+thR);
    XY=[lon lat];
  elseif strcmp(dom,'sqpatch')
    defval('pars',[30 90 10 90]*pi/180);
    thN=pars(1); thS=pars(2); phW=pars(3); phE=pars(4);
    XY=[phW pi/2-thN ; phW pi/2-thS ; phE pi/2-thS ; ...
	phE pi/2-thN ; phW pi/2-thN]*180/pi;
  else
    defval('pars',10);
    XY=eval(sprintf('%s(%i)',dom,pars));
    thN=90-max(XY(:,2)); thN=thN*pi/180;
    thS=90-min(XY(:,2)); thS=thS*pi/180;
  end

  % No more set-up after this point
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Introduce and dimensionalize variables and arrays
  [dems,dels,mz,lmc,mzin,mzout,bigm,bigl]=addmon(Lmax);
  dimK=(Lmax+1)^2; lenm=length(dems);
  Klmlmp=repmat(NaN,dimK,dimK);

  % Calculating different Gauss-Legendre points for all possible product
  % degrees is not a good idea since they get multiplied by more
  % functions of theta 
  intv=cos([thS thN]);
  nGL=max(ngl,2*Lmax);
  % These are going to be the required colatitudes - forget XY
  [w,x,N]=gausslegendrecof(nGL,[],intv);
  disp(sprintf('%i Gauss-Legendre points and weights calculated',N))

  % First calculate the Legendre functions themselves
  % Note that dimK==sum(dubs)
  dubs=repmat(2,lenm,1); dubs(mz)=1; comdi=[];
  % First, calculate all the Legendre functions themselves
  Xlm=repmat(NaN,length(x),lenm);

  % Calculate the Legendre polynomials
  ind=0;
  for l=0:Lmax  
    Xlm(:,ind+1:ind+l+1)=[legendre(l,x(:)','sch')*sqrt(2*l+1)]';
    ind=ind+l+1;
  end

  % Calculate the Legendre products for all combinations of l and m
  XlmXlmp=repmat(NaN,length(x),((lenm^2)+lenm)/2);
  GLint=repmat(NaN,((lenm^2)+lenm)/2,1);
  index=0;
  h=waitbar(0,'KERNELC: Calculating all Legendre products');
  for lm1dex=1:lenm
    l1=dels(lm1dex);
    m1=dems(lm1dex);
    pos1=1/2*(l1)^2+1/2*l1+m1+1;
    comdi=[comdi ; dubs(lm1dex:lenm)];
    % Note that the last index will be ((lenm^2)+lenm)/2
    for lm2dex=lm1dex:lenm
      l2=dels(lm2dex);
      m2=dems(lm2dex);
      index=index+1;
      pos2=1/2*(l2)^2+1/2*l2+m2+1;
      % Calculate products of Legendre functions
      XlmXlmp(:,index)=Xlm(:,pos1).*Xlm(:,pos2);
      waitbar(2*index/((lenm^2)+lenm),h);
    end
  end
  delete(h)
  % disp('Legendre products calculated')

  % In our ordering, the -1 precedes 1 and stands for the cosine term 
  comdex=[1:((lenm^2)+lenm)/2]';
  coss=gamini(comdex,comdi);
  % Need a vector of length "index" that points to the right combination in
  % XlmXlmp for the next array we are designing. First, find the positions
  % we've been missing 
  h=[dimK:-1:1']'; k=find(dems); kk=k+[1:length(k)]';
  % Where to insert other elements
  inpo=[indeks(cumsum(skip(h,kk)),k)+1]';
  % How many elements to insert
  inel=h(kk);
  % Which elements to insert
  beg=inpo-h(k)+[1:length(inel)]';
  ent=inpo-h(k)+inel+[0:length(inel)-1]';
  ins=[];
  for ind=1:length(beg)
    ins=[ins coss(beg(ind):ent(ind))];
  end
  % And how to do it
  bigo=insert(coss,ins,gamini(inpo,inel));

  % Get the integration info for the domain
  switch dom
   case 'patch'
    % Get the parameters of the dom
    phint=dphpatch(acos(x),thR,th0,ph0);
   case 'sqpatch'
    % Always the same longitudinal integration interval
    phint=repmat([phW phE],length(x),1);
   otherwise
    defval('Nk',10);
    % Now we may have multiple pairs
    phint=dphregion(acos(x)*180/pi,Nk,dom);
    phint=phint*pi/180;
  end

  % The number of elements that will be calculated is
  nel=(dimK^2+dimK)/2;
  % Calculate matrix elements
  index=0;
  undex=0;
  andex=1;
  h=waitbar(0,'KERNELC: Evaluating integrals and assembling matrix');
  for lm1dex=1:dimK
    l1=bigl(lm1dex);
    m1=bigm(lm1dex);
    index=index+1;
    ondex=0;
    I=repmat(NaN,length(x),dimK-lm1dex+1);
    for lm2dex=lm1dex:dimK
      l2=bigl(lm2dex);
      m2=bigm(lm2dex);
      ondex=ondex+1;
      undex=undex+1;
      waitbar(undex/nel,h);
      % Now evaluate the longitudinal integrals at the GL points
      if m1>0 & m2>0
	I(:,ondex)=sinsin(acos(x),m1,m2,phint);
      elseif m1<=0 & m2<=0
	I(:,ondex)=coscos(acos(x),m1,m2,phint);
      elseif m1>0 & m2<=0 % Got rid of redundant ,pars below here
	I(:,ondex)=sincos(acos(x),m1,m2,phint);
      elseif m1<=0 & m2>0
	I(:,ondex)=sincos(acos(x),m2,m1,phint);
      end
    end
    % Calculate the entire integral and distribute over the kernel
    Klmlmp(index,lm1dex:dimK)=...
	(w(:)'*(XlmXlmp(:,bigo(andex:undex)).*I));
    
    % And symmetrize them
    Klmlmp(lm1dex+1:dimK,index)=Klmlmp(index,lm1dex+1:dimK)';
    andex=undex+1;

    % Verify right away that the first value correctly gives the area of
    % the patch 
    if lm1dex==1
      if strcmp(dom,'patch')
	parea=2*pi*(1-cos(thR));
	apo=abs(parea-Klmlmp(1))/parea;
	disp(sprintf('Area of the patch approximated to within %5.2f %s',...
		     apo*100,'%'))
	if apo*100>1;
	  error('Something wrong with the area element: radians/degrees ?')
	end
      elseif strcmp(dom,'sqpatch')
	parea=[cos(thN)-cos(thS)]*[phE-phW];
	apo=abs(parea-Klmlmp(1))/parea;
	disp(sprintf('Area of the patch approximated to within %5.2f %s',...
		     apo*100,'%'))
	if apo*100>1;
	  error('Something wrong with the area element: radians/degrees ?')
	end
      else
	disp(sprintf('Area of the domain approximated as %8.3e',...
		     Klmlmp(1)))
      end
    end
  end

  % To make this exactly equivalent to Tony's \ylm:
  Klmlmp=Klmlmp/4/pi;
  
  % This then makes Klmlmp(1) the fractional area on the sphere

  delete(h)
  disp(sprintf('KERNELC (Matrix)  took %8.4f s',etime(clock,t0)))
  % Save this now
  if ~strcmp(dom,'patch')
    eval(sprintf('save %s Lmax Klmlmp dom ngl XY',fnpl)) 
  end
end