o[1] = m[3] * ax + m[4] * ay + m[5];
  return o;
}

//o = vNewTransform (a, m)
//  Array o
//  Array a
//  Array m
function vNewTransform (a, m) {
  let ax = a[0], ay = a[1];  //double
  return [m[0] * ax + m[1] * ay + m[2],
          m[3] * ax + m[4] * ay + m[5]];
}

//------------------------------------------------------------------------
//translate 平行移動
//
//    [1 0 d]   [x]   [x+d]
//    [0 1 e] * [y] = [y+e]
//    [0 0 1]   [1]   [ 1 ]
//
//  a+b

//a = vTranslateV (a, b)
//  Array a
//  Array b
function vTranslateV (a, b) {
  a[0] += b[0];
  a[1] += b[1];
  return a;
}

//o = vLetTranslateV (o, a, b)
//  Array o
//  Array a
//  Array b
function vLetTranslateV (o, a, b) {
  o[0] = a[0] + b[0];
  o[1] = a[1] + b[1];
  return o;
}

//o = vNewTranslateV (a, b)
//  Array o
//  Array a
//  Array b
function vNewTranslateV (a, b) {
  return [a[0] + b[0],
          a[1] + b[1]];
}

//a = vTranslateX (a, d)
//  Array a
//  double d
function vTranslateX (a, d) {
  a[0] += d;
  return a;
}

//o = vLetTranslateX (o, a, d)
//  Array o
//  Array a
//  double d
function vLetTranslateX (o, a, d) {
  o[0] = a[0] + d;
  o[1] = a[1];
  return o;
}

//o = vNewTranslateX (a, d)
//  Array o
//  Array a
//  double d
function vNewTranslateX (a, d) {
  return [a[0] + d,
          a[1]];
}

//a = vTranslateY (a, e)
//  Array a
//  double e
function vTranslateY (a, e) {
  a[1] += e;
  return a;
}

//o = vLetTranslateY (o, a, e)
//  Array o
//  Array a
//  double e
function vLetTranslateY (o, a, e) {
  o[0] = a[0];
  o[1] = a[1] + e;
  return o;
}

//o = vNewTranslateY (a, e)
//  Array o
//  Array a
//  double e
function vNewTranslateY (a, e) {
  return [a[0],
          a[1] + e];
}

//a = vTranslateXY (a, d, e)
//  Array a
//  double d
//  double e
function vTranslateXY (a, d, e) {
  a[0] += d;
  a[1] += e;
  return a;
}

//o = vLetTranslateXY (o, a, d, e)
//  Array o
//  Array a
//  double d
//  double e
function vLetTranslateXY (o, a, d, e) {
  o[0] = a[0] + d;
  o[1] = a[1] + e;
  return o;
}

//o = vNewTranslateXY (a, d, e)
//  Array o
//  Array a
//  double d
//  double e
function vNewTranslateXY (a, d, e) {
  return [a[0] + d,
          a[1] + e];
}

//------------------------------------------------------------------------
//trisect 3分割
//  2/3*a+1/3*b

//a = vTrisect (a, b)
//  Array a
//  Array b
function vTrisect (a, b) {
  a[0] = TWO_3 * a[0] + ONE_3 * b[0];
  a[1] = TWO_3 * a[1] + ONE_3 * b[1];
  return a;
}

//o = vLetTrisect (o, a, b)
//  Array o
//  Array a
//  Array b
function vLetTrisect (o, a, b) {
  o[0] = TWO_3 * a[0] + ONE_3 * b[0];
  o[1] = TWO_3 * a[1] + ONE_3 * b[1];
  return o;
}

//o = vNewTrisect (a, b)
//  Array o
//  Array a
//  Array b
function vNewTrisect (a, b) {
  return [TWO_3 * a[0] + ONE_3 * b[0],
          TWO_3 * a[1] + ONE_3 * b[1]];
}

//------------------------------------------------------------------------
//zero 0ベクトル
//  0

//o = vLetZero (o)
//  Array o
function vLetZero (o) {
  o[0] = 0;
  o[1] = 0;
  return o;
}

