(function(){science = {version: “1.7.0”}; // semver science.ascending = function(a, b) {

return a - b;

}; // Euler's constant. science.EULER = .5772156649015329; // Compute exp(x) - 1 accurately for small x. science.expm1 = function(x) {

return (x < 1e-5 && x > -1e-5) ? x + .5 * x * x : Math.exp(x) - 1;

}; science.functor = function(v) {

return typeof v === "function" ? v : function() { return v; };

}; // Based on: // www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/ science.hypot = function(x, y) {

x = Math.abs(x);
y = Math.abs(y);
var max,
    min;
if (x > y) { max = x; min = y; }
else       { max = y; min = x; }
var r = min / max;
return max * Math.sqrt(1 + r * r);

}; science.quadratic = function() {

var complex = false;

function quadratic(a, b, c) {
  var d = b * b - 4 * a * c;
  if (d > 0) {
    d = Math.sqrt(d) / (2 * a);
    return complex
      ? [{r: -b - d, i: 0}, {r: -b + d, i: 0}]
      : [-b - d, -b + d];
  } else if (d === 0) {
    d = -b / (2 * a);
    return complex ? [{r: d, i: 0}] : [d];
  } else {
    if (complex) {
      d = Math.sqrt(-d) / (2 * a);
      return [
        {r: -b, i: -d},
        {r: -b, i: d}
      ];
    }
    return [];
  }
}

quadratic.complex = function(x) {
  if (!arguments.length) return complex;
  complex = x;
  return quadratic;
};

return quadratic;

}; // Constructs a multi-dimensional array filled with zeroes. science.zeroes = function(n) {

var i = -1,
    a = [];
if (arguments.length === 1)
  while (++i < n)
    a[i] = 0;
else
  while (++i < n)
    a[i] = science.zeroes.apply(
      this, Array.prototype.slice.call(arguments, 1));
return a;

}; science.vector = {}; science.vector.cross = function(a, b) {

// TODO how to handle non-3D vectors?
// TODO handle 7D vectors?
return [
  a[1] * b[2] - a[2] * b[1],
  a[2] * b[0] - a[0] * b[2],
  a[0] * b[1] - a[1] * b[0]
];

}; science.vector.dot = function(a, b) {

var s = 0,
    i = -1,
    n = Math.min(a.length, b.length);
while (++i < n) s += a[i] * b[i];
return s;

}; science.vector.length = function(p) {

return Math.sqrt(science.vector.dot(p, p));

}; science.vector.normalize = function(p) {

var length = science.vector.length(p);
return p.map(function(d) { return d / length; });

}; // 4x4 matrix determinant. science.vector.determinant = function(matrix) {

var m = matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]);
return (
  m[12] * m[9]  * m[6]  * m[3]  - m[8] * m[13] * m[6]  * m[3]  -
  m[12] * m[5]  * m[10] * m[3]  + m[4] * m[13] * m[10] * m[3]  +
  m[8]  * m[5]  * m[14] * m[3]  - m[4] * m[9]  * m[14] * m[3]  -
  m[12] * m[9]  * m[2]  * m[7]  + m[8] * m[13] * m[2]  * m[7]  +
  m[12] * m[1]  * m[10] * m[7]  - m[0] * m[13] * m[10] * m[7]  -
  m[8]  * m[1]  * m[14] * m[7]  + m[0] * m[9]  * m[14] * m[7]  +
  m[12] * m[5]  * m[2]  * m[11] - m[4] * m[13] * m[2]  * m[11] -
  m[12] * m[1]  * m[6]  * m[11] + m[0] * m[13] * m[6]  * m[11] +
  m[4]  * m[1]  * m[14] * m[11] - m[0] * m[5]  * m[14] * m[11] -
  m[8]  * m[5]  * m[2]  * m[15] + m[4] * m[9]  * m[2]  * m[15] +
  m[8]  * m[1]  * m[6]  * m[15] - m[0] * m[9]  * m[6]  * m[15] -
  m[4]  * m[1]  * m[10] * m[15] + m[0] * m[5]  * m[10] * m[15]);

}; // Performs in-place Gauss-Jordan elimination. // // Based on Jarno Elonen's Python version (public domain): // elonen.iki.fi/code/misc-notes/python-gaussj/index.html science.vector.gaussjordan = function(m, eps) {

if (!eps) eps = 1e-10;

var h = m.length,
    w = m[0].length,
    y = -1,
    y2,
    x;

while (++y < h) {
  var maxrow = y;

  // Find max pivot.
  y2 = y; while (++y2 < h) {
    if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y]))
      maxrow = y2;
  }

  // Swap.
  var tmp = m[y];
  m[y] = m[maxrow];
  m[maxrow] = tmp;

  // Singular?
  if (Math.abs(m[y][y]) <= eps) return false;

  // Eliminate column y.
  y2 = y; while (++y2 < h) {
    var c = m[y2][y] / m[y][y];
    x = y - 1; while (++x < w) {
      m[y2][x] -= m[y][x] * c;
    }
  }
}

// Backsubstitute.
y = h; while (--y >= 0) {
  var c = m[y][y];
  y2 = -1; while (++y2 < y) {
    x = w; while (--x >= y) {
      m[y2][x] -=  m[y][x] * m[y2][y] / c;
    }
  }
  m[y][y] /= c;
  // Normalize row y.
  x = h - 1; while (++x < w) {
    m[y][x] /= c;
  }
}
return true;

}; // Find matrix inverse using Gauss-Jordan. science.vector.inverse = function(m) {

var n = m.length
    i = -1;

// Check if the matrix is square.
if (n !== m[0].length) return;

// Augment with identity matrix I to get AI.
m = m.map(function(row, i) {
  var identity = new Array(n),
      j = -1;
  while (++j < n) identity[j] = i === j ? 1 : 0;
  return row.concat(identity);
});

// Compute IA^-1.
science.vector.gaussjordan(m);

// Remove identity matrix I to get A^-1.
while (++i < n) {
  m[i] = m[i].slice(n);
}

return m;

}; science.vector.multiply = function(a, b) {

var m = a.length,
    n = b[0].length,
    p = b.length,
    i = -1,
    j,
    k;
if (p !== a[0].length) throw {"error": "columns(a) != rows(b); " + a[0].length + " != " + p};
var ab = new Array(m);
while (++i < m) {
  ab[i] = new Array(n);
  j = -1; while(++j < n) {
    var s = 0;
    k = -1; while (++k < p) s += a[i][k] * b[k][j];
    ab[i][j] = s;
  }
}
return ab;

}; science.vector.transpose = function(a) {

var m = a.length,
    n = a[0].length,
    i = -1,
    j,
    b = new Array(n);
while (++i < n) {
  b[i] = new Array(m);
  j = -1; while (++j < m) b[i][j] = a[j][i];
}
return b;

}; })()