%  In this example, we will composite surface temperature and 500mb
%  height around the cold tongue index.  This routine should be easily
%  adapted to the homework, where we composite 500mb height around the
%  cold tongue index.  Be sure to replace the appropriate contour
%  intervals, when plotting the results. 
%
%  1.  Load in data
%  2.  Define a basis
%  3.  Perform t-test on spatial data
%  4.  Plot results
%

%  Start by clearing the workspace
clear

%  1.  Load in data
load cti.djf.1950.2005.mat
load tsfc.NA.1950.2005.mat

%  Take a look at what variables we have:
whos

%  Note the variables 'lims', which are the bounds 
%  [minlon maxlon minlat maxlat] of the tsfc data.

%  2.  Define a basis - do this any way you want.  What am I doing?
basisp = find(cti > 1);
basisn = find(cti < -1);
Np = length(basisp);
Nn = length(basisn);

%  Plot results
figure(1);
orient tall;
wysiwyg

%  Set up a subplot in the figure window:
subplot(3,1,1);

%  Plot a line, and +'s and o's for positive and negative ENSO events.
pp = plot(yr, cti, 'k', ...    
          yr(basisp), cti(basisp), '+k', ...
          yr(basisn), cti(basisn), 'ok'); 

%  Now, make the plot pretty:
axis([1949 2006 -2.5 3.5])
set(gca, 'XTick', 1950:5:2005, 'fontsize', 10);
grid on
t1 = title('Dec-Feb Cold Tongue Index (ENSO)', 'fontsize', 14);
yl = ylabel('CTI', 'fontsize', 12);
xl = xlabel('Year (Corresponds to January)', 'fontsize', 12);

%  3.  Perform t-test on spatial data
%  3.1   Identify points with 'equal' variance
%  3.2   Calculate t-statistic with pooled std. dev.
%  3.3   Calculate t-statistic with different std. dev.
%  3.4   Identify regions outside of the 95% confidence levels
%  3.5   Plot results

%  We'll use the difference in means quite a bit, so let's go ahead and
%  calculate some stuff...
tsfcdiff = squeeze(mean(tsfc(basisp,:,:))) - ...
    squeeze(mean(tsfc(basisn,:,:)));
stdtsfcp = squeeze(std(tsfc(basisp,:,:)));
stdtsfcn = squeeze(std(tsfc(basisn,:,:)));

%  3.1:  Identify points with equal variance
%  1.  90% significance level
%  2.  H0:  variance is the same
%      HA:  Variance is different
%  3.  Statistic:  2-tailed F-test
%  4.  Critical region:  finv(0.05,Np-1,Nn-1) < F < finv(0.95,Np-1,Nn-1)
%  5.  Evaluate the statistic (below)

F = stdtsfcp.^2 ./ stdtsfcn.^2;
flow = finv(0.05, Np-1, Nn-1);
fhi = finv(0.95, Np-1, Nn-1);

samevar = find((flow < F) & (F < fhi));
diffvar = find((flow > F) | (F > fhi));

%  3.2  Calculate t-test with pooled std. dev.

%  Set up t-test with equal variance.
%  1.  95% significance level
%  2.  H0:  SST during El Nino events is different than La Nina events
%      HA:  SST is no different during the two events
%  3.  Statistic:  2-tailed t-test with N1+N2-2 DOF's
%  4.  Critical region:  tlow < T < thi
%  5.  Evaluate the statistic ...

%  Start by setting up an empty array for the T variable:
T = repmat(NaN, size(tsfcdiff));

%  Now, calculate t-value for equal variance
spooled = sqrt(...
    ((Np-1)*stdtsfcp.^2 + (Nn-1)*stdtsfcn.^2) ...
    / (Np + Nn - 2) );
T(samevar) = tsfcdiff(samevar) ./ (spooled(samevar) * sqrt(1/Np + 1/Nn));

%  3.3  Calculate t-statistic with different std. dev.
T(diffvar) = tsfcdiff(diffvar) ./ ...
    sqrt(stdtsfcp(diffvar).^2/Np + stdtsfcn(diffvar).^2/Nn);

%  3.4  Identify regions outside the 95% confidence level

%  Start by identifying DOF's from unequal variance, then replace points
%  with equal variance with DOF = Np+Nn-2
DOF = repmat(NaN, size(tsfcdiff));
DOF = ((stdtsfcp.^2/Np) + (stdtsfcn.^2/Nn)).^2 ./ ...
      (((stdtsfcp.^2/Np).^2/(Np-1)) + ((stdtsfcn.^2/Nn).^2/(Nn-1)));
DOF(samevar) = Np + Nn - 2;
DOF = round(DOF);

%  Set up array with 0's where H0 is accepted, and 1's where it is
%  rejected.
H = ones(size(tsfcdiff));
H((tinv(0.025, DOF) < T) & (T < tinv(0.975, DOF))) = 0;

%  3.5  Plot results
subplot(3,1,2); cla;

%  Set up a map axis - this requires the mapping toolbox, which should be
%  on the computers in the 14th floor lounge.
axesm('eqdcylin', ...
      'maplonlimit', lims(1:2), 'maplatlimit', lims(3:4));

%  Produce a filled contour plot, with contours / shading every 0.5 deg C.
hh = contourfm(lat, lon, tsfcdiff, [-5:.5:5]*30);

%  Now, center the color map on 0, and make it saturate at -2.5 and 2.5
caxis([-2.5 2.5]*30);

%  Add a pretty land map - this is a classic example of MATLAB being
%  ridiculously more complicated than is needed.  Don't worry about the
%  details, just plug in the following lines...
lm = worldlo('POline');
for i = 1:length(lm);
  lm(i).otherproperty = {'color', 'k'};
end
displaym(lm);

%  Want states?  Try:  "help usalo".  Have fun.

%  Make it pretty
gridm on;
tightmap;

%  Now, add contours where the above null hypothesis is rejected.
hold on;
[cc2, hh2] = contourm(lat, lon, H, [-0.5 0.5], 'linewidth', 2, 'color', 'k');
hold off;

%  Add a title / labels, etc.
t1 = title('ENSO Composite:  Surface Temperature', 'fontsize', 14);
xl = xlabel('Basis:  CTI \pm 1^oC', 'fontsize', 12);
yl = ylabel('Contour:  0.5^oC', 'fontsize', 12);