//o = vNewZero ()
//  Array o
function vNewZero () {
  return [0,
          0];
}



//========================================================================
//空間ベクトル(w～)
//
//    [x]   [w[0]]
//    [y] = [w[1]]
//    [z]   [w[2]]
//    [1]   [  1 ]
//

//------------------------------------------------------------------------
//angle なす角
//  acos((a.b)/(|a|*|b|))
//  ゼロベクトルは不可

//t = wAngle (a, b)
//  double t
//  Array a
//  Array b
function wAngle (a, b) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  return acos ((ax * bx + ay * by + az * bz) / (hypot3 (ax, ay, az) * hypot3 (bx, by, bz)));
}

//------------------------------------------------------------------------
//angle a 三角形abcの角a
//  辺abと辺acのなす角

//t = wAngleA (a, b, c)
//  double t
//  Array a
//  Array b
//  Array c
function wAngleA (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  bx -= ax;
  by -= ay;
  bz -= az;
  cx -= ax;
  cy -= ay;
  cz -= az;
  return acos ((bx * cx + by * cy + bz * cz) / (hypot3 (bx, by, bz) * hypot3 (cx, cy, cz)));
}

//------------------------------------------------------------------------
//areaoftriangle 三角形の面積
//  三角形abcの面積
//  平面
//    s=1/2*abs((bx-ax)*(cy-ay)-(cx-ax)*(by-ay))
//  空間
//    s=1/2*sqrt(abs(b-a)^2*abs(c-a)^2-((b-a).(c-a))^2)
function wAreaOfTriangle (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];
  let bx = b[0], by = b[1], bz = b[2];
  let cx = c[0], cy = c[1], cz = c[2];
  bx -= ax;
  by -= ay;
  bz -= az;
  cx -= ax;
  cy -= ay;
  cz -= az;
  return 0.5 * sqrt ((bx * bx + by * by + bz * bz) * (cx * cx + cy * cy + cz * cz) - (bx * cx + by * cy + bz * cz) ** 2);
}

//------------------------------------------------------------------------
//azimuthYZ 方位角(+y→+z)

//t = wAzimuthYZ (a)
//  double t
//  Array a
function wAzimuthYZ (a) {
  return atan2 (a[2], a[1]);
}

//------------------------------------------------------------------------
//azimuthZX 方位角(+z→+x)

//t = wAzimuthZX (a)
//  double t
//  Array a
function wAzimuthZX (a) {
  return atan2 (a[0], a[2]);
}

//------------------------------------------------------------------------
//azimuthXY 方位角(+x→+y)

//t = wAzimuthXY (a)
//  double t
//  Array a
function wAzimuthXY (a) {
  return atan2 (a[1], a[0]);
}

//------------------------------------------------------------------------
//bisect 2分割(中点)
//  (a+b)/2

//a = wBisect (a, b)
//  Array a
//  Array b
function wBisect (a, b) {
  a[0] = (a[0] + b[0]) * 0.5;
  a[1] = (a[1] + b[1]) * 0.5;
  a[2] = (a[2] + b[2]) * 0.5;
  return a;
}

//o = wLetBisect (o, a, b)
//  Array o
//  Array a
//  Array b
function wLetBisect (o, a, b) {
  o[0] = (a[0] + b[0]) * 0.5;
  o[1] = (a[1] + b[1]) * 0.5;
  o[2] = (a[2] + b[2]) * 0.5;
  return o;
}

//o = wNewBisect (a, b)
//  Array o
//  Array a
//  Array b
function wNewBisect (a, b) {
  return [(a[0] + b[0]) * 0.5,
          (a[1] + b[1]) * 0.5,
          (a[2] + b[2]) * 0.5];
}

//------------------------------------------------------------------------
//centroid 三角形の重心
//  中線(頂点と辺の中点を結ぶ直線)の交点
//  G=(A+B+C)/3

//a = wCentroid (a, b, c)
//  Array a
//  Array b
//  Array c
function wCentroid (a, b, c) {
  a[0] = (a[0] + b[0] + c[0]) * ONE_3;
  a[1] = (a[1] + b[1] + c[1]) * ONE_3;
  a[2] = (a[2] + b[2] + c[2]) * ONE_3;
  return a;
}

