diff --git a/Tetrahedrons/MVec3d.cs b/Tetrahedrons/MVec3d.cs deleted file mode 100644 index 4191e8a..0000000 --- a/Tetrahedrons/MVec3d.cs +++ /dev/null @@ -1,271 +0,0 @@ -using System; -using OpenTK; - -namespace Tetrahedrons -{ - // ReSharper disable once InconsistentNaming - public struct MVec3d - { - public static readonly MVec3d Zero = new MVec3d(0, 0, 0, 0, 0, 0, 0, 0); - public static readonly MVec3d One = new MVec3d(1, 1, 1, 1, 1, 1, 1, 1); - - public static readonly MVec3d Unit = new MVec3d(1, 0, 0, 0, 0, 0, 0, 0); - public static readonly MVec3d Unit1 = new MVec3d(0, 1, 0, 0, 0, 0, 0, 0); - public static readonly MVec3d Unit2 = new MVec3d(0, 0, 1, 0, 0, 0, 0, 0); - public static readonly MVec3d Unit3 = new MVec3d(0, 0, 0, 1, 0, 0, 0, 0); - public static readonly MVec3d Unit23 = new MVec3d(0, 0, 0, 0, 1, 0, 0, 0); - public static readonly MVec3d Unit31 = new MVec3d(0, 0, 0, 0, 0, 1, 0, 0); - public static readonly MVec3d Unit12 = new MVec3d(0, 0, 0, 0, 0, 0, 1, 0); - public static readonly MVec3d Unit123 = new MVec3d(0, 0, 0, 0, 0, 0, 0, 1); - - public double E; - public double E1; - public double E2; - public double E3; - public double E23; - public double E31; - public double E12; - public double E123; - - public MVec3d(double e, double e1, double e2, double e3, double e23, double e31, double e12, double e123) - { - E = e; - E1 = e1; - E2 = e2; - E3 = e3; - E23 = e23; - E31 = e31; - E12 = e12; - E123 = e123; - } - - public MVec3d(MVec3d s, MVec3d v, MVec3d b, MVec3d t) - : this(s.E, v.E1, v.E2, v.E3, b.E23, b.E31, b.E12, t.E123) - { - } - - public static MVec3d Scalar(double e) => new MVec3d(e, 0, 0, 0, 0, 0, 0, 0); - public static MVec3d Vector(double e1, double e2, double e3) => new MVec3d(0, e1, e2, e3, 0, 0, 0, 0); - public static MVec3d Vector(Vector3d v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static MVec3d Bivector(double e23, double e31, double e12) => new MVec3d(0, 0, 0, 0, e23, e31, e12, 0); - public static MVec3d Trivector(double e123) => new MVec3d(0, 0, 0, 0, 0, 0, 0, e123); - public static MVec3d Complex(double real, double imag) => new MVec3d(real, 0, 0, 0, 0, 0, 0, imag); - public static MVec3d Complex(double w, double i, double j, double k) => new MVec3d(w, 0, 0, 0, i, j, k, 0); - public static MVec3d Complex(double angle, MVec3d plane) => Math.Cos(angle) + Math.Sin(angle) * plane; - - public static MVec3d Add(MVec3d m, MVec3d n) - { - return new MVec3d(m.E + n.E, - m.E1 + n.E1, m.E2 + n.E2, m.E3 + n.E3, - m.E23 + n.E23, m.E31 + n.E31, m.E12 + n.E12, - m.E123 + n.E123); - } - - public static MVec3d Sub(MVec3d m, MVec3d n) - { - return new MVec3d(m.E - n.E, - m.E1 - n.E1, m.E2 - n.E2, m.E3 - n.E3, - m.E23 - n.E23, m.E31 - n.E31, m.E12 - n.E12, - m.E123 - n.E123); - } - - public static MVec3d Mult(MVec3d m, MVec3d n) - { - var r = new MVec3d(); - - if (Math.Abs(m.E) > 1e-8) - { - r.E += m.E * n.E; - r.E1 += m.E * n.E1; - r.E2 += m.E * n.E2; - r.E3 += m.E * n.E3; - r.E23 += m.E * n.E23; - r.E31 += m.E * n.E31; - r.E12 += m.E * n.E12; - r.E123 += m.E * n.E123; - } - - if (Math.Abs(m.E1) > 1e-8) - { - r.E += +m.E1 * n.E1; - r.E1 += +m.E1 * n.E; - r.E2 += +m.E1 * n.E12; - r.E3 += -m.E1 * n.E31; - r.E23 += +m.E1 * n.E123; - r.E31 += -m.E1 * n.E3; - r.E12 += +m.E1 * n.E2; - r.E123 += +m.E1 * n.E23; - } - - if (Math.Abs(m.E2) > 1e-8) - { - r.E += +m.E2 * n.E2; - r.E1 += -m.E2 * n.E12; - r.E2 += +m.E2 * n.E; - r.E3 += +m.E2 * n.E23; - r.E23 += +m.E2 * n.E3; - r.E31 += +m.E2 * n.E123; - r.E12 += -m.E2 * n.E1; - r.E123 += +m.E2 * n.E31; - } - - if (Math.Abs(m.E3) > 1e-8) - { - r.E += +m.E3 * n.E3; - r.E1 += +m.E3 * n.E31; - r.E2 += -m.E3 * n.E23; - r.E3 += +m.E3 * n.E; - r.E23 += -m.E3 * n.E2; - r.E31 += +m.E3 * n.E1; - r.E12 += +m.E3 * n.E123; - r.E123 += +m.E3 * n.E12; - } - - if (Math.Abs(m.E23) > 1e-8) - { - r.E += -m.E23 * n.E23; - r.E1 += -m.E23 * n.E123; - r.E2 += +m.E23 * n.E3; - r.E3 += -m.E23 * n.E2; - r.E23 += +m.E23 * n.E; - r.E31 += +m.E23 * n.E12; - r.E12 += -m.E23 * n.E31; - r.E123 += +m.E23 * n.E1; - } - - if (Math.Abs(m.E31) > 1e-8) - { - r.E += -m.E31 * n.E31; - r.E1 += -m.E31 * n.E3; - r.E2 += -m.E31 * n.E123; - r.E3 += +m.E31 * n.E1; - r.E23 += -m.E31 * n.E12; - r.E31 += +m.E31 * n.E; - r.E12 += +m.E31 * n.E23; - r.E123 += +m.E31 * n.E2; - } - - if (Math.Abs(m.E12) > 1e-8) - { - r.E += -m.E12 * n.E12; - r.E1 += +m.E12 * n.E2; - r.E2 += -m.E12 * n.E1; - r.E3 += -m.E12 * n.E123; - r.E23 += +m.E12 * n.E31; - r.E31 += -m.E12 * n.E23; - r.E12 += +m.E12 * n.E; - r.E123 += +m.E12 * n.E3; - } - - if (Math.Abs(m.E123) > 1e-8) - { - r.E += -m.E123 * n.E123; - r.E1 += -m.E123 * n.E23; - r.E2 += -m.E123 * n.E31; - r.E3 += -m.E123 * n.E12; - r.E23 += +m.E123 * n.E1; - r.E31 += +m.E123 * n.E2; - r.E12 += +m.E123 * n.E3; - r.E123 += +m.E123 * n.E; - } - - return r; - } - - public double Norm2 => E * E + E1 * E1 + E2 * E2 + E3 * E3 + E23 * E23 + E31 * E31 + E12 * E12 + E123 * E123; - public double Norm => Math.Sqrt(Norm2); - - public MVec3d Inv - { - get - { - var n2 = Norm2; - return new MVec3d(E / n2, -E1 / n2, -E2 / n2, -E3 / n2, -E23 / n2, -E31 / n2, -E12 / n2, E123 / n2); -// return new MVec3d(E / n2, E1 / n2, E2 / n2, E3 / n2, E23 / n2, E31 / n2, E12 / n2, E123 / n2); - } - } - - public MVec3d Grade(int k) - { - switch (k) - { - case 0: - return Scalar(E); - case 1: - return Vector(E1, E2, E3); - case 2: - return Bivector(E23, E31, E12); - case 3: - return Trivector(E123); - default: - return new MVec3d(); - } - } - - public static MVec3d Inner(MVec3d m, MVec3d n) - { - var m0 = m.Grade(0); - var m1 = m.Grade(1); - var m2 = m.Grade(2); - var m3 = m.Grade(3); - - var n0 = n.Grade(0); - var n1 = n.Grade(1); - var n2 = n.Grade(2); - var n3 = n.Grade(3); - - return new MVec3d( - m0 * n0 + m1 * n1 + m2 * n2 + m3 * n3, - m0 * n1 + m1 * n2 + m2 * n3, - m0 * n2 + m1 * n3, - m0 * n3 - ); - } - - public static MVec3d Outer(MVec3d m, MVec3d n) - { - var m0 = m.Grade(0); - var m1 = m.Grade(1); - var m2 = m.Grade(2); - var m3 = m.Grade(3); - - var n0 = n.Grade(0); - var n1 = n.Grade(1); - var n2 = n.Grade(2); - var n3 = n.Grade(3); - - return new MVec3d( - m0 * n0, - m1 * n0 + m0 * n1, - m2 * n0 + m1 * n1 + m0 * n2, - m3 * n0 + m2 * n1 + m1 * n2 + m0 * m3 - ); - } - - public static MVec3d operator ~(MVec3d m) => m.Inv; - public static MVec3d operator *(MVec3d m, MVec3d n) => Mult(m, n); - - public static MVec3d operator +(MVec3d m, MVec3d n) => Add(m, n); - public static MVec3d operator -(MVec3d m, MVec3d n) => Sub(m, n); - - public static MVec3d operator &(MVec3d m, MVec3d n) => Inner(m, n); - public static MVec3d operator ^(MVec3d m, MVec3d n) => Outer(m, n); - - public MVec3d ScalarPart => Scalar(E); - public MVec3d VectorPart => Vector(E1, E2, E3); - public MVec3d BivectorPart => Bivector(E23, E31, E12); - public MVec3d TrivectorPart => Trivector(E123); - - public static implicit operator MVec3d(Vector3 v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static implicit operator MVec3d(Vector3d v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static implicit operator MVec3d(double s) => new MVec3d(s, 0, 0, 0, 0, 0, 0, 0); - - public static explicit operator double(MVec3d m) => m.E; - public static explicit operator Vector3d(MVec3d m) => new Vector3d(m.E1, m.E2, m.E3); - - public override string ToString() - { - return $"({E:0.00} {E1:0.00}e1 {E2:0.00}e2 {E3:0.00}e3 {E23:0.00}e23 {E31:0.00}e31 {E12:0.00}e12 {E123:0.00}e123)"; - } - } -} \ No newline at end of file diff --git a/Tetrahedrons/MVec4D.cs b/Tetrahedrons/MVec4D.cs new file mode 100644 index 0000000..9a0d8aa --- /dev/null +++ b/Tetrahedrons/MVec4D.cs @@ -0,0 +1,1121 @@ +using System; +using System.Collections.Generic; +using System.Data.SqlClient; +using System.Dynamic; +using System.Runtime.CompilerServices; +using System.Runtime.Remoting.Messaging; +using System.Text; +using OpenTK; +using OpenTK.Graphics.ES30; + +namespace Tetrahedrons +{ + public struct MVec4D + { + #region Units + + public static readonly MVec4D Zero = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D One = new MVec4D(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1); + + public static readonly MVec4D Unit = new MVec4D(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitX = new MVec4D(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitY = new MVec4D(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitZ = new MVec4D(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitW = new MVec4D(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXy = new MVec4D(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXz = new MVec4D(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitYz = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitYw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitZw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXyz = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0); + public static readonly MVec4D UnitXyw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0); + public static readonly MVec4D UnitXzw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0); + public static readonly MVec4D UnitYzw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0); + public static readonly MVec4D UnitXyzw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1); + + #endregion + + public double S, X, Y, Z, W, Xy, Xz, Xw, Yz, Yw, Zw, Xyz, Xyw, Xzw, Yzw, Xyzw; + + public double Norm2 => S * S + X * X + Y * Y + Z * Z + W * W + Xy * Xy + Xz * Xz + Xw * Xw + Yz * Yz + Yw * Yw + Zw * Zw + Xyz * Xyz + + Xyw * Xyw + Xzw * Xzw + Yzw * Yzw + Xyzw * Xyzw; + + public double Norm => Math.Sqrt(Norm2); + + public MVec4D Normalized => this * (1 / Norm); + + public MVec4D Grade(int k) + { + switch (k) + { + case 0: + return Vec0(S); + case 1: + return Vec1(X, Y, Z, W); + case 2: + return Vec2(Xy, Xz, Xw, Yz, Yw, Zw); + case 3: + return Vec3(Xyz, Xyw, Xzw, Yzw); + case 4: + return Vec4(Xyzw); + default: + return Zero; + } + } + + public Vector4d V4 => new Vector4d(X, Y, Z, W); + public Vector3d V3 => new Vector3d(X, Y, Z); + + #region Swizzle + + public double Yx + { + get => -Xy; + set => Xy = -value; + } + + public double Zx + { + get => -Xz; + set => Xz = -value; + } + + public double Wx + { + get => -Xw; + set => Xw = -value; + } + + public double Zy + { + get => -Yz; + set => Yz = -value; + } + + public double Wy + { + get => -Yw; + set => Yw = -value; + } + + public double Wz + { + get => -Zw; + set => Zw = -value; + } + + public double Xzy + { + get => -Xyz; + set => Xyz = -value; + } + + public double Yxz + { + get => -Xyz; + set => Xyz = -value; + } + + public double Yzx + { + get => Xyz; + set => Xyz = value; + } + + public double Zxy + { + get => Xyz; + set => Xyz = value; + } + + public double Zyx + { + get => -Xyz; + set => Xyz = -value; + } + + public double Xwy + { + get => -Xyw; + set => Xyw = -value; + } + + public double Yxw + { + get => -Xyw; + set => Xyw = -value; + } + + public double Ywx + { + get => Xyw; + set => Xyw = value; + } + + public double Wxy + { + get => Xyw; + set => Xyw = value; + } + + public double Wyx + { + get => -Xyw; + set => Xyw = -value; + } + + public double Xwz + { + get => -Xzw; + set => Xzw = -value; + } + + public double Zxw + { + get => -Xzw; + set => Xzw = -value; + } + + public double Zwx + { + get => Xzw; + set => Xzw = value; + } + + public double Wxz + { + get => Xzw; + set => Xzw = value; + } + + public double Wzx + { + get => -Xzw; + set => Xzw = -value; + } + + public double Ywz + { + get => -Yzw; + set => Yzw = -value; + } + + public double Zyw + { + get => -Yzw; + set => Yzw = -value; + } + + public double Zwy + { + get => Yzw; + set => Yzw = value; + } + + public double Wyz + { + get => Yzw; + set => Yzw = value; + } + + public double Wzy + { + get => -Yzw; + set => Yzw = -value; + } + + public double Xywz + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Xzyw + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Xzwy + { + get => Xyzw; + set => Xyzw = value; + } + + public double Xwyz + { + get => Xyzw; + set => Xyzw = value; + } + + public double Xwzy + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Yxzw + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Yxwz + { + get => Xyzw; + set => Xyzw = value; + } + + public double Yzxw + { + get => Xyzw; + set => Xyzw = value; + } + + public double Yzwx + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Ywxz + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Ywzx + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zxyw + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zxwy + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Zyxw + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Zywx + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zwxy + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zwyx + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wxyz + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wxzy + { + get => Xyzw; + set => Xyzw = value; + } + + public double Wyxz + { + get => Xyzw; + set => Xyzw = value; + } + + public double Wyzx + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wzxy + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wzyx + { + get => Xyzw; + set => Xyzw = value; + } + + #endregion + + #region Inverses + + public MVec4D Neg => new MVec4D(-S, -X, -Y, -Z, -W, -Xy, -Xz, -Xw, -Yz, -Yw, -Zw, -Xyz, -Xyw, -Xzw, -Yzw, -Xyzw); + + public MVec4D? Inv + { + get + { + var v = Mul(this, Rev); + if (v == Vec0(v.S)) + return Mul(1d / v.S, Rev); + + return null; + } + } + + public MVec4D Rev => new MVec4D(S, X, Y, Z, W, Yx, Zx, Wx, Zy, Wy, Wz, Zyx, Wyx, Wzx, Wzy, Wzyx); + + public MVec4D Conj => new MVec4D(Xyzw, -Yzw, Xzw, -Xyw, -Xyz, Zw, Yw, Yz, Xw, Xz, Xy, -W, Z, Y, X, S); + + public MVec4D Dual => new MVec4D(Xyzw, -Yzw, -Xzw, -Xyw, Xyz, -Zw, Yw, -Yz, -Xw, Xz, -Xy, W, -Z, Y, -Xz, S); + + #endregion + + #region Constructor + + public MVec4D(double s, double x, double y, double z, double w, double xy, double xz, double xw, double yz, double yw, double zw, double xyz, + double xyw, double xzw, double yzw, double xyzw) + { + S = s; + X = x; + Y = y; + Z = z; + W = w; + Xy = xy; + Xz = xz; + Xw = xw; + Yz = yz; + Yw = yw; + Zw = zw; + Xyz = xyz; + Xyw = xyw; + Xzw = xzw; + Yzw = yzw; + Xyzw = xyzw; + } + + public MVec4D(MVec4D s, MVec4D v, MVec4D b, MVec4D t, MVec4D q) + { + S = s.S; + X = v.X; + Y = v.Y; + Z = v.Z; + W = v.W; + Xy = b.Xy; + Xz = b.Xz; + Xw = b.Xw; + Yz = b.Yz; + Yw = b.Yw; + Zw = b.Zw; + Xyz = t.Xyz; + Xyw = t.Xyw; + Xzw = t.Xzw; + Yzw = t.Yzw; + Xyzw = q.Xyzw; + } + + #endregion + + #region Factory Methods + + public static MVec4D Vec0(double s) + { + return new MVec4D(s, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + } + + public static MVec4D Vec1(double x, double y, double z, double w) + { + return new MVec4D(0, x, y, z, w, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + } + + public static MVec4D Vec2(double xy, double xz, double xw, double yz, double yw, double zw) + { + return new MVec4D(0, 0, 0, 0, 0, xy, xz, xw, yz, yw, zw, 0, 0, 0, 0, 0); + } + + public static MVec4D Vec3(double xyz, double xyw, double xzw, double yzw) + { + return new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, xyz, xyw, xzw, yzw, 0); + } + + public static MVec4D Vec4(double xyzw) + { + return new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, xyzw); + } + + public static MVec4D Rotor(double angle, MVec4D plane) + { + return Add(Math.Cos(angle), Mul(Math.Sin(angle), plane)); + } + + #endregion + + #region Operations + + public static MVec4D Add(MVec4D m, MVec4D n) + { + return new MVec4D( + m.S + n.S, + m.X + n.X, + m.Y + n.Y, + m.Z + n.Z, + m.W + n.W, + m.Xy + n.Xy, + m.Xz + n.Xz, + m.Xw + n.Xw, + m.Yz + n.Yz, + m.Yw + n.Yw, + m.Zw + n.Zw, + m.Xyz + n.Xyz, + m.Xyw + n.Xyw, + m.Xzw + n.Xzw, + m.Yzw + n.Yzw, + m.Xyzw + n.Xyzw + ); + } + + public static MVec4D Add(MVec4D m, double c) => Add(c, m); + + public static MVec4D Add(double c, MVec4D m) + { + return new MVec4D( + m.S + c, + m.X, + m.Y, + m.Z, + m.W, + m.Xy, + m.Xz, + m.Xw, + m.Yz, + m.Yw, + m.Zw, + m.Xyz, + m.Xyw, + m.Xzw, + m.Yzw, + m.Xyzw + ); + } + + public static MVec4D Sub(MVec4D m, MVec4D n) + { + return new MVec4D( + m.S - n.S, + m.X - n.X, + m.Y - n.Y, + m.Z - n.Z, + m.W - n.W, + m.Xy - n.Xy, + m.Xz - n.Xz, + m.Xw - n.Xw, + m.Yz - n.Yz, + m.Yw - n.Yw, + m.Zw - n.Zw, + m.Xyz - n.Xyz, + m.Xyw - n.Xyw, + m.Xzw - n.Xzw, + m.Yzw - n.Yzw, + m.Xyzw - n.Xyzw + ); + } + + public static MVec4D Sub(MVec4D m, double c) + { + return new MVec4D( + m.S - c, + m.X, + m.Y, + m.Z, + m.W, + m.Xy, + m.Xz, + m.Xw, + m.Yz, + m.Yw, + m.Zw, + m.Xyz, + m.Xyw, + m.Xzw, + m.Yzw, + m.Xyzw + ); + } + + public static MVec4D Sub(double c, MVec4D m) + { + return new MVec4D( + c - m.S, + -m.X, + -m.Y, + -m.Z, + -m.W, + -m.Xy, + -m.Xz, + -m.Xw, + -m.Yz, + -m.Yw, + -m.Zw, + -m.Xyz, + -m.Xyw, + -m.Xzw, + -m.Yzw, + -m.Xyzw + ); + } + + public static MVec4D Mul(MVec4D m, MVec4D n) + { + var r = new MVec4D(); + if (Math.Abs(m.S) > 1e-10) + { + r.S += m.S * n.S; + r.X += m.S * n.X; + r.Y += m.S * n.Y; + r.Z += m.S * n.Z; + r.W += m.S * n.W; + r.Xy += m.S * n.Xy; + r.Xz += m.S * n.Xz; + r.Xw += m.S * n.Xw; + r.Yz += m.S * n.Yz; + r.Yw += m.S * n.Yw; + r.Zw += m.S * n.Zw; + r.Xyz += m.S * n.Xyz; + r.Xyw += m.S * n.Xyw; + r.Xzw += m.S * n.Xzw; + r.Yzw += m.S * n.Yzw; + r.Xyzw += m.S * n.Xyzw; + } + + if (Math.Abs(m.X) > 1e-10) + { + r.S += m.X * n.X; + r.X += m.X * n.S; + r.Y += m.X * n.Xy; + r.Z += m.X * n.Xz; + r.W += m.X * n.Xw; + r.Xy += m.X * n.Y; + r.Xz += m.X * n.Z; + r.Xw += m.X * n.W; + r.Yz += m.X * n.Xyz; + r.Yw += m.X * n.Xyw; + r.Zw += m.X * n.Xzw; + r.Xyz += m.X * n.Yz; + r.Xyw += m.X * n.Yw; + r.Xzw += m.X * n.Zw; + r.Yzw += m.X * n.Xyzw; + r.Xyzw += m.X * n.Yzw; + } + + if (Math.Abs(m.Y) > 1e-10) + { + r.S += m.Y * n.Y; + r.X += m.Y * n.Yx; + r.Y += m.Y * n.S; + r.Z += m.Y * n.Yz; + r.W += m.Y * n.Yw; + r.Yx += m.Y * n.X; + r.Xz += m.Y * n.Yxz; + r.Xw += m.Y * n.Yxw; + r.Yz += m.Y * n.Z; + r.Yw += m.Y * n.W; + r.Zw += m.Y * n.Yzw; + r.Yxz += m.Y * n.Xz; + r.Yxw += m.Y * n.Xw; + r.Xzw += m.Y * n.Yxzw; + r.Yzw += m.Y * n.Zw; + r.Yxzw += m.Y * n.Xzw; + } + + if (Math.Abs(m.Z) > 1e-10) + { + r.S += m.Z * n.Z; + r.X += m.Z * n.Zx; + r.Y += m.Z * n.Zy; + r.Z += m.Z * n.S; + r.W += m.Z * n.Zw; + r.Xy += m.Z * n.Zxy; + r.Zx += m.Z * n.X; + r.Xw += m.Z * n.Zxw; + r.Zy += m.Z * n.Y; + r.Yw += m.Z * n.Zyw; + r.Zw += m.Z * n.W; + r.Zxy += m.Z * n.Xy; + r.Xyw += m.Z * n.Zxyw; + r.Zxw += m.Z * n.Xw; + r.Zyw += m.Z * n.Yw; + r.Zxyw += m.Z * n.Xyw; + } + + if (Math.Abs(m.W) > 1e-10) + { + r.S += m.W * n.W; + r.X += m.W * n.Wx; + r.Y += m.W * n.Wy; + r.Z += m.W * n.Wz; + r.W += m.W * n.S; + r.Xy += m.W * n.Wxy; + r.Xz += m.W * n.Wxz; + r.Wx += m.W * n.X; + r.Yz += m.W * n.Wyz; + r.Wy += m.W * n.Y; + r.Wz += m.W * n.Z; + r.Xyz += m.W * n.Wxyz; + r.Wxy += m.W * n.Xy; + r.Wxz += m.W * n.Xz; + r.Wyz += m.W * n.Yz; + r.Wxyz += m.W * n.Xyz; + } + + if (Math.Abs(m.Xy) > 1e-10) + { + r.S += m.Xy * n.Yx; + r.X += m.Xy * n.Y; + r.Y += m.Yx * n.X; + r.Z += m.Xy * n.Yxz; + r.W += m.Xy * n.Yxw; + r.Xy += m.Xy * n.S; + r.Xz += m.Xy * n.Yz; + r.Xw += m.Xy * n.Yw; + r.Yz += m.Yx * n.Xz; + r.Yw += m.Yx * n.Xw; + r.Zw += m.Xy * n.Yxzw; + r.Xyz += m.Xy * n.Z; + r.Xyw += m.Xy * n.W; + r.Xzw += m.Xy * n.Yzw; + r.Yzw += m.Yx * n.Xzw; + r.Xyzw += m.Xy * n.Zw; + } + + if (Math.Abs(m.Xz) > 1e-10) + { + r.S += m.Xz * n.Zx; + r.X += m.Xz * n.Z; + r.Y += m.Xz * n.Zxy; + r.Z += m.Zx * n.X; + r.W += m.Xz * n.Zxw; + r.Xy += m.Xz * n.Zy; + r.Xz += m.Xz * n.S; + r.Xw += m.Xz * n.Zw; + r.Zy += m.Zx * n.Xy; + r.Yw += m.Xz * n.Zxyw; + r.Zw += m.Zx * n.Xw; + r.Xzy += m.Xz * n.Y; + r.Xyw += m.Xz * n.Zyw; + r.Xzw += m.Xz * n.W; + r.Zyw += m.Zx * n.Xyw; + r.Xzyw += m.Xz * n.Yw; + } + + if (Math.Abs(m.Xw) > 1e-10) + { + r.S += m.Xw * n.Wx; + r.X += m.Xw * n.W; + r.Y += m.Xw * n.Wxy; + r.Z += m.Xw * n.Wxz; + r.W += m.Wx * n.X; + r.Xy += m.Xw * n.Wy; + r.Xz += m.Xw * n.Wz; + r.Xw += m.Xw * n.S; + r.Yz += m.Xw * n.Wxyz; + r.Wy += m.Wx * n.Xy; + r.Wz += m.Wx * n.Xz; + r.Xyz += m.Xw * n.Wyz; + r.Xwy += m.Xw * n.Y; + r.Xwz += m.Xw * n.Z; + r.Wyz += m.Wx * n.Xyz; + r.Xwyz += m.Xw * n.Yz; + } + + if (Math.Abs(m.Yz) > 1e-10) + { + r.S += m.Yz * n.Zy; + r.X += m.Yz * n.Zyx; + r.Y += m.Yz * n.Z; + r.Z += m.Zy * n.Y; + r.W += m.Yz * n.Zyw; + r.Yx += m.Yz * n.Zx; + r.Zx += m.Zy * n.Yx; + r.Xw += m.Yz * n.Zyxw; + r.Yz += m.Yz * n.S; + r.Yw += m.Yz * n.Zw; + r.Zw += m.Zy * n.Yw; + r.Yzx += m.Yz * n.X; + r.Yxw += m.Yz * n.Zxw; + r.Zxw += m.Zy * n.Yxw; + r.Yzw += m.Yz * n.W; + r.Yzxw += m.Yz * n.Xw; + } + + if (Math.Abs(m.Yw) > 1e-10) + { + r.S += m.Yw * n.Wy; + r.X += m.Yw * n.Wyx; + r.Y += m.Yw * n.W; + r.Z += m.Yw * n.Wyz; + r.W += m.Wy * n.Y; + r.Yx += m.Yw * n.Wx; + r.Xz += m.Yw * n.Wyxz; + r.Wx += m.Wy * n.Yx; + r.Yz += m.Yw * n.Wz; + r.Yw += m.Yw * n.S; + r.Wz += m.Wy * n.Yz; + r.Yxz += m.Yw * n.Wxz; + r.Ywx += m.Yw * n.X; + r.Wxz += m.Wy * n.Yxz; + r.Ywz += m.Yw * n.Z; + r.Ywxz += m.Yw * n.Xz; + } + + if (Math.Abs(m.Zw) > 1e-10) + { + r.S += m.Zw * n.Wz; + r.X += m.Zw * n.Wzx; + r.Y += m.Zw * n.Wzy; + r.Z += m.Zw * n.W; + r.W += m.Wz * n.Z; + r.Xy += m.Zw * n.Wzxy; + r.Zx += m.Zw * n.Wx; + r.Wx += m.Wz * n.Zx; + r.Zy += m.Zw * n.Wy; + r.Wy += m.Wz * n.Zy; + r.Zw += m.Zw * n.S; + r.Zxy += m.Zw * n.Wxy; + r.Wxy += m.Wz * n.Zxy; + r.Zwx += m.Zw * n.X; + r.Zwy += m.Zw * n.Y; + r.Zwxy += m.Zw * n.Xy; + } + + if (Math.Abs(m.Xyz) > 1e-10) + { + r.S += m.Xyz * n.Zyx; + r.X += m.Xyz * n.Zy; + r.Y += m.Yxz * n.Zx; + r.Z += m.Zxy * n.Yx; + r.W += m.Xyz * n.Zyxw; + r.Xy += m.Xyz * n.Z; + r.Xz += m.Xzy * n.Y; + r.Xw += m.Xyz * n.Zyw; + r.Yz += m.Yzx * n.X; + r.Yw += m.Yxz * n.Zxw; + r.Zw += m.Zxy * n.Yxw; + r.Xyz += m.Xyz * n.S; + r.Xyw += m.Xyz * n.Zw; + r.Xzw += m.Xzy * n.Yw; + r.Yzw += m.Yzx * n.Xw; + r.Xyzw += m.Xyz * n.W; + } + + if (Math.Abs(m.Xyw) > 1e-10) + { + r.S += m.Xyw * n.Wyx; + r.X += m.Xyw * n.Wy; + r.Y += m.Yxw * n.Wx; + r.Z += m.Xyw * n.Wyxz; + r.W += m.Wxy * n.Yx; + r.Xy += m.Xyw * n.W; + r.Xz += m.Xyw * n.Wyz; + r.Xw += m.Xwy * n.Y; + r.Yz += m.Yxw * n.Wxz; + r.Yw += m.Ywx * n.X; + r.Wz += m.Wxy * n.Yxz; + r.Xyz += m.Xyw * n.Wz; + r.Xyw += m.Xyw * n.S; + r.Xwz += m.Xwy * n.Yz; + r.Ywz += m.Ywx * n.Xz; + r.Xywz += m.Xyw * n.Z; + } + + if (Math.Abs(m.Xzw) > 1e-10) + { + r.S += m.Xzw * n.Wzx; + r.X += m.Xzw * n.Wz; + r.Y += m.Xzw * n.Wzxy; + r.Z += m.Zxw * n.Wx; + r.W += m.Wxz * n.Zx; + r.Xy += m.Xzw * n.Wzy; + r.Xz += m.Xzw * n.W; + r.Xw += m.Xwz * n.Z; + r.Zy += m.Zxw * n.Wxy; + r.Wy += m.Wxz * n.Zxy; + r.Zw += m.Zwx * n.X; + r.Xzy += m.Xzw * n.Wy; + r.Xwy += m.Xwz * n.Zy; + r.Xzw += m.Xzw * n.S; + r.Zwy += m.Zwx * n.Xy; + r.Xzwy += m.Xzw * n.Y; + } + + if (Math.Abs(m.Yzw) > 1e-10) + { + r.S += m.Yzw * n.Wzy; + r.X += m.Yzw * n.Wzyx; + r.Y += m.Yzw * n.Wz; + r.Z += m.Zyw * n.Wy; + r.W += m.Wyz * n.Zy; + r.Yx += m.Yzw * n.Wzx; + r.Zx += m.Zyw * n.Wyx; + r.Wx += m.Wyz * n.Zyx; + r.Yz += m.Yzw * n.W; + r.Yw += m.Ywz * n.Z; + r.Zw += m.Zwy * n.Y; + r.Yzx += m.Yzw * n.Wx; + r.Ywx += m.Ywz * n.Zx; + r.Zwx += m.Zwy * n.Yx; + r.Yzw += m.Yzw * n.S; + r.Yzwx += m.Yzw * n.X; + } + + if (Math.Abs(m.Xyzw) > 1e-10) + { + r.S += m.Xyzw * n.Wzyx; + r.X += m.Xyzw * n.Wzy; + r.Y += m.Yxzw * n.Wzx; + r.Z += m.Zxyw * n.Wyx; + r.W += m.Wxyz * n.Zyx; + r.Xy += m.Xyzw * n.Wz; + r.Xz += m.Xzyw * n.Wy; + r.Xw += m.Xwyz * n.Zy; + r.Yz += m.Yzxw * n.Wx; + r.Yw += m.Ywxz * n.Zx; + r.Zw += m.Zwxy * n.Yx; + r.Xyz += m.Xyzw * n.W; + r.Xyw += m.Xywz * n.Z; + r.Xzw += m.Xzwy * n.Y; + r.Yzw += m.Yzwx * n.X; + r.Xyzw += m.Xyzw * n.S; + } + + return r; + } + + public static MVec4D Mul(MVec4D m, double c) => Mul(c, m); + + public static MVec4D Mul(double c, MVec4D m) + { + return new MVec4D( + c * m.S, + c * m.X, + c * m.Y, + c * m.Z, + c * m.W, + c * m.Xy, + c * m.Xz, + c * m.Xw, + c * m.Yz, + c * m.Yw, + c * m.Zw, + c * m.Xyz, + c * m.Xyw, + c * m.Xzw, + c * m.Yzw, + c * m.Xyzw + ); + } + + public static MVec4D Div(MVec4D m, double c) => Mul(m, 1 / c); + + public static MVec4D Inner(MVec4D m, MVec4D n) + { + var m0 = m.Grade(0); + var m1 = m.Grade(1); + var m2 = m.Grade(2); + var m3 = m.Grade(3); + var m4 = m.Grade(4); + + var n0 = n.Grade(0); + var n1 = n.Grade(1); + var n2 = n.Grade(2); + var n3 = n.Grade(3); + var n4 = n.Grade(4); + + return new MVec4D( + m0 * n0 + m1 * n1 + m2 * n2 + m3 * n3 + n4 * n4, + m0 * n1 + m1 * n2 + m2 * n3 + m3 * n4, + m0 * n2 + m1 * n3 + m2 * n4, + m0 * n3 + m1 * n4, + m0 * n4 + ); + } + + public static MVec4D Inner(MVec4D m, double c) => Mul(m, c); + + public static MVec4D Inner(double c, MVec4D n) => Mul(c, n); + + public static MVec4D Outer(MVec4D m, MVec4D n) + { + var m0 = m.Grade(0); + var m1 = m.Grade(1); + var m2 = m.Grade(2); + var m3 = m.Grade(3); + var m4 = m.Grade(4); + + var n0 = n.Grade(0); + var n1 = n.Grade(1); + var n2 = n.Grade(2); + var n3 = n.Grade(3); + var n4 = n.Grade(4); + + return new MVec4D( + m0 * n0, + m1 * n0 + m0 * n1, + m2 * n0 + m1 * n1 + m0 * n2, + m3 * n0 + m2 * n1 + m1 * n2 + m0 * n3, + m4 * n0 + m3 * n1 + m2 * n2 + m1 * n3 + m0 * n4 + ); + } + + public static MVec4D Outer(MVec4D m, double c) => Mul(m, c); + + public static MVec4D Outer(double c, MVec4D n) => Mul(c, n); + + public static MVec4D Transform(MVec4D m, MVec4D t) => Mul(Mul(t, m), t.Rev); + + public static MVec4D? Add(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; + public static MVec4D? Add(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; + public static MVec4D? Add(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; + + public static MVec4D? Sub(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Sub(m.Value, n.Value) : null; + public static MVec4D? Sub(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Sub(m.Value, n.Value) : null; + public static MVec4D? Sub(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Sub(m.Value, n.Value) : null; + + public static MVec4D? Mul(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Mul(m.Value, n.Value) : null; + public static MVec4D? Mul(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Mul(m.Value, n.Value) : null; + public static MVec4D? Mul(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Mul(m.Value, n.Value) : null; + + public static MVec4D? Div(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Div(m.Value, n.Value) : null; + public static MVec4D? Div(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Div(m.Value, n.Value) : null; + public static MVec4D? Div(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Div(m.Value, n.Value) : null; + + public static MVec4D? Inner(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Inner(m.Value, n.Value) : null; + public static MVec4D? Inner(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Inner(m.Value, n.Value) : null; + public static MVec4D? Inner(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Inner(m.Value, n.Value) : null; + + public static MVec4D? Outer(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + public static MVec4D? Outer(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + public static MVec4D? Outer(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + + public static MVec4D? Transform(MVec4D? m, MVec4D? t) => Mul(Mul(t, m), t?.Rev); + + #endregion + + #region Operators + + public static MVec4D operator +(MVec4D m, MVec4D n) => Add(m, n); + public static MVec4D operator +(double c, MVec4D n) => Add(c, n); + public static MVec4D operator +(MVec4D m, double c) => Add(m, c); + + public static MVec4D operator -(MVec4D m, MVec4D n) => Sub(m, n); + public static MVec4D operator -(double c, MVec4D n) => Sub(c, n); + public static MVec4D operator -(MVec4D m, double c) => Sub(m, c); + + public static MVec4D operator *(MVec4D m, MVec4D n) => Mul(m, n); + public static MVec4D operator *(double c, MVec4D n) => Mul(c, n); + public static MVec4D operator *(MVec4D m, double c) => Mul(m, c); + + public static MVec4D operator &(MVec4D m, MVec4D n) => Inner(m, n); + public static MVec4D operator &(MVec4D m, double c) => Inner(m, c); + public static MVec4D operator &(double c, MVec4D n) => Inner(c, n); + + public static MVec4D operator ^(MVec4D m, MVec4D n) => Outer(m, n); + public static MVec4D operator ^(MVec4D m, double c) => Outer(m, c); + public static MVec4D operator ^(double c, MVec4D n) => Outer(c, n); + + public static MVec4D operator |(MVec4D m, MVec4D t) => Transform(m, t); + + + public static MVec4D? operator +(MVec4D? m, MVec4D? n) => Add(m, n); + public static MVec4D? operator +(double? c, MVec4D? n) => Add(c, n); + public static MVec4D? operator +(MVec4D? m, double? c) => Add(m, c); + + public static MVec4D? operator -(MVec4D? m, MVec4D? n) => Sub(m, n); + public static MVec4D? operator -(double? c, MVec4D? n) => Sub(c, n); + public static MVec4D? operator -(MVec4D? m, double? c) => Sub(m, c); + + public static MVec4D? operator *(MVec4D? m, MVec4D? n) => Mul(m, n); + public static MVec4D? operator *(double? c, MVec4D? n) => Mul(c, n); + public static MVec4D? operator *(MVec4D? m, double? c) => Mul(m, c); + + public static MVec4D? operator &(MVec4D? m, MVec4D? n) => Inner(m, n); + public static MVec4D? operator &(MVec4D? m, double? c) => Inner(m, c); + public static MVec4D? operator &(double? c, MVec4D? n) => Inner(c, n); + + public static MVec4D? operator ^(MVec4D? m, MVec4D? n) => Outer(m, n); + public static MVec4D? operator ^(MVec4D? m, double? c) => Outer(m, c); + public static MVec4D? operator ^(double? c, MVec4D? n) => Outer(c, n); + + public static MVec4D? operator |(MVec4D? m, MVec4D? t) => Transform(m, t); + + + public static MVec4D operator -(MVec4D m) => m.Neg; + public static MVec4D? operator ~(MVec4D m) => m.Inv; + public static MVec4D operator !(MVec4D m) => m.Rev; + + public static implicit operator MVec4D(double v) => Vec0(v); + public static implicit operator MVec4D(Vector3d v) => Vec1(v.X, v.Y, v.Z, 0); + public static implicit operator MVec4D(Vector4d v) => Vec1(v.X, v.Y, v.Z, v.W); + + #endregion + + #region Equality + + public static bool operator ==(MVec4D m, MVec4D n) + { + var diff = Sub(m, n).Norm2; + return Math.Abs(diff) <= 10e-8; + } + + public static bool operator !=(MVec4D m, MVec4D n) + { + var diff = Sub(m, n).Norm2; + return Math.Abs(diff) > 10e-8; + } + + #endregion + + #region Formatting + + public override string ToString() + { + var sb = new List(); + if (Math.Abs(S) > 10e-2) sb.Add($"{S:F2}"); + if (Math.Abs(X) > 10e-2) sb.Add($"{X:F2}e1"); + if (Math.Abs(Y) > 10e-2) sb.Add($"{Y:F2}e2"); + if (Math.Abs(Z) > 10e-2) sb.Add($"{Z:F2}e3"); + if (Math.Abs(W) > 10e-2) sb.Add($"{W:F2}e4"); + if (Math.Abs(Xy) > 10e-2) sb.Add($"{Xy:F2}e12"); + if (Math.Abs(Yz) > 10e-2) sb.Add($"{Yz:F2}e23"); + if (Math.Abs(Zw) > 10e-2) sb.Add($"{Zw:F2}e34"); + if (Math.Abs(Wx) > 10e-2) sb.Add($"{Wx:F2}e41"); + if (Math.Abs(Xz) > 10e-2) sb.Add($"{Xz:F2}e13"); + if (Math.Abs(Yw) > 10e-2) sb.Add($"{Yw:F2}e24"); + if (Math.Abs(Xyz) > 10e-2) sb.Add($"{Xyz:F2}e123"); + if (Math.Abs(Yzw) > 10e-2) sb.Add($"{Yzw:F2}e234"); + if (Math.Abs(Zwx) > 10e-2) sb.Add($"{Zwx:F2}e341"); + if (Math.Abs(Wxy) > 10e-2) sb.Add($"{Wxy:F2}e412"); + if (Math.Abs(Xyzw) > 10e-2) sb.Add($"{Xyzw:F2}e1234"); + return $"({string.Join(" + ", sb)})"; + } + + #endregion + } +} \ No newline at end of file diff --git a/Tetrahedrons/Program.cs b/Tetrahedrons/Program.cs index 6175764..e2f16f6 100644 --- a/Tetrahedrons/Program.cs +++ b/Tetrahedrons/Program.cs @@ -1,13 +1,14 @@ -using OpenTK; +using System; +using OpenTK; namespace Tetrahedrons { - internal class Program : GameWindow - { - public static void Main(string[] args) - { - using(var p = new TetrahedronWindow()) - p.Run(); - } - } + internal class Program : GameWindow + { + public static void Main(string[] args) + { + using (var p = new TetrahedronWindow()) + p.Run(); + } + } } \ No newline at end of file diff --git a/Tetrahedrons/Simplex.cs b/Tetrahedrons/Simplex.cs new file mode 100644 index 0000000..13d9cc1 --- /dev/null +++ b/Tetrahedrons/Simplex.cs @@ -0,0 +1,148 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Net.Mail; + +namespace Tetrahedrons +{ + public class Simplex + { + public MVec4D[] Verts; + + public int[] Faces; + public int[] Edjes; + + public Simplex(IEnumerable verts) + : this(verts.ToArray()) + { + } + + public Simplex(params MVec4D[] verts) + { + Verts = verts.ToArray(); + + var e = new List(); + for (var i = 0; i < Verts.Length; i++) + for (var j = i + 1; j < Verts.Length; j++) + { + e.Add(i); + e.Add(j); + } + + Edjes = e.ToArray(); + + Faces = new int[0]; + } + + public Simplex(IEnumerable verts, IEnumerable faces) + : this(verts) + { + Faces = faces.ToArray(); + } + + public Simplex(IEnumerable verts, IEnumerable faces, IEnumerable edges) + { + Verts = verts.ToArray(); + Faces = faces.ToArray(); + Edjes = edges.ToArray(); + } + + public bool[] Sides(MVec4D blade, MVec4D pivot) + { + var sides = new bool[Verts.Length]; + + var perp = blade.Dual; + + for (var i = 0; i < Verts.Length; i++) + { + var m = Verts[i]; + var p = m - pivot; + var rej = (p ^ blade) * ~blade; + if (!rej.HasValue) continue; + var dot = rej & perp; + + sides[i] = dot.Value.S > 0; + } + + return sides; + } + + public Simplex Intersect(MVec4D blade, MVec4D pivot) + { + var sides = Sides(blade, pivot); + var verts = new List(); + var edges = new List(); + + var face_verts = new List[Verts.Length]; + for (var i = 0; i < face_verts.Length; i++) + face_verts[i] = new List(); + + { + var k = 0; + for (var i = 0; i < Verts.Length; i++) + for (var j = i + 1; j < Verts.Length; j++) + { + if (sides[i] == sides[j]) continue; + var m = Verts[i] - pivot; + var n = Verts[j] - pivot; + + var v = m - n; + var alph = (blade ^ m) * ~(blade ^ v); + if (!alph.HasValue) continue; + var p = m - v * alph.Value; + + verts.Add(p); + face_verts[i].Add(k); + face_verts[j].Add(k); + + k++; + } + } + + var faces = new List(); + foreach (var vlst in face_verts) + { + for (var i = 0; i < vlst.Count - 2; i++) + for (var j = i + 1; j < vlst.Count - 1; j++) + for (var k = j + 1; k < vlst.Count; k++) + { + faces.Add(vlst[i]); + faces.Add(vlst[j]); + faces.Add(vlst[k]); + + edges.Add(vlst[i]); + edges.Add(vlst[j]); + edges.Add(vlst[j]); + edges.Add(vlst[k]); + edges.Add(vlst[k]); + edges.Add(vlst[i]); + } + } + + if (faces.Count == 6) + { + for (var i = 0; i < 3; i++) + { + var a = faces[i]; + var b = faces[(1 + i) % 3]; + var c = faces[3 + i]; + var d = faces[3 + (1 + i) % 3]; + + faces.Add(c); + faces.Add(a); + faces.Add(b); + faces.Add(b); + faces.Add(d); + faces.Add(c); + + edges.Add(a); + edges.Add(c); + edges.Add(b); + edges.Add(d); + } + } + + return new Simplex(verts, faces, edges); + } + } +} \ No newline at end of file diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index 31ea6a8..c772f40 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -1,332 +1,153 @@ -using System; -using System.Collections.Generic; +using System; +using System.Diagnostics.PerformanceData; using System.Drawing; +using System.Drawing.Printing; +using System.Net; using System.Net.Mail; +using System.Runtime.InteropServices.ComTypes; using System.Security; +using System.Xml.XPath; using OpenTK; using OpenTK.Audio.OpenAL; using OpenTK.Graphics.OpenGL; using OpenTK.Input; +using OpenTK.Platform.Windows; namespace Tetrahedrons { - public struct Polygon - { - public Vector3d[] Points; - public Color Color; + public class TetrahedronWindow : GameWindow + { + private Matrix4 _proj3d; + private Matrix4 _view; - public Polygon(Vector3d[] points) : this(points, Color.DodgerBlue) - { - } + private MVec4D _a = MVec4D.UnitX; + private MVec4D _b = MVec4D.UnitY; - public Polygon(Vector3d[] points, Color color) - { - Points = points; - Color = color; - } + private Simplex[] _simplexes; + private Simplex[] _intersections; + + private double _t; + private bool _pause; + private bool _wire; + + protected override void OnLoad(EventArgs e) + { + base.OnLoad(e); + + X = (DisplayDevice.Default.Width - Width) / 2; + Y = (DisplayDevice.Default.Height - Height) / 2; + + _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(.86f, .5f, 1), Vector3.UnitZ); + + const int n = 5; + _simplexes = new Simplex[n * n * n * n]; + _intersections = new Simplex[_simplexes.Length]; + + var off = (new Vector4d(n - 1, n - 1, n - 1, n - 1)) / 2; + int i = 0; + for (var x = 0; x < n; x++) + for (var y = 0; y < n; y++) + for (var z = 0; z < n; z++) + for (var w = 0; w < n; w++) + { + var d = new Vector4d(x, y, z, w); + _simplexes[i] = new Simplex( // not-quite-regular pentatope + new Vector4d(+.5, -.5, -.5, -.5) + d - off, + new Vector4d(-.5, +.5, -.5, -.5) + d - off, + new Vector4d(-.5, -.5, +.5, -.5) + d - off, + new Vector4d(-.5, -.5, -.5, +.5) + d - off, + new Vector4d(+.5, +.5, +.5, +.5) + d - off); + _intersections[i] = new Simplex(); + i++; + } + } + + protected override void OnRenderFrame(FrameEventArgs e) + { + base.OnRenderFrame(e); + + GL.Viewport(ClientRectangle); + + GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit); + GL.PointSize(10f); + GL.LineWidth(3f); + + GL.PushMatrix(); + + GL.MatrixMode(MatrixMode.Projection); + GL.LoadMatrix(ref _proj3d); + GL.MatrixMode(MatrixMode.Modelview); + GL.LoadMatrix(ref _view); + + GL.Enable(EnableCap.DepthTest); + + GL.Begin(PrimitiveType.Lines); + GL.Color3(0f, 0, 0); + foreach (var intr in _intersections) + foreach (var f in intr.Edjes) + Util.Vertex3(intr.Verts[f]); + GL.End(); + + GL.Begin(PrimitiveType.Triangles); + foreach (var intr in _intersections) + foreach (var f in intr.Faces) + { + Util.Color3(intr.Verts[f]); + Util.Vertex3(intr.Verts[f]); + } + + GL.End(); + GL.Disable(EnableCap.DepthTest); + + if (_wire) + { + GL.Begin(PrimitiveType.Lines); + foreach (var sim in _simplexes) + foreach (var f in sim.Edjes) + Util.Vertex4(sim.Verts[f]); - public void Draw() - { - GL.Begin(PrimitiveType.TriangleStrip); - GL.Color3(Color); - foreach (var v in Points) - GL.Vertex3(v); GL.End(); - } + } - public static Polygon operator +(Polygon p, MVec3d m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (Vector3d) (p.Points[i] + m); - return new Polygon(pts, p.Color); - } + SwapBuffers(); + } - public static Polygon operator -(Polygon p, MVec3d m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (Vector3d) (p.Points[i] - m); - return new Polygon(pts, p.Color); - } + protected override void OnUpdateFrame(FrameEventArgs e) + { + base.OnUpdateFrame(e); - public static Polygon operator *(Polygon p, MVec3d m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (Vector3d) (p.Points[i] * m); - return new Polygon(pts, p.Color); - } + var w = 10f; + var d = w; + _proj3d = Matrix4.CreateOrthographic(w, w * Height / Width, -d / 2, d / 2); - public static Polygon operator ^(Polygon p, MVec3d m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (Vector3d) (p.Points[i] ^ m); - return new Polygon(pts, p.Color); - } + if (_pause) return; - public static Polygon operator &(Polygon p, MVec3d m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (Vector3d) (p.Points[i] & m); - return new Polygon(pts, p.Color); - } - } + _t += e.Time; - public class Tetrahedron - { - public readonly Vector3d[] Points; - public Color FillColor = Color.DodgerBlue; - public Color SurfaceColor = Color.GhostWhite; - public Color WireColor = Color.Maroon; + var pln = MVec4D.UnitXy + MVec4D.UnitXz + MVec4D.UnitYw; + var r = MVec4D.Rotor(e.Time / 10, pln.Normalized); + foreach (var pent in _simplexes) + for (var i = 0; i < pent.Verts.Length; i++) + pent.Verts[i] |= r; - public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4) - { - Points = new[] {p1, p2, p3, p4}; - } + var blade = MVec4D.UnitXyz; + var pivot = MVec4D.Zero; - public void Draw(bool wire = true) - { - if (wire) - { - GL.Begin(PrimitiveType.Lines); - GL.Color3(WireColor); + for (var i = 0; i < _simplexes.Length; i++) + { + _intersections[i] = _simplexes[i].Intersect(blade, pivot); + } + } - for (var i = 0; i < 4; i++) - for (var j = i + 1; j < 4; j++) - { - GL.Vertex3(Points[i]); - GL.Vertex3(Points[j]); - } + protected override void OnKeyDown(KeyboardKeyEventArgs e) + { + base.OnKeyDown(e); - GL.End(); - } - else - { - GL.Begin(PrimitiveType.Triangles); - GL.Color3(SurfaceColor); + if (e.Key == Key.Space) + _pause = !_pause; - for (var i = 0; i < 4; i++) - for (var j = i + 1; j < 4; j++) - for (var k = j + 1; k < 4; k++) - { - GL.Vertex3(Points[i]); - GL.Vertex3(Points[j]); - GL.Vertex3(Points[k]); - } - - GL.End(); - } - } - - public Polygon Intersection(MVec3d blade, MVec3d pivot) - { - var pts = new Vector3d[4]; - for (var i = 0; i < 4; i++) - { - pts[i] = (Vector3d) (Points[i] - pivot); - } - - var a = blade * MVec3d.Unit123; - var sides = new double[4]; - for (var i = 0; i < 4; i++) - { - sides[i] = (a & ((pts[i] ^ blade) * ~blade)).E; // normal dot rejection - } - - var vecs = new List(4); - - for (var i = 0; i < 4; i++) - { - for (var j = i + 1; j < 4; j++) - { - if (sides[i] * sides[j] >= 0) continue; // both pts on same side - - var t = ((blade ^ pts[i]) * ~(blade ^ (pts[j] - pts[i]))).E; - var point = pts[i] + t * (pts[j] - pts[i]); - vecs.Add((Vector3d) (point + pivot)); - } - } - - return new Polygon(vecs.ToArray()) {Color = FillColor}; - } - } - - public class TetrahedronWindow : GameWindow - { - private Matrix4 _proj3D; - private Matrix4 _proj2D; - private Matrix4 _view; - - private double _t; - private bool _pause; - private bool _div; - private Tetrahedron _tetra; - private Polygon _polyg; - - private MVec3d _b1 = MVec3d.Unit1; - private MVec3d _b2 = MVec3d.Unit3; - private MVec3d _pivot = MVec3d.Zero; - - private Polygon _square = new Polygon(new[] {new Vector3d(-1, -1, 0), new Vector3d(-1, 1, 0), new Vector3d(1, -1, 0), new Vector3d(1, 1, 0),}, Color.Red); - - protected override void OnLoad(EventArgs e) - { - base.OnLoad(e); - - X = (DisplayDevice.Default.Width - Width) / 2; - Y = (DisplayDevice.Default.Height - Height) / 2; - - _tetra = new Tetrahedron(new Vector3d(1, 1, 1), new Vector3d(-1, -1, 1), new Vector3d(-1, 1, -1), new Vector3d(1, -1, -1)); - _polyg = new Polygon(new[] {new Vector3d(-1, -1, 0), new Vector3d(-1, 1, 0), new Vector3d(1, -1, 0), new Vector3d(1, 1, 0)}); - - _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(1, .6f, 1f), Vector3.UnitZ); - } - - protected override void OnRenderFrame(FrameEventArgs e) - { - base.OnRenderFrame(e); - - GL.Viewport(ClientRectangle); - - GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit); - GL.PointSize(10f); - - GL.PushMatrix(); - - if (_div) - { - GL.MatrixMode(MatrixMode.Projection); - GL.LoadMatrix(ref _proj2D); - - _polyg.Draw(); - } - else - { - GL.MatrixMode(MatrixMode.Projection); - GL.LoadMatrix(ref _proj3D); - GL.MatrixMode(MatrixMode.Modelview); - GL.LoadMatrix(ref _view); - - - GL.Enable(EnableCap.DepthTest); - _tetra.Draw(wire: false); - GL.Disable(EnableCap.DepthTest); - _polyg.Draw(); - - GL.Begin(PrimitiveType.Points); - GL.Color3(Color.DarkBlue); - GL.Vertex3((Vector3d) _pivot); - GL.End(); - - GL.Begin(PrimitiveType.Lines); - GL.Vertex3((Vector3d) _pivot); - GL.Vertex3((Vector3d) (_pivot + (_b2 ^ _b1) * MVec3d.Unit123 * .5)); - GL.End(); - - GL.Enable(EnableCap.DepthTest); - _tetra.Draw(wire: true); - } - - GL.PopMatrix(); - - SwapBuffers(); - } - - protected override void OnUpdateFrame(FrameEventArgs e) - { - base.OnUpdateFrame(e); - - _proj3D = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2); - _proj2D = Matrix4.CreateOrthographic(4, 4f * Height / Width, -1, 1); - - _t += e.Time; - - if (Keyboard[Key.Left]) - { - var rot = MVec3d.Complex(-e.Time, MVec3d.Unit12); - _b1 = rot * _b1 * ~rot; // todo make blade classes - _b2 = rot * _b2 * ~rot; - } - - if (Keyboard[Key.Right]) - { - var rot = MVec3d.Complex(e.Time, MVec3d.Unit12); - _b1 = rot * _b1 * ~rot; - _b2 = rot * _b2 * ~rot; - } - - if (Keyboard[Key.Up]) - { - _pivot -= (_b1 ^ _b2) * MVec3d.Unit123 * e.Time; - } - - if (Keyboard[Key.Down]) - { - _pivot += (_b1 ^ _b2) * MVec3d.Unit123 * e.Time; - } - - if (_pause) - return; - -// var b1 = MVec3d.Unit2; -// var b2 = MVec3d.Unit3; -// var pivot = new Vector3d(Math.Cos(_t), 0, 0); - -// var b1 = MVec3d.Vector(0, Math.Cos(_t * .2), Math.Sin(_t * .5)); -// var b2 = MVec3d.Vector(Math.Cos(_t * .6), 0, Math.Sin(_t * .4)); -// var pivot = Vector3d.Zero; - -// var b1 = MVec3d.Vector(Math.Cos(_t/2), Math.Sin(_t/2), 0); -// var b2 = MVec3d.Unit3; -// var pivot = Vector3d.Zero; - -// var blade = MVec3d.Outer(b1, b2); - -// var rot = MVec3d.Complex(_t / 2, MVec3d.Unit12); -// var b1 = rot * MVec3d.Unit2 * ~rot; -// var b2 = rot * MVec3d.Unit3 * ~rot; -// var pln = MVec3d.Outer(b1, b2); -// var blade = pln; -// var pivot = Vector3d.Zero; - - _polyg = _tetra.Intersection(_b1 ^ _b2, _pivot); - - if (_div) - { - var px = ((_pivot & _b1) * ~_b1) & _b1; - var py = ((_pivot ^ _b1) * ~_b1) & ((_b2 ^ _b1) * ~_b1); - - for (var i = 0; i < _polyg.Points.Length; i++) - { - MVec3d p = _polyg.Points[i]; - var x = ((p & _b1) * ~_b1) & _b1; - var y = ((p ^ _b1) * ~_b1) & ((_b2 ^ _b1) * ~_b1); - _polyg.Points[i] = new Vector3d(x.E - px.E, y.E - py.E, 0); - } - } - } - - protected override void OnKeyDown(KeyboardKeyEventArgs e) - { - base.OnKeyDown(e); - - switch (e.Key) - { - case Key.Space: - _pause = !_pause; - break; - case Key.A: - _div = !_div; - break; - case Key.R: - _b1 = MVec3d.Unit1; - _b2 = MVec3d.Unit3; - _pivot = MVec3d.Zero; - break; - } - } - } + if (e.Key == Key.Comma) + _wire = !_wire; + } + } } \ No newline at end of file diff --git a/Tetrahedrons/Tetrahedrons.csproj b/Tetrahedrons/Tetrahedrons.csproj index a76912a..546af1a 100644 --- a/Tetrahedrons/Tetrahedrons.csproj +++ b/Tetrahedrons/Tetrahedrons.csproj @@ -42,10 +42,12 @@ - + + + diff --git a/Tetrahedrons/Util.cs b/Tetrahedrons/Util.cs new file mode 100644 index 0000000..6f3aa54 --- /dev/null +++ b/Tetrahedrons/Util.cs @@ -0,0 +1,25 @@ +using OpenTK.Graphics.OpenGL; + +namespace Tetrahedrons +{ + public static class Util + { + public static void Color3(MVec4D m) + { + m *= .5; + GL.Color3(.5 + m.X, .5 + m.Y, .5 * m.Z); + } + + public static void Vertex4(MVec4D m) + { + var i = m.W / 2; + GL.Color3(.5 + i, 0, .5 - i); + GL.Vertex3(m.X, m.Y, m.Z); + } + + public static void Vertex3(MVec4D m) + { + GL.Vertex3(m.X, m.Y, m.Z); + } + } +} \ No newline at end of file