No big news yet, I've been studying a little bit analytic methods for shading and occlusion, I can't report anything really now, even because I'm not yet satisfied with what I've done.

But I'd like to share this link: http://www.me.utexas.edu/~howell/tablecon.html, differential element to finite area is what you’ll need.

Also, you might find this useful, if you're starting to play around with spherical harmonics, it's a small snipped of what I've been doing, with Mathematica:

(* analytic solution for real spherical harmonics test *)

shIndices[level_] := (Range[-#1, #1] & ) /@ Range[0, level]

shGetNormFn[l_, m_] := Sqrt[((2*l + 1)*(l - m)!)/(4*Pi*(l + m)!)]

shGetFn[l_, m_] :=

Piecewise[{{shGetNormFn[l, 0]*LegendreP[l, 0, Cos[\[Theta]]],

m == 0}, {Sqrt[2]*shGetNormFn[l, m]*Cos[m*\[Phi]]*

LegendreP[l, m, Cos[\[Theta]]],

m > 0}, {Sqrt[2]*shGetNormFn[l, -m]*Sin[(-m)*\[Phi]]*

LegendreP[l, -m, Cos[\[Theta]]], m <>

shFunctions[level_] :=

MapIndexed[

Function[{list, currlevel}, (shGetFn[currlevel - 1, #1] & ) /@

list], shIndices[level]]

shGenCoeffs[shfns_, fn_] :=

Map[Integrate[#1*fn[\[Theta], \[Phi]]*Sin[\[Theta]], {\[Theta], 0,

Pi}, {\[Phi], 0, 2*Pi}] & , shfns, {2}]

shReconstruct[shfns_, shcoeffs_] :=

Simplify[Plus @@ (Flatten[shcoeffs]*Flatten[shfns]),

Assumptions -> {Element[\[Theta], Reals],

Element[\[Phi], Reals], \[Theta] >= 0, \[Phi] >= 0, \[Theta] <=

Pi, \[Phi] <= 2*Pi}]

shIsZonal[shcoeffs_, level_] :=

Plus @@ (Flatten[shIndices[level]] Flatten[shcoeffs]) == 0

shGetSymConvolveNorm[level_] :=

MapIndexed[

Function[{list, currlevel},

Table[Sqrt[(4 \[Pi])/(2 currlevel + 1)], {Length[list]}]],

shIndices[level]]

shGetSymCoeffs[shcoeffs_] :=

Table[#1[[Ceiling[Length[#1]/2]]], {Length[#1]}] & /@ shcoeffs

shSymConvolve[shcoeffs_, shsymkerncoeffs_,

level_] := (Check[shIsZonal[shsymkerncoeffs], err];

shGetSymConvolveNorm[level] shcoeffs shGetSymCoeffs[

shsymkerncoeffs])

(* tests.... *)

testnumlevels = 2

testfn[a_, b_] :=

Cos[a]^10*UnitStep[Cos[a]] (*symmetric on the z axis*)

(*testfn[a_,b_]:= (a/Pi)^4*)

shfns = shFunctions[testnumlevels]

testfncoeffs = shGenCoeffs[shfns, testfn]

shIsZonal[testfncoeffs, testnumlevels]

testfnrec = {\[Theta], \[Phi]} \[Function]

Evaluate[shReconstruct[shfns, testfncoeffs]]

SphericalPlot3D[{testfn[\[Theta], \[Phi]],

testfnrec[\[Theta], \[Phi]]}, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},

Mesh -> False, PlotRange -> Full]

testfn2[a_, b_] := UnitStep[Cos[a] Sin[b]](*asymmetric*)

testfn2coeffs = shGenCoeffs[shfns, testfn2]

testfn3coeffs =

shSymConvolve[testfn2coeffs, testfncoeffs, testnumlevels]

testfn2rec = {\[Theta], \[Phi]} \[Function]

Evaluate[shReconstruct[shfns, testfn2coeffs]]

testfn3rec = {\[Theta], \[Phi]} \[Function]

Evaluate[shReconstruct[shfns, testfn3coeffs]]

SphericalPlot3D[{testfn2[\[Theta], \[Phi]],(*testfn2rec[\[Theta],\

\[Phi]],*)testfn3rec[\[Theta], \[Phi]]}, {\[Theta], 0, Pi}, {\[Phi],

0, 2 Pi}, Mesh -> False, PlotRange -> Full]

But I'd like to share this link: http://www.me.utexas.edu/~howell/tablecon.html, differential element to finite area is what you’ll need.

Also, you might find this useful, if you're starting to play around with spherical harmonics, it's a small snipped of what I've been doing, with Mathematica:

(* analytic solution for real spherical harmonics test *)

shIndices[level_] := (Range[-#1, #1] & ) /@ Range[0, level]

shGetNormFn[l_, m_] := Sqrt[((2*l + 1)*(l - m)!)/(4*Pi*(l + m)!)]

shGetFn[l_, m_] :=

Piecewise[{{shGetNormFn[l, 0]*LegendreP[l, 0, Cos[\[Theta]]],

m == 0}, {Sqrt[2]*shGetNormFn[l, m]*Cos[m*\[Phi]]*

LegendreP[l, m, Cos[\[Theta]]],

m > 0}, {Sqrt[2]*shGetNormFn[l, -m]*Sin[(-m)*\[Phi]]*

LegendreP[l, -m, Cos[\[Theta]]], m <>

shFunctions[level_] :=

MapIndexed[

Function[{list, currlevel}, (shGetFn[currlevel - 1, #1] & ) /@

list], shIndices[level]]

shGenCoeffs[shfns_, fn_] :=

Map[Integrate[#1*fn[\[Theta], \[Phi]]*Sin[\[Theta]], {\[Theta], 0,

Pi}, {\[Phi], 0, 2*Pi}] & , shfns, {2}]

shReconstruct[shfns_, shcoeffs_] :=

Simplify[Plus @@ (Flatten[shcoeffs]*Flatten[shfns]),

Assumptions -> {Element[\[Theta], Reals],

Element[\[Phi], Reals], \[Theta] >= 0, \[Phi] >= 0, \[Theta] <=

Pi, \[Phi] <= 2*Pi}]

shIsZonal[shcoeffs_, level_] :=

Plus @@ (Flatten[shIndices[level]] Flatten[shcoeffs]) == 0

shGetSymConvolveNorm[level_] :=

MapIndexed[

Function[{list, currlevel},

Table[Sqrt[(4 \[Pi])/(2 currlevel + 1)], {Length[list]}]],

shIndices[level]]

shGetSymCoeffs[shcoeffs_] :=

Table[#1[[Ceiling[Length[#1]/2]]], {Length[#1]}] & /@ shcoeffs

shSymConvolve[shcoeffs_, shsymkerncoeffs_,

level_] := (Check[shIsZonal[shsymkerncoeffs], err];

shGetSymConvolveNorm[level] shcoeffs shGetSymCoeffs[

shsymkerncoeffs])

(* tests.... *)

testnumlevels = 2

testfn[a_, b_] :=

Cos[a]^10*UnitStep[Cos[a]] (*symmetric on the z axis*)

(*testfn[a_,b_]:= (a/Pi)^4*)

shfns = shFunctions[testnumlevels]

testfncoeffs = shGenCoeffs[shfns, testfn]

shIsZonal[testfncoeffs, testnumlevels]

testfnrec = {\[Theta], \[Phi]} \[Function]

Evaluate[shReconstruct[shfns, testfncoeffs]]

SphericalPlot3D[{testfn[\[Theta], \[Phi]],

testfnrec[\[Theta], \[Phi]]}, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},

Mesh -> False, PlotRange -> Full]

testfn2[a_, b_] := UnitStep[Cos[a] Sin[b]](*asymmetric*)

testfn2coeffs = shGenCoeffs[shfns, testfn2]

testfn3coeffs =

shSymConvolve[testfn2coeffs, testfncoeffs, testnumlevels]

testfn2rec = {\[Theta], \[Phi]} \[Function]

Evaluate[shReconstruct[shfns, testfn2coeffs]]

testfn3rec = {\[Theta], \[Phi]} \[Function]

Evaluate[shReconstruct[shfns, testfn3coeffs]]

SphericalPlot3D[{testfn2[\[Theta], \[Phi]],(*testfn2rec[\[Theta],\

\[Phi]],*)testfn3rec[\[Theta], \[Phi]]}, {\[Theta], 0, Pi}, {\[Phi],

0, 2 Pi}, Mesh -> False, PlotRange -> Full]