//o = wLetCentroid (o, a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wLetCentroid (o, a, b, c) {
  o[0] = (a[0] + b[0] + c[0]) * ONE_3;
  o[1] = (a[1] + b[1] + c[1]) * ONE_3;
  o[2] = (a[2] + b[2] + c[2]) * ONE_3;
  return o;
}

//o = wNewCentroid (a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wNewCentroid (a, b, c) {
  return [(a[0] + b[0] + c[0]) * ONE_3,
          (a[1] + b[1] + c[1]) * ONE_3,
          (a[2] + b[2] + c[2]) * ONE_3];
}

//------------------------------------------------------------------------
//circumcenter 三角形の外心
//  外接円の中心。辺の垂直二等分線の交点
//  O=(sin(2*A)*A+sin(2*B)*B+sin(2*C)*C)/(sin(2*A)+sin(2*B)+sin(2*C))

//a = wCircumcenter (a, b, c)
//  Array a
//  Array b
//  Array c
function wCircumcenter (a, b, c) {
  let cosa = wCosAngleA (a, b, c);  //double
  let cosb = wCosAngleA (b, c, a);  //double
  let cosc = wCosAngleA (c, a, b);  //double
  let sin2a = 2 * sqrt (1 - cosa * cosa) * cosa;  //double
  let sin2b = 2 * sqrt (1 - cosb * cosb) * cosb;  //double
  let sin2c = 2 * sqrt (1 - cosc * cosc) * cosc;  //double
  let d = sin2a + sin2b + sin2c;  //double
  a[0] = (sin2a * a[0] + sin2b * b[0] + sin2c * c[0]) / d;
  a[1] = (sin2a * a[1] + sin2b * b[1] + sin2c * c[1]) / d;
  a[2] = (sin2a * a[2] + sin2b * b[2] + sin2c * c[2]) / d;
  return a;
}

//o = wLetCircumcenter (o, a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wLetCircumcenter (o, a, b, c) {
  let cosa = wCosAngleA (a, b, c);  //double
  let cosb = wCosAngleA (b, c, a);  //double
  let cosc = wCosAngleA (c, a, b);  //double
  let sin2a = 2 * sqrt (1 - cosa * cosa) * cosa;  //double
  let sin2b = 2 * sqrt (1 - cosb * cosb) * cosb;  //double
  let sin2c = 2 * sqrt (1 - cosc * cosc) * cosc;  //double
  let d = sin2a + sin2b + sin2c;  //double
  o[0] = (sin2a * a[0] + sin2b * b[0] + sin2c * c[0]) / d;
  o[1] = (sin2a * a[1] + sin2b * b[1] + sin2c * c[1]) / d;
  o[2] = (sin2a * a[2] + sin2b * b[2] + sin2c * c[2]) / d;
  return o;
}

//o = wNewCircumcenter (a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wNewCircumcenter (a, b, c) {
  let cosa = wCosAngleA (a, b, c);  //double
  let cosb = wCosAngleA (b, c, a);  //double
  let cosc = wCosAngleA (c, a, b);  //double
  let sin2a = 2 * sqrt (1 - cosa * cosa) * cosa;  //double
  let sin2b = 2 * sqrt (1 - cosb * cosb) * cosb;  //double
  let sin2c = 2 * sqrt (1 - cosc * cosc) * cosc;  //double
  let d = sin2a + sin2b + sin2c;  //double
  return [(sin2a * a[0] + sin2b * b[0] + sin2c * c[0]) / d,
          (sin2a * a[1] + sin2b * b[1] + sin2c * c[1]) / d,
          (sin2a * a[2] + sin2b * b[2] + sin2c * c[2]) / d];
}

//------------------------------------------------------------------------
//copy コピー
//  a

//o = wLetCopy (o, a)
//  Array o
//  Array a
function wLetCopy (o, a) {
  o[0] = a[0];
  o[1] = a[1];
  o[2] = a[2];
  return o;
}

