function R = ser(s,varargin);
% SER - An Algorithm for Computing Spherical Error Radius
%
% USAGE:
%
% R = ser(s)
% R = ser(s,prob)
% R = ser(s,prob,flag)
%
% where s = [sx], [sx,sy], or [sx,sy,sz] is a 1, 2, or 3 dimensional
% vector of standard deviations; 0 < prob < 1 is the probability; and
% flag ~= 0 invokes a diagnostic subroutine.
%
% This function computes the Spherical Error Radius (SER) from inputs
% consisting of the square roots of the eigenvalues of a covariance matrix
% (equivalently, from sigma-x, sigma-y, and sigma-z of a trivariate normal
% distribution in a coordinate system where there is no cross-correlation
% between variables.) This means that if you have a covariance matrix and
% wish to compute the SER, simply obtain the square roots of the eigenvalues
% and use these as inputs. For example, list them via "sqrt(eig(C))" where C
% is your covariance matrix.
%
% The SER is the radius of a sphere which contains a fraction of probability
% equal to the input "prob," which is asumed to be 0.95 if omitted.
%
% REFERENCES:
%
%   Kleder, Michael (2004). "An algorithm for converting covariance to
%     spherical error probable." Mathworks Central File Exchange.
%     (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
%     objectId=5688&objectType=file)
%   Koopman, Ray (2004). "Joint normal distribution integral over the unit
%     sphere." The Math Forum.
%     (http://mathforum.org/kb/message.jspa?messageID=1552353&tstart=0)

%   Copyright(c)2006 Version 1.00
%     Tom Davis   (tdavis@eng.usf.edu)
%     Mike Kleder (mkleder@hotmail.com)

msg=['SER requires a vector of standard deviations ',...
  'with length less than or equal to three.'];
if nargin<1
  error(msg)
end
if length(s)>3 | length(s)<1
  error(msg)
end
if nargin>1, prob=varargin{1}; else, prob=0.95; end
if nargin>2, flag=varargin{2}; else, flag=0; end
options=optimset(@fzero);
lowend=eps; highend=5*max(s);

warning('off','MATLAB:quad:ImproperFcnValue')
try
  R=fzero(@enclosederr,[lowend,highend],options,s,prob);
catch
  R=abs(fzero(@enclosederr,lowend,options,s,prob));
end
warning('on','MATLAB:quad:ImproperFcnValue')

if flag
  diagnostic(R,s) % set flag to verify results against RANDN
end

%--------------------------------------------------------------------
function enclerr=enclosederr(R,s,prob)
n=length(s);
if n==1
  enclerr=erf(R/(sqrt(2)*s));
else
  enclerr=quad(@dens,0,R,[],[],R,s);
end
enclerr=prob-enclerr;

%--------------------------------------------------------------------
function c=dens(r,R,s)
n=length(s);
a= 0.25*(1/s(1)^2-1/s(2)^2);
b=-0.25*(1/s(1)^2+1/s(2)^2);
c=(r.*exp(b*r.^2).*besseli(0,a*r.^2))/(s(1)*s(2));
if n==3
  c=c.*erf(sqrt(0.5*(R^2-r.^2))/s(3));
end

%--------------------------------------------------------------------
function diagnostic(R,s)
disp('Running test of SER.')
n=length(s);
for i=10:-1:1
  fprintf('%g ',i);
  x=randn(1e6,1)*s(1);
  if n>1, y=randn(1e6,1)*s(2); else y=zeros(1e6,1); end
  if n>2, z=randn(1e6,1)*s(3); else z=zeros(1e6,1); end
  r=sqrt(sum([x y z].^2,2));
  m(i)=sum(r<R)/length(r);
end
disp(' ');
disp([num2str(mean(m)*100),...
  '% of random points were located within the computed radius.'])