-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathswspace2d.m
150 lines (134 loc) · 4 KB
/
swspace2d.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
function swspace2d(method,scalem)
% SWSPACE2D(method,scalem)
%
% FIGURE 1 of SIMONS & WANG, doi: 10.1007/s13137-011-0016-z
%
% Spatial rendition of the radial part of the 2D Cartesiandisk tapers scaled
% to the unit interval. Checks that finite-length effects are not so grave
% by ensuring that the "small" Fourier components could be "zero" after all.
%
% INPUT:
%
% method 'SE' by Slepian "extension" [default]
% 'GL' by direct Gauss-Legendre integration
% scalem 1 Scales solution for weightless orthogonality
% 0 Scales solution for orthogonality with weight x [default]
%
% SEE ALSO:
%
% SWDISK, SWSPECTRAL2D, and SDWSPACE, SDWSPECTRAL
%
% Last modified by fjsimons-at-alum.mit.edu, 07/27/2022
% The method of computation
defval('method','SE');
% The method of scaling, if 0 don't do sqrt(x) scaling of Slepian
defval('scalem',0);
% Number of tapers per order
NM=4;
% Maximum order considered
M=4;
% Shannon number in question
N=42;
% The abscissas - as long as you like to display though see the
% wavenumber-domain representation and SWSPECTRAL2D for the implications
% of changing these values
x=linspace(0,5,2^12);
% This only if you do anything in the wavenumber domain
[fax,selekt]=fftaxis1D(x,length(x),max(x));
kax=2*pi*fax;
% Limiting wavenumber in question
K=2*sqrt(N);
% Bandlimiting array on the full set of wavenumbers
selK=abs([kax ; -flipud(kax(2:end-1))])>K;
% Start the computation and the figure
clf
[ah,ha,H]=krijetem(subnum(M+1,NM));
for m=0:M
% The above are the unscaled functions "phi" or thus "g"
[E1,V1{m+1},Nm,c,x,E2]=swdisk(m,N,NM,[],x,method,scalem);
% This is the integral over the restricted interval
for inds=1:length(V1{m+1})
werisit=x(:)<=1 & ~isnan(E1(:,inds));
difer(2*pi*trapeze(x(werisit),...
x(werisit)'.*E1(werisit,inds).^2)-V1{m+1}(inds))
end
E1m(m+1)=max(abs(E1(:)));
if strcmp(method,'GL') && scalem==1
% Calculate the truly bandlimited spectrum
FE1=fft(E1); FE1(selK,:)=0;
% And calculate the truly bandlimited eigenfunctions
E1F=ifft(FE1);
% If the input was Hermitian you're all set
difer(imag(E1F))
% Turns out in this case you almost get the input back - but only
% when the normalization is unweighted, if not you'd probably have to
% put the sqrt(x) in there with the ifft somehow.
end
% Plot the space-domain solution
for ondex=1:NM
axes(ah(m*NM+ondex))
pm(m+1,ondex)=plot(x,E1(:,ondex),'-','Color',grey(6));
hold on
pg(m+1,ondex)=plot(x(x<=1),E1(x<=1,ondex),'-','Color','k');
if strcmp(method,'GL') && scalem==1
% Plot the ifft with the small components removed
pgg(m+1,ondex)=plot(x,E1F(:,ondex),'-','Color','r');
end
drawnow
end
end
% Cosmetics
set(ah,'xlim',[0 max(x)],'xtick',[0 1 max(x)],'xgrid','on','ygrid','on',...
'ylim',[-1.2 1.2]*max(E1m))
nolabels(ah(1:end-NM),1)
nolabels(ha(M+2:end),2)
longticks(ah,1/2)
% Eigenvalue labels
nf=9;
lox={'ll','ur','lr','ul',...
'll','ur','lr','ul',...
'll','ur','lr','ul',...
'll','ur','lr','ul',...
'll','ur','lr','ul'};
for m=0:M
for ondex=1:NM
axes(ah(m*NM+ondex))
t{m+1,ondex}=sprintf('%s = %9.6f','\lambda',V1{m+1}(ondex));
if V1{m+1}(ondex)<0.975
lox='lr';
else
lox='ur';
end
[bh(m+1,ondex),th(m+1,ondex)]=...
boxtex(lox,ah(m*NM+ondex),t{m+1,ondex},nf,[],0.75);
end
end
set(th,'FontSize',nf-1)
for index=NM*M+1:NM*(M+1)
axes(ah(index))
xl(index)=xlabel(sprintf('radius %s=r/R','\xi'));
end
for index=1:M+1
axes(ha(index))
yl(index)=ylabel(sprintf('m = %i',index-1));
end
serre(H,1/3,'across')
% So here this goes wrong in the second version
for index=1:NM
axes(ah(index))
tlb(index)=title(sprintf('%s = %i','\alpha',index));
end
set([pm(:) ; pg(:)],'LineW',1)
set(ah,'Fontsize',nf)
fig2print(gcf,'portrait')
% set(gcf,'Color','w','Inv','off')
shrink(ah,1,1/1.1)
if strcmp(method,'GL') && scalem==1
set([pgg(:)],'LineW',0.5)
q=input('Delete red lines? [y] ','s');
defval('q','y')
if strcmp(q,'y')
delete(pgg)
end
end
figdisp([],[],[],0)