-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathdsbpboxcap.m
131 lines (117 loc) · 3.31 KB
/
dsbpboxcap.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
function dsbpboxcap
% DSBPBOXCAP
%
% Makes a plot of the boxcar power spectrum
% Simons & Dahlen (2007), Figure 2
%
% Last modified by fjsimons-at-alum.mit.edu, 12/19/2006
% For some reason, the second time round, the plot looks right, but not
% the first... rerun!
TH1=[10 20 30];
[Bl1,dels]=bpboxcap(TH1,100,[],0,1);
TH2=[30 20 10];
[Bl2,dels]=bpboxcap(TH2,100,[],0,2);
Bl2(2:2:end,:)=NaN;
for index=1:length(TH1)
% Note that this is COLUMN per COLUMN !!
Bl1(:,index)=decibel(Bl1(:,index));
Bl2(:,index)=decibel(Bl2(:,index));
% Make sure the max is the B0 term
difer(max(Bl1(:,index))-Bl1(1,index))
difer(max(Bl2(:,index))-Bl2(1,index))
end
% Make the figure
clf
[ah,ha]=krijetem(subnum(3,2));
minb=min([Bl1(:) ; Bl2(:)]);
maxb=-minb;
tix=[0:20:maxb];
widr=[0 5];
% Single cap
for index=1:3
axes(ha(index))
b1(index)=bar(dels,Bl1(:,index)-minb,1);
set(gca,'ytick',sort(maxb-tix),'ytickl',fliplr(-tix))
set(gca,'ylim',minmax(Bl1(:,index)-minb)+widr)
xl(index)=xlabel('degree p');
yl(index)=ylabel('power B_p/B_0 (dB)');
end
% Double cap
for index=1:3
axes(ha(index+3))
b2(index)=bar(dels,Bl2(:,index)-minb,1);
set(gca,'ytick',sort(maxb-tix),'ytickl',fliplr(-tix))
set(gca,'ylim',minmax(Bl2(:,index)-minb)+widr)
xl(index+3)=xlabel('degree p');
yl(index+3)=ylabel('power B_p/B_0 (dB)');
end
% Cosmetics
set(ah,'xlim',[-1 101],'ylim',[30 maxb]+widr)
set([b1 b2],'FaceC',grey)
longticks(ah)
serre(ah(1:2),1/2,'across')
serre(ah(3:4),1/2,'across')
serre(ah(5:6),1/2,'across')
% For some reason this doesn't work
serre(ha(1:3),1/3,'down')
serre(ha(4:6),1/3,'down')
% Leave to last
nolabels(ah(1:4),1)
nolabels(ha(4:6),2)
delete(xl([1 2 4 5]))
delete(yl([4:6]))
% Make miraculous adjustment
set(ah(6),'position',...
[getpos(ah(6),1) getpos(ah(5),2) ...
getpos(ah(6),3) getpos(ah(5),4)])
% Logo scale and location
logsk=20;
cloc='ll';
% Put in equivalent wavelength labels
for index=1:3
% Get equivalent wavelengths for the single cap
for N=1:5
skel(N)=max(roots([1 1 -(N*180./TH1(index))^2]));
end
% Put in cap size labels
legsi{index}=sprintf(' %s = %i%s','\Theta',TH1(index),str2mat(176));
[bh(index),th(index)]=boxtex('ur',ha(index),legsi{index},12);
% Put on top labels
[ax(index),axl(index)]=...
xtraxis(ha(index),skel,1:5,'number of wavelengths in \itR');
% Put the cap logo on
lah(index)=caplogo(ha(index),1,cloc,1/logsk,1/logsk,20);
end
clear skel
EN=1:2:5;
for index=1:3
% Get equivalent wavelengths for the double cap
for N=1:length(EN)
skel(N)=max(roots([1 1 -(EN(N)*180./(90-TH2(index)))^2]));
end
% Put in cap size labels
legsi{index+3}=sprintf(' %s = %i%s','\Theta',90-TH2(index),str2mat(176));
[bh(index+3),th(index+3)]=boxtex('ur',ha(index+3),legsi{index+3},12);
% Put on top labels
[ax(index+3),axl(index+3)]=...
xtraxis(ha(index+3),skel,EN,'number of wavelengths in \itR');
% Put the cap logo on
[lah(index+3),xal(index+3)]=caplogo(ha(index+3),3,cloc,1/logsk,1/logsk,...
20);
end
set(ax,'xgrid','on')
set(ah,'ygrid','on')
delete(axl([2 3 5 6]))
longticks(ax)
% Make the extra axis come last
for index=1:length(ah)
axes(ah(index))
set(ah(index),'color','none')
end
for index=1:length(lah)
axes(lah(index))
end
set(gcf,'color','w','inverthardcopy','off')
fig2print(gcf,'portrait')
figdisp
disp('Rerun this for the correct label dimensions!')