Global Index (short | long) | Local contents | Local Index (short | long)
tfield1 = tfield1_10000;
| This script calls | |
|---|---|
clean
cd ~/model/enso
load tfield1;
yr = tfield1(:,1)./360;
ct = tfield1(:,2);
ct = annave(ct);
ctm = ct;
yrm = yr;
cd ~/model/enso
save ct_monthly.mat ctm yrm
ct = ctm;
if 1;
ntim = length(ctm);
ct = zeros((ntim/3), 1);
for i = 1:(ntim/3);
ind = 3*(i-1)+[1:3];
ct(i) = mean(ctm(ind));
end
yr2 = [1/4:1/4:(100000)]';
end
cd ~/model/enso
save tfield_seas.mat ct yr2
% Calculate spectrum
clean
cd ~/model/enso
load tfield_seas.mat
ctm = ct;
nfft = 128;
nx = length(ctm); noverlap = 0.875*nfft; n2 = nfft/2;
dofx = round(nx/n2);
lev5 = finv(0.95, dofx, 1200000);
lev1 = finv(0.99, dofx, 1200000);
[p,f] = spectrum(ctm, nfft, noverlap, hanning(nfft), 1, 'none');
f2 = f;
cd ~/model/enso
load tfield_ann.mat
% AR1 process
n2b = 128/2;
rho = corr(ct, ct, 1);
freq = (pi/(n2b))*[0:(n2b)];
rn = (1-rho^2) ./ (1-2*rho*cos(freq)+rho^2);
rn = rn / mean(rn);
rn = rn * mean(p(:,1));
maxrn = find((f.*rn')==max(f.*rn'));
% AR2 process
n2b = 128/2;
rho1 = corr(ct, ct, 1);
rho2 = corr(ct, ct, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
% R = sqrt(-1*phi2);
% f = acos(phi1/(2*sqrt(-1*phi2)))/(2*pi);
freq = (pi/(n2b))*(0:n2b);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq)-...
2*phi2*cos(2*freq)+2*phi1*phi2*cos(freq));
rn2 = rn2 / mean(rn2);
rn2 = rn2 * mean(p(1:65,1));
maxrn2 = find((f.*rn2')==max(f.*rn2'));
plot((f2), (p(:,1)), 'k', (freq/(2*pi)), (rn2), 'k')
plot(log(f2), log(p(:,1)), 'k', log(freq/(2*pi)), log(rn2), 'k', ...
log(f2), log(lev5*rn2), '--k')
% Plot the data
figure(1); fo(1); clf;
subplot(4,1,1);
plot(yr2(1:100000), ct, 'k');
axis([-1 100001 -3.1 3.1]);
set(gca, 'fontsize', 8);
t1 = title('100Ky Annually resolved CT time series, Battisti model');
x1 = xlabel('time (Kyr)');
set(t1, 'fontsize', 10);
ylabel('^\circ C');
set(gca, 'YTick', [-5:1:5], 'XTick', 0:1e4:1e5, ...
'XTickLabel', [0:10:100]);
subplot(4,2,3);
x = [-5:.1:5];
N = hist(ct,x);
h = bar(x, N);
axis([-2.25 2.25 -0.2e3 7.2e3])
set(gca, 'XTick', [-5:1:5], 'YTick', [0:1e3:15e3], ...
'YTickLabel', [0:1:15], 'fontsize', 8);
t1 = title('CT Histogram (Annual)')
set(t1, 'fontsize', 9);
x1 = xlabel('^\circ C');
y1 = ylabel('Count (in thousands)')
tem = standardize(ct);
moms = [moment(ct, 2) moment(tem, 3) moment(tem, 4)];
subplot(4,2,4); cla;
set(gca, 'Visible', 'off');
axis([0 1 0 1]);
text1 = ['Mean '; 'Variance';'Skewness'; 'Kurtosis'];
tx1 = text(0*ones(4,1),[0.7:-.2:0.1],text1);
tx1b = text(0.35, 0.9, 'CT Index');
tx2b = text(0.35, 0.7, ['0.0 ^\circC']);
tx3b = text(0.35, 0.5, [num2str(moms(1)) ' ^\circC^2']);
tx4b = text(0.35, 0.3, [num2str(moms(2))]);
tx5b = text(0.35, 0.1, [num2str(moms(3))]);
tx1c = text(0.8, 0.9, 'Normal');
tx2c = text(0.8, 0.7, '-');
tx3c = text(0.8, 0.5, '-');
tx4c = text(0.8, 0.3, '0');
tx5c = text(0.8, 0.1, '3');
%figure(3); fo(1); clf;
% Plot spectra
subplot(4,2,6);
pos = get(gca, 'Position'); pos(1) = pos(1)+0.3*pos(3);
pos(3) = 0.7*pos(3); set(gca, 'Position', pos);
set(gca, 'Visible', 'off');
axis([0 0.5 0 1]);
text(0.1,0.65,'AR1');
text(0.1,0.45,['Efold: ' ...
num2str(round(100./(2*pi*f2(maxrn)))/100) ...
' yrs']);
subplot(4,2,5); cla;
pos = get(gca, 'Position');
pos(3) = 1.5*pos(3); set(gca, 'Position', pos);
hh = plot(log(f2), f2.*p(:,1), 'k', ...
log(f2), f2'.*rn, '-k', ...
log(f2), lev5*f2'.*rn, '--k', ...
log(f2), lev1*f2'.*rn, '-.k');
set(hh(1), 'linewidth', 2);
% hh = plot(log(f2), log(p(:,1)), 'k', ...
% log(f2), log(rn), '--k');
axis([-log(128)-0.1 log(0.5)+0.1 -0.01 .16]);
set(gca, 'XTick', -log(2.^[10:-1:-6]), ...
'XTickLabel', 2.^[10:-1:-6], 'YTick', [0:.025:3]);
grid on
set(gca, 'fontsize', 8);
% Plot spectra
subplot(4,2,8);
pos = get(gca, 'Position'); pos(1) = pos(1)+0.3*pos(3);
pos(3) = 0.7*pos(3); set(gca, 'Position', pos);
set(gca, 'Visible', 'off');
axis([0 0.5 0 1]);
text(0.1,0.65,'AR2');
text(0.1,0.45,['Period: ' ...
num2str(round(100./(f2(maxrn2)))/100) ...
' yrs']);
subplot(4,2,7); cla;
pos = get(gca, 'Position');
pos(3) = 1.5*pos(3); set(gca, 'Position', pos);
hh = plot(log(f2), f2.*p(:,1), 'k', ...
log(f2), f2'.*rn2, '-k', ...
log(f2), lev5*f2'.*rn2, '--k', ...
log(f2), lev1*f2'.*rn2, '-.k');
set(hh(1), 'linewidth', 2);
axis([-log(128)-0.1 log(0.5)+0.1 -0.01 .16]);
set(gca, 'XTick', -log(2.^[10:-1:-6]), ...
'XTickLabel', 2.^[10:-1:-6], 'YTick', [0:.025:3]);
grid on
set(gca, 'fontsize', 8);
cd ~/matlab/CSIRO/Longrun/Figs
print -dps2 bat_annual_spectrum.ps