#!BPY

"""
Name: 'Nebula Generator'
Blender: 248
Group: 'Misc'
Tooltip: 'Generate point clouds that resemble nebulae, clouds or even explosions'
"""

__author__ = "Alex 'CubOfJudahsLion' Feterman"
__url__ = (	"blender",\
			"Discussion @Elysiun, http://www.elysiun.com/forum/viewtopic.php?t=54758&sid=c04e554aa38f179a879e834155e443cf http://blenderartists.org/forum/showthread?t=53356")
__version__ = "237"

__bpydoc__ = """\
A tool to generate meshes that resemble clouds of various forms, mainly those produced
from a central source. Use to produce anything, from nebulae to explosions.

Usage:<br/>

Run the script in the script editing window.

Notes:<br/>

While this script generates a cloud and a texture when it's run, it can be used to
generate various kinds of clouds. Play with it to find out!
"""


# Copyright (C) 2005 Alexander Feterman
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.




import Blender, math, types, datetime
from Blender import BGL, Draw, Material, NMesh, Noise, Object, Scene, Window
from Blender.BGL import *


# ====== constants
invNorQ		= (2*math.pi)**0.5	# inverse of normal sigma=1 mu=1 quofficient
norQ		= 1.0 / invNorQ		# actual quofficient
G			= 6.67*10**-11		# universal gravitational constant


# ====== probability distributons

# A comment on the normal distribution:
# It is one of the most common probability distributions in everyday phenomena;
# it has been demonstrated that any probability distribution of means is a normal
# distribution, so it naturally applies to moving clouds of gas.

def Normal(x):
	"""
		This is the standard normal distribution
	"""
	global norQ
	return norQ*math.exp(-x*x/2.0)


def NormalP(x, mu, sigma):
	"""
		Parametrical normal distribution
	"""
	global norQ
	return norQ/sigma * math.exp((x-mu)*(mu-x)/(2*sigma*sigma))


def InvNormal(y, mu, sigma):
	"""
		Inverse of parametric normal distribution
	"""
	return sigma * math.sqrt(-2 * math.log(y)) + mu


# ====== vector class for simple point manipulation
class Vector:
	def __init__(self, x=0, y=0, z=0):
		self.x = x
		self.y = y
		self.z = z

	def __neg__(self):
		return Vector(-self.x, -self.y, -self.z)

	def __add__(self, v):
		return Vector(self.x + v.x, self.y + v.y, self.z + v.z)

	#__radd__ = __add__

	def __sub__(self, v):
		return Vector(self.x - v.x, self.y - v.y, self.z - v.z)

	def __mul__(self, p):
		typ = type(p)
		if typ == types.IntType or typ == types.LongType or typ == types.FloatType:
			return Vector(self.x * p, self.y * p, self.z * p)
		else:
			return Vector(self.x * p.x, self.y * p.y, self.z * p.z)

	__rmul__ = __mul__

	def __div__(self, p):
		typ = type(p)
		if typ == types.IntType or typ == types.LongType or typ == types.FloatType:
			return Vector(self.x / p, self.y / p, self.z / p)
		else:
			return Vector(self.x / p.x, self.y / p.y, self.z / p.z)

	def __rdiv__(self, p):
		return Vector(p / self.x, p / self.y, p / self.z)

	def __repr__(self):
		return "(" + str(self.x) + "," + str(self.y) + "," + str(self.z) + ")"

	__str__ = __repr__

	def radius(self):
		return (self.x * self.x + self.y * self.y + self.z * self.z) ** 0.5

	def distance(self, v):
		return (self-v).radius()


# a class to handle material colors as one item
class RGB:
	def __init__(self, r=0, g=0, b=0):
		self.r = r
		self.g = g
		self.b = b

	def __repr__(self):
		return "(" + str(self.r) + "," + str(self.g) + "," + str(self.b) + ")"

	__str__ = __repr__