//o = wNewCopy (a)
//  Array o
//  Array a
function wNewCopy (a) {
  return a.concat ();
}

//------------------------------------------------------------------------
//cos angle なす角のコサイン
//  cos(t)=(a.b)/(|a|*|b|)
//  ゼロベクトルは不可

//cost = wCosAngle (a, b)
//  double cost
//  Array a
//  Array b
function wCosAngle (a, b) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  return (ax * bx + ay * by + az * bz) / (hypot3 (ax, ay, az) * hypot3 (bx, by, bz));
}

//------------------------------------------------------------------------
//cos angle a 三角形abcの角aのコサイン
//  辺abと辺acのなす角のコサイン

//cost = wCosAngleA (a, b, c)
//  double cost
//  Array a
//  Array b
//  Array c
function wCosAngleA (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  bx -= ax;
  by -= ay;
  bz -= az;
  cx -= ax;
  cy -= ay;
  cz -= az;
  return (bx * cx + by * cy + bz * cz) / (hypot3 (bx, by, bz) * hypot3 (cx, cy, cz));
}

//------------------------------------------------------------------------
//cross クロス積(外積,ベクトル積)
//  |a×b|=|a|*|b|*sin(t)

//a = wCross (a, b)
//  Array a
//  Array b
function wCross (a, b) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  a[0] = ay * bz - az * by;
  a[1] = az * bx - ax * bz;
  a[2] = ax * by - ay * bx;
  return a;
}

//o = wLetCross (o, a, b)
//  Array o
//  Array a
//  Array b
function wLetCross (o, a, b) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  o[0] = ay * bz - az * by;
  o[1] = az * bx - ax * bz;
  o[2] = ax * by - ay * bx;
  return o;
}

//o = wNewCross (a, b)
//  Array o
//  Array a
//  Array b
function wNewCross (a, b) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  return [ay * bz - az * by,
          az * bx - ax * bz,
          ax * by - ay * bx];
}

//------------------------------------------------------------------------
//distance 点と点の距離
//  |b-a|

//s = wDistance (a, b)
//  double s
//  Array a
//  Array b
function wDistance (a, b) {
  return hypot3 (b[0] - a[0], b[1] - a[1], b[2] - a[2]);
}

//------------------------------------------------------------------------
//distance2 点と点の距離の2乗
//  |b-a|^2

//s = wDistance2 (a, b)
//  double s
//  Array a
//  Array b
function wDistance2 (a, b) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  bx -= ax;
  by -= ay;
  bz -= az;
  return bx * bx + by * by + bz * bz;
}

//------------------------------------------------------------------------
//dot ドット積(内積,スカラー積)
//  a.b=|a|*|b|*cos(angle(a,b))

//s = wDot (a, b)
//  double s
//  Array a
//  Array b
function wDot (a, b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

//------------------------------------------------------------------------
//hat 単位ベクトル化
//  a/|a|
//  ゼロベクトルは不可

//a = wHat (a)
//  Array a
function wHat (a) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let s = hypot3 (ax, ay, az);  //double
  a[0] = ax / s;
  a[1] = ay / s;
  a[2] = az / s;
  return a;
}

//o = wLetHat (o, a)
//  Array o
//  Array a
function wLetHat (o, a) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let s = hypot3 (ax, ay, az);  //double
  o[0] = ax / s;
  o[1] = ay / s;
  o[2] = az / s;
  return o;
}

//o = wNewHat (a)
//  Array o
//  Array a
function wNewHat (a) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let s = hypot3 (ax, ay, az);  //double
  return [ax / s,
          ay / s,
          az / s];
}

//------------------------------------------------------------------------
//height 三角形abcの高さ
//  点aと直線bcの距離。点aから直線bcへ下ろした垂線の長さ

//s = wHeight (a, b, c)
//  double s
//  Array a
//  Array b
//  Array c
function wHeight (a, b, c) {
  return wDistance (a, wNewPerpendicularFoot (a, b, c));
}

//------------------------------------------------------------------------
//incenter 三角形の内心
//  内接円の中心。角の二等分線の交点
//  I=(a*A+b*B+c*C)/(a+b+c)

