% PS09 answer. This answer m-file is based on the answers provided by Emily
%   Van Ark and Rachel Stanley.
%
% Started 2008:11:17 D. Glover, WHOI assembled from parts of student's answers
% Modif'd 2008:12:02 DMG for further embellishments
% Modif'd 2008:12:06 DMG added 64-color anomaly color map
% Modif'd 2008:12:22 DMG to move the text into a PDF file and to place a,
%                        b, c, or d on each subplot

load bdalev.mat
load purp_map4
load anom64cmap

Adiff=zeros(4,5); % initialize matrix for the average differences
Asqdev=zeros(4,5); % initialize matrix for the squared deviation
V=[-3:.5:3]';   % vector to use to draw contours and contour labels
Kv=[0.5e-4, 0.75e-4, 1.0e-4, 1.25e-4, 1.5e-4];
wf=[0 0.5 1 2];

% We will completely use a datafile and then go onto next one so don't have to rename
% every matrix... Also we will use a for loop to load in datafiles, plot up differences,  
% and calculate average difference and average square deviation 

for j=1:5       % do loop with various kz first for ease in plotting -- i want to make one figure for each k
   figure       % start a new figure for each Kz value
   for i=1:4
      eval(sprintf('load d%d%d',i,j)); % eg load in data files, for example d11
      subplot(4,1,i);
      imagesc(tf,z,(Ta-bdat),[-2 +2]);colorbar;colormap(cmap);hold;
      [C,Hc]=contour(tf,z,Ta-bdat,[-2:0.5:2],'k');
      set(Hc,'showtext','on');
      if i==1
          title([sprintf('Run d%d%d',i,j) ' Climatic']);
          text(40,-60,'{\Delta}T ({\circ}C)');
          text(0,-75,['K_{v}=' sprintf('%1.2e',Kv(j))]);
          text(2.5,250,'{\bf a}','fontsize',12);
      elseif i==2
          title([sprintf('Run d%d%d',i,j) ' Half Stochastic']);
          text(2.5,250,'{\bf b}','fontsize',12);
      elseif i==3
          title([sprintf('Run d%d%d',i,j) ' Full Stochastic']);
          text(2.5,250,'{\bf c}','fontsize',12);
      else
          title([sprintf('Run d%d%d',i,j) ' Double Stochastic']);
          xlabel('Month');
          text(2.5,250,'{\bf d}','fontsize',12);
      end
      ylabel('depth (m)');
      hold;
      Adiff(i,j)=mean(mean(Ta-bdat));   % calculate and save the average difference   
      Asqdev(i,j)= mean(mean((Ta-bdat).^2));    % calc and save the average square deviation
   end
end

disp('Hit any key to continue')
disp(' ')
pause

figure;
subplot(121);
clabel(contourf(Kv,wf,Adiff));axis square;colormap(purp_map4)
%imagesc(Kv,wf,Adiff);colorbar;colormap(purp_map4);hold;axis square
%[C,Hc]=contour(Kv,wf,Adiff,[-2:0.5:2],'k');
%set(Hc,'showtext','on');
title('Average Difference'); xlabel('K_{v}'); ylabel('Stochastic Forcing');
subplot(122);
clabel(contourf(Kv,wf,Asqdev));axis square;colormap(purp_map4)
%imagesc(Kv,wf,Asqdev);colorbar;colormap(purp_map4);hold;axis square
%[C,Hc]=contour(Kv,wf,Asqdev,[-2:0.5:2],'k');
%set(Hc,'showtext','on');
title('Average Square Deviation'); xlabel('K_{v}'); ylabel('Stochastic Forcing');