# ====== material properties
class Gas:
	"""
		Stores properties of a gaseos substance in the nebula
	"""
	
	# molecular gas is generally stable, so it doesn't ionize and it's opaque
	MOLECULAR	= 0
	OPAQUE		= 0
	# whereas non-grouped elements are sensitive to ionization and thus irradiate
	IONIZED		= 1
	LUMINOUS	= 1

	def __init__(self, gasName, rgbColor, blendMode, relativeProportion, relativeWeight=1):
		"""
			gasName:			It's a name for the material, preferrably easily identifiable.
			rgbColor:			Median spectral color (translated to RGB) of the material.
								This color is supposed to be influenced by
			blendMode:			whether the material is molecular or ionized or not.
			relativeProportion:	abundance of the element in linear proportion of the others;
								it doesn't matter what unit you use as long as you use the
								same one for all others
			relativeWeight:		the weight of the material; since we're representing particles,
								it's not out of the question to represent the sum of atomic
								weights in the molecule or the atomic weight of single (hot and
								ionized) atoms.
		"""
		self.name		= gasName
		self.mode		= blendMode
		self.rgb		= rgbColor
		self.massMul	= relativeWeight
		self.proportion	= relativeProportion

	def __repr__(self):
		if(self.mode == Gas.LUMINOUS):
			blending = "* "
		else:
			blending = " "
		return self.name + blending + "x" + str(self.massMul) + ", c=" + str(self.rgb) + ", 1/" + str(1/self.proportion)

	__str__ = __repr__


# ====== material catalog
class NebulaMix:
	"""
		Description of a mixture of gases in a nebula
	"""

	def normalizeProportions(self):
		"""
			Changes the proportions of materials so the sum of the
			proportions equals 1
		"""
		proportionSum = 0
		for matIdx in range(len(self.materials)):
			proportionSum += self.materials[matIdx].proportion
		for matIdx in range(len(self.materials)):
			self.materials[matIdx].proportion /= proportionSum

	def __init__(self, *gasMix):
		"""
			The gas mix is a variant array of Gas elements.
		"""
		self.materials = gasMix
		# we normalize their proportions for random material generation
		self.normalizeProportions()

	def getRandomMaterialIndex(self):
		"""
			Gets a random element, following the probability distribution
			of the material.
		"""
		r, i, sum = Noise.random(), 0, self.materials[0].proportion
		while r > sum:
			i	+= 1;
			sum	+= self.materials[i].proportion
		return i


	def __len__(self): return len(self.materials)

	def __getitem__(self, key): return self.materials[key]

	def __repr__(self):
		matLen	= len(self.materials)
		sRep	= "[\n\t"
		for i in range(matLen):
			sRep += repr(self.materials[i])
			if i+1 < matLen:
				sRep += ",\n\t"
		sRep += "\n]"
		return sRep


# ====== nebulae creation+evolution functions
def add_fn(a,b): return a+b


