function RGB = XYZ_to_RGB(XYZ, Primary)
% function RGB = Color_Convert_XYZ_to_RGB(XYZ, Primary)
%
% Converts from Judd-Vos XYZ to Monitor RGB (for monitor with Primaries defined in Judd-Vos XYZ)
%
% Takes: XYZ (N, 3) - XYZ coordinates
%        Primary (3, 3) - XYZ coordinates of red, green and blue guns / filters of monitor
% Returns: RGB (N, 3) - RGB coordinates (scaled from 0 to 1)
% Dependencies: Variable_Report.m
%
% Created 2015-06-26 by KCM
%
% Updated 2015-07-10 by KCM

%% Arguments
if ~exist('XYZ', 'var') || isempty(XYZ)
    fprintf(char(strcat({'\nColor_Convert_XYZ_to_RGB: '}, {'''XYZ'' must be provided!\n\n'})))
    RGB = []; return
end
if ~ismatrix(XYZ) || size(XYZ, 2) ~= 3
    fprintf(char(strcat({'\nColor_Convert_XYZ_to_RGB: '}, {'''XYZ'' must have size (N, 3)!\n\n'})))
    RGB = []; return
end
if sum(sum(~isfinite(XYZ))) || sum(sum(~isreal(XYZ)))
    fprintf(char(strcat({'\nColor_Convert_XYZ_to_RGB: '}, ...
        {'''XYZ'' must be finite and real!\n\n'})))
    RGB = []; return
end
if ~exist('Primary', 'var') || isempty(Primary)
    fprintf(char(strcat({'\nColor_Convert_XYZ_to_RGB: '}, {'''Primary'' must be provided!\n\n'})))
    RGB = []; return
end
if ~ismatrix(Primary) || size(Primary, 1) ~= 3 || size(Primary, 2) ~= 3
    fprintf(char(strcat({'\nColor_Convert_XYZ_to_RGB: '}, ...
        {'''Primary'' must have size (3, 3)!\n\n'})))
    RGB = []; return
end
if sum(sum(~isfinite(Primary))) || sum(sum(~isreal(Primary)))
    fprintf(char(strcat({'\nColor_Convert_XYZ_to_RGB: '}, ...
        {'''Primary'' must be finite and real!\n\n'})))
    RGB = []; return
end

%% Denominator
denom = Primary(1, 1) * Primary(2, 2) * Primary(3, 3) + ...
    Primary(2, 1) * Primary(3, 2) * Primary(1, 3) + ...
    Primary(3, 1) * Primary(1, 2) * Primary(2, 3) - ...
    Primary(1, 1) * Primary(3, 2) * Primary(2, 3) - ...
    Primary(2, 1) * Primary(1, 2) * Primary(3, 3) - ...
    Primary(3, 1) * Primary(2, 2) * Primary(1, 3);

%% Components
rX = (Primary(2, 2) * Primary(3, 3) - Primary(2, 3) * Primary(3, 2)) / denom;
rY = (Primary(2, 3) * Primary(3, 1) - Primary(2, 1) * Primary(3, 3)) / denom;
rZ = (Primary(2, 1) * Primary(3, 2) - Primary(2, 2) * Primary(3, 1)) / denom;
gX = (Primary(3, 2) * Primary(1, 3) - Primary(3, 3) * Primary(1, 2)) / denom;
gY = (Primary(3, 3) * Primary(1, 1) - Primary(3, 1) * Primary(1, 3)) / denom;
gZ = (Primary(3, 1) * Primary(1, 2) - Primary(3, 2) * Primary(1, 1)) / denom;
bX = (Primary(1, 2) * Primary(2, 3) - Primary(1, 3) * Primary(2, 2)) / denom;
bY = (Primary(1, 3) * Primary(2, 1) - Primary(1, 1) * Primary(2, 3)) / denom;
bZ = (Primary(1, 1) * Primary(2, 2) - Primary(1, 2) * Primary(2, 1)) / denom; clear Primary denom

%% RGB Values
tR = rX .* XYZ(:, 1) + rY .* XYZ(:, 2) + rZ .* XYZ(:, 3); clear rX rY rZ
tG = gX .* XYZ(:, 1) + gY .* XYZ(:, 2) + gZ .* XYZ(:, 3); clear gX gY gZ
tB = bX .* XYZ(:, 1) + bY .* XYZ(:, 2) + bZ .* XYZ(:, 3); clear XYZ bX bY bZ
RGB = [tR, tG, tB]; clear tR tG tB
RGB = round(RGB .* 1000) ./ 1000;

% %% Variable Report (Housecleaning: for spotting uncleared variables.  Omit if desired)
% Variable_Report(whos, 'RGB')

end
