function swfried2d(method) % swfried2d(method) % % FIGURE 4 of SIMONS & WANG, GEM 2011 % % INPUT: % % method 'SE' by Slepian "extension" [faster!] % 'GL' by direct Gauss-Legendre integration [slow!] % 'RS' by direct Riemann summation [slow!] % 'compare' by comparison of the methods % % Makes a fried-egg plot for +- angular orders for concentration to a % disk-shaped domain. If 'method' is 'GL' The eigenvalues will be ever so % slightly different from the "truth" that we know so well - they may not % occur exactly in pairs, etc. But overall this is pretty well behaved to % a part in a million or so starting from reasonable choice. % % SEE ALSO: % % SWVALS2D % % Last modified by fjsimons-at-alum.mit.edu, 03/03/2011 % eps2png -width 850 swfried2d.eps % convert swfried2d.png swfried2d.gif % Somehow this performs best with a KEYBOARD before the print statement % Here we both simultaneously check and compute the eigenfunctions defval('method','SE'); % Make a circle of some description circn=256; XY=[cos(linspace(0,2*pi,circn)) ; sin(linspace(0,2*pi,circn))]'; % What is the mean distance between the points mdis=mean(sum(sqrt(diff(XY).^2),2)); % Perform the localization N=42; % The Shannon number J=30; % The number of Slepian functions that will be calculated % Define the (resolution of the) output grid j=7; xp=linspace(-2,2,2^j); yp=linspace(-2,2,2^j); [XP,YP]=meshgrid(xp,yp); XYP=[XP(:) YP(:)]; % Now do the actual calculation switch method case {'GL','RS'} % This returns all of the numerical eigenvalus [G,H,V,K,XYP,XY,A]=localization2D(XY,N,J,method,[],[],[],XYP); G=reshape(G,2^j,2^j,J); % Make the figure - with G/E or GV/EV for comparison theG=G; theV=V; case 'SE' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Now we need to compare this to the other way of calculating it % The upper of number of orders should "cover" all the cases above M=8; % Calculate at least all of those that have eigenvalues of larger than 0.5 [Nm,Nsum]=nsubm(N,M); % Form the spatial eigenfunction R=sqrt(XP.^2+YP.^2); TH=atan2(YP,XP); [Rx,iR,jR]=unique(R); % Calculate the radial functions at the unique radial points that will be needed for m=0:M [E{m+1},EV{m+1}]=swdisk(m,N,ceil(Nm(m+1)),[],Rx,'SE',0); EM{m+1}=repmat(m,1,ceil(Nm(m+1))); end % Sort them all and compare eigenfunctions with the same eigenvalue E=[E{:}]; EV=[EV{:}]; EM=[EM{:}]; [EV,i]=sort(EV,'descend'); E=E(:,i); EM=EM(i); % Repeat the nonzero orders twice and take as many as you had dbl=~~EM+1; % The indexing sequence in the non-repeated vectors seq=indeks(gamini(1:length(EV),dbl),1:J); % This will change the sine into the cosine later % If there is no difference, it is a pair addph=~diff([0 seq]); % Make the eigenvalue sequence with the repeats EV=indeks(gamini(EV,dbl),1:J); % Reassemble and combine with the sine/cosine functions GV=nan(size(R,1),size(R,2),J); for i=1:length(EV) % Just leave it here - it's the GL that is in trouble GV(:,:,i)=reshape(E(jR,seq(i)),size(R)).*... cos(EM(seq(i)).*TH-addph(i)*pi/2); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Make the figure - with G/E or GV/EV for comparison theG=GV; theV=EV; case 'compare' keyboard % Not properly written but largely irrelevant now % Check how close these eigenvalues are now difer(EV-V,4) % Also compare the eigenfunctions and make sure that the normalization % is identical in both cases. end % They should differ slightly - in the sign - in the phase clf [ah,ha,H]=krijetem(subnum(6,5)); c11=[min(XYP(:,1)) max(XYP(:,2))]; cmn=[max(XYP(:,1)) min(XYP(:,1))]; defval('bw',1) for index=1:length(ah) axes(ah(index)) % Check if the color saturation is appropriate if bw==1 imagefnan(c11,cmn,setnans(theG(:,:,index)),gray(21),[-1 1],grey(5)) else imagefnan(c11,cmn,setnans(theG(:,:,index)),kelicol,[-1 1]) end hold on % Plot the circle plot(XY(:,1),XY(:,2),'k--') % Plot the square plot(c11([1 2 2 1 1]),cmn([1 1 2 2 1]),'k-') switch method case {'GL','RS'} tl(index)=title(sprintf('%s_{%i} = %9.6f','\lambda',... index,theV(index))); case 'SE' tl(index)=title(sprintf('%s_{%i} = %9.6f ; m = %i',... '\lambda',... index,theV(index),... (-1).^[addph(index)-1]*EM(seq(index)))); end axis off drawnow end % Cosmetics (best after KEYBOARD and manually line by line) movev(tl,-0.3) set(tl,'FontS',8) serre(H',[],'down') fig2print(gcf,'tall') figdisp([],[],[],0) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function data=setnans(data) data(abs(data)