def Evolve_Cuadratic(vertexes):
	"""
		This variant is the O(n*(n-1)/2) version, every point contributing
		against every other point--should run at least twice as fast
	"""
	# parameters used
	global maxD, totalPoints, evolveIter, useLineAlgo, turbDistr, probBranch
	# init
	speeds		= [Vector(0,0,0) for v in vertexes]
	cumVertexes	= []
	for e in range(evolveIter):
		if e % 10 == 0:
			#print "Iteration #", e
			Window.DrawProgressBar(0.99*e/evolveIter, "Iteration %d of %d" % (e,evolveIter))
		newVertexes	= []
		accel		= [Vector(0,0,0) for w in vertexes]
		i = 0
		vlen = len(vertexes)
		while i < vlen:
			for j in range(i+1, vlen):
				# this interaction allows 
				dist		=	vertexes[j].distance(vertexes[i])
				factor		=	(vertexes[j] - vertexes[i]) / (9000 * dist * dist * dist)
				accel[i]	+=	factor # * mass[j]
				accel[j]	-=	factor # * mass[i]
			# use previous speed and calculated accelereation to move particle
			iters	= 0
			cont	= True
			# now for the branches, having at least one main branch
			while cont:
				# get the turbulence
				turb	=	Noise.vTurbulence((vertexes[i].x,vertexes[i].y,vertexes[i].z,), 1, 0)
				# turn it into a vector and scale it according to the distribution
				# we don't want a zero here
				rDone	=	False
				while not rDone:
					norm	=	Vector(	InvNormal(Noise.random(), turbScale['mu'], turbScale['sigma']),\
										InvNormal(Noise.random(), turbScale['mu'], turbScale['sigma']),\
										InvNormal(Noise.random(), turbScale['mu'], turbScale['sigma']) )
					if norm.x and norm.y and norm.z:
						rDone = True
				# turn it into a vector and scale it according to the distribution
				vturb	=	Vector(turb[0],turb[1],turb[2])*norm
				# during one discretized interval (this)
				newVertexes.append(vertexes[i]+speeds[i]+accel[i]/2.0+vturb)
				# now for the next branch, which is a bernoulli experiment of
				# up to eight tries
				iters	+=	1
				if iters == 8 or iters + i == vlen or Noise.random() > probBranch:
					cont = False
			speeds[i]	+= accel[i]
			i			+= iters
		cumVertexes	+=	newVertexes
		vertexes	=	newVertexes
	return cumVertexes


def Evolve_Linear(vertexes):
	"""
		This function uses vectorial gravitational acceleration
		to evolve the system through mutual attraction between particles:

					G m1 (^v1-^v2) / |^v1-^v2|**3

		Where G ~= 6.67x10**-11. Assuming that the
		mass of each particle is as large as xG**-1, the final equation is:

					x * (^v1 - ^v2) / |^v1 - ^v2| ** 3

		An exhaustive method is O(n**2), but we won't calculate attraction
		from each particle to every other particle. Instead, we'll simplify
		and calculate and approximation through cumulative centers of mass
		which will solve our problem in O(kn) time.

		The formula for the center of mass is:

					SUM m[i] * pos[i] for all i/SUM m[i] for all i
	"""
	# parameters used
	global maxD, totalPoints, evolveIter, useLineAlgo, turbDistr, probBranch
	# init
	vlen		= len(vertexes)
	speeds		= [Vector(0,0,0) for v in vertexes]
	cumVertexes = []
	for e in range(evolveIter):
		if e % 10 == 0:
			#print "Iteration #", e
			Window.DrawProgressBar(0.99*e/evolveIter, "Iteration %d of %d" % (e,evolveIter))
		# first: calculate join mass and center of mass
		center	= Vector(0,0,0)
		mass	= 0
		for i in range(vlen):
			mass	+= 1.0 # mass[i]
			center	+= vertexes[i] # * mass[i]
		# we'll save the final division (center /= mass) for later
		# now let's go through vertexes to change them
		newVertexes = []
		i = 0
		while i < vlen:
			# the plan is to exclude this particle from the center of mass calculations
			massEx		= mass - 1.0 # * x * mass[i]
			# observe that it's easier to calculate the center
			# that excludes this mass without the division
			centerEx	= (center - vertexes[i]) / massEx
			dist		= centerEx.distance(vertexes[i])
			# now for the branches, having at least one main branch			
			accel			= massEx * (centerEx - vertexes[i]) / (9000 * dist * dist * dist)
			#print "accel =", accel
			# get the turbulence
			turb		=	Noise.vTurbulence((vertexes[i].x,vertexes[i].y,vertexes[i].z,), 1, 0)
			# turn it into a vector and scale it according to the distribution
			# we don't want a zero here
			rDone	=	False
			while not rDone:
				norm	=	Vector(	InvNormal(Noise.random(), turbDistr['mu'], turbDistr['sigma']),\
									InvNormal(Noise.random(), turbDistr['mu'], turbDistr['sigma']),\
									InvNormal(Noise.random(), turbDistr['mu'], turbDistr['sigma']) )
				if norm.x and norm.y and norm.z:
					rDone = True
			# turn it into a vector and scale it according to the distribution
			vturb	=	Vector(turb[0],turb[1],turb[2])*norm
			#accel	+=	Vector(turb[0],turb[1],turb[2])*norm
			# use previous speed and calculated accelereation to move particle
			# during one discretized interval (this)
			newVertexes.append(vertexes[i]+speeds[i]+accel/2.0+vturb)
			speeds[i]	+= accel
			i			+= 1
		cumVertexes	+=	newVertexes
		vertexes	=	newVertexes
	return cumVertexes


