@@ 0,0 1,194 @@
+
+A = 4 * pi
+SH1_C0 = sqrt(3 / (4 * pi))
+SH1_C1 = sqrt(1 / (4 * pi))
+SH1_Basis = (SH1_C0, SH1_C0, SH1_C0, SH1_C1)
+
+# clean solid angle computation for bcc grid
+
+def vxor3 (a, b):
+ return (int(a[0]) ^^ int(b[0]), int(a[1]) ^^ int(b[1]), int(a[2]) ^^ int(b[2]))
+def vand3 (a, b):
+ return (a[0] & b[0], a[1] & b[1], a[2] & b[2])
+def vor3 (a, b):
+ return (a[0] | b[0], a[1] | b[1], a[2] | b[2])
+def vadd3 (a, b):
+ return (a[0] + b[0], a[1] + b[1], a[2] + b[2])
+def vadd3s (a, b):
+ return (a[0] + b, a[1] + b, a[2] + b)
+def vsub3 (a, b):
+ return (a[0] - b[0], a[1] - b[1], a[2] - b[2])
+def vsub3s (a, b):
+ return (a[0] - b, a[1] - b, a[2] - b)
+def vmul3 (a, b):
+ return (a[0] * b[0], a[1] * b[1], a[2] * b[2])
+def vmul4 (a, b):
+ return (a[0] * b[0], a[1] * b[1], a[2] * b[2], a[3] * b[3])
+def vmul4s (a, b):
+ return (a[0] * b, a[1] * b, a[2] * b, a[3] * b)
+def vmul3s (a, b):
+ return (a[0] * b, a[1] * b, a[2] * b)
+
+def dot3 (a, b):
+ return (a[0] * b[0] + a[1] * b[1] + a[2] * b[2])
+def length (v):
+ return sqrt(dot3(v,v))
+def normalize (v):
+ return vmul3s(v, 1 / length(v))
+
+def center3 (a,b,c):
+ return vmul3s(vadd3(vadd3(a,b),c),1/3)
+def center4 (a,b,c,d):
+ return vmul3s(vadd3(vadd3(a,b),vadd3(c,d)),1/4)
+
+def bcc (v):
+ s = (v[0] + v[1] + v[2]) / 6
+ return vsub3s(v, s)
+
+def cc (v):
+ s = (v[0] + v[1] + v[2]) / 3
+ return vadd3s(v, s)
+
+def anglebasis(a, b):
+ ab = dot3(a, b)
+ aabb = dot3(a,a) * dot3(b, b)
+ return (ab / sqrt(aabb), sqrt(1 - ab * ab / aabb))
+
+def solid_angle (u, v, w):
+ return \
+ acos ((u[0] - v[0] * w[0]) / (v[1] * w[1])) + \
+ acos ((v[0] - w[0] * u[0]) / (w[1] * u[1])) + \
+ acos ((w[0] - u[0] * v[0]) / (u[1] * v[1])) - pi
+
+def simplex_solid_angle(U, V, W):
+ u = anglebasis(V, W)
+ v = anglebasis(W, U)
+ w = anglebasis(U, V)
+ return solid_angle(u, v, w)
+
+def edge (code):
+ return (code & 1, (code >> 1) & 1, (code >> 2) & 1)
+
+def simplex (n):
+ a = (0b101100110010011001 >> (n * 3)) & 7
+ b = (0b100110010011001101 >> (n * 3)) & 7
+ return (bcc(edge(a)), bcc(edge(b)), bcc(edge(0b000)), bcc(edge(0b111)))
+
+def simplex_code (n):
+ a = (0b101100110010011001 >> (n * 3)) & 7
+ b = (0b100110010011001101 >> (n * 3)) & 7
+ return (edge(a), edge(b), edge(0b000), edge(0b111))
+
+def face_vertices (s, n):
+ return (
+ s[(0b00011011 >> (n * 2)) & 3],
+ s[(0b01001110 >> (n * 2)) & 3],
+ s[(0b10110001 >> (n * 2)) & 3])
+
+def adj (s, f):
+ o = [0,0,0]
+ if (f == 0):
+ simplex = s + 5
+ elif (f == 1):
+ simplex = s + 1
+ elif (f == 2):
+ o[s >> 1] = 1
+ simplex = s + 4 - ((s & 1) << 1)
+ elif (f == 3): # (f == 3)
+ o[((s + 3) % 6) >> 1] = -1
+ simplex = s + 2 + ((s & 1) << 1)
+ else:
+ assert False
+ return ((simplex % 6), (f ^^ 1), bcc(tuple(o)))
+
+syms = []
+symtable = {}
+symtable[0] = 0
+usedsyms = {}
+rsymtable = {}
+def symlify(x, s = None):
+ global nextsym
+ if (not (x in symtable)):
+ assert s
+ k = 1
+ while (s + str(k)) in usedsyms:
+ k += 1
+ usedsyms[s + str(k)] = True
+ d = var(s + str(k))
+ sx = x.simplify()
+ if sx.substitute(rsymtable) < 0:
+ sx = -sx
+ symtable[-sx] = -d
+ symtable[sx] = d
+ symtable[d] = d
+ symtable[-d] = -d
+ rsymtable[d] = sx
+ syms.append((d,sx))
+ return symtable[x]
+
+table = [[[] for i in range(4)] for i in range(6)]
+
+# gather pattern
+for i in range(6):
+ #print("simplex",i)
+ s = simplex(i)
+ for f in range(4):
+ #print("face",f)
+ ai,af,ofs = adj(i,f)
+ aa = simplex(ai)
+ aa = vadd3(aa[0],ofs),vadd3(aa[1],ofs),vadd3(aa[2],ofs),vadd3(aa[3],ofs)
+ c = center4(*aa)
+ for gg in range(3):
+ g = (f + 1 + gg) % 4
+ vs = face_vertices(s, g)
+ v0 = vsub3(vs[0], c)
+ v1 = vsub3(vs[1], c)
+ v2 = vsub3(vs[2], c)
+ vc = center3(vs[0],vs[1],vs[2])
+ n = vsub3(vc,c)
+ n = (*normalize(n), 1)
+ sa = symlify(simplex_solid_angle(v0,v1,v2).simplify(), 'sa')
+ shn = vmul4(n,SH1_Basis)
+ shn = ([symlify(x,'t') for x in shn])
+ table[i][f].append(vmul4s(shn, sa))
+
+ #print()
+
+for i in range(6):
+ for f in range(4):
+ m = table[i][f]
+ m.sort(key = lambda k: str(k[3]))
+
+for q in range(3):
+ for f in range(4):
+ for i in range(6):
+ W = table[i][f][q]
+ x,y,z,w = W
+ if (q < 2):
+ b = 'u'
+ else:
+ b = 'v'
+ x = symlify(x, b)
+ y = symlify(y, b)
+ z = symlify(z, b)
+ w = symlify(w, 'q')
+
+print("Tx3x4x6")
+for i in range(6):
+ print(" Tx3x4")
+ for f in range(4):
+ s = " Tx3 "
+ for q in range(3):
+ W = table[i][f][q]
+ x,y,z,w = W
+ x = symlify(x)
+ y = symlify(y)
+ z = symlify(z)
+ w = symlify(w)
+ s += "(T " + repr(x) + " " + repr(y) + " " + repr(z) + " " + repr(w) + ") "
+ print(s)
+print()
+print("symbols:")
+for k,v in syms:
+ print(k,":=",v)
+
@@ 17,6 17,10 @@
* plane along +YZ splits simplices 1 2 3 | 4 5 0
* plane along +XZ splits simplices 0 1 2 | 3 4 5
* plane along +XY splits simplices 5 0 1 | 2 3 4
+ * the solid angles of a simplices neighbors visible faces are, approximately
+ a := 2*acos(3/sqrt(14)) - acos(6/7); for the smaller one
+ and
+ b := (pi - a)/2 = pi/2 - acos(3/sqrt(14)) + acos(6/7)/2
using import glm
@@ 38,14 42,14 @@ fn simplex-volume (A B C D)
fn anglebasis (a b)
ab := (dot a b)
aabb := (dot a a) * (dot b b)
- vec2
+ dvec2
ab / (sqrt aabb)
- sqrt (1.0 - ab * ab / aabb)
+ sqrt (1.0:f64 - ab * ab / aabb)
# for when both arguments are already normalized
fn nanglebasis (a b)
ab := (dot a b)
- vec2 ab (sqrt (1.0 - ab * ab))
+ dvec2 ab (sqrt (1.0:f64 - ab * ab))
# solid angle of basis angles, specified as sin/cos pairs
fn simplex-solid-angle (u v w)
@@ 54,8 58,8 @@ fn simplex-solid-angle (u v w)
acos ((u.x - v.x * w.x) / (v.y * w.y))
acos ((v.x - w.x * u.x) / (w.y * u.y))
acos ((w.x - u.x * v.x) / (u.y * v.y))
- -pi
- ? (A != A) 0.0 A
+ -pi:f64
+ ? (A != A) 0.0:f64 A
# solid angle of basis from normalized vectors
fn simplex-points-solid-angle (U V W)
@@ 76,6 80,10 @@ fn pyramid-points-solid-angle (A B C D)
simplex-solid-angle bc ac ab
simplex-solid-angle da ac cd
+A:f64 := 4.0:f64 * pi:f64
+SH1_Basis:f64 := (sqrt ((dvec4 (dvec3 3.0:f64) 1.0:f64) / A:f64)) #* (dvec4 1 1 1 1)
+SH1_Basis := (vec4 SH1_Basis:f64)
+
inline cc->bcc (p)
p - (p.x + p.y + p.z) / 6.0
@@ 89,6 97,9 @@ fn cos_angle (a b)
# constructive edges are either unit length (orthogonal) or (sqrt 3) / 2 (body diagonal)
+fn iedge (code)
+ (ivec3 code (code >> 1) (code >> 2)) & 1
+
fn edge (code)
vec3 ((ivec3 code (code >> 1) (code >> 2)) & 1)
@@ 178,6 189,59 @@ fn neighbor (s f)
trap;
_ (simplex % 6) (f ^ 1) (copy o)
+define sh1-constants
+ sa2 := 2 * (acos (3 / (sqrt 14.0))) - (acos (6 / 7)); sa1 := (pi - sa2) / 2
+ t4 := 1 / (2 * (sqrt pi))
+ q1 := sa1 * t4; q2 := sa2 * t4
+ t11 := 1 / (sqrt (87 * pi))
+ t1 := 7 / 2 * t11; t2 := 2 * t11; t3 := 7 * t11
+ t9 := 8 * t11; t10 := 1 / 2 * t11; t12 := 4 * t11
+ t6 := (sqrt (3 / (29 * pi)))
+ t5 := 5 / 2 * t6; t7 := 3 / 2 * t6; t8 := 2 * t6
+ u1 := sa1 * t1; u2 := sa1 * t2; u3 := sa1 * t3
+ u4 := sa1 * t5; u5 := sa1 * t6; u6 := sa1 * t8
+ u7 := sa1 * t7; u8 := sa1 * t10; u9 := sa1 * t9
+ u10 := sa1 * t11
+ v1 := sa2 * t7; v2 := sa2 * t6; v3 := sa2 * t8
+ v4 := sa2 * t3; v5 := sa2 * t10; v6 := sa2 * t12
+ T := vec4; Tx3 := (array T 3); Tx3x4 := (array Tx3 4)
+ (array Tx3x4 6)
+ Tx3x4
+ Tx3 (T u1 u2 -u3 q1) (T u4 0 -u5 q1) (T v1 -v2 -v3 q2)
+ Tx3 (T 0 -u4 u5 q1) (T -u2 -u1 u3 q1) (T v2 -v1 v3 q2)
+ Tx3 (T -u6 -u7 -u5 q1) (T -u9 -u8 u10 q1) (T -v4 v5 -v6 q2)
+ Tx3 (T u8 u9 -u10 q1) (T u7 u6 u5 q1) (T -v5 v4 v6 q2)
+ Tx3x4
+ Tx3 (T -u2 u3 -u1 q1) (T 0 u5 -u4 q1) (T v2 v3 -v1 q2)
+ Tx3 (T u4 -u5 0 q1) (T u1 -u3 u2 q1) (T v1 -v3 -v2 q2)
+ Tx3 (T -u6 -u5 -u7 q1) (T -u9 u10 -u8 q1) (T -v4 -v6 v5 q2)
+ Tx3 (T u8 -u10 u9 q1) (T u7 u5 u6 q1) (T -v5 v6 v4 q2)
+ Tx3x4
+ Tx3 (T -u3 u1 u2 q1) (T -u5 u4 0 q1) (T -v3 v1 -v2 q2)
+ Tx3 (T u5 0 -u4 q1) (T u3 -u2 -u1 q1) (T v3 v2 -v1 q2)
+ Tx3 (T -u5 -u6 -u7 q1) (T u10 -u9 -u8 q1) (T -v6 -v4 v5 q2)
+ Tx3 (T -u10 u8 u9 q1) (T u5 u7 u6 q1) (T v6 -v5 v4 q2)
+ Tx3x4
+ Tx3 (T -u1 -u2 u3 q1) (T -u4 0 u5 q1) (T -v1 v2 v3 q2)
+ Tx3 (T 0 u4 -u5 q1) (T u2 u1 -u3 q1) (T -v2 v1 -v3 q2)
+ Tx3 (T -u7 -u6 -u5 q1) (T -u8 -u9 u10 q1) (T v5 -v4 -v6 q2)
+ Tx3 (T u9 u8 -u10 q1) (T u6 u7 u5 q1) (T v4 -v5 v6 q2)
+ Tx3x4
+ Tx3 (T u2 -u3 u1 q1) (T 0 -u5 u4 q1) (T -v2 -v3 v1 q2)
+ Tx3 (T -u4 u5 0 q1) (T -u1 u3 -u2 q1) (T -v1 v3 v2 q2)
+ Tx3 (T -u7 -u5 -u6 q1) (T -u8 u10 -u9 q1) (T v5 -v6 -v4 q2)
+ Tx3 (T u9 -u10 u8 q1) (T u6 u5 u7 q1) (T v4 v6 -v5 q2)
+ Tx3x4
+ Tx3 (T u3 -u1 -u2 q1) (T u5 -u4 0 q1) (T v3 -v1 v2 q2)
+ Tx3 (T -u5 0 u4 q1) (T -u3 u2 u1 q1) (T -v3 -v2 v1 q2)
+ Tx3 (T -u5 -u7 -u6 q1) (T u10 -u8 -u9 q1) (T -v6 v5 -v4 q2)
+ Tx3 (T -u10 u9 u8 q1) (T u5 u6 u7 q1) (T v6 v4 -v5 q2)
+
+global sh1-constants = sh1-constants
+
+inline gather-solid-angle (s f)
+ unpack (sh1-constants @ s @ f)
+
inline sort3 (pred a b c)
let a b =
if (pred a b) (_ b a)
@@ 217,6 281,23 @@ inline simplexfacecoords (s f ofs)
print "simplices" i j "faces" f g
#print va vb
+fn dv4xvT (a)
+ dmat4
+ a.x * a
+ a.y * a
+ a.z * a
+ a.w * a
+
+type+ dmat4x4
+ inline __+ (cls T)
+ static-if (cls == T)
+ inline (a b)
+ dmat4
+ (a @ 0) + (b @ 0)
+ (a @ 1) + (b @ 1)
+ (a @ 2) + (b @ 2)
+ (a @ 3) + (b @ 3)
+
for i in (range 6)
s := (simplex i)
print "simplex" i
@@ 227,45 308,59 @@ for i in (range 6)
inline (i) (cc->bcc (edge i))
#inline (i) (edge i)
unpack s
- center :=
- /
- +
- v @ 0
- v @ 1
- v @ 2
- v @ 3
- 4
- print
- simplex-volume (v @ 0) (v @ 1) (v @ 2) (v @ 3)
+ local sa = 0.0:f64
+ probe := (dvec4 (normalize (dvec3 1 1 1)) 1) * pi:f64 * SH1_Basis:f64
+ #probe := (dvec4 1 0 0 1) * 4 * pi:f64 * SH1_Basis:f64
+ local sh = (dvec4 0)
for f in (range 4)
- let v1 v2 v3 =
+ print "face" f
+ #let v1 v2 v3 =
va-map
inline (i) (v @ i)
unpack (simplex-local-face-vertices f)
- c := ((+ v1 v2 v3) / 3)
let i1 f1 d1 = (neighbor i f)
local vb =
arrayof vec3
va-map
- inline (i) (cc->bcc (edge i))
+ inline (i) (cc->bcc ((edge i) + (vec3 d1)))
#inline (i) (edge i)
unpack (simplex i1)
- let vb1 vb2 vb3 =
- va-map
- inline (i) (vb @ i)
- unpack (simplex-local-face-vertices f1)
- c1 := (((+ vb1 vb2 vb3) / 3) + (cc->bcc (vec3 d1)))
- #c1 := ((+ vb1 vb2 vb3) / 3)
- #A := (triangle-area v1 v2 v3)
- #print v1 v2 v3 A
- print (triangle-normal v1 v2 v3)
- print (c - c1)
- print (dot (triangle-normal v1 v2 v3) (triangle-normal vb1 vb2 vb3))
- #print
- simplex-points-solid-angle
- normalize (v1 - center)
- normalize (v2 - center)
- normalize (v3 - center)
+ center :=
+ /
+ +
+ vb @ 0
+ vb @ 1
+ vb @ 2
+ vb @ 3
+ 4
+ for gg in (range 3)
+ let g = ((f + 1 + gg) % 4)
+ let vb1 vb2 vb3 =
+ va-map
+ inline (i) (v @ i)
+ unpack (simplex-local-face-vertices g)
+ ct := (((vb1 + vb2 + vb3) / 3) - center)
+ ct := (normalize ct)
+ gsa :=
+ simplex-points-solid-angle
+ normalize (dvec3 (vb1 - center))
+ normalize (dvec3 (vb2 - center))
+ normalize (dvec3 (vb3 - center))
+ shb := (dvec4 ct 1) * SH1_Basis:f64 * gsa
+ let n0 n1 n2 = (gather-solid-angle i f)
+ let l1 l2 l3 =
+ length ((vec4 shb) - n0)
+ length ((vec4 shb) - n1)
+ length ((vec4 shb) - n2)
+ print (min l1 l2 l3)
+ assert ((min l1 l2 l3) < 0.00001)
+ sh += (dot shb probe) * shb
+ print "neighbor face" i1 g ct gsa
+ sa += gsa
+ #sh += shc * probe
+ print "total" (sa / (4 * pi:f64))
+ print
+ (dot sh sh) / (4 * pi)
print;