/* Authors: Alexei Boronine, Bastien Dejean */

/* Matrices generated with the following algorithm:                 */
/* http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html */
/* Which makes use of this transformation:                          */
/* http://www.brucelindbloom.com/index.html?Eqn_xyY_to_XYZ.html     */

xyY_to_XYZ(x, y, Y) := [x * Y / y, Y, (1 - x - y) * Y / y];

/* sRGB primaries with D65 white point               */
/* http://en.wikipedia.org/wiki/SRGB#The_sRGB_gamut  */
/* http://en.wikipedia.org/wiki/Illuminant_D65       */

x_r: rat(0.64);
y_r: rat(0.33);
x_g: rat(0.30);
y_g: rat(0.60);
x_b: rat(0.15);
y_b: rat(0.06);
x_w: rat(0.3127);
y_w: rat(0.3290);

[X_r, Y_r, Z_r]: xyY_to_XYZ(x_r, y_r, 1);
[X_g, Y_g, Z_g]: xyY_to_XYZ(x_g, y_g, 1);
[X_b, Y_b, Z_b]: xyY_to_XYZ(x_b, y_b, 1);
[X_w, Y_w, Z_w]: xyY_to_XYZ(x_w, y_w, 1);
 
M_P: matrix([X_r, X_g, X_b],
            [Y_r, Y_g, Y_b],
            [Z_r, Z_g, Z_b]);
 
M_S: invert(M_P) . [X_w, Y_w, Z_w];

M_RGB_XYZ: M_P * transpose(addcol(M_S, M_S, M_S));
M_XYZ_RGB: invert(M_RGB_XYZ);

fpprec: 16;
bfloat(M_RGB_XYZ);
bfloat(M_XYZ_RGB);

/* XYZ tristimulus values for D65, used as reference white point here:
http://www.brucelindbloom.com/index.html?Eqn_xyY_to_XYZ.html */

xyY_to_XYZ(rat(0.3127), rat(0.3290), 1);