#
# create lookup textures for Perez allweather sky model
#
# Author: Alex V. Boreskoff <steps3d@narod.ru>, <alexboreskoff@mtu-net.ru>
#

import PIL, Image, math

# compute zenith Yxy values

def getZenith ( t, thetaSun ):
	t2  = t*t						# turbidity squared
	chi = (4.0 / 9.0 - t / 120.0 ) * (math.pi - 2.0 * thetaSun )
	theta2 = thetaSun * thetaSun
	theta3 = thetaSun * theta2
	Y   = (4.0453 * t - 4.9710) * math.tan ( chi ) - 0.2155 * t + 2.4192

	x   = ( 0.00165 * theta3 - 0.00375 * theta2 + 0.00209 * thetaSun + 0.0)  * t2 + \
	      (-0.02903 * theta3 + 0.06377 * theta2 - 0.03202 * thetaSun + 0.00394) * t +\
	      ( 0.11693 * theta3 - 0.21196 * theta2 + 0.06052 * thetaSun + 0.25886)

	y = ( 0.00275 * theta3 - 0.00610 * theta2 + 0.00317 * thetaSun + 0.0)     * t2 +\
	    (-0.04214 * theta3 + 0.08970 * theta2 - 0.04153 * thetaSun + 0.00516) * t +\
	    ( 0.15346 * theta3 - 0.26756 * theta2 + 0.06670 * thetaSun + 0.26688)

	return [Y, y, x]

def perezFunc ( a, b, c, d, e, theta, gamma ):
	ct = math.cos ( theta )
	cg = math.cos ( gamma )
	return (1.0 + a * math.exp(b/ct)) * (1.0 + c * math.exp(d * gamma) + e*cg*cg)

									# now write texture for Yxy-zenith ( turbidity, cos(thetaSun) )
def writeZenithTable ( name, size ):
	image   = Image.new ( "RGB", (size, size) )
	invSize = 1.0 / float ( size )

	for i in range ( size ):
		for j in range ( size ):
			turbidity = 1.0 + 30.0 * i * invSize
			cosTheta  = i * invSize
			thetaSun  = math.acos ( cosTheta )
			value     = getZenith ( turbidity, thetaSun )

									# rescale Y
			value [0] = 1.0 - math.exp ( -math.fabs ( value [0] ) / 25.0 )

									# encode as RGB
			color = ( int ( 255*value [0] ), int ( 255*value [1] ), int ( 255 * value [2] ) )

			image.putpixel ( (i,j), color )

			print value

	image.save ( name, "PNG" )

									# write table perezFunc ( cosTheta, cosGamma, turbidity )
def writePerezFunc ( name, size ):
	image   = Image.new ( "RGB", (size, size) )
	invSize = 1.0 / float ( size )
	value   = range ( 3 )

	print "*************Perez func******************"
	for i in range ( size ):
		t  = 1.0 + 30 * invSize

		aY =  0.17872 * t - 1.46303
		bY = -0.35540 * t + 0.42749
		cY = -0.02266 * t + 5.32505
		dY =  0.12064 * t - 2.57705
		eY = -0.06696 * t + 0.37027

		ax = -0.01925 * t - 0.25922
		bx = -0.06651 * t + 0.00081
		cx = -0.00041 * t + 0.21247
		dx = -0.06409 * t - 0.89887
		ex = -0.00325 * t + 0.04517

		ay = -0.01669 * t - 0.26078
		by = -0.09495 * t + 0.00921
		cy = -0.00792 * t + 0.21023
		dy = -0.04405 * t - 1.65369
		ey = -0.01092 * t + 0.05291

		for j in range ( size ):
			for k in range ( size ):
				cosTheta  = float ( j + 1 )  / float ( size + 1 )
				cosGamma  = float ( k + 1 )  / float ( size + 1 )
				theta     = math.acos ( cosTheta )
				gamma     = math.acos ( cosGamma )
				value [0] = perezFunc ( aY, bY, cY, dY, eY, theta, gamma )
				value [1] = perezFunc ( ax, bx, cx, dx, ex, theta, gamma )
				value [2] = perezFunc ( ay, by, cy, dy, ey, theta, gamma )

									# rescale Y
#				value [0] = 1.0 - math.exp ( -math.fabs ( value [0] ) / 25.0 )
#
									# fix xy overflow
#				if value [1] > 1.0:
#					value [1] = 1.0

#				if value [2] > 1.0:
#					value [2] = 1.0

									# encode as RGB
				color = ( int ( 255*value [0] ), int ( 255*value [1] ), int ( 255 * value [2] ) )

				image.putpixel ( (j,k), color )

				print value

		image.save ( name % i, "PNG" )

writeZenithTable ( 'zenith.png',     64 )
writePerezFunc   ( 'perez-%02d.png', 32 )