Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
86 lines (79 sloc) 2.25 KB
function [E,V,th,C]=sdwcapt2(TH,L,m,nth,vcut,grd)
% [E,V,th,C]=SDWCAPT2(TH,L,m,nth,vcut,grd)
%
% Spherical harmonic localization to DOUBLE 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
%
% Last modified by fjsimons-at-alum.mit.edu, 10/03/2008
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);
if(m>L)
error('Order cannot exceed degree')
end
% Filename of saved data
fnpl=sprintf('%s/SDW-%i-%i-%i-%i.mat',...
fullfile(getenv('IFILES'),'SDWCAPT2'),TH,L,nth,m);
if exist(fnpl,'file')==2
eval(sprintf('load %s',fnpl))
disp(sprintf('%s loaded by SDWCAPT2',fnpl))
else
[E,V,th,C,ngl1,ngl2,K]=sdwcap2(TH,L,m,nth,vcut,grd);
% Use the "trick": expand truncated basis functions
% Put on the black lines
belt=[th>TH]&[th<180-TH];
E(belt,:)=0;
if m==0
% Expand to Nyquist; watch for the factors!
C2=th2pl(E,size(E,1)-1)*sqrt(4*pi);
else
C2=th2plm(E,size(E,1)-1,m);
end
% Don't need to do this, but it should return the zeroed input E!
% Reexpand
% E2=pl2th(C2,[],1); % Watch for the factors!
% Now E with the zeroed portions and E2 should be identical
% and the C2's are still orthogonal
C=C2;
% Output in degrees
eval(sprintf('save %s E V L TH C th',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