#!/usr/bin/env python 

import math,random
 
class SHSample():
	sph=(0.0,0.0)
	vec=(0.0,0.0,0.0)
	coeff={}
	pass
 
def P(l,m,x):
	#Associated Legendre Polynomial P(l,m,x) at x
	pmm = 1.0
	if m > 0:
		somx2=math.sqrt(1.0-(x*x))
		fact = 1.0
		for i in xrange(1,m+1,1):
			pmm *= (-fact)*somx2
			fact += 2.0
	if l == m:
		return pmm
 
	pmmp1 = x * ((2.0*m)+1.0)*pmm
 
	if l == m+1:
		return pmmp1
 
	plm = 0.0
 
	for ll in xrange(m+2,l+1,1):
		plm = ((2.0*ll - 1.0) * x * pmmp1 - (ll + m - 1.0) * pmm) / (ll - m);
		pmm = pmmp1
		pmmp1 = plm
 
	return plm
 
def K(l,m):
	temp = float((((2.0*l)+1.0)*math.factorial(l-m))/(4.0*math.pi*math.factorial(l+m)))
	return math.sqrt(temp)
 
def SH(l,m,theta,phi):
	sqrt2 = math.sqrt(2.0)
	if m==0:
		return K(l,0)*P(l,0,math.cos(theta))
	elif m > 0:
		return sqrt2*K(l,m)*math.cos(m*phi)*P(l,m,math.cos(theta))
	else:
		return sqrt2*K(l,-m)*math.sin(-m*phi)*P(l,-m,math.cos(theta))
 
 
def setupSamples(sqrtNumSamples=16,numBands=4):
	counter = 0
	oneOverN = 1.0/float(sqrtNumSamples)
	samples = [SHSample() for i in range(sqrtNumSamples*sqrtNumSamples)]
 
	for i in range(sqrtNumSamples):
		for j in range(sqrtNumSamples):
			x = (i+ random.random())*oneOverN
			theta = 2.0*math.acos(math.sqrt(1-x))
 
			y = (j+ random.random())*oneOverN
			phi = 2.0*math.pi*y
 
			samples[counter].sph=(theta,phi)
 
			vec = (math.sin(theta)*math.cos(phi),\
			math.sin(theta)*math.sin(phi),\
			math.cos(theta))
			samples[counter].vec = vec
 
			tmpDict = {}
			for l in range(numBands):
				for m in xrange(-l,l+1,1):
					index = l*(l+1)+m
					sh= SH(l,m,theta,phi)
					tmpDict[index]=sh
			samples[counter].coeff=tmpDict
			counter +=1
	return samples
 
sampleCounter = 0
 
totalSamples = 1024
numBands = 4

sqrtSamples = int(math.sqrt(totalSamples))

labels = {\
"vector":"uniform vector samplesVector[%s]" % totalSamples,
"sph0":"uniform float samplesSph0[%s]" % totalSamples,
"sph1":"uniform float samplesSph1[%s]" % totalSamples,
"coeff0":"uniform float samplesCoeffs0[%s]" % totalSamples,
"coeff1":"uniform float samplesCoeffs1[%s]" % totalSamples,
"coeff2":"uniform float samplesCoeffs2[%s]" % totalSamples,
"coeff3":"uniform float samplesCoeffs3[%s]" % totalSamples,
"coeff4":"uniform float samplesCoeffs4[%s]" % totalSamples,
"coeff5":"uniform float samplesCoeffs5[%s]" % totalSamples,
"coeff6":"uniform float samplesCoeffs6[%s]" % totalSamples,
"coeff7":"uniform float samplesCoeffs7[%s]" % totalSamples,
"coeff8":"uniform float samplesCoeffs8[%s]" % totalSamples,
}

samples = setupSamples(sqrtNumSamples=sqrtSamples,numBands=numBands)

print "#define SHVECTOR %s = %s;\n\n" % (labels["vector"],"{%s}" % ",".join([str("vector \"world\" " + str(i.vec)) for i in samples]))
print "#define SHSPH0 %s = %s;\n\n" % (labels["sph0"],"{%s}" % ",".join([str(i.sph[0]) for i in samples]))
print "#define SHSPH1 %s = %s;\n\n" % (labels["sph1"],"{%s}" % ",".join([str(i.sph[1]) for i in samples]))
print "#define SHCOEFF0 %s = %s;\n\n" % (labels["coeff0"],"{%s}" % ",".join([str(i.coeff[0]) for i in samples]))
print "#define SHCOEFF1 %s = %s;\n\n" % (labels["coeff1"],"{%s}" % ",".join([str(i.coeff[1]) for i in samples]))
print "#define SHCOEFF2 %s = %s;\n\n" % (labels["coeff2"],"{%s}" % ",".join([str(i.coeff[2]) for i in samples]))
print "#define SHCOEFF3 %s = %s;\n\n" % (labels["coeff3"],"{%s}" % ",".join([str(i.coeff[3]) for i in samples]))
print "#define SHCOEFF4 %s = %s;\n\n" % (labels["coeff4"],"{%s}" % ",".join([str(i.coeff[4]) for i in samples]))
print "#define SHCOEFF5 %s = %s;\n\n" % (labels["coeff5"],"{%s}" % ",".join([str(i.coeff[5]) for i in samples]))
print "#define SHCOEFF6 %s = %s;\n\n" % (labels["coeff6"],"{%s}" % ",".join([str(i.coeff[6]) for i in samples]))
print "#define SHCOEFF7 %s = %s;\n\n" % (labels["coeff7"],"{%s}" % ",".join([str(i.coeff[7]) for i in samples]))
print "#define SHCOEFF8 %s = %s;\n\n" % (labels["coeff8"],"{%s}" % ",".join([str(i.coeff[8]) for i in samples]))