function [Klmlmp,XY,K1,K]=kernelc(Lmax,dom,pars,ngl,rotb)
% [Klmlmp,XY,K1,K]=KERNELC(Lmax,dom,pars,ngl,rotb)
%
% Calculation of the localization matrix for some domain on the sphere.
% NOT FOR POLAR PATCHES! AND NOT GOOD FOR NEAR-POLAR PATCHES! (See GRUNBAUM)
% NOT WITHOUT MODIFICATIONS FOR REGIONS CONTAINING THE NORTH POLE OR THE
% SOUTH POLE! (See, e.g. 'antarctica'). Watch normalization!
%
% 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', 'gpsnamerica',
%            'antarctica', with specs in 'pars'
%            OR: [lon lat] an ordered list defining a closed curve [degrees]
% 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]
%            OR: the string with the name you want the result saved as
% ngl        The degree of the Gauss-Legendre integration [default: 200] OR
%            'alternative' for an alternative calculation method
% rotb       1 rotate back as required, e.g. for Antarctica [default: 0]
%
% 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, and
%            see, e.g. PLOTSLEP and KLMLMP2ROT for some implementations
% XY         The outlines of the region into which you are localizing
% K1         An intermediate result useful when rotb=1, see KLMLMP2ROT
% K          An verification result useful when rotb=1, see KLMLMP2ROT
% 
% EXAMPLE:
%
% L=19;
% [Klmlmp,XY]=kernelc(L,'australia'); 
% Klmlmp2=kernelc(L,'australia',[],'alternative'); 
% and then DIFER, EIG, PLOTSLEP, etc, to evaluate the difference
%
% kernelc('demo1') % For an illustration of the Antarctica matrix
% kernelc('demo2') % For an illustration of the Antarctica functions
% kernelc('demo3') % For a show of Australia
%
% See also LOCALIZATION, SDWREGIONS, GLMALPHA, DLMLMP, KERNELC2D, PLOTSLEP
%
% Last modified by fjsimons-at-alum.mit.edu, 11/22/2009

t0=clock;
defval('Lmax',18); 
defval('dom','patch')
defval('ngl',200)
defval('rotb',0)
defval('K1',NaN)

