-
Notifications
You must be signed in to change notification settings - Fork 3
/
mvarratios.m
185 lines (172 loc) · 6.85 KB
/
mvarratios.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
function [mt2ws,l,mt2wsinf]=mvarratios(L,TH,sord,l,xver)
% [mt2ws,l,mt2wsinf]=MVARRATIOS(L,TH,sord,l,xver)
%
% Calculates the ratio of the variance of the eigenvalue-weighted
% multitaper estimates to the untapered whole-sphere estimation variance
% which we can derive from first-principle statistics...
% Dahlen and Simons (2008) eq. (178)
%
% INPUT:
%
% L Bandwidth of the multitaper windows, always a scalar
% TH Colatitudinal radius of the cap, in degrees <=180, may be vector
% Colatitudinal halfwidth of the cut, degrees <=90, may be vector
% sord 1 Single cap of diameter 2TH [default], always a scalar
% 2 Double cap left by subtracting belt of width 2TH, always scalar
% l Degrees at which the calculation is to be carried out
%
% OUTPUT:
%
% mt2ws Ratio of the estimation variance of the eigenvalue-weighted
% multitaper estimate to the untapered whole-sphere estimate
% This is a matrix of dimensions length(TH) x length(l)
% l The degrees at which this variance ratio was computed [0->50]
% mt2wsinf Asymptotic approximation to the above for degrees >> bandwidth
%
% SEE ALSO: PERIODOVAR, MULTIVAR, MTVAR
%
% Last modified by fjsimons-at-alum.mit.edu, 04/29/2007
% SHOULD USE THIS FUNCTION IN DSVARRATIO3
% Supply default values
defval('L',20)
defval('TH',30)
defval('sord',1)
defval('l',0:50)
defval('xver',1)
% Highest degree of a zero-j database, after that, switch to approximation
zjmax=500;
% Highest degree of a six-j database, after that, switch to scaled
% full-sphere result - actually, with the improvements, best switch to
% method 3 if you can do this
sjmax=40;
% The maximum degree
Lmax=min(max(l),zjmax);
% Get all the ZEROJ coefficients at the same time
[allW,C0,S0,Leff]=zeroj(repmat(0:2:2*L,1,Lmax+1),...
gamini(0:Lmax,L+1),gamini(0:Lmax,L+1));
% NEED TO BUILD IN WHAT HAPPENS WHEN A IS EQUAL TO ZERO
% Calculate the area under the taper
A=4*pi*spharea(TH,sord);
if sord==3
error('This choice is not allowed')
end
if 2*L<=sjmax
% Calculate the matrix that goes into this - this takes most of the time
% Always only get the evens since we're studying l=l, the variance
[Gp,p,K]=gammap(L,TH,sord,1,1);
disp(sprintf('We are getting a Shannon number of %5.2f',K))
needscaling=0;
elseif 2*L>sjmax
% Cannot do the full works for lack of the 6j coefficients
% But instead simply scale the A=4pi multitaper result
% Initialize some important arrays
bigS=gamini([0:L],(L+1))';
bigSp=repmat([0:L]',(L+1),1);
cG=repmat(NaN,2*L+1,1);
% Calculate the curly Gamma instead of the regular gamma
for e=0:2*L
cG(e+1,1)=sum([2*bigS'+1].*[2*bigSp'+1].*...
zeroj(bigS,e,bigSp,Leff,[],C0,S0).^2);
end
% Turns out that this "curly Gamma" is what regular Gamma goes to when A
% goes to the whole sphere (A=4pi) --- compare only the even degrees
% Under the single cap thing: single cap of entire globe
Gp=cG(1:2:end)'*4*pi/(L+1)^4;
p=0:2:2*L;
% After this, no sense in verifying since the verification wants to
% know the truth, which you know you can't know
xver=0;
% This part of the loop was verified to give the same results even for
% small L as intended
needscaling=1;
end
% Initialize the mt2ws variance ratio
mt2ws=repmat(NaN,length(TH),length(l));
% Better get all of the wigner0j symbols at once here
for ixl=1:length(l)
if l(ixl)<=Lmax
% Stick with the one-blow precalculated ones if it's below the bandwidth
W=allW((L+1)*l(ixl)+1:(L+1)*(l(ixl)+1));
else
% Calculate on the fly rather than from the loaded database
W=wigner0j(2*L,l(ixl),l(ixl),xver);
W=W(1:2:end);
% Compare to the asymptotic representation if you want
md=sum(abs(W-(-1)^l(ixl)*[plm(0:2:2*L,0,0)./sqrt((2*l(ixl)+1))]')/(L+1));
% disp(sprintf(...
% 'With the asymptotic form you''d have errors of %5.3e but luckily' ...
% ' you didn''t'],md))
end
% Only select the evens since we're doing VARIANCE at equal l=l'
mt2ws(:,ixl)=(2*l(ixl)+1)/(4*pi)*[repmat(2*p+1,length(TH),1).*Gp]*[W.^2]';
end
% Now use scaling if you need to to rectify
sexp=-0.88;
if needscaling==1
mt2ws=mt2ws*(A/4/pi)^sexp;
disp(sprintf('MVARRATIOS Using scaling (A/4%s)^%5.3f','\pi',sexp))
end
% FOR SMALL L VERIFY THAT THIS APPROACH GIVES REASONABLE RESULTS BOTH FOR
% THE SINGLE AND THE DOUBLE CAP NOTE THAT K MUST BE VERY BIG
% mvarratios(20,10,2,0:50) % is the double cap with cut 10
% Now come in here and hand-edit to see that it works.
% What should the large-l limit be? Use the Legendre functions instead
% Initialize the large-l limit of this
mt2wsinf=[1/(4*pi)*[repmat(2*p+1,length(TH),1).*Gp]*...
plm(0:2:2*L,0,0).^2]';
%disp(sprintf(['MVARRATIOS: Error of l=inf approximation at'...
% ' l/L=%3i is %4.1e%s'],l(end)/L,...
% abs(mt2ws(end)-mt2wsinf*((A/4/pi)^(sexp*(needscaling==1))))/...
% mt2ws(end)*100,'%'))
if xver==1 && l(1)==0 && ~[sord==1 & TH==0]
% What should the variance ratio be at l=0 for each of these TH's?
% We have three ways to calculate this
mt2ws01=repmat(NaN,1,length(TH));
mt2ws02=repmat(NaN,1,length(TH));
mt2ws03=repmat(NaN,1,length(TH));
% Initialize some important arrays
bigS=gamini([0:L],(L+1))';
bigSp=repmat([0:L]',(L+1),1);
cG=repmat(NaN,2*L+1,1);
% Calculate the "curly Gamma"
for e=0:2*L
cG(e+1,1)=sum([2*bigS'+1].*[2*bigSp'+1].*...
zeroj(bigS,e,bigSp,Leff,[],C0,S0).^2);
end
% Turns out that this curly Gamma is what regular Gamma goes to when A
% goes to the whole sphere (A=4pi) --- compare only the even degrees
% Under the single cap thing: single cap of entire globe
difer(cG(1:2:end)'*4*pi/(L+1)^4-gammap(L,180,1,1,1),[],[],...
'Check for A=4pi from single cap passed')
% Under the double cap thing: subtract belt of nothing
difer(cG(1:2:end)'*4*pi/(L+1)^4-gammap(L,0,2,1,1),[],[],...
'Check for A=4pi from double cap passed')
for inx=1:length(TH)
% Use the power spectrum of the single/double cap
[BeL(inx,:),dels]=bpboxcap(TH(inx),2*L,[],0,sord);
% Compute this alternative
mt2ws01(inx)=1/4/pi/K(inx)^2*[(2*[0:2*L]+1).*BeL(inx,:)]*cG;
end
% Compare just for good measure
difer(mt2ws01'-mt2ws(:,1),[],[],...
'MVARRATIOS: First check for zero-l passed')
% Another check for the l=0 ratio
for inx=1:length(TH)
[G,V,EL,EM,KN]=glmalpha(90*(sord==2)+(-1)^(sord+1)*TH(inx),L,sord);
mt2ws02(inx)=sum(V.^2)/sum(V)^2;
end
difer(mt2ws02'-mt2ws(:,1),[],[],...
'MVARRATIOS: Secnd check for zero-l passed')
% There is a third check which we should put in there as well
% See also Wieczorek and Simons 2005!
if sord==1
for inx=1:length(TH)
% Get the square localization matrix for axisymmetric domains
D=dlmlmp(TH(inx),L,L,xver);
% Have since rewritten this for the double cap, as well
mt2ws03(inx)=sum(abs(D(:).^2))/K(inx).^2;
end
difer(mt2ws03'-mt2ws(:,1),[],[],...
'MVARRATIOS: Third check for zero-l passed')
end
end