function [bias,rmsd,urmsd,dat_std]=target(data,models);
% A utility for making a Target diagram.
% Code by Laurence Anderson, Woods Hole Oceanographic Institution, 9 April 2009.
% Reference:
%   Jolliff, J.K., J.C. Kindle, I. Shulman, B. Penta, M.A.M. Friedrichs, R.
%   Helber and R.A. Arnone, 2009. Summary diagrams for coupled hydrodynamic-
%   ecosystem model skill assessment. Journal of Marine Systems 76, p 64-82.
%
% 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:
%   bias(1,1:m) = bias, normalized by std(data) if std(data) is not zero.
%   rmsd(1,1:m) = Root Mean Square Difference, normalized by std(data) if 
%                 std(data) is not zero.
%  urmsd(1,1:m) = unbiased RMSD, normalized by std(data) if std(data) is
%                 not zero.
%   dat_std(1)  = standard deviation of the data.
%   Target diagram. "1" will be plotted to identify 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='Target diagram'; % title

% check if the input is in the correct format

[sd1,sd2]=size(data);
[sm1,sm2]=size(models);
if sd2>1,
   error(['target.m: Error: data array has wrong dimensions.'])
elseif sm1~=sd1,
   error(['target.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 (actually first moment).

% compute mod_std, bias, rmsd, urmsd
for j=1:m;
    mod_std(1,j)=std(models(:,j),1); % model standard deviations
    bias(1,j)=mod_mean(1,j)-dat_mean;
    xxx=models(:,j)-data;
    rmsd(1,j)=sqrt(mean(xxx.*xxx));
    xxx=(models(:,j)-mod_mean(1,j))-(data-dat_mean);
    urmsd(1,j)=sqrt(mean(xxx.*xxx));
    % apply sign as in Eq. 11 in Jolliff et al. (2009)
    if mod_std(1,j) < dat_std, 
       urmsd(1,j) = -urmsd(1,j);
    end;
end;

% normalize bias, rmsd and urmsd by std(data)
wflag=0; 
if dat_std==0,
   warning(['target.m: Warning: std(data)=0. ' ...
   'Not normalizing units by std(data).'])
   wflag=1;
else,
   bias=bias/dat_std;
   rmsd=rmsd/dat_std;
   urmsd=urmsd/dat_std;
end;
  
% make Target diagram
%   Note cartesian coordinates are used for easier control of axes, etc.

figure
x=urmsd; % x location in Cartesian coordinates
y=bias;  % 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=1.05*max([max(abs(x)) max(abs(y))]);
  if wflag==0, rmax=max(rmax,1); end;
  axis([-rmax rmax -rmax rmax])
  plot(0,0,'k.','markersize',15) % plot optimal point
% add grid lines
  ang=[0:1:360]*pi/180;
  xx=cos(ang); yy=sin(ang);
  if wflag==0,
     plot(xx,yy,'k-') % plot unit circle
  end;
  plot([0 0],[-rmax rmax],'k-')
  plot([-rmax rmax],[0 0],'k-')
  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;
  set(gca,'xtick',[(-nticks*rtick):rtick:(nticks*rtick)],...
          'ytick',[(-nticks*rtick):rtick:(nticks*rtick)])
set(gca,'box','on','tickdir','out')
set(gca,'fontsize',fsize)
if wflag==0,
  xlabel('unbiased RMSD/\sigma_d','fontsize',fsize)
  ylabel('Bias/\sigma_d','fontsize',fsize)
else,
  xlabel('unbiased RMSD','fontsize',fsize)
  ylabel('Bias','fontsize',fsize)
end;
title(titl,'fontsize',fsize)

return
