Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
257 lines (239 sloc) 8.96 KB
function varargout=ptoslep(phi,theta,omega,TH,L,m,amax)
% [lm,cosi,glma,EL,EM,V]=PTOSLEP(phi,theta,omega,TH,L,m,amax)
%
% Computes SINGLE-CAP FIXED-ORDER azimuthally rotated bandlimited or
% passband Slepian tapers centered at a desired geographical location
% from an initial North-Polar location, that is.
%
% INPUT:
%
% phi Longitude of the new center (degrees)
% theta Colatitude of the new center (degrees)
% omega Clockwise azimuthal rotation (degrees)
% This is consistent with PLM2ROT(lmcosi,omega,-theta,-phi)
% TH Radius of the concentration region (degrees)
% L Bandwidth of the Slepian basis, or passband (two degrees)
% m Single angular order of the Slepian basis function, -l=<m<=l
% amax The number of tapers you want [default: all available]
%
% OUTPUT:
%
% lm The degrees and orders of the expansion, an (L+1)^2 X 2 array
% cosi Real spherical harmonic expansion coefficients suited for
% plotting by PLM2XYZ; this is a (L+1)^2 X 2 X amax array
% glma The same as if they came out of GLMALPHA, (L+1)^2 X amax
% After rotation of course these aren't single-order anymore,
% so block-ordering is not an option here
% EL,EM An array with degrees and orders belonging to the latter
% V The eigenvalues or quadratic concentration fractions
%
% EXAMPLE:
%
% ptoslep('demo1') % Some random rotations
% ptoslep('demo2') % Rotations by pi/2/m
% ptoslep('demo3') % Various rotations over omega
% ptoslep('demo4') % Check that unrotated polar caps are same as GLMALPHA
% ptoslep('demo5') % Check that rotations over pi/2/m are as GLMALPHA
% ptoslep('demo6') % Check the results for changed order signs
%
% See also: SDWTARGET, GLMALPHAPTO
%
% Last modified by fjsimons-at-alum.mit.edu, 06/20/2011
% Should I be able to predict the relation between -m (cos) and +m (sin)
% expansion coefficients - AFTER MAKING THESE COMPLEX? Remember, GRUNBAUM
% and SDWCAP are the same for +/- m. But then they go in at the right
% orders into lmcosi. Then rotated. If then made complex, the relation
% should be easy. But should look at UMMP again to see exactly in what
% order we need to present the coefficients.
defval('phi',78)
if ~isstr(phi)
defval('theta',78)
defval('omega',10)
defval('TH',20)
defval('L',18)
defval('m',0)
% Figure out if it's lowpass or bandpass
lp=length(L)==1;
bp=length(L)==2;
maxL=max(L);
% The spherical harmonic dimension
ldim=(L(2-lp)+1)^2-bp*L(1)^2;
defval('amax',maxL-max(abs(m),bp*min(L))+1);
% Amax can't be more than however many are available
amax=min(amax,maxL-max(abs(m),bp*min(L))+1);
% Calculate the taper coefficients centered on the pole
if lp
[E,Vg,th,C,T,V]=grunbaum(TH,L,m,0);
elseif bp
% Note that the small-eigenvalue eigenfunctions might be
% numerically degenerate and thus not as using Grunbaum - if
% you need to compare, compare where the eigenvalues are "big"
[E,V,Np,th,C]=sdwcap(TH,L,m,0,-1);
else
error('The degree range should be either one or two numbers')
end
% Make a blank array of indices and coefficients
[em,el,mzero,jk1,jk2,jk3,jk4,jk5,jk6,ronm]=addmon(maxL);
lm=[el em];
% The final arrays
cosi=repmat(0,[length(em) 2 amax]);
glma=repmat(0,[(maxL+1)^2 amax]);
% The intermediate array
lmcosi=[el em repmat(0,length(em),2)];
% Do it!
for index=1:amax
% Put the zonal coefficients in the right place
% Negative orders should be the cosine tapers
lmcosi(mzero(abs(m)+1:end)+abs(m),3+(m>0))=C(:,index);
% Rotate the tapers to a position on the sphere
tempry=kindeks(plm2rot(lmcosi,omega,-theta,-phi),3:4);
% Put them into the lmcosi type matrix
cosi(:,:,index)=tempry;
% Put them into the not-block-ordered GLMALPHA type matrix
glma(:,index)=tempry(ronm);
end
% Supply the indices for the GLMALPHA type matrix
[EM,EL]=addmout(maxL);
% Assign output
varns={lm,cosi,glma,EL,EM,V};
varargout=varns(1:nargout);
elseif strcmp(phi,'demo1')
phi=round(rand*180)+90;
theta=round(rand*90)+45;
omega=round(rand*360);
TH=30; L=18; m=1;
[lm,cosi]=ptoslep(phi,theta,omega,TH,L,m,4);
clf
% Calculate cap coordinates
[lon1,lat1]=caploc([phi 90-theta],TH);
% Calculate azimuth coordinates, define EAST from SOUTH
[lon2,lat2]=grcazim([phi 90-theta],...
linspace(0,TH*fralmanac('DegDis')/1000,100),180+omega);
for index=1:4
ah(index)=subplot(2,2,index);
r=plm2xyz([lm cosi(:,:,index)],1);
imagef([],[],r)
set(gca,'xtick',[0 phi 360],'xtickl',[0 phi 360])
set(gca,'ytick',[-90 90-theta 90],'ytickl',[-90 90-theta 90])
hold on
p1(index)=plot(lon1,lat1,'k','LineW',1);
p2(index)=plot(lon2,lat2,'k','LineW',1);
end
supertit(ah(1:2),sprintf('lon= %i ; lat= %i ; rot = %i',phi,90-theta,omega));
movev(ah(1:2),-0.025)
longticks(ah)
set(ah,'xgrid','on','ygrid','on')
elseif strcmp(phi,'demo2')
phi=round(rand*180)+90; phi=180;
theta=round(rand*90)+45; theta=90;
% Pick parameters; positive order m is
TH=30; L=12; m=-2;
% Make a "sine taper" from a "cosine taper"
omega=180/2/abs(m);
[lm,cosi1]=ptoslep(phi,theta,omega,TH,L,m,2);
% Make the "sine taper" by flipping the sign of the order
[lm,cosi2]=ptoslep(phi,theta,0,TH,L,-m,2);
% Calculate cap coordinates
[lon1,lat1]=caploc([phi 90-theta],TH);
% Calculate azimuth coordinates, define EAST from SOUTH
[lon2,lat2]=grcazim([phi 90-theta],...
linspace(0,TH*fralmanac('DegDis')/1000,100),180-omega);
% Start plot subdivision
clf
[ah,ha]=krijetem(subnum(2,2));
% Plot the rotated cosine taper
for index=1:2
axes(ah(index))
r=plm2xyz([lm cosi1(:,:,index)],1);
imagef([],[],r)
set(gca,'xtick',[0 phi 360],'xtickl',[0 phi 360])
set(gca,'ytick',[-90 90-theta 90],'ytickl',[-90 90-theta 90])
hold on
p1(index)=plot(lon1,lat1,'k','LineW',1);
p2(index)=plot(lon2,lat2,'k','LineW',1);
end
hold off
% Plot the sine taper
for index=1:2
axes(ah(index+2))
r=plm2xyz([lm cosi2(:,:,index)],1);
imagef([],[],r)
set(gca,'xtick',[0 phi 360],'xtickl',[0 phi 360])
set(gca,'ytick',[-90 90-theta 90],'ytickl',[-90 90-theta 90])
hold on
p1(index)=plot(lon1,lat1,'k','LineW',1);
end
hold off
supertit(ah(1:2),sprintf('lon= %i ; lat= %i ; rot = %i',phi,90-theta,omega));
movev(ah(1:2),-0.025)
longticks(ah)
set(ah,'xgrid','on','ygrid','on')
elseif strcmp(phi,'demo3')
phi=round(rand*180)+90; phi=180;
theta=round(rand*90)+45; theta=90;
% Pick parameters; positive order m is the cosine taper
TH=30; L=12; m=-1;
% Make a "sine taper" from a "cosine taper"
omega=linspace(0,360,8);
clf
% Calculate cap coordinates
[lon1,lat1]=caploc([phi 90-theta],TH);
% Start plot subdivision
[ah,ha]=krijetem(subnum(2,4));
% Plot the rotated cosine taper
for index=1:length(omega)
axes(ah(index))
% Make the taper
[lm,cosi]=ptoslep(phi,theta,omega(index),TH,L,m,1);
% Calculate azimuth coordinates, define EAST from SOUTH
[lon2,lat2]=grcazim([phi 90-theta],...
linspace(0,TH*fralmanac('DegDis')/1000,100),...
180-omega(index));
% Perform the expansion
r=plm2xyz([lm cosi],1);
imagef([],[],r)
set(gca,'xtick',[0 phi 360],'xtickl',[0 phi 360])
set(gca,'ytick',[-90 90-theta 90],'ytickl',[-90 90-theta 90])
hold on
p1(index)=plot(lon1,lat1,'k','LineW',1);
p2(index)=plot(lon2,lat2,'k','LineW',1);
x(index)=xlabel(sprintf('rot = %i',round(omega(index))));
end
supertit(ah(1:length(omega)/2),sprintf('lon= %i ; lat= %i',phi,90-theta));
movev(ah(1:length(omega)/2),-0.025)
longticks(ah)
set(ah,'xgrid','on','ygrid','on')
fig2print(gcf,'landscape')
elseif strcmp(phi,'demo4')
% Check that output from GLMALPHA equals PTOSLEP at North Pole and omega=0
TH=30; L=18; m=-L:L;
% Figure out where the equivalent sits in the output of GLMALPHA
alpha=cumsum([1 L+1 gamini(L:-1:1,2)]);
% See GLMALPHA for a bandpass update
for index=1:length(m)
[lm,cosi,glma,EL,EM]=ptoslep(0,0,0,TH,L,m(index),1);
[G,V,EL2,EM2]=glmalpha(TH,L,1);
difer(glma-G(:,alpha(2*abs(m(index))+(m(index)>=0))));
end
elseif strcmp(phi,'demo5')
% The first cosine coefficient, rotated by omega+180/(2|m|)...
m=-1; omega=30; L=18; phi=78; theta=78; TH=30;
[a,b,G1]=ptoslep(phi,theta,omega+180/2/abs(m),TH,L,m,1);
% ... should be the same as the sine coefficient, rotated by just omega
m=1;
[a,b,G2]=ptoslep(phi,theta,omega,TH,L,m,1);
difer(G1-G2)
% Next check this with GLMALPHAPTO
elseif strcmp(phi,'demo6')
% What should be the relation between +/-m?
% Well, the result is NOT simply a rotation around the z-axis, so much
% is clear, but rather an azimuthal rotation around the center
m=3;
[lm1,cosi1,glma1,EL1,EM1,V1]=ptoslep(40,70,10,30,12,m,1);
[lm2,cosi2,glma2,EL2,EM2,V2]=ptoslep(40,70,10+180/2/abs(m),30,12,-m,1);
difer(lm1-lm2); difer(EL1-EL2); difer(EM1-EM2); difer(V1-V2)
difer(glma1-glma2); difer(cosi1-cosi2)
clf
subplot(211); plotplm([lm1 cosi1(:,:,1)],[],[],4,1)
subplot(212); plotplm([lm2 cosi2(:,:,1)],[],[],4,1)
end
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.