/**
 * Cube color-contrast functions for Google Sheets.
 *
 * INSTALL:  Extensions > Apps Script > paste this whole file > Save (disk icon).
 *           Go back to the sheet; the functions below work like any formula.
 *
 * USAGE (six hex codes in A2:F2):
 *   =MINCONTRAST(A2:F2)              worst-pair CAM16-UCS, normal vision
 *   =MINCONTRAST(A2:F2, "protan")    ... for a protanope
 *   =MINCONTRAST(A2:F2, "deutan")
 *   =MINCONTRAST(A2:F2, "tritan")
 *   =MINCONTRASTPAIR(A2:F2)          which two colors are the bottleneck
 *   =CONTRAST("#ff0000", "#00ff00")  distance between any two colors
 *   =CONTRASTTABLE(A2:F2)            all 15 pairs (spills down 15 rows x 3 cols)
 *
 * Any number of colors (2 or more) is fine; the range may be a row or column.
 * Hex codes may be with or without '#'.  Larger numbers = more contrast.
 *
 * MODEL: sRGB -> XYZ (D65) -> CAM16 -> CAM16-UCS.  Viewing condition:
 *   L_A = 64 cd/m2, Y_b = 20, "average" surround.
 * CVD simulation: Machado, Oliveira & Fernandes (2009), severity 1.0,
 * applied in linear sRGB before CAM16.
 * Matches the Python cube_contrast.py pipeline to ~0.01.
 */

// ----------------------------------------------------------------- constants
var VC_LA = 64.0, VC_YB = 20.0, VC_F = 1.0, VC_C = 0.69, VC_NC = 1.0;
var XYZ_W = [95.047, 100.0, 108.883];               // D65, 2 degree

var M16 = [[ 0.401288, 0.650173, -0.051461],
           [-0.250268, 1.204414,  0.045854],
           [-0.002079, 0.048952,  0.953127]];

var M_sRGB_XYZ = [[0.4124, 0.3576, 0.1805],
                  [0.2126, 0.7152, 0.0722],
                  [0.0193, 0.1192, 0.9505]];

var CVD = {
  protan: [[ 0.152286, 1.052583, -0.204868],
           [ 0.114503, 0.786281,  0.099216],
           [-0.003882,-0.048116,  1.051998]],
  deutan: [[ 0.367322, 0.860646, -0.227968],
           [ 0.280085, 0.672501,  0.047413],
           [-0.011820, 0.042940,  0.968881]],
  tritan: [[ 1.255528,-0.076749, -0.178779],
           [-0.078411, 0.930809,  0.147602],
           [ 0.004733, 0.691367,  0.303900]]
};

// ------------------------------------------------------------------- helpers
function mul_(M, v) {
  return [M[0][0]*v[0]+M[0][1]*v[1]+M[0][2]*v[2],
          M[1][0]*v[0]+M[1][1]*v[1]+M[1][2]*v[2],
          M[2][0]*v[0]+M[2][1]*v[1]+M[2][2]*v[2]];
}