if ~isstr(Lmax)
  if isstr(dom)
    switch dom  
     case 'sqpatch'
      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));
     case 'patch'
      fnpl=sprintf('%s/%s-%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));
     otherwise
      fnpl=sprintf('%s/WREG-%s-%i.mat',...
		   fullfile(getenv('IFILES'),'KERNELC'),dom,Lmax);
    end
  else
    fnpl='neveravailable';
  end

  if exist(fnpl,'file')==2 & ~isstr(ngl)
    load(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
	disp('Really, should be putting in the GRUNBAUM call here')
	error('Not for polar caps! Use GRUNBAUM or SDWCAP instead')
	% BUT IN COMPARING, NOTE THAT THE SIGN MAY BE OFF
      end
      if thR>th0
	disp('Really, should be putting in the GRUNBAUM call here')
	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
      % If it's a named geographical region or a coordinate boundary
      if isstr(dom)
	defval('pars',10);
	% Run the named function to return the coordinates
	if strcmp(dom,'antarctica') && rotb==1
	  % Return the rotation parameters also, to undo later
	  [XY,lonc,latc]=eval(sprintf('%s(%i)',dom,pars));
	else
	  % Don't, the result will be the kernel for the rotated dom
	  XY=eval(sprintf('%s(%i)',dom,pars));
	end
      elseif isstr(pars)
	defval('pars','shouldntihaveadefaulthere');
	XY=dom;
	% Use the input to define the file name that will be created
	fnpl=sprintf('%s/WREG-%s-%i.mat',...
		     fullfile(getenv('IFILES'),'KERNELC'),pars,Lmax);
      else
	XY=dom;
      end
      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,mzi,mzo,bigm,bigl]=addmon(Lmax);
    dimK=(Lmax+1)^2; lenm=length(dems);
    Klmlmp=repmat(NaN,dimK,dimK);

    if isstr(ngl) && strcmp(ngl,'alternative')
      % Hold on and see if we can go more quickly here using GEOBOXCAP
      fax=2^5; 
      % Try oversampling the Nyquist degree by a certain factor
      degres=180/sqrt(Lmax*(Lmax+1))/fax;
      % Calculate the mask function
      [Bl,dels,r,lon,lat]=geoboxcap(Lmax,dom,[],degres);
      % Prepare the reindexing arrays
      % We're not fully using the recursion here, so there is wastage  
      % Perform the masked spherical harmonic transform
      h=waitbar(0,'KERNELC: Loop over all degrees and orders');
      % With the recursions as they are, we are not yet taking full
      % advantage of this method. See Mark Wieczorek's Fortran code which
      % presumably works better for this case.
      for l=0:Lmax
	% Remember the normalization conventions
	theYplus=2*sqrt(pi)*ylm(l,0:l,(90-lat)/180*pi,lon/180*pi);
	theYmins=2*sqrt(pi)*ylm(l,-1:-1:-l,(90-lat)/180*pi,lon/180*pi);
	for m=0:l
	  waitbar((addmup(l)+m)/addmup(Lmax),h);
	  % Return the expansion coefficients in "standard" real
          % harmonics order 
	  lmcosiplus=xyz2plm((-1)^m*theYplus(:,:,m+1).*r,Lmax);
	  % Reindex the coefficients to "standard" localization kernel order
	  % and put them into where the positive m sits
	  posm=addmoff(l-1)+2*m+1;
	  % Redundant check that we are at the right order and degree
	  % difer(bigl(posm)-l)
	  % difer(bigm(posm)-m)
	  Klmlmp(posm,:)=lmcosiplus(2*size(lmcosiplus,1)+mzo)';
	  if m>0
	    % Return the expansion coefficients in "standard" real
            % harmonics order 
	    lmcosimins=xyz2plm((-1)^m*theYmins(:,:,m).*r,Lmax);
	    % Also do the negatives which come right before the positives
	    Klmlmp(posm-1,:)=lmcosimins(2*size(lmcosimins,1)+mzo)';
	  end
	end
      end
      delete(h)
    else
      % 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 longitudinal integration info for the domain
      if isstr(dom)
	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
      else
	% Now we may have multiple pairs
	phint=dphregion(acos(x)*180/pi,[],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)
      % Save this now
      if isstr(pars)
	dom=pars;
      end
      % This is only saved when it's not the alternative calculation method
      % nor when it is still required to unrotate the result, as expected
      % by LOCALIZATION, but note we may want to change this convention later
      if ~strcmp(fnpl,'neveravailable') && rotb==0
	eval(sprintf('save %s Lmax Klmlmp dom ngl XY',fnpl))
      end
    end
    disp(sprintf('KERNELC (Matrix)  took %8.4f s',etime(clock,t0)))
  end

  % By whichever way you obtained the kernel, now check if you might want
  % to rotate it back so its eigenvectors are "right" for Antarctica
  if rotb==1
    % Get the rotation parameters for this particular region
    [XY,lonc,latc]=eval(sprintf('%s(%i)',dom,pars));
    if nargout<4
      % Rotate the kernel, properly
      [Klmlmp,K1]=klmlmp2rot(Klmlmp,lonc,latc);
    else
      % Some extra verification in here
      [Klmlmp,K1,K]=klmlmp2rot(Klmlmp,lonc,latc);
    end
  end
elseif strcmp(Lmax,'demo1')
  L=7;
  % Construct the kernel for Antarctica in rotated space
  K0=kernelc(L,'antarctica',[],[],0);
  % Excessive verification
  xver=0;
  if xver==1
    % Construct the kernel for Antarctica in physical space
    [K2,XY,K1,K]=kernelc(L,'antarctica',[],[],1);
    % This better be very nearly the same:
    difer(K-K0)
  else
    [K2,XY,K1]=kernelc(L,'antarctica',[],[],1);
  end

  % Compare the kernels visually
  clf
  strx='rank | degree and order';
  stry1='degree | degree and order';
  stry2='degree | degree and order';
  stry3='order | degree and order';
  [ah,ha,H]=krijetem(subnum(1,3));
  
  c11=[1 1];
  cmn=[(L+1)^2 (L+1)^2];
  axes(ah(1))
  % Regular view, degrees indicated
  imagefnan(c11,cmn,K0); axis ij
  xl(1)=xlabel(strx);
  yl(1)=ylabel(stry1); 
  t(1)=title('The original kernel');
    
  axes(ah(2))
  % Regular view, degrees indicated
  % Should look "zonal", almost all on m=0
  imagefnan(c11,cmn,K1); axis ij
  xl(2)=xlabel(strx);
  yl(2)=ylabel(stry2); 
  t(2)=title('Once rotated');

  axes(ah(3))
  % Sort orders and degrees from KERNELC to ADDMOUT
  [dems,dels,mz,lmc,mzin,mzo,Km,Kl,rinm]=addmon(L);
  % Now [Kl(rinm) Km(rinm)] is like ADDMOUT
  [EM,EL,mz,blkm,dblk]=addmout(L);
  % Now [Kl(rinm(blkm)),Km(rinm(blkm))] is block sorted
  imagefnan(c11,cmn,K2(rinm(blkm),rinm(blkm))); axis ij
  blox=Km(rinm(blkm));
  ordshew=[1; find(diff(Km(rinm(blkm))))+1];
  xl(3)=xlabel(strx);
  yl(3)=ylabel(stry3); 
  t(3)=title('Fully rotated');

  % Cosmetics
  longticks(ah)
  set(ah(1:2),'ytick',addmoff([0:L]-1)+1,'ytickl',[0:L])
  set(ah(3),'ytick',ordshew,'ytickl',num2cell(Km(ordshew)))
  set(ah,'xtick',[1 (L+1)^2],'xtickl',[1 (L+1)^2])
  fig2print(gcf,'landscape')
  figdisp([],sprintf('antarctica_%i',L))
  
elseif strcmp(Lmax,'demo2')
  L=18;
  % Construct the kernel for Antarctica in rotated space
  K0=kernelc(L,'antarctica',[],[],0);
  % Construct the kernel for Antarctica in physical space
  [K2,XY,K1]=kernelc(L,'antarctica',[],[],1);

  % Now check the diagonalization of the first kernel 
  [C0,V0]=eig(K0);
  [V0,isrt]=sort(sum(real(V0),1),'descend');
  C0=C0(:,isrt(1:addmoff(L)));
 
  % Now check the diagonalization of the second kernel
  [C2,V2]=eig(K2);
  [V2,isrt]=sort(sum(real(V2),1),'descend');
  C2=C2(:,isrt(1:addmoff(L)));

  % Compare the eigenvalues
  difer(V0-V2)
  
  % Compare the eigenfunctions!
  % Clean this up!
  J=5;
  [dems,dels,mz,lmc,mzin]=addmon(L);
  [XY,lonc,latc]=antarctica;

  CC=cellnan(J,length(dels),2);
  for index=1:J
    CC{index}=plm2rot(...
	[dels dems reshape(insert(C0(:,index),0,mzin),2,length(dems))'],...
	lonc,-latc,0);
  end
    
  % Do the plotting
  clf
  [ah,ha,H]=krijetem(subnum(J,3));

  for index=1:J
    axes(ha(index))
    plotslep(C0,index,2);
    axes(ha(index+J))
    p=plotslep(C2,index,2);
    axes(ha(index+2*J))
    d=plotplm(CC{index},[],[],4,1);
    set(ah,'clim',halverange(d,100))
    % Let us appraise the comparison also - knowing that the sign remains
    % arbitrary
    a=minmax(abs(d-p)./d);
    b=minmax(abs(d+p)./d);
    disp(sprintf('The min/max relative difference is %8.3e | %8.3e',...
	 [max([a(1) b(1)]) min([a(2) b(2)])]))
  end  

  % Cosmetics
  nolabels(ha(J+1:end),2)
  nolabels(ah(1:end-3),1)
  longticks(ah)
  fig2print(gcf,'tall')
elseif strcmp(Lmax,'demo3')
  % Reconcile the indexing of KERNELC | KLMLMP and GLMALPHA
  Klmlmp=kernelc(18,'australia');
  % Actually, just see inside GLMALPHA for an example, that will do for now.
end