def MakeNebula():
	"""
		We use a voronoi distance to displace the point and
		use normalized multiplicative turbulence to generate unique,
		sort of wispy, uniformly distributed points.
		Then we evolve the system using gravity and turbulence
		
		maxD		:	Maximum average distance from the origin.
		nebulaMats	:	A NebulaMix containing the description of gases in the nebula.
		totalPoints	:	An approximate point count. Please note that this process is
						highly dependent on probabilistic models and therefore this is
						only a relative goal; the ending cloud can be twice as dense.
		evolveIter	:	Number of iterations to evolve the system (through turbulence and
						gravitational attraction.)
		useLineAlgo	:	There are two algorithms to calculate particle motion: one based on
						centers of mass, which runs in linear (O(kn)) time and another that
						calculates the attraction to every other particle and runs in
						O(kn(n-1)/2). By default, this is true.
		turbDistr	:	Turbulence distribution. Even turbulence operates
						probabilistically here (using a Normal distribution, since we
						model statistical behavior of gas.) It's a dictionary with two
						values: average force of the turbulence (mu) and range of
						variation around the mean (sigma.)
		probBranch	:	We model dense clouds of gases that move about as the result of
						turbulent solar winds, naturally leaving traces behind itself.
						This is the probability that the cloud will branch.
	"""
	# parameters used
	global maxD, totalPoints, evolveIter, useLineAlgo, turbDistr, probBranch
	# start status reports
	print "\n\n\n*** Cloud Generator***\\nn"
	Window.DrawProgressBar(0.0, "Generating cloud")
	print "Generating cloud..."
	# a branch can become up to 8 branches depending on a 'toss of a coin' type
	# of repeated experiment, each iteration allowing another branch
	# let's calculate the median of branches
	probPoints = 0
	for iToss in range(1,9):
		probPoints += iToss * probBranch**(iToss-1)
	probPoints *= 1-probBranch
	# we divide by 4 because the voronoi noise function returns four
	# vectors which we use, then by evolveIter because the evolver
	# creates new points the way a flame fractal would
	pointsPerDimension = (((totalPoints+2.9999)/3.0+evolveIter-1)/evolveIter/probPoints)**(1/3.0)
	#print "ppd:", pointsPerDimension
	# prepare randomization
	Noise.setRandomSeed(0)
	# create object and mesh
	meshObj	= Object.New("Mesh", "Nebula")
	mesh	= NMesh.New("Nebula.mesh")
	meshObj.link(mesh)
	# add random vertexes
	step, cnt	= (maxD*1.01)/pointsPerDimension, 0
	vertexes	= []
	speeds		= []
	x = -maxD/2.0
	#statV = Vector(0,0,0)
	while x < maxD/2.0:
		y = -maxD/2.0
		while y < maxD/2.0:
			z = -maxD/2.0
			while z < maxD/2.0:
				#cnt += 1
				xt, yt, zt	= Noise.vTurbulence((x,y,z,), 2, 0, Noise.NoiseTypes.VORONOI_CRACKLE)
				thread		= Noise.voronoi((x,y,z,))[1]
				sliding		= 0.0
				for iT in range(3):
					xv, yv, zv	=	[abs(vc) for vc in thread[iT]]
					finalV		=	Vector(x*xv*xt, y*yv*yt, z*zv*zt)
					vertexes.append(finalV)
				z += step
			y += step
		x += step
	timeStart = datetime.datetime.now()
	Window.DrawProgressBar(0.001, "Evolving cloud")
	print "Evolving cloud..."
	if useLineAlgo:
		vertexes = Evolve_Linear(vertexes)
	else:
		vertexes = Evolve_Cuadratic(vertexes)
	timeEnd		= datetime.datetime.now()
	timeDiff	= timeEnd - timeStart
	print 'Time taken for evolution:', timeDiff
	for iVert in range(len(vertexes)):
		v = vertexes[iVert]
		mesh.verts.append(NMesh.Vert(v.x, v.y, v.z))
	mesh.update()
	# once created, evolve them
	#mat			=   Material.New("Nebula.mat")
	#mat.mode	=   Material.Modes["HALO"]		| Material.Modes["ZTRANSP"]      \
	#			|   Material.Modes["TRACEABLE"]	| Material.Modes["HALOSHADE"]   \
	#			|   Material.Modes["SHADOW"]	| Material.Modes["VCOL_PAINT"]
	#mat.alpha      = 0.0625
	#mat.add         = 0.25
	#mat.haloSize   = step*3.5
	# and put it in current scene
	Scene.GetCurrent().link(meshObj)
	Window.DrawProgressBar(1.0, "Done.")