function hexToLinearRGB_(hex) {
  var h = String(hex).trim().replace(/^#/, "");
  if (h.length === 3) h = h[0]+h[0]+h[1]+h[1]+h[2]+h[2];
  if (!/^[0-9a-fA-F]{6}$/.test(h)) throw new Error("Bad hex code: " + hex);
  var out = [];
  for (var i = 0; i < 3; i++) {
    var c = parseInt(h.substr(i*2, 2), 16) / 255;
    out.push(c <= 0.04045 ? c/12.92 : Math.pow((c+0.055)/1.055, 2.4));
  }
  return out;
}

function adapt_(x, FL) {                       // CAM16 post-adaptation cone response
  var p = Math.pow(FL * Math.abs(x) / 100, 0.42);
  return (x < 0 ? -1 : 1) * 400 * p / (p + 27.13) + 0.1;
}

// Precompute white-point terms once.
var VC_ = (function () {
  var k = 1 / (5*VC_LA + 1), k4 = k*k*k*k;
  var FL = 0.2*k4*5*VC_LA + 0.1*(1-k4)*(1-k4)*Math.pow(5*VC_LA, 1/3);
  var n = VC_YB / XYZ_W[1];
  var z = 1.48 + Math.sqrt(n);
  var Nbb = 0.725 * Math.pow(1/n, 0.2);
  var D = VC_F * (1 - (1/3.6) * Math.exp((-VC_LA - 42) / 92));
  D = Math.max(0, Math.min(1, D));
  var RGBw = mul_(M16, XYZ_W);
  var Dc = [D*XYZ_W[1]/RGBw[0] + 1 - D,
            D*XYZ_W[1]/RGBw[1] + 1 - D,
            D*XYZ_W[1]/RGBw[2] + 1 - D];
  var RGBaw = [adapt_(Dc[0]*RGBw[0], FL), adapt_(Dc[1]*RGBw[1], FL), adapt_(Dc[2]*RGBw[2], FL)];
  var Aw = (2*RGBaw[0] + RGBaw[1] + RGBaw[2]/20 - 0.305) * Nbb;
  return {FL: FL, n: n, z: z, Nbb: Nbb, Dc: Dc, Aw: Aw};
})();

// XYZ (0-100 scale, D65) -> CAM16-UCS [J', a', b']
function xyzToUCS_(XYZ) {
  var v = VC_;
  var RGB = mul_(M16, XYZ);
  var Ra = adapt_(v.Dc[0]*RGB[0], v.FL);
  var Ga = adapt_(v.Dc[1]*RGB[1], v.FL);
  var Ba = adapt_(v.Dc[2]*RGB[2], v.FL);
  var a = Ra - 12*Ga/11 + Ba/11;
  var b = (Ra + Ga - 2*Ba) / 9;
  var h = Math.atan2(b, a); if (h < 0) h += 2*Math.PI;
  var et = 0.25 * (Math.cos(h + 2) + 3.8);
  var A = (2*Ra + Ga + Ba/20 - 0.305) * v.Nbb;
  var J = 100 * Math.pow(A / v.Aw, VC_C * v.z);
  var t = (50000/13) * VC_NC * v.Nbb * et * Math.sqrt(a*a + b*b) / (Ra + Ga + 21*Ba/20);
  var C = Math.pow(t, 0.9) * Math.sqrt(J/100) * Math.pow(1.64 - Math.pow(0.29, v.n), 0.73);
  var M = C * Math.pow(v.FL, 0.25);
  var Jp = 1.7*J / (1 + 0.007*J);
  var Mp = Math.log(1 + 0.0228*M) / 0.0228;
  return [Jp, Mp*Math.cos(h), Mp*Math.sin(h)];
}

function hexToUCS_(hex, observer) {
  var rgb = hexToLinearRGB_(hex);
  var obs = String(observer || "normal").toLowerCase();
  if (obs !== "normal") {
    if (!CVD[obs]) throw new Error('observer must be normal, protan, deutan or tritan');
    rgb = mul_(CVD[obs], rgb);              // no clipping: transform is linear
  }
  var XYZ = mul_(M_sRGB_XYZ, rgb);
  return xyzToUCS_([XYZ[0]*100, XYZ[1]*100, XYZ[2]*100]);
}

function flatten_(range) {
  var out = [];
  if (!Array.isArray(range)) return [range];
  for (var i = 0; i < range.length; i++) {
    var row = Array.isArray(range[i]) ? range[i] : [range[i]];
    for (var j = 0; j < row.length; j++) {
      if (row[j] !== "" && row[j] !== null && row[j] !== undefined) out.push(row[j]);
    }
  }
  return out;
}

function allPairs_(range, observer) {
  var hexes = flatten_(range);
  if (hexes.length < 2) throw new Error("Need at least two colors");
  var ucs = hexes.map(function (h) { return hexToUCS_(h, observer); });
  var pairs = [];
  for (var i = 0; i < ucs.length; i++) for (var j = i+1; j < ucs.length; j++) {
    var d = Math.sqrt(Math.pow(ucs[i][0]-ucs[j][0], 2) +
                      Math.pow(ucs[i][1]-ucs[j][1], 2) +
                      Math.pow(ucs[i][2]-ucs[j][2], 2));
    pairs.push({a: hexes[i], b: hexes[j], d: d});
  }
  pairs.sort(function (p, q) { return p.d - q.d; });
  return pairs;
}

// ---------------------------------------------------------- sheet functions

/**
 * Minimum pairwise CAM16-UCS contrast across a set of hex colors.
 * @param {Range} range  Cells containing hex codes (e.g. A2:F2).
 * @param {string} observer  "normal" (default), "protan", "deutan" or "tritan".
 * @return {number} Worst-pair contrast; larger is better.
 * @customfunction
 */
function MINCONTRAST(range, observer) {
  return allPairs_(range, observer)[0].d;
}

/**
 * Names the two colors that form the bottleneck pair.
 * @param {Range} range  Cells containing hex codes.
 * @param {string} observer  "normal" (default), "protan", "deutan" or "tritan".
 * @return {string} e.g. "#ffe500 / #f4f4f0"
 * @customfunction
 */
function MINCONTRASTPAIR(range, observer) {
  var p = allPairs_(range, observer)[0];
  return p.a + " / " + p.b;
}

/**
 * CAM16-UCS contrast between two hex colors.
 * @param {string} hex1
 * @param {string} hex2
 * @param {string} observer  "normal" (default), "protan", "deutan" or "tritan".
 * @return {number}
 * @customfunction
 */
function CONTRAST(hex1, hex2, observer) {
  return allPairs_([[hex1, hex2]], observer)[0].d;
}

/**
 * Every pair, sorted worst-first.  Spills into rows below.
 * @param {Range} range  Cells containing hex codes.
 * @param {string} observer  "normal" (default), "protan", "deutan" or "tritan".
 * @return {Array} rows of [color A, color B, contrast]
 * @customfunction
 */
function CONTRASTTABLE(range, observer) {
  return allPairs_(range, observer).map(function (p) { return [p.a, p.b, p.d]; });
}

/**
 * Mean pairwise contrast (secondary statistic; the minimum is the score).
 * @param {Range} range
 * @param {string} observer
 * @return {number}
 * @customfunction
 */
function MEANCONTRAST(range, observer) {
  var p = allPairs_(range, observer), s = 0;
  for (var i = 0; i < p.length; i++) s += p[i].d;
  return s / p.length;
}