a866b82c90e4 — Leonard Ritter 9 days ago
* check-in BCC test cube
4 files changed, 740 insertions(+), 7 deletions(-)

M lib/tukan/dsp/utils.sc
M lib/tukan/random.sc
A => testing/bcc_shearing.sc
A => testing/solid_angle.sc
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)
+