I got this equation I want to integrate over a sphere:

pow( clamp( dot(-vSurfaceToCameraDir, vLight + vNormal * C1), 0, 1 ), C2 )

C1 & C2 are constants and vSurfaceToCameraDir, vLight & vNormal are unit vectors that can point in any direction. If C1=0 this would add up to π/(C2 + 1). However now we have the 2nd spherical integration with vNormal.

Is there some analytical way to compute this? If not I suppose I can do a numeric integration by making lots of random samples but I thought there might be some more elegant & faster way.