//a = wIncenter (a, b, c)
//  Array a
//  Array b
//  Array c
function wIncenter (a, b, c) {
  let la = wDistance (b, c);  //double
  let lb = wDistance (c, a);  //double
  let lc = wDistance (a, b);  //double
  let d = la + lb + lc;  //double
  a[0] = (la * a[0] + lb * b[0] + lc * c[0]) / d;
  a[1] = (la * a[1] + lb * b[1] + lc * c[1]) / d;
  a[2] = (la * a[2] + lb * b[2] + lc * c[2]) / d;
  return a;
}

//o = wLetIncenter (o, a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wLetIncenter (o, a, b, c) {
  let la = wDistance (b, c);  //double
  let lb = wDistance (c, a);  //double
  let lc = wDistance (a, b);  //double
  let d = la + lb + lc;  //double
  o[0] = (la * a[0] + lb * b[0] + lc * c[0]) / d;
  o[1] = (la * a[1] + lb * b[1] + lc * c[1]) / d;
  o[2] = (la * a[2] + lb * b[2] + lc * c[2]) / d;
  return o;
}

//o = wNewIncenter (a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wNewIncenter (a, b, c) {
  let la = wDistance (b, c);  //double
  let lb = wDistance (c, a);  //double
  let lc = wDistance (a, b);  //double
  let d = la + lb + lc;  //double
  return [(la * a[0] + lb * b[0] + lc * c[0]) / d,
          (la * a[1] + lb * b[1] + lc * c[1]) / d,
          (la * a[2] + lb * b[2] + lc * c[2]) / d];
}

//------------------------------------------------------------------------
//norm 長さ(ノルム)
//  |a|=sqrt(a.a)

//s = wNorm (a)
//  double s
//  Array a
function wNorm (a) {
  return hypot3 (a[0], a[1], a[2]);
}

//------------------------------------------------------------------------
//norm2 長さ(ノルム)の2乗
//  |a|^2=a.a

//s = wNorm2 (a)
//  double s
//  Array a
function wNorm2 (a) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  return ax * ax + ay * ay + az * az;
}

//------------------------------------------------------------------------
//normal 三角形(a,b,c)の法線ベクトル
//  |(b-a)×(c-a)|=|b-a|*|c-a|*sin(t)
//  右手座標系のとき、頂点の順序は法線ベクトルの方向に進む右ねじを回転させる方向
//  (多面体の内側から見ると右回り、外側から見ると左回り)
//  法線ベクトルは単位ベクトルになっていないので必要に応じて単位ベクトルに変換すること

//a = wNormal (a, b, c)
//  Array a
//  Array b
//  Array c
function wNormal (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  bx -= ax;
  by -= ay;
  bz -= az;
  cx -= ax;
  cy -= ay;
  cz -= az;
  a[0] = by * cz - bz * cy;
  a[1] = bz * cx - bx * cz;
  a[2] = bx * cy - by * cx;
  return a;
}

//o = wLetNormal (o, a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wLetNormal (o, a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  bx -= ax;
  by -= ay;
  bz -= az;
  cx -= ax;
  cy -= ay;
  cz -= az;
  o[0] = by * cz - bz * cy;
  o[1] = bz * cx - bx * cz;
  o[2] = bx * cy - by * cx;
  return o;
}

//o = wNewNormal (a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wNewNormal (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  bx -= ax;
  by -= ay;
  bz -= az;
  cx -= ax;
  cy -= ay;
  cz -= az;
  return [by * cz - bz * cy,
          bz * cx - bx * cz,
          bx * cy - by * cx];
}

//------------------------------------------------------------------------
//one 1ベクトル
//  1

//o = wLetOneX (o)
//  Array o
function wLetOneX (o) {
  o[0] = 1;
  o[1] = 0;
  o[2] = 0;
  return o;
}

//o = wNewOneX ()
//  Array o
function wNewOneX () {
  return [1,
          0,
          0];
}

//o = wLetOneY (o)
//  Array o
function wLetOneY (o) {
  o[0] = 0;
  o[1] = 1;
  o[2] = 0;
  return o;
}

