Documentation of seasonal_spect

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

save seas_ct.mat ct1 ct5;
save seas_ct2.mat ct ct2;

Cross-Reference Information

Listing of script seasonal_spect

clear
lims = [180 270 -6 6];
tim = 1200:12000;
cd /home/disk/hayes2/dvimont/csiro/data/Individual_levels
filin = 'temp_M_L1_1000_years_new.nc';
nc = netcdf(filin, 'nowrite');
  lat = nc{'latitude'}(:);
  lon = nc{'longitude'}(:);
  [xk, yk] = keep_var(lims, lon, lat);
  ct = nc{'temp'}(tim, 1, yk, xk);
nc = close(nc);
ct = squeeze(ct);
ct = mean(mean(ct, 2), 3);
[ct, clim] = annave(ct);
ntim = length(ct);
for i = 1:ntim/3;
  ind = 3*(i-1) + [1:3];
  ct2(i) = mean(ct(ind));
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Thesis/Data
clear

cd /home/disk/tao/dvimont/matlab/CSIRO/Thesis/Data
load seas_ct2.mat
ntim = length(ct2);
ct1 = ct2;

ct2 = detrend(ct1);
for i = 1:3;
  if i == 2;
    ctann = ct2(1:(ntim+1)/2);
    tit = 'Years 101 - 550';
  elseif i == 3;
    ctann = ct2(((ntim+1)/2+1):ntim);
    tit = 'Years 550 - 1000';
  elseif or(i == 1, i == 4) ;
    ctann = ct2;
    tit = 'Years 101 - 1000';
  end

  nx = length(ctann);

  nfft = 128*4; 
  noverlap = 0.75*nfft;

  [p, f] = spectrum(ctann, nfft, noverlap);
  n2 = nfft/2;
  f2 = 2 * (0:n2)/n2;
  rho = (corr(ctann, ctann, 1) + sqrt(corr(ctann, ctann, 2))) / 2 ;% + ...
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn);
  rn = rn * mean(p(:,1));
  if nfft == 128*4;
    if i ~= 1;
      rner1 = rn * 2.64; rner5 = rn * 2.01;
    else
      rner1 = rn * 2.05; rner5 = rn * 1.68;
    end
  elseif nfft == 256*4;
    if i ~= 1;
      rner1 = rn * 3.55; rner5 = rn * 2.50;
    else
      rner1 = rn * 2.64; rner5 = rn * 2.01;
    end
  end
  dofx = nx / n2; 
  figure(2); figure_orient(1);
%  subplot(4,1,i)
%  pos = get(gca, 'Position');  set(gca, 'Position', [0.25 pos(2) 0.5350 pos(4)]);
  sptalk(4,1,i)
     semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ...
	      f2, rner5, 'b--', f2, rner1, 'b-.')
   if i ~= 4;
     set(gca, 'YTick', [.1 1 10 100], 'XTick', [0:.05:.5]);
     set(gca, 'YTickLabel', [0.1 1 10])
     axis([f2(3) 0.5 .01 5])
   else
     axis([f2(3) 2 .001 5])
     set(gca, 'YTick', [.01 .1 1 10 100], 'XTick', [0:.2:2]);
     set(gca, 'YTickLabel', [0.01 0.1 1 10])
   end
     grid
     ylabel(tit);
     if i == 1;
       title(['Power Spectra for CT index (seasonal data)']);
     end
   if i == 4;
     xlabel(['Frequency:  yr^-^1']);%;  Degrees of Freedom:  ' ...
   end
end
%              num2str(round(dofx))]);%' ...


cd /home/disk/tao/dvimont/matlab/CSIRO/Thesis/Chap1/Plots

cd ~/Thesis/Chap2




%  Plot the seasonal CT index

for i = 1:4
  subplot(4,1,i);
  ind = 200*(i-1)+[1:200] + 1600;
  plot(ind/4, ct1(ind))
  axis([min(ind/4)-0.25 max(ind/4) -1. 1.]);
  hold on;
    vline((min(ind/4)+[1:100]), ':k')
  hold off
%  grid on;
end


%  Plot high and lowpass filtered together

[b1, a1] = butter(9, 0.5/10);
[b2, a2] = butter(9, 2/36);

cthp = ct1 - filtfilt(b1, a1, ct1);
ctlp = filtfilt(b2, a2, ct1);

for i = 1:4
  subplot(4,1,i);
  ind = 200*(i-1)+[1:200] + 1600;
  plot(ind/4, cthp(ind), '-k', ind/4, ctlp(ind), '-r')
  axis([min(ind/4)-0.25 max(ind/4) -1. 1.]);
  hold on;
    vline((min(ind/4)+[1:100]), ':k')
  hold off
%  grid on;
end