M lib/tukan/dsp/utils.sc +5 -1
@@ 29,7 29,11 @@ let
fn hz (c)
""""converts from midi note to frequency in hz
conversion according to https://en.wikipedia.org/wiki/MIDI_tuning_standard
- 440.0 * (2.0 ** ((c - 69.0) / 12.0))
+ 440.0 * (exp2 ((c - 69.0) / 12.0))
+
+fn midinote (hz)
+ """"inverse of hz()
+ (log2 (hz / 440.0)) * 12.0 + 69.0
# OPL oscillators
M lib/tukan/random.sc +4 -6
@@ 264,12 264,10 @@ inline Random (bits opts...)
fn triangle-area (self)
let u =
random-vec2 self
- let sux =
- sqrt u.x
- vec3
- 1.0 - sux
- sux * (1.0 - u.y)
- sux * u.y
+ let u =
+ if (u.x + u.y > 1.0) (1.0 - u)
+ else u
+ vec3 u.x u.y (1.0 - u.x - u.y)
inline triangular (self low high mode)
let T = (typeof low)
A => testing/bcc_shearing.sc +271 -0
@@ 0,0 1,271 @@
+
+# some statistics about the BCC truncated octahedron mesh
+
+ * two types of edges, 3 x length 1, 4 x length sqrt(3)/2
+ * all simplices have the volume 1/12
+ * all simplex faces have the same area 1/sqrt(8) and solid angle pi relative
+ to their centers
+ * simplices are numbered counterclockwise from origin, and faces of
+ simplices of neighboring cubic cells pair up in even and odd
+ permutations, i.e. 1-3 3-5 5-1 are connected, as well as 0-2 2-4 4-0
+ * simplex faces of index 0 and 1 are connected to the previous / next
+ simplex.
+ * simplex faces of index 2 point towards +XYZ, faces of index 3 point
+ towards -XYZ
+ * simplex numbering is so that face[2] memory planes are +XXYYZZ, and
+ face[3] planes are -YZZXXY
+ * 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
+
+using import glm
+
+fn triangle-normal (A B C)
+ normalize (cross (B - A) (C - A))
+fn triangle-plane (A B C)
+ n := (triangle-normal A B C)
+ vec4 n (- (dot n A))
+fn triangle-area (A B C)
+ (length (cross (B - A) (C - A))) / 2.0
+fn triangle-area-normal (A B C)
+ (cross (B - A) (C - A)) / 2.0
+
+# assuming D is at origin
+fn simplex-volume (A B C D)
+ (dot (A - D) (cross (B - D) (C - D))) / 6.0
+
+# returns cos and sin of angle between a and b
+fn anglebasis (a b)
+ ab := (dot a b)
+ aabb := (dot a a) * (dot b b)
+ vec2
+ ab / (sqrt aabb)
+ sqrt (1.0 - 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))
+
+# solid angle of basis angles, specified as sin/cos pairs
+fn simplex-solid-angle (u v w)
+ A :=
+ +
+ 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
+
+# solid angle of basis from normalized vectors
+fn simplex-points-solid-angle (U V W)
+ u := (nanglebasis V W)
+ v := (nanglebasis W U)
+ w := (nanglebasis U V)
+
+ simplex-solid-angle u v w
+
+# solid angle of pyramid from normalized vectors
+fn pyramid-points-solid-angle (A B C D)
+ ab := (nanglebasis A B)
+ ac := (nanglebasis A C)
+ da := (nanglebasis A D)
+ bc := (nanglebasis B C)
+ cd := (nanglebasis C D)
+ +
+ simplex-solid-angle bc ac ab
+ simplex-solid-angle da ac cd
+
+inline cc->bcc (p)
+ p - (p.x + p.y + p.z) / 6.0
+
+inline bcc->cc (p)
+ p + (p.x + p.y + p.z) / 3.0
+
+fn cos_angle (a b)
+ degrees
+ acos
+ dot (normalize a) (normalize b)
+
+# constructive edges are either unit length (orthogonal) or (sqrt 3) / 2 (body diagonal)
+
+fn edge (code)
+ vec3 ((ivec3 code (code >> 1) (code >> 2)) & 1)
+
+# corner edge
+e000 := 0b000
+
+# orthogonal basis edges of unit length; permutations of +-2/3 +-2/3 -+1/3
+e011 := 0b011
+e101 := 0b101
+e110 := 0b110
+
+# diagonal edges of length (sqrt 3) / 2; permutations of +-5/6 -+1/6 -+1/6
+e100 := 0b100
+e010 := 0b010
+e001 := 0b001
+# fourth diagonal of length (sqrt 3) / 2; coordinate 1/2 1/2 1/2
+e111 := 0b111
+
+#
+ V0
+ +
+ F1 /|\ F3
+ \/ | /
+ /\ | \
+ V2 +-------+ V1
+ \ | \/
+ / | /\
+ F0 \|/ F2
+ +
+ V3
+
+ each cell has four vertices, four half-faces, and twelve half-edges
+ we only store the vertices, all other information is implicit
+ per convention, each tet face i starts with vertex i,
+ and each group of 3 half-edges maps to each tet face,
+ in ccw order as seen from the inside; Tets are inverted, that is,
+ their faces point inwards.
+ vertices are enumerated so that a triangle strip of indices
+ 0-1-2-3-0-1 draws the complete tet.
+ ->opposite maps each first triangle edge (0,3,6,9) to the adjacent
+ flipped triangle edge of the other half-face.
+ ->star maps each vertex to the previous and next edge in an adjacent
+ cell with the same starting vertex, describing a linked list of
+ all cells connected to this vertex. The list is circular.
+
+# all simplices within the cube have edge 0b000 - 0b111
+# the other two points are, counterclockwise
+ 0b001 0b011 0b010 0b110 0b100 0b101
+inline simplex (n)
+ # returns indices of nth simplex (0..5)
+ q := ((ivec2 0b101100110010011001 0b100110010011001101) >> (n * 3)) & 7
+ ivec4 q 0b000 0b111
+inline simplex-local-face-vertices (n)
+ # returns local vertex indices of nth face (0..3)
+ # 0 -> 3 2 1
+ # 1 -> 2 3 0
+ # 2 -> 1 0 3
+ # 3 -> 0 1 2
+ ((ivec3
+ 0b00011011
+ 0b01001110
+ 0b10110001) >> (n * 2)) & 3
+#
+ simplex pairs between cells (face 2 to face 3):
+ +X: 0-4, 1-3 (a-b = -+2)
+ +Y: 2-0, 3-5 (a+b = 6)
+ +Z: 4-2, 5-1 (a+b = 4)
+ all even and odd numbers pair up
+
+fn neighbor (s f)
+ # given simplex index and face number, produces simplex
+ index, face number and cell offset of the neighboring face
+ local o = (ivec3 0)
+ let simplex =
+ switch f
+ case 0
+ s + 5
+ case 1
+ s + 1
+ case 2
+ o @ (s >> 1) = 1
+ s + 4 - ((s & 1) << 1)
+ case 3
+ o @ (((s + 3) % 6) >> 1) = -1
+ s + 2 + ((s & 1) << 1)
+ default
+ trap;
+ _ (simplex % 6) (f ^ 1) (copy o)
+
+inline sort3 (pred a b c)
+ let a b =
+ if (pred a b) (_ b a)
+ else (_ a b)
+ let b c =
+ if (pred b c) (_ c b)
+ else (_ b c)
+ let a b =
+ if (pred a b) (_ b a)
+ else (_ a b)
+ _ a b c
+
+inline simplexfacecoords (s f ofs)
+ triangle-plane
+ va-map
+ inline (i)
+ (edge i) + ofs
+ va-map
+ inline (i) (s @ i)
+ unpack (simplex-local-face-vertices f)
+
+#for i in (range 6)
+ s := (simplex i)
+ # vertices
+ for f in (range 4)
+ let va =
+ simplexfacecoords s f (vec3 0)
+ for j in (range 6)
+ t := (simplex j)
+ for g in (range 4)
+ let vb =
+ simplexfacecoords t g (vec3 0 0 1)
+ #print i (unpack va)
+ #print j (unpack vb)
+ if (i != j)
+ if (va == vb)
+ print "simplices" i j "faces" f g
+ #print va vb
+
+for i in (range 6)
+ s := (simplex i)
+ print "simplex" i
+ # vertices
+ local v =
+ arrayof vec3
+ va-map
+ 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)
+ for f in (range 4)
+ 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) (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)
+
+ print;
+
A => testing/solid_angle.sc +460 -0
@@ 0,0 1,460 @@
+
+import ..lib.tukan.use
+using import glm
+
+# returns cos and sin of angle between a and b
+fn anglebasis (a b)
+ ab := (dot a b)
+ aabb := (dot a a) * (dot b b)
+ vec2
+ ab / (sqrt aabb)
+ sqrt (1.0 - 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))
+
+# solid angle of basis angles, specified as sin/cos pairs
+fn simplex-solid-angle (u v w)
+ A :=
+ +
+ 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
+
+# solid angle of basis from normalized vectors
+fn simplex-points-solid-angle (U V W)
+ u := (nanglebasis V W)
+ v := (nanglebasis W U)
+ w := (nanglebasis U V)
+
+ simplex-solid-angle u v w
+
+# solid angle of pyramid from normalized vectors
+fn pyramid-points-solid-angle (A B C D)
+ ab := (nanglebasis A B)
+ ac := (nanglebasis A C)
+ da := (nanglebasis A D)
+ bc := (nanglebasis B C)
+ cd := (nanglebasis C D)
+ +
+ simplex-solid-angle bc ac ab
+ simplex-solid-angle da ac cd
+
+do
+ # octahedral corner
+ C := (vec3 1 1 1)
+ a := 0.0
+ U := (mix (vec3 1 0 0) C a)
+ V := (mix (vec3 0 1 0) C a)
+ W := (mix (vec3 0 0 1) C a)
+
+ print
+ /
+ simplex-points-solid-angle U V W
+ 2 * tau
+
+do
+ # cube face
+ w := 2.0
+ A := (normalize (vec3 -1 -1 1))
+ B := (normalize (vec3 (-1 + w) -1 1))
+ C := (normalize (vec3 (-1 + w) (-1 + w) 1))
+ D := (normalize (vec3 -1 (-1 + w) 1))
+ print
+ /
+ pyramid-points-solid-angle A B C D
+ 2 * tau
+
+do
+ # 0.400669685
+ # 0.42347
+
+ # neighboring cube face
+ A := (normalize (vec3 -0.5 -0.5 1.5))
+ B := (normalize (vec3 0.5 -0.5 1.5))
+ C := (normalize (vec3 0.5 0.5 1.5))
+ D := (normalize (vec3 -0.5 0.5 1.5))
+ u :=
+ pyramid-points-solid-angle A B C D
+ # side cube face
+ A := (normalize (vec3 -0.5 0.5 0.5))
+ B := (normalize (vec3 0.5 0.5 0.5))
+ C := (normalize (vec3 0.5 0.5 1.5))
+ D := (normalize (vec3 -0.5 0.5 1.5))
+ v :=
+ pyramid-points-solid-angle A B C D
+ print u
+ print v
+ print ((u + 4 * v) / (2 * tau))
+
+# factorial
+fn ! (n)
+ fold (x = 1) for i in (range 1 (n + 1))
+ x * i
+
+# from https://basesandframes.wordpress.com/2016/05/11/spherical-harmonic-lighting-the-gritty-details/
+fn P (l m x)
+ x as:= f64
+ # evaluate an Associated Legendre Polynomial P(l,m,x) at x
+ let pmm =
+ if (m > 0)
+ somx2 := (sqrt ((1.0:f64 - x) * (1.0:f64 + x)))
+ fact := 1.0:f64
+ pmm := 1.0:f64
+ _
+ fold (pmm fact) for i in (range 0 m)
+ _ (pmm * -fact * somx2)
+ fact + 2.0:f64
+ ()
+ else 1.0:f64
+ if (l == m) pmm
+ else
+ pmmp1 := x * (2.0:f64 * m as f64 + 1.0:f64) * pmm
+ if (l == (m + 1)) pmmp1
+ else
+ pll := 0.0:f64
+ ll := m + 2
+ _
+ fold (pll pmm pmmp1) for ll in (range (m + 2) l)
+ _
+ ((2.0:f64 * ll as f64 - 1.0) * x * pmmp1
+ - (ll + m - 1) as f64 * pmm) / (ll - m) as f64
+ pmmp1
+ pll
+ ()
+
+for i in (range 11)
+ i := i / 10 * 2 - 1
+ print i (P 1 1 i) (- (sqrt ((1.0:f64 - i) * (1.0:f64 + i))))
+
+fn K (l m)
+ m := (abs m)
+ sqrt
+ ((2 * l + 1) * (! (l - m))) as f64 / (2 * tau:f64 * (! (l + m)) as f64)
+
+fn Y (l m theta phi)
+ if (m > 0)
+ (sqrt 2.0:f64) * (K l m) * (cos (m as f64 * phi)) * (P l m (cos theta))
+ elseif (m < 0)
+ (sqrt 2.0:f64) * (K l m) * (sin (-m as f64 * phi)) * (P l -m (cos theta))
+ else # m == 0
+ (K l 0) * (P l 0 (cos theta))
+
+# theta and phi from normal vector
+fn sh1_normal_angle (dir)
+ dir2 := dir.xy
+ l := (length dir2)
+ vec2 (acos dir.z) (atan2 dir2.y dir2.x)
+
+# normal vector from theta and phi
+fn sh1_angle_normal (angles)
+ let ct cp = (unpack (cos angles))
+ let st sp = (unpack (sin angles))
+ vec3 (st * cp) (st * sp) ct
+
+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)
+print "SH1 Basis:" SH1_Basis
+
+do
+ # theta = vertical; phi = horizontal
+ let tp =
+ sh1_normal_angle (normalize (vec3 1 0 1))
+ theta := tp.x as f64
+ phi := tp.y as f64
+ print theta phi
+ print (sh1_angle_normal (vec2 theta phi))
+ print "Y 1 1 = " ((Y 1 1 theta phi) / SH1_Basis.x)
+ print "Y 1 -1 = " ((Y 1 -1 theta phi) / SH1_Basis.y)
+ print "Y 1 0 = " ((Y 1 0 theta phi) / SH1_Basis.z)
+ print "Y 0 0 = " ((Y 0 0 theta phi) / SH1_Basis.w)
+
+print "(K 0 0) =" (K 0 0)
+print "(K 1 0) =" (K 1 0)
+print "(K 1 1) =" (K 1 1)
+
+do
+ # uniform dip
+ u := SH1_Basis * (vec4 0 0 0 1)
+ v := SH1_Basis * (vec4 (/ (sqrt 3.0)) 0 0 (/ 4))
+
+ # K²*u.x*v.x + L²*u.w*v.w = L
+ # K²*u.x*v.x + L²*u.w*v.w = L
+
+ print (4.0 * pi * (dot u u))
+ print (4.0 * pi * (dot v v))
+ print (4.0 * pi * (dot u v))
+
+# sh1 unit vector for normalized direction
+inline sh1_unit (dir)
+ vec4 dir 1.0
+
+# sh1 unit vector for normalized direction
+#fn sh1_unit (dir)
+ dir2 := dir.xy
+ l := (length dir2)
+ let cp sp = (unpack (? (l == 0.0) (vec2 0) (dir2 / l)))
+ let ct st = dir.z (sqrt (1 - dir.z * dir.z))
+ vec4 (st * cp) (st * sp) ct 1.0
+
+print
+ sh1_unit (vec3 1 0 0)
+ sh1_unit (vec3 0 1 0)
+ sh1_unit (vec3 0 0 1)
+print
+ sh1_unit (normalize (vec3 0 1 1))
+ sh1_unit (normalize (vec3 1 0 1))
+ sh1_unit (normalize (vec3 1 1 0))
+print
+ sh1_unit (normalize (vec3 1 1 1))
+
+fn sh_rotate (dir coeffs)
+ * (vec4 dir 1.0)
+ vec4 coeffs.x coeffs.x -coeffs.x coeffs.y
+
+fn sh_project_cone (dir angle)
+ let ca sa = (cos angle) (sin angle)
+ sh_rotate dir
+ vec2
+ 0.75 * sa * sa
+ 0.5 * (1 - ca)
+
+fn sh1_clamped_cosine_lobe (dir)
+ k0 := (sqrt (pi / 2.0))
+ k1 := (sqrt (pi / 3.0))
+ * (vec4 dir 1)
+ vec4 -k1 -k1 k1 k0
+
+fn sh_project_lobe (dir)
+ sh_rotate dir (vec2 0.5 0.25)
+
+#print
+ dot
+ sh_project_lobe
+ vec3 1 0 0
+ sh_project_cone
+ vec3 1 0 0
+ 2 * tau / 2
+
+do
+ # integrate cosine lobe 2d
+ fn triangle_normal (A B C)
+ normalize (cross (B - A) (C - A))
+ fn triangle_area (A B C)
+ (length (cross (B - A) (C - A))) / 2.0
+ fn triangle_area_normal (A B C)
+ (cross (B - A) (C - A)) / 2.0
+
+ fn fixnan (x)
+ ? (x != x) 0.0 x
+
+ L := (normalize (vec2 1 0))
+
+ N := 1000
+ stp_phi := (/ N) * tau
+ local area = 0.0
+ local rad = 0.0
+ r := 0.5 * pi
+ k := (cos (r / 2.0))
+ # k = -1, n = 0
+ # k = -0.5, n = sqrt 3
+ # k = 0.0, n = sqrt 3
+ # for k = 0.707107, n = 2
+ print "k=" k
+ local nagg = (vec2 0)
+ for phi in (range N)
+ phi0 := phi as f32 * stp_phi
+ phi1 := phi0 + stp_phi
+ n0 := (vec2 (cos phi0) (sin phi0))
+ n1 := (vec2 (cos phi1) (sin phi1))
+ d01 := (n1 - n0)
+ A := (length d01)
+ n := (normalize (vec2 -d01.y d01.x))
+ #l := (max 0.0 (dot n L))
+ l := (dot n L) * 0.5 + 0.5
+ #l := (step 0.0 ((dot n L) - k))
+ nagg += n * l * A
+ rad += l * A
+ area += A
+ print "2d area" area
+ print "2d radiance" rad
+ print "2d normal" nagg
+ print "est" ((acos k) * 2.0) ((sin (r / 2.0)) * 2.0)
+
+do
+ # integrate cosine lobe 3d
+
+ fn triangle_normal (A B C)
+ normalize (cross (B - A) (C - A))
+ fn triangle_area (A B C)
+ (length (cross (B - A) (C - A))) / 2.0
+ fn triangle_area_normal (A B C)
+ (cross (B - A) (C - A)) / 2.0
+
+ fn fixnan (x)
+ ? (x != x) 0.0 x
+
+ L := (normalize (vec3 1 1 1))
+
+ N := 1000
+ stp_theta := (/ N) * pi
+ stp_phi := (/ N) * tau
+ local area = 0.0
+ local rad = 0.0
+ for theta in (range N)
+ theta0 := theta as f32 * stp_theta
+ theta1 := theta0 + stp_theta
+ for phi in (range N)
+ phi0 := phi as f32 * stp_phi
+ phi1 := phi0 + stp_phi
+ n00 := (sh1_angle_normal (vec2 theta0 phi0))
+ n10 := (sh1_angle_normal (vec2 theta1 phi0))
+ n11 := (sh1_angle_normal (vec2 theta1 phi1))
+ n01 := (sh1_angle_normal (vec2 theta0 phi1))
+ n1 := (triangle_area_normal n00 n10 n11)
+ n2 := (triangle_area_normal n11 n01 n00)
+ A1 := (length n1)
+ A2 := (length n2)
+ l1 := (max 0.0 (dot n1 L))
+ l2 := (max 0.0 (dot n2 L))
+ rad += l1 + l2
+ area += A1 + A2
+ print "3d area" (area / (2 * tau))
+ print "3d radiance" rad
+
+do
+ # monte carlo integration
+
+ fn sphere-area (u)
+ let phi =
+ * 2.0 pi u.x
+ let rho_c =
+ 2.0 * u.y - 1.0
+ let rho_s =
+ sqrt
+ 1.0 - (rho_c * rho_c)
+ vec3
+ rho_s * (cos phi)
+ rho_s * (sin phi)
+ rho_c
+
+ fn roberts-R2-quasirandom (n)
+ g := 1.32471795724474602596
+ Q := (vec2 (/ g) (/ (g * g)))
+ (0.5 + Q * n as f32) % 1.0
+
+ using import tukan.random
+
+ local rng : (Random)
+
+ fn triangle_normal (A B C)
+ normalize (cross (B - A) (C - A))
+ fn triangle_area (A B C)
+ (length (cross (B - A) (C - A))) / 2.0
+ fn triangle_area_normal (A B C)
+ (cross (B - A) (C - A)) / 2.0
+
+ fn fixnan (x)
+ ? (x != x) 0.0 x
+
+ L := (normalize (vec3 1 1 1))
+
+ N := 1000
+ local area = 0.0
+ local rad = 0.0
+ local SH = (vec4 0.0)
+
+ # integral of step(0, dot(n, L) - k)
+ # is A*(1-k)/2
+ # so a disc light source with radius cos(x) has power (1-cos(x))/2
+ # the function is symmetric around k=1
+ # a light source of solid angle R integrates to (1-cos(R/4))/2 units of light
+
+ r := 2 * pi
+ #(2 * pi * (1 - (cos (r/4))))
+ k := (cos (r / 4))
+
+ # -1.0 ~ 0
+ # -0.5 ~ 0.251112
+ # 0.0 ~ 0.5
+ # 0.5 ~ 0.75
+ # 1 ~ 1
+
+ print "k=" k
+ for s in (range N)
+ u := (roberts-R2-quasirandom s)
+ n := (sphere-area u)
+ #n := ('sphere-area rng)
+ #l := (max 0.0 (dot n L))
+ #l := (dot n L) * 0.5 + 0.5
+ l := (step 0.0 ((dot n L) - k))
+ rad += l
+ SH += (vec4 n 1) * l
+ factor := (2.0 * tau / N as f32)
+ SH0 := SH / N as f32
+ SH = SH * factor
+ print "3d area" (2 * tau)
+ print "3d radiance" (rad * factor) ((-k + 1) / 2)
+ est := (2 * pi * (1 - k))
+ est_n := pi * (1 - k * k)
+ print "estimated" est_n est
+ print "SH" SH (vec3 SH.xyz) (length (vec3 SH.xyz))
+ SH = (SH * SH1_Basis)
+ print "SH" SH (dot SH SH)
+
+fn ambient ()
+ vec4 0 0 0 (2 * tau)
+fn lambert (n)
+ vec4 (n * tau / 3.0) pi
+fn half-lambert (n)
+ vec4 (n * tau / 3.0) tau
+fn hemispheric (n)
+ vec4 (n * pi) tau
+fn light (n r)
+ """"r is the solid angle
+ k := (cos (r / 4))
+ vec4 (n * (pi * (1 - k * k))) (tau * (1 - k))
+
+inline steradians (a)
+ a / 180.0 * tau
+
+run-stage;
+
+# 0.0 is buttery, zero directional energy, good for
+# 0.7 is pretty good
+# 1.0
+SHSharpness := 1.0 # 2.0
+fn irradiance_probe (v)
+ sh_c0 := (2.0 - SHSharpness) * 1.0
+ sh_c1 := SHSharpness * 2.0 / 3.0
+ vec4 (v.xyz * sh_c1) (v.w * sh_c0)
+
+do
+ L := (light (vec3 0 0 1) (steradians 180.0))
+ LA := (ambient) #/ 0.2
+ fn test_probe (L n)
+ n := (normalize n)
+ S := (irradiance_probe (vec4 n 1))
+
+ print "L" L
+ print "S" S
+ s := (dot (L * SH1_Basis) (S * SH1_Basis))
+ print s
+
+ print "ambient light:"
+ test_probe LA (vec3 0 0 -1)
+ test_probe LA (vec3 0 1 0)
+ test_probe LA (vec3 0 0 1)
+
+ print "hemispheric light:"
+ test_probe L (vec3 0 0 -1)
+ test_probe L (vec3 0 0.5 -0.5)
+ test_probe L (vec3 0 1 0)
+ test_probe L (vec3 0 0.5 0.5)
+ test_probe L (vec3 0 0 1)
+