function [corr,nstd]=taylor(data,models);
% A utility for making a Taylor diagram.
% Code by Laurence Anderson, Woods Hole Oceanographic Institution, 2 July 2008.
% Reference:
%   Taylor, K.E., 2001. Summarizing multiple aspects of model performance in
%   a single diagram. Journal of Geophysical Research 106(D7), p 7183-7192.
%
% Input:
%   data(1:n,1) = vector containing "n" data values.
%                 Note the data must be all in the same units.
%   models(1:n,1:m) = array containing "n" values from "m" models.
%                 Note these "n" model values are assumed to correspond
%                 in time and space with the data values.
% Output:
%   corr(1,1:m) = model-data correlation coefficients
%   nstd(1,1:m)  = data-normalized standard deviations
%   Taylor diagram. "1" will be plotted for the first model, 
%                   "2" for the second model, etc.

% common user-tunable parameters

tsize=15; % fontsize for model identifiers
fsize=15; % fontsize for axis labels
titl='Taylor diagram'; % title

% check if the input is in the correct format

[sd1,sd2]=size(data);
[sm1,sm2]=size(models);
if sd2>1,
   error(['plot_taylor.m: Error: data array has wrong dimensions.'])
elseif sd1==1,
   error(['plot_taylor.m: Error: Only one data point.'...
          ' Correlation coefficient is ill-defined.'])
elseif sm1~=sd1,
   error(['plot_taylor.m: Error: models array has wrong dimensions.'])
end;
n=sd1; % number of data points
m=sm2; % number of models

% compute statistics

dat_mean=mean(data); % data mean
mod_mean=mean(models); % model means
dat_std=std(data,1); % data standard deviation
                     % Note to get the correct value of the correlation
                     % coefficient, the standard deviation must use division
                     % by "n" not "n-1".
if dat_std==0,
   error(['plot_taylor.m: Error: std(data)=0.' ...
          ' Correlation coefficient is ill-defined.'])
end;
for j=1:m;
    mod_std(1,j)=std(models(:,j),1); % model standard deviations
    if mod_std(1,j)==0,
       % In this case, the correlation coefficient = 0/0, ill-defined.
       % Yet nstd = 0, making the placement on the Taylor diagram
       % unambigously at the origin.
       % So for plotting temporarily set corr = 0.
       % This permits plotting of "persistence" models.
       corr(1,j)=0; % correlation coefficient
    else,
      % correlation coefficient
      corr(1,j)=mean((models(:,j)-mod_mean(1,j)).*(data-dat_mean))/...
                (mod_std(1,j)*dat_std);
    end;
end;
% data-normalized standard deviation 
nstd=mod_std/dat_std;
   
% make Taylor diagram
%   Note cartesian coordinates are used for easier control of axes, etc.

figure
angle=acos(corr); % radians
x=nstd.*cos(angle); % x location in Cartesian coordinates
y=nstd.*sin(angle); % y location in Cartesian coordinates
for j=1:m;
  text(x(j),y(j),num2str(round(j)),'horizontalalignment','center',...
      'verticalalignment','middle','fontsize',tsize); hold on
end;
% set axes
  axis square
  rmax=max([1 1.05*max(abs(x)) 1.05*max(y)]);
  if min(x)>=0, % plot only first quadrant, if possible
    axis([0 rmax 0 rmax])
  else
    set(gca,'PlotBoxAspectRatio',[1 .5 1])
    axis([-rmax rmax 0 rmax])
  end;
% add grid lines
  ang=[0:1:180]*pi/180;
  xx=cos(ang); yy=sin(ang);
  plot(xx,yy,'k-') % plot unit circle
  plot(1,0,'k.','markersize',15) % plot optimal point
  rtick1=rmax/4; % first-guess grid line increment
  % pick the nearest increment of 0.2 0.5 1 2 5 10 20 50
  rticks=[0.2 0.5 1 2 5 10 20 50];
  [yr,ir]=sort(abs(rtick1-rticks));
  rtick=rticks(ir(1)); % grid line increment
  nticks=floor(sqrt(5)*rmax/rtick); % number of grid lines
  for i=1:nticks;
      plot(i*rtick*xx,i*rtick*yy,'k:')
  end;
  rsqrt3=rmax/sqrt(3);
  plot([0 rmax],[0 rsqrt3],'k:')  % 30 degrees
  plot([0 rsqrt3],[0 rmax],'k:')  % 60 degrees
  plot([0 0],[0 rmax],'k:')       % 90 degrees
  plot([0 -rsqrt3],[0 rmax],'k:') % 120 degrees
  plot([0 -rmax],[0 rsqrt3],'k:') % 150 degrees
  set(gca,'xtick',[(-nticks*rtick):rtick:(nticks*rtick)],...
          'ytick',[0:rtick:(nticks*rtick)])
set(gca,'box','on','tickdir','out')
set(gca,'fontsize',fsize)
xlabel('(\sigma_m/\sigma_d) cos(acos(r))','fontsize',fsize)
ylabel('(\sigma_m/\sigma_d) sin(acos(r))','fontsize',fsize)
title(titl,'fontsize',fsize)

% for output, reset corr to 0/0 = NaN if mod_std=0
for j=1:m; if mod_std(1,j)==0, corr(1,j)=NaN; end; end;

return