//o = wNewOneY ()
//  Array o
function wNewOneY () {
  return [0,
          1,
          0];
}

//o = wLetOneZ (o)
//  Array o
function wLetOneZ (o) {
  o[0] = 0;
  o[1] = 0;
  o[2] = 1;
  return o;
}

//o = wNewOneZ ()
//  Array o
function wNewOneZ () {
  return [0,
          0,
          1];
}

//------------------------------------------------------------------------
//orthocenter 三角形の垂心
//  頂点から辺に下ろした垂線の交点
//  H=(tan(A)*A+tan(B)*B+tan(C)*C)/(tan(A)+tan(B)+tan(C))

//a = wOrthocenter (a, b, c)
//  Array a
//  Array b
//  Array c
function wOrthocenter (a, b, c) {
  let cosa = wCosAngleA (a, b, c);  //double
  let cosb = wCosAngleA (b, c, a);  //double
  let cosc = wCosAngleA (c, a, b);  //double
  let tana = sqrt (1 - cosa * cosa) / cosa;  //double
  let tanb = sqrt (1 - cosb * cosb) / cosb;  //double
  let tanc = sqrt (1 - cosc * cosc) / cosc;  //double
  let d = tana + tanb + tanc;  //double
  a[0] = (tana * a[0] + tanb * b[0] + tanc * c[0]) / d;
  a[1] = (tana * a[1] + tanb * b[1] + tanc * c[1]) / d;
  a[2] = (tana * a[2] + tanb * b[2] + tanc * c[2]) / d;
  return a;
}

//o = wLetOrthocenter (o, a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wLetOrthocenter (o, a, b, c) {
  let cosa = wCosAngleA (a, b, c);  //double
  let cosb = wCosAngleA (b, c, a);  //double
  let cosc = wCosAngleA (c, a, b);  //double
  let tana = sqrt (1 - cosa * cosa) / cosa;  //double
  let tanb = sqrt (1 - cosb * cosb) / cosb;  //double
  let tanc = sqrt (1 - cosc * cosc) / cosc;  //double
  let d = tana + tanb + tanc;  //double
  o[0] = (tana * a[0] + tanb * b[0] + tanc * c[0]) / d;
  o[1] = (tana * a[1] + tanb * b[1] + tanc * c[1]) / d;
  o[2] = (tana * a[2] + tanb * b[2] + tanc * c[2]) / d;
  return o;
}

//o = wNewOrthocenter (a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wNewOrthocenter (a, b, c) {
  let cosa = wCosAngleA (a, b, c);  //double
  let cosb = wCosAngleA (b, c, a);  //double
  let cosc = wCosAngleA (c, a, b);  //double
  let tana = sqrt (1 - cosa * cosa) / cosa;  //double
  let tanb = sqrt (1 - cosb * cosb) / cosb;  //double
  let tanc = sqrt (1 - cosc * cosc) / cosc;  //double
  let d = tana + tanb + tanc;  //double
  return [(tana * a[0] + tanb * b[0] + tanc * c[0]) / d,
          (tana * a[1] + tanb * b[1] + tanc * c[1]) / d,
          (tana * a[2] + tanb * b[2] + tanc * c[2]) / d];
}

//------------------------------------------------------------------------
//perpendicular foot 垂線の足
//  点aから直線bcに下ろした垂線の足d
//  平面
//    直線bc上の点dを通って直線bcと直交する直線上の点a
//    (dx,dy)=(bx+(cx-bx)*s,by+(cy-by)*s)
//    (cx-bx)*(ax-dx)+(cy-by)*(ay-dy)=0
//    (cx-bx)*(ax-(bx+(cx-bx)*s))+(cy-by)*(ay-(by+(cy-by)*s))=0
//    s=((ax-bx)*(cx-bx)+(ay-by)*(cy-by))/((cx-bx)^2+(cy-by)^2)
//  空間
//    直線bc上の点dを通って直線bcと直交する平面上の点a
//    (dx,dy,dz)=(bx+(cx-bx)*s,by+(cy-by)*s,bz+(cz-bz)*s)
//    (cx-bx)*(ax-dx)+(cy-by)*(ay-dy)+(cz-bz)*(az-dz)=0
//    (cx-bx)*(ax-(bx+(cx-bx)*s))+(cy-by)*(ay-(by+(cy-by)*s))+(cz-bz)*(az-(bz+(cz-bz)*s))=0
//    s=((ax-bx)*(cx-bx)+(ay-by)*(cy-by)+(az-bz)*(cz-bz))/((cx-bx)^2+(cy-by)^2+(cz-bz)^2)