# our material catalogs should imitate the relative abundance of the color-yielding or absorbing elements in emissive and reflective nebulae
realSpectrum = NebulaMix(	Gas("H+", RGB(0x30, 0, 0), Gas.IONIZED, 600, 1), \
							Gas("H2", RGB(0xC0, 0x27, 0x12), Gas.MOLECULAR, 200, 2), \
							Gas("He", RGB(0,0,0), Gas.MOLECULAR, 100, 4), \
							Gas("OIII", RGB(0, 0xFF, 0x20), Gas.IONIZED, 100, 15), \
							Gas("NII", RGB(0, 0, 0xFF), Gas.IONIZED, 5, 14), \
							Gas("SII", RGB(0x75, 0, 0), Gas.IONIZED, 5, 32)	)


# ====== parameters

maxD			= 1.0
useLineAlgo		= True
totalPoints		= 10000
turbMu			= 0.0
turbSigma		= 1.0
turbDistr		= {'mu':turbMu, 'sigma':turbSigma}
probBranch		= 0.5
evolveIter		= 100

id_maxD			= 1001
id_useLineAlgo	= 1002
id_totalPoints	= 1003
id_turbMu		= 1004
id_turbSigma	= 1005
id_probBranch	= 1006
id_evolveIter	= 1007
id_Go			= 1010
id_Exit			= 1101


# ====== main functions
def ButtonEventHandler(buttonId):
	global maxD, useLineAlgo, totalPoints, turbMu, turbSigma, turbDistr, probBranch, evolveIter
	global id_maxD, id_totalPoints, id_evolveIter, id_turbMu, id_turbSigma, id_useLineAlgo, id_probBranch, id_Exit, id_Go
	global btn_evolveIter, btn_useLineAlgo, btn_maxD, btn_totalPoints, btn_turbMu, btn_turbSigma, btn_probBranch
	if buttonId == id_Exit:
		Draw.Exit()
	elif buttonId == id_Go:
		MakeNebula()
	elif buttonId == id_maxD:
		maxD = btn_maxD.val
	elif buttonId == id_totalPoints:
		totalPoints = btn_totalPoints.val
	elif buttonId == id_evolveIter:
		evolveIter = btn_evolveIter.val
	elif buttonId == id_turbMu:
		turbMu = btn_turbMu.val
		turbDistr = {'mu':turbMu, 'sigma':turbSigma}
	elif buttonId == id_turbSigma:
		turbSigma = btn_turbSigma.val
		turbDistr = {'mu':turbMu, 'sigma':turbSigma}
	elif buttonId == id_useLineAlgo:
		useLineAlgo = btn_useLineAlgo.val == 1
	elif buttonId == id_probBranch:
		probBranch = btn_probBranch.val


