function [X,Y,Z]=spec2xyz(lambda,spectrum);
%% Computes X, Y, Z color coordinates from spectral data
%  using a trapeziodal integration of interpolated xbar,ybar,zbar data
%
%  usage:
%  [X,Y,Z]=spec2xyz(lambda,spectrum);
%
%  lambda is a vector of wavelengths in nanometers
%  spectrum is a data cube in any units. The first two dimensions
%  are spatial and the third is wavelength.  The wavelengths are
%  given by lambda from 1 to length of lambda in the third
%  dimension of spectrum.
%  X,Y,Z are the color coordinates in the (units of
%  spectrum)*nanometers
%  They are 2-D arrays with size given by the first two dimensions
%  of spectrum.
%
%  The user can, if desired, compute x,y using x=X/(X+Y+Z) and
%  y=Y/(X+Y+Z) 
%
load('xyzdata');
test=find(lambda>=1000*a(1,1)&lambda<=1000*a(end,1));
xbar=zeros(length(lambda),1,1);ybar=xbar;zbar=xbar;temp=xbar;
xbar(ind2sub(size(xbar),test))=linterp(1000*a(:,1),a(:,2),lambda(test));
ybar(ind2sub(size(xbar),test))=linterp(1000*b(:,1),b(:,2),lambda(test));
zbar(ind2sub(size(xbar),test))=linterp(1000*c(:,1),c(:,2),lambda(test));
temp(ind2sub(size(xbar),test))=lambda(test);
% permute xbar, ybar, and zbar to 1x1x(size(lambda))
xbarper = permute(xbar, [2,3,1]);
ybarper  = permute(ybar, [2,3,1]);
zbarper  = permute(zbar, [2,3,1]);
lambdaper= permute(temp, [2,3,1]);
% replicate xbarper, ybarper, and zbarper to be size(r)
xbar3d = repmat(xbarper, [size(spectrum,1),size(spectrum,2),1]);
ybar3d = repmat(ybarper, [size(spectrum,1),size(spectrum,2),1]);
zbar3d = repmat(zbarper, [size(spectrum,1),size(spectrum,2),1]);
lambda3d = repmat(lambdaper, [size(spectrum,1),size(spectrum,2),1]);

xintegrand=xbar3d.*spectrum;
yintegrand=ybar3d.*spectrum;
zintegrand=zbar3d.*spectrum;

% X=trapz(lambda,permute(xintegrand, [3,1,2]));
% Y=trapz(lambda,permute(yintegrand, [3,1,2]));
% Z=trapz(lambda,permute(zintegrand, [3,1,2]));

X=(sum(xintegrand(:,:,2:end).*lambda3d(:,:,1:end-1)+...
       xintegrand(:,:,1:end-1).*lambda3d(:,:,2:end),3)...
   +xintegrand(:,:,end).*lambda3d(:,:,end)...
   -xintegrand(:,:,1).*lambda3d(:,:,1))/2;


Y=(sum(yintegrand(:,:,2:end).*lambda3d(:,:,1:end-1)+...
       yintegrand(:,:,1:end-1).*lambda3d(:,:,2:end),3)...
   +yintegrand(:,:,end).*lambda3d(:,:,end)...
   -yintegrand(:,:,1).*lambda3d(:,:,1))/2;


Z=(sum(zintegrand(:,:,2:end).*lambda3d(:,:,1:end-1)+...
       zintegrand(:,:,1:end-1).*lambda3d(:,:,2:end),3)...
   +zintegrand(:,:,end).*lambda3d(:,:,end)...
   -zintegrand(:,:,1).*lambda3d(:,:,1))/2;