//a = wPerpendicularFoot (a, b, c)
//  Array a
//  Array b
//  Array c
function wPerpendicularFoot (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  let s = (((ax - bx) * (cx - bx) + (ay - by) * (cy - by) + (az - bz) * (cz - bz)) /
           ((cx - bx) * (cx - bx) + (cy - by) * (cy - by) + (cz - bz) * (cz - bz)));  //double
  a[0] = bx + (cx - bx) * s;
  a[1] = by + (cy - by) * s;
  a[2] = bz + (cz - bz) * s;
  return a;
}

//o = wLetPerpendicularFoot (o, a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wLetPerpendicularFoot (o, a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  let s = (((ax - bx) * (cx - bx) + (ay - by) * (cy - by) + (az - bz) * (cz - bz)) /
           ((cx - bx) * (cx - bx) + (cy - by) * (cy - by) + (cz - bz) * (cz - bz)));  //double
  o[0] = bx + (cx - bx) * s;
  o[1] = by + (cy - by) * s;
  o[2] = bz + (cz - bz) * s;
  return o;
}

//o = wNewPerpendicularFoot (a, b, c)
//  Array o
//  Array a
//  Array b
//  Array c
function wNewPerpendicularFoot (a, b, c) {
  let ax = a[0], ay = a[1], az = a[2];  //double
  let bx = b[0], by = b[1], bz = b[2];  //double
  let cx = c[0], cy = c[1], cz = c[2];  //double
  let s = (((ax - bx) * (cx - bx) + (ay - by) * (cy - by) + (az - bz) * (cz - bz)) /
           ((cx - bx) * (cx - bx) + (cy - by) * (cy - by) + (cz - bz) * (cz - bz)));  //double
  return [bx + (cx - bx) * s,
          by + (cy - by) * s,
          bz + (cz - bz) * s];
}

//------------------------------------------------------------------------
//reverse 逆向きベクトル
//  -a

//a = wReverse (a)
//  Array a
function wReverse (a) {
  a[0] = -a[0];
  a[1] = -a[1];
  a[2] = -a[2];
  return a;
}

//o = wLetReverse (o, a)
//  Array o
//  Array a
function wLetReverse (o, a) {
  o[0] = -a[0];
  o[1] = -a[1];
  o[2] = -a[2];
  return o;
}

//o = wNewReverse (a)
//  Array o
//  Array a
function wNewReverse (a) {
  return [-a[0],
          -a[1],
          -a[2]];
}

//------------------------------------------------------------------------
//rotateXY 回転(+x→+y)
//
//    [cos(t) -sin(t) 0 0]   [x]   [cos(t)*x-sin(t)*y]
//    [sin(t)  cos(t) 0 0] * [y] = [sin(t)*x+cos(t)*y]
//    [  0       0    1 0]   [z]   [        z        ]
//    [  0       0    0 1]   [1]   [        1        ]
//

//a = wRotateXY (a, t)
//  Array a
//  double t
function wRotateXY (a, t) {
  let c = cos (t), s = sin (t);  //double
  let ax = a[0], ay = a[1];  //double
  a[0] = c * ax - s * ay;
  a[1] = s * ax + c * ay;
  return a;
}

//o = wLetRotateXY (o, a, t)
//  Array o
//  Array a
//  double t
function wLetRotateXY (o, a, t) {
  let c = cos (t), s = sin (t);  //double
  let ax = a[0], ay = a[1];  //double
  o[0] = c * ax - s * ay;
  o[1] = s * ax + c * ay;
  o[2] = a[2];
  return o;
}

//o = wNewRotateXY (a, t)
//  Array o