function CCT = SPDToCCT(SPD,S_SPD)
% CCT = SPDToCCT(SPD,S_SPD)
%
% Calculate Correlated Colour Temperature (CCT) from a spectral power
% distribution.
%
% Works by computing 1960 uv coordinates of black body radiators and then
% finding closest to the uv coordinates of the passed spectral power
% distribution.
%
% The source code contains examples.  Edit and then select and execute.

% History:
%
%                Written by Danny Garside 
% 05/04/19  dhb  Try to fix line breaks, tune up comments.
% 05/10/19  dhb  Use bsxfun rather than -, because dimension expansion does
%                not work in older versions of Matlab.

% Examples:
%{ 
    % Demonstrate proper inversion on black body radiators
    S_SPD = [380,5,81];
    theCCT = 6500;
    SPD = GenerateBlackBody(theCCT,SToWls(S_SPD));
    foundCCT = SPDToCCT(SPD,S_SPD);
    if (theCCT == foundCCT)
        fprintf('Self-inversion works\n');
    else
        fprintf('Input CCT is %d, found is %d\n',theCCT,foundCCT);
    end
%}
%{
    % Demonstrate approximate inversion on CIE daylights
    % This comes out at 6501 rather than 6500 for unknown reasons.
    load B_cieday
    theCCT = 6500;
    SPD = GenerateCIEDay(theCCT,B_cieday);
    foundCCT = SPDToCCT(SPD,S_cieday);
    if (theCCT == foundCCT)
        fprintf('Self-inversion works\n');
    else
        fprintf('Input CCT is %d, found is %d\n',theCCT,foundCCT);
    end
%}
%{
    % Calculate a bunch
    load spd_houser.mat
    CCT = SPDToCCT(spd_houser,S_houser);
%}


%% Calculate colorimetry
load T_xyz1931.mat T_xyz1931 S_xyz1931
if any(S_xyz1931 ~= S_SPD)
    T_xyz1931 = SplineCmf(S_xyz1931,T_xyz1931,S_SPD); %set to same wavelength range and sampling interval
end
SPD_XYZ = T_xyz1931*SPD;
SPD_uv = XYZTouv(SPD_XYZ,'Compute1960');

%% Calculate look ups
range = 1000:10000;
lookup_uv = zeros(2,length(range));
for i=1:length(range)
    lookup_uv(:,i) = XYZTouv(T_xyz1931*GenerateBlackBody(range(i),SToWls(S_SPD)),'Compute1960');
end

%% Find closest
for i=1:size(SPD,2)
    [~,minloc(i)] = min(sqrt(sum((bsxfun(@minus,lookup_uv,SPD_uv(:,i))).^2)));
    CCT(i) = range(minloc(i));
end

%% Check
if any(CCT == range(1) | CCT == range(end))
    warning(['You have exceeded the boundary of the temps searched for CCT. Current range is ',num2str(range(1)),'K:',num2str(range(end)),'K'])
end