%  This file goes through the steps involved in generating the maps in
%  question 2 of HW2.  This exercise involves libraries in
%  ~dvimont/matlab.  If you do not already have a matlab directory in
%  your home directory, run the following commands in an open terminal
%  window:
%
%  cp /f5/raid5/class2/fall03/dvimont/matlab.tar .
%  tar -xf matlab.tar
%
%  If you already have a matlab directory there, and have your own
%  startup files and such, then you may want to be more selective about
%  which files you copy.  In either case, you should copy the files
%  before starting MATLAB.
%

%  Start by clearing the workspace.
clear

%  Next, load data for the monthly-average regression maps

%load /Users/dvimont/Class/aos575/Homework/Data/HW2_prob2_monthly.mat
load /f5/raid5/class2/fall03/dvimont/Class/Data/HW2_prob2_monthly.mat

%  Take a look at what variables are in the workspace by typing 'who', or
%  'whos': 

who    %  Short list of variables in the workspace

whos   %  Long list of variables in the workspace

%  Note that the variable 'description' is in the workspace.  This is a
%  string variable that contains a description of all the other
%  variables.  To see the variable 'description', type 'description' on
%  the command line:

description

%  Let's start by generating regression maps of SST onto the CTI.  Start
%  by standardizing the CTI (place a semi-colon on the end of the line to
%  avoid seeing the results of the operation):

cti = (cti_mon - mean(cti_mon)) / std(cti_mon);

%  Now, let's get the regression map of of SST onto the CTI.  First, we
%  need to get the SST data into matrix format.  Currently, it's a
%  3-dimensional array, where the first, second and third dimensions are
%  time, latitude, and longitude, respecitvely (in general, NetCDF
%  climate data have dimension time, level, latitude, longitude):

[ntim, nlat, nlon] = size(skt_mon);

%  Now, reshape into a two-dimensional matrix (time X space):

skt_mon = reshape(skt_mon, ntim, nlat*nlon);
size(skt_mon)  %  Just to double-check

%  I've already removed the monthly climatology (hence the mean) of the
%  data, so we can go ahead and get the regression maps as the covariance
%  of the cti with each column of skt_mon.  Note that the variance of
%  'cti' is one, so we don't need to divide by this in the regression
%  equation:  a1_map = cti' * skt_mon / (cti' * cti);

a1_map = cti' * skt_mon / (ntim - 1);

%  The correlation coefficient involves another step:  First, we need to
%  divide each column of skt_mon by its standard deviation, then get the
%  covariance.  To divide each element in skt_mon, use the element
%  division (which has a '.' in front of the '/':  ./):

corr_map = cti' * (skt_mon ./ (ones(ntim,1) * std(skt_mon))) / (ntim - 1);

%  Finally, reshape a1_map and corr_map so that it has dimensions (lat X
%  lon):

a1_map = reshape(a1_map, nlat, nlon);
corr_map = reshape(corr_map, nlat, nlon);

%  Now, we're ready to generate a contour plot of this.  First, open a
%  figure window - I'm going to use the 'tall' orientation, and place two
%  maps in each window (using the subplot command):

figure_tall(1);
subplot(2,1,1); cla;

%  We'll use the routine 'make_contour_plot' to plot the map.  This is
%  just a whole slew of commands I've bunched into one file.  For
%  information on how to use the routine, type:

help make_contour_plot

%  Now, contour the map, following directions as you go:

make_contour_plot(lat_skt, lon_skt, a1_map);

%  Do the same for the correlation map, in a new subplot:

subplot(2,1,2); cla;
make_contour_plot(lat_skt, lon_skt, corr_map);

%  Finally, print this

%  Now, let's do the same with 500mb height
%  Reshape hgt_mon
[ntim, nlat, nlon] = size(hgt_mon);
hgt_mon = reshape(hgt_mon, ntim, nlat*nlon);

%  Get regression and correlation maps
a1_map = cti' * hgt_mon / (ntim - 1);
corr_map = cti' * (hgt_mon ./ (ones(ntim,1) * std(hgt_mon))) ...
             ./ (ntim - 1);

%  And plot these - let's plot the correlation and regression on the same
%  map.

figure_tall(2); clf;  %  Clears the current figure window

subplot(2,1,1); cla;
  a1_map = reshape(a1_map, nlat, nlon);
  make_contour_plot(lat_hgt, lon_hgt, a1_map);
  hold on;
    corr_map = reshape(corr_map, nlat, nlon);
    make_surface_plot(lat_hgt, lon_hgt, corr_map);
  hold off;

%  Add a colorbar:
cb = colorbar3('bottom');

%  Now, plot the same plot, except using annually resolved data (the winter
%  average data):

%load /Users/dvimont/Class/aos575/Homework/Data/HW2_prob2_winave.mat
load /f5/raid5/class2/fall03/dvimont/Class/Data/HW2_prob2_winave.mat

%  Standardize the cti for winter this time:
cti = (cti_win - mean(cti_win)) / std(cti_win);

%  Reshape into a data matrix:
[nyr, nlat, nlon] = size(hgt_win);
hgt_win = reshape(hgt_win, nyr, nlat*nlon);

%  And get the regression maps:
a1_map = cti' * hgt_win / (nyr - 1);
corr_map = cti' * (hgt_win ./ (ones(nyr,1) * std(hgt_win))) / (nyr -1);

subplot(2,1,2); cla;
  a1_map = reshape(a1_map, nlat, nlon);
  make_contour_plot(lat_hgt, lon_hgt, a1_map);
  hold on;
    corr_map = reshape(corr_map, nlat, nlon);
    make_surface_plot(lat_hgt, lon_hgt, corr_map);
  hold off;