def DrawHandler():
	"""
		Draws interface
	"""
	# ids required for each button
	global maxD, useLineAlgo, totalPoints, turbMu, turbSigma, turbDistr, probBranch, evolveIter
	global id_maxD, id_totalPoints, id_evolveIter, id_turbMu, id_turbSigma, id_useLineAlgo, id_probBranch, id_Exit, id_Go
	global btn_evolveIter, btn_useLineAlgo, btn_maxD, btn_totalPoints, btn_turbMu, btn_turbSigma, btn_probBranch
	# get the size of the area
	aW,aH = Window.GetAreaSize()
	# black background
	glClearColor(0.0, 0.0, 0.0, 0.0)
	glClear(GL_COLOR_BUFFER_BIT)
	# orange title bar with bold black text
	glColor3f(1,0.4,0)
	glRecti(0,aH,250,aH-16)
	aH -= 16
	nebGenTitle = "Nebula Generator"
	glColor3f(0,0,0)
	glRasterPos2i(0, aH+2)
	Draw.Text(nebGenTitle,'large')
	glRasterPos2i(1, aH+2)
	Draw.Text(nebGenTitle,'large')
	# and a close button
	Xwidth =  Draw.GetStringWidth("X")+7
	Draw.PushButton("X", id_Exit, 250-Xwidth, aH, Xwidth, 14, "Exit")
	# now for the controls
	# first the maximum distance
	aH -= 16
	btn_maxD = Draw.Slider("InitialSpread", id_maxD, 0, aH, 250, 14, maxD, 0.1, 20.0, 0, "Approximate dimension of cube area for the initial super-dense cloud configuration")
	# fast/slow method
	aH -= 16
	btn_useLineAlgo = Draw.Menu("Algorithm%t|Fast O(n)%x1|Thorough O(n*(n-1)/2)%x2",\
								id_useLineAlgo, 0, aH, 250, 14, 1,\
								"'Fast' uses approximated centers of mass. "+
								"'Thorough' calculates the effect of every particle on every other one. "+\
								"Normally, the fast algorithm is sufficient.")
	# total point count
	aH -= 16
	btn_totalPoints = Draw.Number(	"Particles", id_totalPoints, 0, aH, 250, 14, totalPoints, 100, 1000000,\
								 	"Approximate final number of particles in the cloud. "+
									"Since this is a probabilisitc model, an exact "+\
									"result cannot be guaranteed.")
	# branching probability
	aH -= 16
	btn_probBranch = Draw.Slider(	"ProbBranch", id_probBranch, 0, aH, 250, 14, probBranch, 0.0, 1.0, 0,\
									"Probability that a dense cloud divides in two")
	# number of iterations
	aH -= 16
	btn_evolveIter = Draw.Number("Iterations", id_evolveIter, 0, aH, 250, 14, evolveIter, 10, 500000, "Number of cycles of cloud evolution")
	# turbulence parameters
	glColor3f(0.5,0.2,0)
	glRecti(0, aH, 250, aH-16)
	aH -= 16
	glColor3f(0.75,0.75,0.75)
	glRasterPos2i(0,aH+2)
	Draw.Text("Turbulence")
	aH -= 16
	btn_turbMu = Draw.Number("Mean", id_turbMu, 0, aH, 125, 14, turbMu, 0.0, 10.0, "Typical strength of turbulence of the solar winds")
	btn_turbSigma = Draw.Number("Variation", id_turbSigma, 126, aH, 125, 14, turbSigma, 0.01, 10.0, "Average range of turbulence of the solar winds")
	# of course we need a 'go' button
	glColor3f(0.5,0.2,0)
	glRecti(0, aH, 250, aH-16)
	aH -= 16
	Draw.PushButton("Make Cloud", id_Go, 20, aH, 210, 14, "Exits and starts the simulation to create the cloud")


# go!
Draw.Register(DrawHandler, None, ButtonEventHandler)