Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
167 lines (140 sloc) 4.55 KB
function [E,V,N,th,C,E0]=sdwcapt(TH,L,m,nth,vcut,grd)
% [E,V,N,th,C,E0]=SDWCAPT(TH,L,m,nth,vcut,grd)
%
% Spherical harmonic localization to a spherical polar cap:
% spatially concentrated and optimally bandlimited solutions.
%
% INPUT:
%
% TH Angular extent of the spherical cap, in degrees
% L Bandwidth
% m Angular order of the required data window, -l<m<l
% nth Number of colatitudes to evaluate
% vcut Cut-off eigenvalue [default= eps*10]
% grd 1 Colatitudes only; returns matrix E
% 2 Colatitude/Longitude; returns cell E
%
% OUTPUT:
%
% E Optimally concentrated tapers, expanded to space domain
% V Eigenvalues, sorted
% N Sum of all the eigenvalues
% th Colatitudes at which the functions are evaluated, in degrees
% C Optimally concentrated tapers; coefficients
% E0 Tapers with zeros added to make them of length nth
%
% SEE ALSO:
%
% CDLKERNEL
%
% Last modified by fjsimons-at-alum.mit.edu, 05/21/2009
defval('TH',40)
defval('L',18)
defval('m',0)
defval('nth',720)
defval('vcut',eps*10)
defval('grd',1)
% Work with the absolute value of m
mor=m;
m=abs(m);
fnpl=sprintf('%s/SDW-%i-%i-%i-%i.mat',...
fullfile(getenv('IFILES'),'SDWCAPT'),TH,L,nth,m);
if exist(fnpl,'file')==2 % & 1==3
eval(sprintf('load %s',fnpl))
disp(sprintf('%s loaded by SDWCAPT',fnpl))
else
N=(L+1)^2*(1-cos(TH/180*pi))/2;
% Approximate Nyquist degree is $\pi/d\theta$
if nth < (L+1) & nth~=0
error('Sample finer to avoid aliasing')
end
% Regular sampling at number of requested points
nTH=ceil(TH/180*nth);
% But you sample only until TH
th=linspace(0,TH/180*pi,nTH);
% Convert to radians
thor=th*180/pi;
TH=TH*pi/180;
% The degree of the integrand will be 2*Lmax
% Gauss-Legendre integration points and weights scaled on the interval
% [cos(TH) 1] after transformation of variables.
ngl=max(2*L,200*(m~=0));
% If you want to be fancy do this
% [jk1,V2,jk2,jk3,jk4,ngl1,ngl2]=sdwcap(TH*180/pi,L,m);
% ngl=ngl2
[w,xGL,nsel]=gausslegendrecof(ngl,[],[cos(TH) 1]);
disp(sprintf('Gauss-Legendre with %i points',nsel))
% Later on, we will need points over the entire interval from cos([0,pi])
% Interpolation points at the requested number of TH's
xint=cos(th);
% Construct the Christoffel-Darboux kernel
[KGL,Kint]=cdlkernel(xGL,xint,L,m);
% The interval rolls out for Gauss-Legendre quadrature
% but you need the scaling; contains longitudinal integral
normz=(1+(m==0))/4;
KGL=normz*KGL;
Kint=normz*Kint;
% Calculate eigenvectors and values of symmetrized kernel
[EGL,V]=eig(diag(sqrt(w))*KGL*diag(sqrt(w)));
% Need to unscale E
EGL=EGL./repmat(sqrt(w(:)),1,size(EGL,2));
% With this, inner product over domain becomes V; no dot!
EGL=EGL*sqrt(V)/sqrt((1+(m==0))*pi);
% Should you ever want to plot E at the GL values.
th=acos(xGL);
% Order eigenvalues and eigenfunctions downgoing
[V,isrt]=sort(sum(V,1),'descend');
% V=fliplr(V); EGL=EGL(:,fliplr(isrt));
EGL=EGL(:,isrt);
% We know how many there should be
EGL=EGL(:,1:L-m+1);
V=V(:,1:L-m+1);
% Only return nonzero and real eigenvalues (order matters)
EGL=EGL(:,V>vcut & imag(V)==0);
V=V(V>vcut & imag(V)==0);
% But usually we "interpolate" them using the same kernel at the desired
% values.
E=(Kint*diag(w))*EGL*diag(1./V);
% Make E start with a positive lobe
for index=1:size(E,2)
E(:,index)=E(:,index)*sign(E(2,index));
end
% Verify automatic orthonormality over restricted interval
% orv1=w(:)'*[EGL.^2]*(1+(m==0))*pi;
% Calculate the latitudinal and longitudinal integral
orv2=[EGL'*diag(w)*EGL]*(1+(m==0))*pi;
disp(sprintf('Orthonormality: Mean absolute error %8.3e',...
mean(mean(abs(orv2-diag(V))))))
% Output coefficients C
E0=[E ; zeros(nth-size(E,1),size(E,2))];
if m==0
% These are normalized to the AREA of the unit sphere
[C,rnk]=th2pl(E0,[],'im');
% Now they're normalized to unit area for the unit sphere
C=C*sqrt(4*pi);
else
% Make it -2 today, August 11th, 2004
[C,rnk]=th2plm(E0,nth-2,m);
C=[C ; repmat(NaN,1,size(C,2))];
end
nlon=ceil(2*nth-1);
% Output in degrees
TH=TH*180/pi;
th=thor;
eval(sprintf('save %s E V L N TH C th nlon E0',fnpl))
end
if grd==2
% Output on full grid
if mor<=0
EE=E; clear E
for index=1:size(EE,2)
E{index}=EE(:,index)*cos(m*linspace(0,2*pi,nlon));
end
end
if mor>0
EE=E; clear E
for index=1:size(EE,2)
E{index}=EE(:,index)*sin(m*linspace(0,2*pi,nlon));
end
end
end