Google Summer of Code project: "Smoke Simulator Improvements & Fire".

Documentation & Test blend files:
------------------
http://wiki.blender.org/index.php/User:MiikaH/GSoC-2012-Smoke-Simulator-Improvements

Credits:
------------------
Miika Hamalainen (MiikaH): Student / Main programmer

Daniel Genrich (Genscher): Mentor / Programmer of merged patches from Smoke2 branch
Google: For Google Summer of Code 2012

11 files changed, 1231 insertions, 242 deletions

diff --git a/intern/smoke/CMakeLists.txt b/intern/smoke/CMakeLists.txt
index 9524f7b2935..11f173a6835 100644
--- a/intern/smoke/CMakeLists.txt
+++ b/intern/smoke/CMakeLists.txt
@@ -40,6 +40,7 @@ set(SRC
 	intern/FLUID_3D_SOLVERS.cpp
 	intern/FLUID_3D_STATIC.cpp
 	intern/LU_HELPER.cpp
+	intern/spectrum.cpp
 	intern/SPHERE.cpp
 	intern/WTURBULENCE.cpp
 	intern/smoke_API.cpp
@@ -53,6 +54,7 @@ set(SRC
 	intern/LU_HELPER.h
 	intern/MERSENNETWISTER.h
 	intern/OBSTACLE.h
+	intern/spectrum.h
 	intern/SPHERE.h
 	intern/VEC3.h
 	intern/WAVELET_NOISE.h
diff --git a/intern/smoke/extern/smoke_API.h b/intern/smoke/extern/smoke_API.h
index a0eb1bf38e0..98c63f08fbd 100644
--- a/intern/smoke/extern/smoke_API.h
+++ b/intern/smoke/extern/smoke_API.h
@@ -37,17 +37,23 @@ extern "C" {
 
 struct FLUID_3D;
 
-// export
-void smoke_export(struct FLUID_3D *fluid, float *dt, float *dx, float **dens, float **densold, float **heat, float **heatold, float **vx, float **vy, float **vz, float **vxold, float **vyold, float **vzold, unsigned char **obstacles);
-
 // low res
-struct FLUID_3D *smoke_init(int *res, float *p0, float dtdef);
+struct FLUID_3D *smoke_init(int *res, float dx, float dtdef, int use_heat, int use_fire, int use_colors);
 void smoke_free(struct FLUID_3D *fluid);
 
-void smoke_initBlenderRNA(struct FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli);
-void smoke_step(struct FLUID_3D *fluid, float dtSubdiv);
+void smoke_initBlenderRNA(struct FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli, float *burning_rate,
+						  float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp);
+void smoke_step(struct FLUID_3D *fluid, float gravity[3], float dtSubdiv);
 
 float *smoke_get_density(struct FLUID_3D *fluid);
+float *smoke_get_flame(struct FLUID_3D *fluid);
+float *smoke_get_fuel(struct FLUID_3D *fluid);
+float *smoke_get_react(struct FLUID_3D *fluid);
+float *smoke_get_color_r(struct FLUID_3D *fluid);
+float *smoke_get_color_g(struct FLUID_3D *fluid);
+float *smoke_get_color_b(struct FLUID_3D *fluid);
+void smoke_get_rgba(struct FLUID_3D *fluid, float *data, int sequential);
+void smoke_get_rgba_from_density(struct FLUID_3D *fluid, float color[3], float *data, int sequential);
 float *smoke_get_heat(struct FLUID_3D *fluid);
 float *smoke_get_velocity_x(struct FLUID_3D *fluid);
 float *smoke_get_velocity_y(struct FLUID_3D *fluid);
@@ -68,19 +74,44 @@ size_t smoke_get_index2d(int x, int max_x, int y);
 void smoke_dissolve(struct FLUID_3D *fluid, int speed, int log);
 
 // wavelet turbulence functions
-struct WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype);
+struct WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype, int use_fire, int use_colors);
 void smoke_turbulence_free(struct WTURBULENCE *wt);
 void smoke_turbulence_step(struct WTURBULENCE *wt, struct FLUID_3D *fluid);
 
 float *smoke_turbulence_get_density(struct WTURBULENCE *wt);
+float *smoke_turbulence_get_color_r(struct WTURBULENCE *wt);
+float *smoke_turbulence_get_color_g(struct WTURBULENCE *wt);
+float *smoke_turbulence_get_color_b(struct WTURBULENCE *wt);
+void smoke_turbulence_get_rgba(struct WTURBULENCE *wt, float *data, int sequential);
+void smoke_turbulence_get_rgba_from_density(struct WTURBULENCE *wt, float color[3], float *data, int sequential);
+float *smoke_turbulence_get_flame(struct WTURBULENCE *wt);
+float *smoke_turbulence_get_fuel(struct WTURBULENCE *wt);
+float *smoke_turbulence_get_react(struct WTURBULENCE *wt);
 void smoke_turbulence_get_res(struct WTURBULENCE *wt, int *res);
+int smoke_turbulence_get_cells(struct WTURBULENCE *wt);
 void smoke_turbulence_set_noise(struct WTURBULENCE *wt, int type);
 void smoke_initWaveletBlenderRNA(struct WTURBULENCE *wt, float *strength);
-
 void smoke_dissolve_wavelet(struct WTURBULENCE *wt, int speed, int log);
 
-// export
-void smoke_turbulence_export(struct WTURBULENCE *wt, float **dens, float **densold, float **tcu, float **tcv, float **tcw);
+/* export */
+void smoke_export(struct FLUID_3D *fluid, float *dt, float *dx, float **dens, float **react, float **flame, float **fuel, float **heat, float **heatold,
+				  float **vx, float **vy, float **vz, float **r, float **g, float **b, unsigned char **obstacles);
+void smoke_turbulence_export(struct WTURBULENCE *wt, float **dens, float **react, float **flame, float **fuel,
+							 float **r, float **g, float **b, float **tcu, float **tcv, float **tcw);
+
+/* flame spectrum */
+void flame_get_spectrum(unsigned char *spec, int width, float t1, float t2);
+
+/* data fields */
+int smoke_has_heat(struct FLUID_3D *fluid);
+int smoke_has_fuel(struct FLUID_3D *fluid);
+int smoke_has_colors(struct FLUID_3D *fluid);
+int smoke_turbulence_has_fuel(struct WTURBULENCE *wt);
+int smoke_turbulence_has_colors(struct WTURBULENCE *wt);
+
+void smoke_ensure_heat(struct FLUID_3D *fluid);
+void smoke_ensure_fire(struct FLUID_3D *fluid, struct WTURBULENCE *wt);
+void smoke_ensure_colors(struct FLUID_3D *fluid, struct WTURBULENCE *wt, float init_r, float init_g, float init_b);
 
 #ifdef __cplusplus
 }
diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp
index e006132ea8f..b7b353352e2 100644
--- a/intern/smoke/intern/FLUID_3D.cpp
+++ b/intern/smoke/intern/FLUID_3D.cpp
@@ -44,16 +44,11 @@
 // Construction/Destruction
 //////////////////////////////////////////////////////////////////////
 
-FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
+FLUID_3D::FLUID_3D(int *res, float dx, float dtdef, int init_heat, int init_fire, int init_colors) :
 	_xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f)
 {
 	// set simulation consts
 	_dt = dtdef;	// just in case. set in step from a RNA factor
-	
-	// start point of array
-	_p0[0] = p0[0];
-	_p0[1] = p0[1];
-	_p0[2] = p0[2];
 
 	_iterations = 100;
 	_tempAmb = 0; 
@@ -72,7 +67,10 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 	*/
 	
 	// scale the constants according to the refinement of the grid
-	_dx = 1.0f / (float)_maxRes;
+	if (!dx)
+		_dx = 1.0f / (float)_maxRes;
+	else
+		_dx = dx;
 	_constantScaling = 64.0f / _maxRes;
 	_constantScaling = (_constantScaling < 1.0f) ? 1.0f : _constantScaling;
 	_vorticityEps = 2.0f / _constantScaling; // Just in case set a default value
@@ -94,8 +92,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 	_zForce       = new float[_totalCells];
 	_density      = new float[_totalCells];
 	_densityOld   = new float[_totalCells];
-	_heat         = new float[_totalCells];
-	_heatOld      = new float[_totalCells];
 	_obstacles    = new unsigned char[_totalCells]; // set 0 at end of step
 
 	// For threaded version:
@@ -103,7 +99,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 	_yVelocityTemp = new float[_totalCells];
 	_zVelocityTemp = new float[_totalCells];
 	_densityTemp   = new float[_totalCells];
-	_heatTemp      = new float[_totalCells];
 
 	// DG TODO: check if alloc went fine
 
@@ -111,8 +106,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 	{
 		_density[x]      = 0.0f;
 		_densityOld[x]   = 0.0f;
-		_heat[x]         = 0.0f;
-		_heatOld[x]      = 0.0f;
 		_xVelocity[x]    = 0.0f;
 		_yVelocity[x]    = 0.0f;
 		_zVelocity[x]    = 0.0f;
@@ -128,6 +121,25 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 		_obstacles[x]    = false;
 	}
 
+	/* heat */
+	_heat = _heatOld = _heatTemp = NULL;
+	if (init_heat) {
+		initHeat();
+	}
+	// Fire simulation
+	_flame = _fuel = _fuelTemp = _fuelOld = NULL;
+	_react = _reactTemp = _reactOld = NULL;
+	if (init_fire) {
+		initFire();
+	}
+	// Smoke color
+	_color_r = _color_rOld = _color_rTemp = NULL;
+	_color_g = _color_gOld = _color_gTemp = NULL;
+	_color_b = _color_bOld = _color_bTemp = NULL;
+	if (init_colors) {
+		initColors(0.0f, 0.0f, 0.0f);
+	}
+
 	// boundary conditions of the fluid domain
 	// set default values -> vertically non-colliding
 	_domainBcFront = true;
@@ -138,9 +150,70 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 	_domainBcRight	= _domainBcLeft;
 
 	_colloPrev = 1;	// default value
+}
 
-	setBorderObstacles(); // walls
+void FLUID_3D::initHeat()
+{
+	if (!_heat) {
+		_heat         = new float[_totalCells];
+		_heatOld      = new float[_totalCells];
+		_heatTemp      = new float[_totalCells];
 
+		for (int x = 0; x < _totalCells; x++)
+		{
+			_heat[x]         = 0.0f;
+			_heatOld[x]      = 0.0f;
+		}
+	}
+}
+
+void FLUID_3D::initFire()
+{
+	if (!_flame) {
+		_flame		= new float[_totalCells];
+		_fuel		= new float[_totalCells];
+		_fuelTemp	= new float[_totalCells];
+		_fuelOld	= new float[_totalCells];
+		_react		= new float[_totalCells];
+		_reactTemp	= new float[_totalCells];
+		_reactOld	= new float[_totalCells];
+
+		for (int x = 0; x < _totalCells; x++)
+		{
+			_flame[x]		= 0.0f;
+			_fuel[x]		= 0.0f;
+			_fuelTemp[x]	= 0.0f;
+			_fuelOld[x]		= 0.0f;
+			_react[x]		= 0.0f;
+			_reactTemp[x]	= 0.0f;
+			_reactOld[x]	= 0.0f;
+		}
+	}
+}
+
+void FLUID_3D::initColors(float init_r, float init_g, float init_b)
+{
+	if (!_color_r) {
+		_color_r		= new float[_totalCells];
+		_color_rOld		= new float[_totalCells];
+		_color_rTemp	= new float[_totalCells];
+		_color_g		= new float[_totalCells];
+		_color_gOld		= new float[_totalCells];
+		_color_gTemp	= new float[_totalCells];
+		_color_b		= new float[_totalCells];
+		_color_bOld		= new float[_totalCells];
+		_color_bTemp	= new float[_totalCells];
+
+		for (int x = 0; x < _totalCells; x++)
+		{
+			_color_r[x]		= _density[x] * init_r;
+			_color_rOld[x]	= 0.0f;
+			_color_g[x]		= _density[x] * init_g;
+			_color_gOld[x]	= 0.0f;
+			_color_b[x]		= _density[x] * init_b;
+			_color_bOld[x]	= 0.0f;
+		}
+	}
 }
 
 void FLUID_3D::setBorderObstacles()
@@ -204,7 +277,6 @@ FLUID_3D::~FLUID_3D()
 	if (_heat) delete[] _heat;
 	if (_heatOld) delete[] _heatOld;
 	if (_obstacles) delete[] _obstacles;
-    // if (_wTurbulence) delete _wTurbulence;
 
 	if (_xVelocityTemp) delete[] _xVelocityTemp;
 	if (_yVelocityTemp) delete[] _yVelocityTemp;
@@ -212,23 +284,48 @@ FLUID_3D::~FLUID_3D()
 	if (_densityTemp) delete[] _densityTemp;
 	if (_heatTemp) delete[] _heatTemp;
 
+	if (_flame) delete[] _flame;
+	if (_fuel) delete[] _fuel;
+	if (_fuelTemp) delete[] _fuelTemp;
+	if (_fuelOld) delete[] _fuelOld;
+	if (_react) delete[] _react;
+	if (_reactTemp) delete[] _reactTemp;
+	if (_reactOld) delete[] _reactOld;
+
+	if (_color_r) delete[] _color_r;
+	if (_color_rOld) delete[] _color_rOld;
+	if (_color_rTemp) delete[] _color_rTemp;
+	if (_color_g) delete[] _color_g;
+	if (_color_gOld) delete[] _color_gOld;
+	if (_color_gTemp) delete[] _color_gTemp;
+	if (_color_b) delete[] _color_b;
+	if (_color_bOld) delete[] _color_bOld;
+	if (_color_bTemp) delete[] _color_bTemp;
+
     // printf("deleted fluid\n");
 }
 
 // init direct access functions from blender
-void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision)
+void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision, float *burning_rate,
+							  float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp)
 {
 	_alpha = alpha;
 	_beta = beta;
 	_dtFactor = dt_factor;
 	_vorticityRNA = vorticity;
 	_borderColli = borderCollision;
+	_burning_rate = burning_rate;
+	_flame_smoke = flame_smoke;
+	_flame_smoke_color = flame_smoke_color;
+	_flame_vorticity = flame_vorticity;
+	_ignition_temp = flame_ignition_temp;
+	_max_temp = flame_max_temp;
 }
 
 //////////////////////////////////////////////////////////////////////
 // step simulation once
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::step(float dt)
+void FLUID_3D::step(float dt, float gravity[3])
 {
 #if 0
 	// If border rules have been changed
@@ -281,7 +378,7 @@ void FLUID_3D::step(float dt)
 
 		wipeBoundariesSL(zBegin, zEnd);
 		addVorticity(zBegin, zEnd);
-		addBuoyancy(_heat, _density, zBegin, zEnd);
+		addBuoyancy(_heat, _density, gravity, zBegin, zEnd);
 		addForce(zBegin, zEnd);
 
 #if PARALLEL==1
@@ -312,13 +409,15 @@ void FLUID_3D::step(float dt)
 		if (i==0)
 		{
 #endif
-		project();
+			project();
 #if PARALLEL==1
 		}
-		else
+		else if (i==1)
 		{
 #endif
-		diffuseHeat();
+			if (_heat) {
+				diffuseHeat();
+			}
 #if PARALLEL==1
 		}
 	}
@@ -338,6 +437,13 @@ void FLUID_3D::step(float dt)
 	SWAP_POINTERS(_density, _densityOld);
 	SWAP_POINTERS(_heat, _heatOld);
 
+	SWAP_POINTERS(_fuel, _fuelOld);
+	SWAP_POINTERS(_react, _reactOld);
+
+	SWAP_POINTERS(_color_r, _color_rOld);
+	SWAP_POINTERS(_color_g, _color_gOld);
+	SWAP_POINTERS(_color_b, _color_bOld);
+
 	advectMacCormackBegin(0, _zRes);
 
 #if PARALLEL==1
@@ -398,11 +504,8 @@ void FLUID_3D::step(float dt)
 	SWAP_POINTERS(_yVelocity, _yForce);
 	SWAP_POINTERS(_zVelocity, _zForce);
 
-
-
-
 	_totalTime += _dt;
-	_totalSteps++;		
+	_totalSteps++;
 
 	for (int i = 0; i < _totalCells; i++)
 	{
@@ -643,6 +746,15 @@ void FLUID_3D::wipeBoundaries(int zBegin, int zEnd)
 	setZeroBorder(_yVelocity, _res, zBegin, zEnd);
 	setZeroBorder(_zVelocity, _res, zBegin, zEnd);
 	setZeroBorder(_density, _res, zBegin, zEnd);
+	if (_fuel) {
+		setZeroBorder(_fuel, _res, zBegin, zEnd);
+		setZeroBorder(_react, _res, zBegin, zEnd);
+	}
+	if (_color_r) {
+		setZeroBorder(_color_r, _res, zBegin, zEnd);
+		setZeroBorder(_color_g, _res, zBegin, zEnd);
+		setZeroBorder(_color_b, _res, zBegin, zEnd);
+	}
 }
 
 void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
@@ -668,6 +780,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 			_yVelocity[index] = 0.0f;
 			_zVelocity[index] = 0.0f;
 			_density[index] = 0.0f;
+			if (_fuel) {
+				_fuel[index] = 0.0f;
+				_react[index] = 0.0f;
+			}
+			if (_color_r) {
+				_color_r[index] = 0.0f;
+				_color_g[index] = 0.0f;
+				_color_b[index] = 0.0f;
+			}
 
 			// right slab
 			index += _xRes - 1;
@@ -675,6 +796,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 			_yVelocity[index] = 0.0f;
 			_zVelocity[index] = 0.0f;
 			_density[index] = 0.0f;
+			if (_fuel) {
+				_fuel[index] = 0.0f;
+				_react[index] = 0.0f;
+			}
+			if (_color_r) {
+				_color_r[index] = 0.0f;
+				_color_g[index] = 0.0f;
+				_color_b[index] = 0.0f;
+			}
 		}
 
 	/////////////////////////////////////
@@ -690,6 +820,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 			_yVelocity[index] = 0.0f;
 			_zVelocity[index] = 0.0f;
 			_density[index] = 0.0f;
+			if (_fuel) {
+				_fuel[index] = 0.0f;
+				_react[index] = 0.0f;
+			}
+			if (_color_r) {
+				_color_r[index] = 0.0f;
+				_color_g[index] = 0.0f;
+				_color_b[index] = 0.0f;
+			}
 
 			// top slab
 			index += slabSize - _xRes;
@@ -697,6 +836,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 			_yVelocity[index] = 0.0f;
 			_zVelocity[index] = 0.0f;
 			_density[index] = 0.0f;
+			if (_fuel) {
+				_fuel[index] = 0.0f;
+				_react[index] = 0.0f;
+			}
+			if (_color_r) {
+				_color_r[index] = 0.0f;
+				_color_g[index] = 0.0f;
+				_color_b[index] = 0.0f;
+			}
 
 		}
 
@@ -717,6 +865,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 			_yVelocity[index] = 0.0f;
 			_zVelocity[index] = 0.0f;
 			_density[index] = 0.0f;
+			if (_fuel) {
+				_fuel[index] = 0.0f;
+				_react[index] = 0.0f;
+			}
+			if (_color_r) {
+				_color_r[index] = 0.0f;
+				_color_g[index] = 0.0f;
+				_color_b[index] = 0.0f;
+			}
     }
 
 	if (zEnd == _zRes)
@@ -735,6 +892,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 				_yVelocity[indexx] = 0.0f;
 				_zVelocity[indexx] = 0.0f;
 				_density[indexx] = 0.0f;
+				if (_fuel) {
+					_fuel[index] = 0.0f;
+					_react[index] = 0.0f;
+				}
+				if (_color_r) {
+					_color_r[index] = 0.0f;
+					_color_g[index] = 0.0f;
+					_color_b[index] = 0.0f;
+				}
 			}
 	}
 
@@ -781,35 +947,6 @@ void FLUID_3D::project()
 	if(!_domainBcTop) setNeumannZ(_zVelocity, _res, 0, _zRes);
 	else setZeroZ(_zVelocity, _res, 0, _zRes);
 
-	/*
-	{
-		float maxx = 0, maxy = 0, maxz = 0;
-		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
-		{
-			if(_xVelocity[i] > maxx)
-				maxx = _xVelocity[i];
-			if(_yVelocity[i] > maxy)
-				maxy = _yVelocity[i];
-			if(_zVelocity[i] > maxz)
-				maxz = _zVelocity[i];
-		}
-		printf("Max velx: %f, vely: %f, velz: %f\n", maxx, maxy, maxz);
-	}
-	*/
-
-	/*
-	{
-		float maxvalue = 0;
-		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
-		{
-			if(_heat[i] > maxvalue)
-				maxvalue = _heat[i];
-
-		}
-		printf("Max heat: %f\n", maxvalue);
-	}
-	*/
-
 	// calculate divergence
 	index = _slabSize + _xRes + 1;
 	for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
@@ -1007,7 +1144,6 @@ void FLUID_3D::diffuseHeat()
 	for (int x = 0; x < _totalCells; x++)
 		if (_obstacles[x])
 			_heat[x] = 0.0f;
-
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -1175,7 +1311,7 @@ void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
 //////////////////////////////////////////////////////////////////////
 // add buoyancy forces
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
+void FLUID_3D::addBuoyancy(float *heat, float *density, float gravity[3], int zBegin, int zEnd)
 {
 	int index = zBegin*_slabSize;
 
@@ -1183,7 +1319,10 @@ void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
 		for (int y = 0; y < _yRes; y++)
 			for (int x = 0; x < _xRes; x++, index++)
 			{
-				_zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender
+				float buoyancy = *_alpha * density[index] + (*_beta * (((heat) ? heat[index] : 0.0f) - _tempAmb));
+				_xForce[index] -= gravity[0] * buoyancy;
+				_yForce[index] -= gravity[1] * buoyancy;
+				_zForce[index] -= gravity[2] * buoyancy;
 			}
 }
 
@@ -1192,8 +1331,10 @@ void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::addVorticity(int zBegin, int zEnd)
 {
+	// set flame vorticity from RNA value
+	float flame_vorticity = (*_flame_vorticity)/_constantScaling;
 	//int x,y,z,index;
-	if(_vorticityEps<=0.) return;
+	if(_vorticityEps+flame_vorticity<=0.) return;
 
 	int _blockSize=zEnd-zBegin;
 	int _blockTotalCells = _slabSize * (_blockSize+2);
@@ -1300,14 +1441,15 @@ void FLUID_3D::addVorticity(int zBegin, int zEnd)
 					float magnitude = sqrtf(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
 					if (magnitude > FLT_EPSILON)
 					{
+						float flame_vort = (_fuel) ? _fuel[index]*flame_vorticity : 0.0f;
 						magnitude = 1.0f / magnitude;
 						N[0] *= magnitude;
 						N[1] *= magnitude;
 						N[2] *= magnitude;
 
-						_xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * eps;
-						_yForce[index] += (N[2] * _xVorticity[vIndex] - N[0] * _zVorticity[vIndex]) * _dx * eps;
-						_zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
+						_xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * (eps + flame_vort);
+						_yForce[index] += (N[2] * _xVorticity[vIndex] - N[0] * _zVorticity[vIndex]) * _dx * (eps + flame_vort);
+						_zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * (eps + flame_vort);
 					}
 					}	// if
 					vIndex++;
@@ -1328,14 +1470,9 @@ void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
 {
 	Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
 
-	if(!_domainBcLeft) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
-	else setZeroX(_xVelocityOld, res, zBegin, zEnd);
-
-	if(!_domainBcFront) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
-	else setZeroY(_yVelocityOld, res, zBegin, zEnd); 
-
-	if(!_domainBcTop) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
-	else setZeroZ(_zVelocityOld, res, zBegin, zEnd);
+	setZeroX(_xVelocityOld, res, zBegin, zEnd);
+	setZeroY(_yVelocityOld, res, zBegin, zEnd);
+	setZeroZ(_zVelocityOld, res, zBegin, zEnd);
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -1355,7 +1492,18 @@ void FLUID_3D::advectMacCormackEnd1(int zBegin, int zEnd)
 	// advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res)
 
 	advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, res, zBegin, zEnd);
-	advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
+	if (_heat) {
+		advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
+	}
+	if (_fuel) {
+		advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _fuelOld, _fuelTemp, res, zBegin, zEnd);
+		advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _reactOld, _reactTemp, res, zBegin, zEnd);
+	}
+	if (_color_r) {
+		advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_rOld, _color_rTemp, res, zBegin, zEnd);
+		advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_gOld, _color_gTemp, res, zBegin, zEnd);
+		advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_bOld, _color_bTemp, res, zBegin, zEnd);
+	}
 	advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, res, zBegin, zEnd);
 	advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, res, zBegin, zEnd);
 	advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, res, zBegin, zEnd);
@@ -1371,17 +1519,30 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
 	const float dt0 = _dt / _dx;
 	Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
 
-	// use force array as temp arrays
+	// use force array as temp array
 	float* t1 = _xForce;
 
 	// advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles)
 
+	/* finish advection */
 	advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd);
-	advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
+	if (_heat) {
+		advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
+	}
+	if (_fuel) {
+		advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _fuelOld, _fuel, _fuelTemp, t1, res, _obstacles, zBegin, zEnd);
+		advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _reactOld, _react, _reactTemp, t1, res, _obstacles, zBegin, zEnd);
+	}
+	if (_color_r) {
+		advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_rOld, _color_r, _color_rTemp, t1, res, _obstacles, zBegin, zEnd);
+		advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_gOld, _color_g, _color_gTemp, t1, res, _obstacles, zBegin, zEnd);
+		advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_bOld, _color_b, _color_bTemp, t1, res, _obstacles, zBegin, zEnd);
+	}
 	advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocityTemp, _xVelocity, t1, res, _obstacles, zBegin, zEnd);
 	advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd);
 	advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd);
 
+	/* set boundary conditions for velocity */
 	if(!_domainBcLeft) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
 	else setZeroX(_xVelocityTemp, res, zBegin, zEnd);				
 
@@ -1391,40 +1552,71 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
 	if(!_domainBcTop) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
 	else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
 
+	/* clear data boundaries */
 	setZeroBorder(_density, res, zBegin, zEnd);
-	setZeroBorder(_heat, res, zBegin, zEnd);
-#if 0
-	{
-	const size_t index_ = _slabSize + _xRes + 1;
-	int bb=0;
-	int bt=0;
+	if (_fuel) {
+		setZeroBorder(_fuel, res, zBegin, zEnd);
+		setZeroBorder(_react, res, zBegin, zEnd);
+	}
+	if (_color_r) {
+		setZeroBorder(_color_r, res, zBegin, zEnd);
+		setZeroBorder(_color_g, res, zBegin, zEnd);
+		setZeroBorder(_color_b, res, zBegin, zEnd);
+	}
+}
 
-	if (zBegin == 0) {bb = 1;}
-	if (zEnd == _zRes) {bt = 1;}
-	
-	for (int z = zBegin + bb; z < zEnd - bt; z++)
+
+void FLUID_3D::processBurn(float *fuel, float *smoke, float *react, float *flame, float *heat,
+						   float *r, float *g, float *b, int total_cells, float dt)
+{
+	float burning_rate = *_burning_rate;
+	float flame_smoke = *_flame_smoke;
+	float ignition_point = *_ignition_temp;
+	float temp_max = *_max_temp;
+
+	for (int index = 0; index < total_cells; index++)
 	{
-		size_t index = index_ +(z-1)*_slabSize;
+		float orig_fuel = fuel[index];
+		float orig_smoke = smoke[index];
+		float smoke_emit = 0.0f;
+		float react_coord = 0.0f;
+
+		/* process fuel */
+		fuel[index] -= burning_rate * dt;
+		if (fuel[index] < 0.0f) fuel[index] = 0.0f;
+		/* process reaction coordinate */
+		if (orig_fuel) {
+			react[index] *= fuel[index]/orig_fuel;
+			react_coord = react[index];
+		}
+		else {
+			react[index] = 0.0f;
+		}
 
-		for (int y = 1; y < _yRes - 1; y++, index += 2)
-		{
-			for (int x = 1; x < _xRes - 1; x++, index++)
-			{
-				// clean custom velocities from moving obstacles again
-				if (_obstacles[index])
-				{
-					_xVelocity[index] =
-					_yVelocity[index] =
-					_zVelocity[index] = 0.0f;
-				}
-			}
+		/* emit smoke based on fuel burn rate and "flame_smoke" factor */
+		smoke_emit = (orig_fuel < 1.0f) ? (1.0f - orig_fuel)*0.5f : 0.0f;
+		smoke_emit = (smoke_emit + 0.5f) * (orig_fuel-fuel[index]) * 0.1f * flame_smoke;
+		smoke[index] += smoke_emit;
+		CLAMP(smoke[index], 0.0f, 1.0f);
+
+		/* model flame temperature curve from the reaction coordinate (fuel) */
+		if (react_coord>0.0f) {
+			/* do a smooth falloff for rest of the values */
+			flame[index] = pow(react_coord, 0.5f);
+		}
+		else
+			flame[index] = 0.0f;
+
+		/* set fluid temperature from the flame temperature profile */
+		if (heat && flame[index])
+			heat[index] = (1.0f-flame[index])*ignition_point + flame[index]*temp_max;
+
+		/* mix new color */
+		if (r && smoke_emit) {
+			float smoke_factor = smoke[index]/(orig_smoke+smoke_emit);
+			r[index] = (r[index] + _flame_smoke_color[0] * smoke_emit) * smoke_factor;
+			g[index] = (g[index] + _flame_smoke_color[1] * smoke_emit) * smoke_factor;
+			b[index] = (b[index] + _flame_smoke_color[2] * smoke_emit) * smoke_factor;
 		}
 	}
-	}
-#endif
-
-	/*int begin=zBegin * _slabSize;
-	int end=begin + (zEnd - zBegin) * _slabSize;
-  for (int x = begin; x < end; x++)
-    _xForce[x] = _yForce[x] = 0.0f;*/
 }
diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h
index ac77148ce84..8cadf3bc989 100644
--- a/intern/smoke/intern/FLUID_3D.h
+++ b/intern/smoke/intern/FLUID_3D.h
@@ -46,11 +46,16 @@ class WTURBULENCE;
 class FLUID_3D  
 {
 	public:
-		FLUID_3D(int *res, /* int amplify, */ float *p0, float dtdef);
+		FLUID_3D(int *res, float dx, float dtdef, int init_heat, int init_fire, int init_colors);
 		FLUID_3D() {};
 		virtual ~FLUID_3D();
 
-		void initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli);
+		void initHeat();
+		void initFire();
+		void initColors(float init_r, float init_g, float init_b);
+
+		void initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli, float *burning_rate,
+							float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *ignition_temp, float *max_temp);
 		
 		// create & allocate vector noise advection 
 		void initVectorNoise(int amplify);
@@ -58,7 +63,7 @@ class FLUID_3D
 		void addSmokeColumn();
 		static void addSmokeTestCase(float* field, Vec3Int res);
 
-		void step(float dt);
+		void step(float dt, float gravity[3]);
 		void addObstacle(OBSTACLE* obstacle);
 
 		const float* xVelocity() { return _xVelocity; }; 
@@ -115,6 +120,27 @@ class FLUID_3D
 		float* _heatTemp;
 		float* _densityTemp;
 
+		// fire simulation
+		float *_flame;
+		float *_fuel;
+		float *_fuelTemp;
+		float *_fuelOld;
+		float *_react;
+		float *_reactTemp;
+		float *_reactOld;
+
+		// smoke color
+		float *_color_r;
+		float *_color_rOld;
+		float *_color_rTemp;
+		float *_color_g;
+		float *_color_gOld;
+		float *_color_gTemp;
+		float *_color_b;
+		float *_color_bOld;
+		float *_color_bTemp;
+
+
 		// CG fields
 		int _iterations;
 
@@ -153,14 +179,16 @@ class FLUID_3D
 		void wipeBoundariesSL(int zBegin, int zEnd);
 		void addForce(int zBegin, int zEnd);
 		void addVorticity(int zBegin, int zEnd);
-		void addBuoyancy(float *heat, float *density, int zBegin, int zEnd);
+		void addBuoyancy(float *heat, float *density, float gravity[3], int zBegin, int zEnd);
 
 		// solver stuff
 		void project();
 		void diffuseHeat();
+		void diffuseColor();
 		void solvePressure(float* field, float* b, unsigned char* skip);
 		void solvePressurePre(float* field, float* b, unsigned char* skip);
 		void solveHeat(float* field, float* b, unsigned char* skip);
+		void solveDiffusion(float* field, float* b, float* factor);
 
 
 		// handle obstacle boundaries
@@ -174,6 +202,16 @@ class FLUID_3D
 		void advectMacCormackEnd1(int zBegin, int zEnd);
 		void advectMacCormackEnd2(int zBegin, int zEnd);
 
+		/* burning */
+		float *_burning_rate; // RNA pointer
+		float *_flame_smoke; // RNA pointer
+		float *_flame_smoke_color; // RNA pointer
+		float *_flame_vorticity; // RNA pointer
+		float *_ignition_temp; // RNA pointer
+		float *_max_temp; // RNA pointer
+		void processBurn(float *fuel, float *smoke, float *react, float *flame, float *heat,
+						 float *r, float *g, float *b, int total_cells, float dt);
+
 		// boundary setting functions
 		static void copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd);
 		static void copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd);
diff --git a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
index 3cf94eb00a9..42a8b2d6923 100644
--- a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
+++ b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
@@ -165,7 +165,6 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
 	if (_Acenter)  delete[] _Acenter;
 }
 
-
 void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
 {
 	int x, y, z;
diff --git a/intern/smoke/intern/FLUID_3D_STATIC.cpp b/intern/smoke/intern/FLUID_3D_STATIC.cpp
index c7a0ddcd04c..ac485ad983a 100644
--- a/intern/smoke/intern/FLUID_3D_STATIC.cpp
+++ b/intern/smoke/intern/FLUID_3D_STATIC.cpp
@@ -92,18 +92,10 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd)
 			// left slab
 			index = y * res[0] + z * slabSize;
 			field[index] = field[index + 2];
-			/* only allow outwards flux */
-			if(field[index]>0.) field[index] = 0.;
-			index += 1;
-			if(field[index]>0.) field[index] = 0.;
 
 			// right slab
 			index = y * res[0] + z * slabSize + res[0] - 1;
 			field[index] = field[index - 2];
-			/* only allow outwards flux */
-			if(field[index]<0.) field[index] = 0.;
-			index -= 1;
-			if(field[index]<0.) field[index] = 0.;
 		}
  }
 
@@ -120,18 +112,10 @@ void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd)
 			// front slab
 			index = x + z * slabSize;
 			field[index] = field[index + 2 * res[0]];
-			/* only allow outwards flux */
-			if(field[index]>0.) field[index] = 0.;
-			index += res[0];
-			if(field[index]>0.) field[index] = 0.;
 
 			// back slab
 			index = x + z * slabSize + slabSize - res[0];
 			field[index] = field[index - 2 * res[0]];
-			/* only allow outwards flux */
-			if(field[index]<0.) field[index] = 0.;
-			index -= res[0];
-			if(field[index]<0.) field[index] = 0.;
 		}
 }
 
@@ -152,14 +136,6 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd)
 				// front slab
 				index = x + y * res[0];
 				field[index] = field[index + 2 * slabSize];
-				/* only allow outwards flux */
-
-				// DG: Disable this for z-axis.
-				// The problem is that smoke somehow gets sucked in again
-				// from the TOP slab when this is enabled
-				// if(field[index]>0.) field[index] = 0.;
-				// index += slabSize;
-				// if(field[index]>0.) field[index] = 0.;
 			}
 	}
 
@@ -170,10 +146,6 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd)
 				// back slab
 				index = x + y * res[0] + cellsslab;
 				field[index] = field[index - 2 * slabSize];
-				/* only allow outwards flux */
-				if(field[index]<0.) field[index] = 0.;
-				index -= slabSize;
-				if(field[index]<0.) field[index] = 0.;
 			}
 	}
 		
diff --git a/intern/smoke/intern/WTURBULENCE.cpp b/intern/smoke/intern/WTURBULENCE.cpp
index 671198065e8..1c89f5d681b 100644
--- a/intern/smoke/intern/WTURBULENCE.cpp
+++ b/intern/smoke/intern/WTURBULENCE.cpp
@@ -51,7 +51,7 @@ static const float persistence = 0.56123f;
 //////////////////////////////////////////////////////////////////////
 // constructor
 //////////////////////////////////////////////////////////////////////
-WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype)
+WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype, int init_fire, int init_colors)
 {
 	// if noise magnitude is below this threshold, its contribution
 	// is negilgible, so stop evaluating new octaves
@@ -87,12 +87,26 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no
 	// allocate high resolution density field
 	_totalStepsBig = 0;
 	_densityBig = new float[_totalCellsBig];
-	_densityBigOld = new float[_totalCellsBig]; 
+	_densityBigOld = new float[_totalCellsBig];
 	
 	for(int i = 0; i < _totalCellsBig; i++) {
 		_densityBig[i] = 
 		_densityBigOld[i] = 0.;
 	}
+
+	/* fire */
+	_flameBig = _fuelBig = _fuelBigOld = NULL;
+	_reactBig = _reactBigOld = NULL;
+	if (init_fire) {
+		initFire();
+	}
+	/* colors */
+	_color_rBig = _color_rBigOld = NULL;
+	_color_gBig = _color_gBigOld = NULL;
+	_color_bBig = _color_bBigOld = NULL;
+	if (init_colors) {
+		initColors(0.0f, 0.0f, 0.0f);
+	}
 	
 	// allocate & init texture coordinates
 	_tcU = new float[_totalCellsSm];
@@ -128,12 +142,64 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no
 	*/
 }
 
+void WTURBULENCE::initFire()
+{
+	if (!_fuelBig) {
+		_flameBig = new float[_totalCellsBig];
+		_fuelBig = new float[_totalCellsBig];
+		_fuelBigOld = new float[_totalCellsBig];
+		_reactBig = new float[_totalCellsBig];
+		_reactBigOld = new float[_totalCellsBig];
+
+		for(int i = 0; i < _totalCellsBig; i++) {
+			_flameBig[i] = 
+			_fuelBig[i] = 
+			_fuelBigOld[i] = 0.;
+			_reactBig[i] = 
+			_reactBigOld[i] = 0.;
+		}
+	}
+}
+
+void WTURBULENCE::initColors(float init_r, float init_g, float init_b)
+{
+	if (!_color_rBig) {
+		_color_rBig = new float[_totalCellsBig];
+		_color_rBigOld = new float[_totalCellsBig];
+		_color_gBig = new float[_totalCellsBig];
+		_color_gBigOld = new float[_totalCellsBig];
+		_color_bBig = new float[_totalCellsBig];
+		_color_bBigOld = new float[_totalCellsBig];
+
+		for(int i = 0; i < _totalCellsBig; i++) {
+			_color_rBig[i] = _densityBig[i] * init_r;
+			_color_rBigOld[i] = 0.0f;
+			_color_gBig[i] = _densityBig[i] * init_g;
+			_color_gBigOld[i] = 0.0f;
+			_color_bBig[i] = _densityBig[i] * init_b;
+			_color_bBigOld[i] = 0.0f;
+		}
+	}
+}
+
 //////////////////////////////////////////////////////////////////////
 // destructor
 //////////////////////////////////////////////////////////////////////
 WTURBULENCE::~WTURBULENCE() {
   delete[] _densityBig;
   delete[] _densityBigOld;
+  if (_flameBig) delete[] _flameBig;
+  if (_fuelBig) delete[] _fuelBig;
+  if (_fuelBigOld) delete[] _fuelBigOld;
+  if (_reactBig) delete[] _reactBig;
+  if (_reactBigOld) delete[] _reactBigOld;
+
+  if (_color_rBig) delete[] _color_rBig;
+  if (_color_rBigOld) delete[] _color_rBigOld;
+  if (_color_gBig) delete[] _color_gBig;
+  if (_color_gBigOld) delete[] _color_gBigOld;
+  if (_color_bBig) delete[] _color_bBig;
+  if (_color_bBigOld) delete[] _color_bBigOld;
 
   delete[] _tcU;
   delete[] _tcV;
@@ -757,8 +823,10 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 	// enlarge timestep to match grid
 	const float dt = dtOrg * _amplify;
 	const float invAmp = 1.0f / _amplify;
-	float *tempBig1 = (float *)calloc(_totalCellsBig, sizeof(float));
-	float *tempBig2 = (float *)calloc(_totalCellsBig, sizeof(float));
+	float *tempFuelBig = NULL, *tempReactBig = NULL;
+	float *tempColor_rBig = NULL, *tempColor_gBig = NULL, *tempColor_bBig = NULL;
+	float *tempDensityBig = (float *)calloc(_totalCellsBig, sizeof(float));
+	float *tempBig = (float *)calloc(_totalCellsBig, sizeof(float));
 	float *bigUx = (float *)calloc(_totalCellsBig, sizeof(float));
 	float *bigUy = (float *)calloc(_totalCellsBig, sizeof(float));
 	float *bigUz = (float *)calloc(_totalCellsBig, sizeof(float)); 
@@ -767,11 +835,21 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 	float *eigMin  = (float *)calloc(_totalCellsSm, sizeof(float));
 	float *eigMax  = (float *)calloc(_totalCellsSm, sizeof(float));
 
+	if (_fuelBig) {
+		tempFuelBig = (float *)calloc(_totalCellsBig, sizeof(float));
+		tempReactBig = (float *)calloc(_totalCellsBig, sizeof(float));
+	}
+	if (_color_rBig) {
+		tempColor_rBig = (float *)calloc(_totalCellsBig, sizeof(float));
+		tempColor_gBig = (float *)calloc(_totalCellsBig, sizeof(float));
+		tempColor_bBig = (float *)calloc(_totalCellsBig, sizeof(float));
+	}
+
 	memset(_tcTemp, 0, sizeof(float)*_totalCellsSm);
 
 
 	// prepare textures
-	advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
+	advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempDensityBig, tempBig);
 
 	// do wavelet decomposition of energy
 	computeEnergy(_energy, xvel, yvel, zvel, obstacles);
@@ -972,6 +1050,11 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 
   // prepare density for an advection
   SWAP_POINTERS(_densityBig, _densityBigOld);
+  SWAP_POINTERS(_fuelBig, _fuelBigOld);
+  SWAP_POINTERS(_reactBig, _reactBigOld);
+  SWAP_POINTERS(_color_rBig, _color_rBigOld);
+  SWAP_POINTERS(_color_gBig, _color_gBigOld);
+  SWAP_POINTERS(_color_bBig, _color_bBigOld);
 
   // based on the maximum velocity present, see if we need to substep,
   // but cap the maximum number of substeps to 5
@@ -1017,7 +1100,21 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 		int zEnd = (int)((float)(i+1)*partSize + 0.5f);
 #endif
 		FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
-		    _densityBigOld, tempBig1, _resBig, zBegin, zEnd);
+		    _densityBigOld, tempDensityBig, _resBig, zBegin, zEnd);
+		if (_fuelBig) {
+			FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
+				_fuelBigOld, tempFuelBig, _resBig, zBegin, zEnd);
+			FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
+				_reactBigOld, tempReactBig, _resBig, zBegin, zEnd);
+		}
+		if (_color_rBig) {
+			FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
+				_color_rBigOld, tempColor_rBig, _resBig, zBegin, zEnd);
+			FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
+				_color_gBigOld, tempColor_gBig, _resBig, zBegin, zEnd);
+			FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
+				_color_bBigOld, tempColor_bBig, _resBig, zBegin, zEnd);
+		}
 #if PARALLEL==1
 	}
 
@@ -1030,18 +1127,43 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 		int zEnd = (int)((float)(i+1)*partSize + 0.5f);
 #endif
 		FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
-		    _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL, zBegin, zEnd);
+		    _densityBigOld, _densityBig, tempDensityBig, tempBig, _resBig, NULL, zBegin, zEnd);
+		if (_fuelBig) {
+			FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
+				_fuelBigOld, _fuelBig, tempFuelBig, tempBig, _resBig, NULL, zBegin, zEnd);
+			FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
+				_reactBigOld, _reactBig, tempReactBig, tempBig, _resBig, NULL, zBegin, zEnd);
+		}
+		if (_color_rBig) {
+			FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
+				_color_rBigOld, _color_rBig, tempColor_rBig, tempBig, _resBig, NULL, zBegin, zEnd);
+			FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
+				_color_gBigOld, _color_gBig, tempColor_gBig, tempBig, _resBig, NULL, zBegin, zEnd);
+			FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
+				_color_bBigOld, _color_bBig, tempColor_bBig, tempBig, _resBig, NULL, zBegin, zEnd);
+		}
 #if PARALLEL==1
 	}
 	}
 #endif
 
-    if (substep < totalSubsteps - 1) 
+	if (substep < totalSubsteps - 1) {
       SWAP_POINTERS(_densityBig, _densityBigOld);
+	  SWAP_POINTERS(_fuelBig, _fuelBigOld);
+	  SWAP_POINTERS(_reactBig, _reactBigOld);
+	  SWAP_POINTERS(_color_rBig, _color_rBigOld);
+	  SWAP_POINTERS(_color_gBig, _color_gBigOld);
+	  SWAP_POINTERS(_color_bBig, _color_bBigOld);
+	}
   } // substep
 
-  free(tempBig1);
-  free(tempBig2);
+  free(tempDensityBig);
+  if (tempFuelBig) free(tempFuelBig);
+  if (tempReactBig) free(tempReactBig);
+  if (tempColor_rBig) free(tempColor_rBig);
+  if (tempColor_gBig) free(tempColor_gBig);
+  if (tempColor_bBig) free(tempColor_bBig);
+  free(tempBig);
   free(bigUx);
   free(bigUy);
   free(bigUz);
@@ -1050,6 +1172,15 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
   
   // wipe the density borders
   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0 , _resBig[2]);
+  if (_fuelBig) {
+	FLUID_3D::setZeroBorder(_fuelBig, _resBig, 0 , _resBig[2]);
+	FLUID_3D::setZeroBorder(_reactBig, _resBig, 0 , _resBig[2]);
+  }
+  if (_color_rBig) {
+	  FLUID_3D::setZeroBorder(_color_rBig, _resBig, 0 , _resBig[2]);
+	  FLUID_3D::setZeroBorder(_color_gBig, _resBig, 0 , _resBig[2]);
+	  FLUID_3D::setZeroBorder(_color_bBig, _resBig, 0 , _resBig[2]);
+  }
     
   // reset texture coordinates now in preparation for next timestep
   // Shouldn't do this before generating the noise because then the 
diff --git a/intern/smoke/intern/WTURBULENCE.h b/intern/smoke/intern/WTURBULENCE.h
index f31ca100fdf..1655bd95d32 100644
--- a/intern/smoke/intern/WTURBULENCE.h
+++ b/intern/smoke/intern/WTURBULENCE.h
@@ -36,10 +36,13 @@ class WTURBULENCE
 {
 	public:
 		// both config files can be NULL, altCfg might override values from noiseCfg
-		WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype);
+		WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype, int init_fire, int init_colors);
 
 		/// destructor
 		virtual ~WTURBULENCE();
+
+		void initFire();
+		void initColors(float init_r, float init_g, float init_b);
 		
 		void setNoise(int type);
 		void initBlenderRNA(float *strength);
@@ -63,6 +66,8 @@ class WTURBULENCE
 
 		// access functions
 		inline float* getDensityBig() { return _densityBig; }
+		inline float* getFlameBig() { return _flameBig; }
+		inline float* getFuelBig() { return _fuelBig; }
 		inline float* getArrayTcU() { return _tcU; }
 		inline float* getArrayTcV() { return _tcV; }
 		inline float* getArrayTcW() { return _tcW; }
@@ -111,6 +116,18 @@ class WTURBULENCE
 
 		float* _densityBig;
 		float* _densityBigOld;
+		float* _flameBig;
+		float* _fuelBig;
+		float* _fuelBigOld;
+		float* _reactBig;
+		float* _reactBigOld;
+
+		float* _color_rBig;
+		float* _color_rBigOld;
+		float* _color_gBig;
+		float* _color_gBigOld;
+		float* _color_bBig;
+		float* _color_bBigOld;
 
 		// texture coordinates for noise
 		float* _tcU;
diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp
index 4bbf8e0a82b..e51c3176699 100644
--- a/intern/smoke/intern/smoke_API.cpp
+++ b/intern/smoke/intern/smoke_API.cpp
@@ -30,6 +30,7 @@
 
 #include "FLUID_3D.h"
 #include "WTURBULENCE.h"
+#include "spectrum.h"
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -37,22 +38,16 @@
 
 #include "../extern/smoke_API.h"  /* to ensure valid prototypes */
 
-// y in smoke is z in blender
-extern "C" FLUID_3D *smoke_init(int *res, float *p0, float dtdef)
+extern "C" FLUID_3D *smoke_init(int *res, float dx, float dtdef, int use_heat, int use_fire, int use_colors)
 {
-	// smoke lib uses y as top-bottom/vertical axis where blender uses z
-	FLUID_3D *fluid = new FLUID_3D(res, p0, dtdef);
-
-	// printf("xres: %d, yres: %d, zres: %d\n", res[0], res[1], res[2]);
-
+	FLUID_3D *fluid = new FLUID_3D(res, dx, dtdef, use_heat, use_fire, use_colors);
 	return fluid;
 }
 
-extern "C" WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype)
+extern "C" WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype, int use_fire, int use_colors)
 {
-	// initialize wavelet turbulence
 	if(amplify)
-		return new WTURBULENCE(res[0],res[1],res[2], amplify, noisetype);
+		return new WTURBULENCE(res[0],res[1],res[2], amplify, noisetype, use_fire, use_colors);
 	else 
 		return NULL;
 }
@@ -71,7 +66,6 @@ extern "C" void smoke_turbulence_free(WTURBULENCE *wt)
 
 extern "C" size_t smoke_get_index(int x, int max_x, int y, int max_y, int z /*, int max_z */)
 {
-	// // const int index = x + y * smd->res[0] + z * smd->res[0]*smd->res[1];
 	return x + y * max_x + z * max_x*max_y;
 }
 
@@ -80,137 +74,134 @@ extern "C" size_t smoke_get_index2d(int x, int max_x, int y /*, int max_y, int z
 	return x + y * max_x;
 }
 
-extern "C" void smoke_step(FLUID_3D *fluid, float dtSubdiv)
+extern "C" void smoke_step(FLUID_3D *fluid, float gravity[3], float dtSubdiv)
 {
-	fluid->step(dtSubdiv);
+	if (fluid->_fuel) {
+		fluid->processBurn(fluid->_fuel, fluid->_density, fluid->_react, fluid->_flame, fluid->_heat,
+						   fluid->_color_r, fluid->_color_g, fluid->_color_b, fluid->_totalCells, (*fluid->_dtFactor)*dtSubdiv);
+	}
+	fluid->step(dtSubdiv, gravity);
 }
 
 extern "C" void smoke_turbulence_step(WTURBULENCE *wt, FLUID_3D *fluid)
 {
+	if (wt->_fuelBig) {
+		fluid->processBurn(wt->_fuelBig, wt->_densityBig, wt->_reactBig, wt->_flameBig, 0,
+						   wt->_color_rBig, wt->_color_gBig, wt->_color_bBig, wt->_totalCellsBig, fluid->_dt);
+	}
 	wt->stepTurbulenceFull(fluid->_dt/fluid->_dx, fluid->_xVelocity, fluid->_yVelocity, fluid->_zVelocity, fluid->_obstacles); 
 }
 
-extern "C" void smoke_initBlenderRNA(FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli)
+extern "C" void smoke_initBlenderRNA(FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli, float *burning_rate,
+									 float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp)
 {
-	fluid->initBlenderRNA(alpha, beta, dt_factor, vorticity, border_colli);
+	fluid->initBlenderRNA(alpha, beta, dt_factor, vorticity, border_colli, burning_rate, flame_smoke, flame_smoke_color, flame_vorticity, flame_ignition_temp, flame_max_temp);
 }
 
-extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log)
+extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength)
 {
-	float *density = fluid->_density;
-	//float *densityOld = fluid->_densityOld;
-	float *heat = fluid->_heat;
+	wt->initBlenderRNA(strength);
+}
 
+static void data_dissolve(float *density, float *heat, float *r, float *g, float *b, int total_cells, int speed, int log)
+{
 	if(log)
 	{
 		/* max density/speed = dydx */
-		float dydx = 1.0 / (float)speed;
-		size_t size= fluid->_xRes * fluid->_yRes * fluid->_zRes;
+		float fac = 1.0f - (1.0f / (float)speed);
 
-		for(size_t i = 0; i < size; i++)
+		for(size_t i = 0; i < total_cells; i++)
 		{
-			density[i] *= (1.0 - dydx);
-
-			if(density[i] < 0.0f)
-				density[i] = 0.0f;
-
-			heat[i] *= (1.0 - dydx);
-
-			/*if(heat[i] < 0.0f)
-				heat[i] = 0.0f;*/
+			/* density */
+			density[i] *= fac;
+
+			/* heat */
+			if (heat) {
+				heat[i] *= fac;
+			}
+
+			/* color */
+			if (r) {
+				r[i] *= fac;
+				g[i] *= fac;
+				b[i] *= fac;
+			}
 		}
 	}
 	else // linear falloff
 	{
 		/* max density/speed = dydx */
-		float dydx = 1.0 / (float)speed;
-		size_t size= fluid->_xRes * fluid->_yRes * fluid->_zRes;
+		float dydx = 1.0f / (float)speed;
 
-		for(size_t i = 0; i < size; i++)
+		for(size_t i = 0; i < total_cells; i++)
 		{
+			float d = density[i];
+			/* density */
 			density[i] -= dydx;
-
 			if(density[i] < 0.0f)
 				density[i] = 0.0f;
 
-			if(abs(heat[i]) < dydx) heat[i] = 0.0f;
-			else if (heat[i]>0.0f) heat[i] -= dydx;
-			else if (heat[i]<0.0f) heat[i] += dydx;
+			/* heat */
+			if (heat) {
+				if(abs(heat[i]) < dydx) heat[i] = 0.0f;
+				else if (heat[i]>0.0f) heat[i] -= dydx;
+				else if (heat[i]<0.0f) heat[i] += dydx;
+			}
+
+			/* color */
+			if (r && d) {
+				r[i] *= (density[i]/d);
+				g[i] *= (density[i]/d);
+				b[i] *= (density[i]/d);
+			}
 				
 		}
 	}
 }
 
-extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log)
-{
-	float *density = wt->getDensityBig();
-	Vec3Int r = wt->getResBig();
-
-	if(log)
-	{
-		/* max density/speed = dydx */
-		float dydx = 1.0 / (float)speed;
-		size_t size= r[0] * r[1] * r[2];
-
-		for(size_t i = 0; i < size; i++)
-		{
-			density[i] *= (1.0 - dydx);
-
-			if(density[i] < 0.0f)
-				density[i] = 0.0f;
-		}
-	}
-	else // linear falloff
-	{
-		/* max density/speed = dydx */
-		float dydx = 1.0 / (float)speed;
-		size_t size= r[0] * r[1] * r[2];
-
-		for(size_t i = 0; i < size; i++)
-		{
-			density[i] -= dydx;
-
-			if(density[i] < 0.0f)
-				density[i] = 0.0f;				
-		}
-	}
-}
-
-extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength)
+extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log)
 {
-	wt->initBlenderRNA(strength);
+	data_dissolve(fluid->_density, fluid->_heat, fluid->_color_r, fluid->_color_g, fluid->_color_b, fluid->_totalCells, speed, log);
 }
 
-template < class T > inline T ABS( T a )
+extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log)
 {
-	return (0 < a) ? a : -a ;
+	data_dissolve(wt->_densityBig, 0, wt->_color_rBig, wt->_color_gBig, wt->_color_bBig, wt->_totalCellsBig, speed, log);
 }
 
-extern "C" void smoke_export(FLUID_3D *fluid, float *dt, float *dx, float **dens, float **densold, float **heat, float **heatold, float **vx, float **vy, float **vz, float **vxold, float **vyold, float **vzold, unsigned char **obstacles)
+extern "C" void smoke_export(FLUID_3D *fluid, float *dt, float *dx, float **dens, float **react, float **flame, float **fuel, float **heat, 
+							 float **heatold, float **vx, float **vy, float **vz, float **r, float **g, float **b, unsigned char **obstacles)
 {
 	*dens = fluid->_density;
-	*densold = fluid->_densityOld;
+	*fuel = fluid->_fuel;
+	*react = fluid->_react;
+	*flame = fluid->_flame;
 	*heat = fluid->_heat;
 	*heatold = fluid->_heatOld;
 	*vx = fluid->_xVelocity;
 	*vy = fluid->_yVelocity;
 	*vz = fluid->_zVelocity;
-	*vxold = fluid->_xVelocityOld;
-	*vyold = fluid->_yVelocityOld;
-	*vzold = fluid->_zVelocityOld;
+	*r = fluid->_color_r;
+	*g = fluid->_color_g;
+	*b = fluid->_color_b;
 	*obstacles = fluid->_obstacles;
 	*dt = fluid->_dt;
 	*dx = fluid->_dx;
-
 }
 
-extern "C" void smoke_turbulence_export(WTURBULENCE *wt, float **dens, float **densold, float **tcu, float **tcv, float **tcw)
+extern "C" void smoke_turbulence_export(WTURBULENCE *wt, float **dens, float **react, float **flame, float **fuel,
+										float **r, float **g, float **b , float **tcu, float **tcv, float **tcw)
 {
 	if(!wt)
 		return;
 
 	*dens = wt->_densityBig;
-	*densold = wt->_densityBigOld;
+	*fuel = wt->_fuelBig;
+	*react = wt->_reactBig;
+	*flame = wt->_flameBig;
+	*r = wt->_color_rBig;
+	*g = wt->_color_gBig;
+	*b = wt->_color_bBig;
 	*tcu = wt->_tcU;
 	*tcv = wt->_tcV;
 	*tcw = wt->_tcW;
@@ -221,6 +212,16 @@ extern "C" float *smoke_get_density(FLUID_3D *fluid)
 	return fluid->_density;
 }
 
+extern "C" float *smoke_get_fuel(FLUID_3D *fluid)
+{
+	return fluid->_fuel;
+}
+
+extern "C" float *smoke_get_react(FLUID_3D *fluid)
+{
+	return fluid->_react;
+}
+
 extern "C" float *smoke_get_heat(FLUID_3D *fluid)
 {
 	return fluid->_heat;
@@ -256,15 +257,137 @@ extern "C" float *smoke_get_force_z(FLUID_3D *fluid)
 	return fluid->_zForce;
 }
 
+extern "C" float *smoke_get_flame(FLUID_3D *fluid)
+{
+	return fluid->_flame;
+}
+
+extern "C" float *smoke_get_color_r(FLUID_3D *fluid)
+{
+	return fluid->_color_r;
+}
+
+extern "C" float *smoke_get_color_g(FLUID_3D *fluid)
+{
+	return fluid->_color_g;
+}
+
+extern "C" float *smoke_get_color_b(FLUID_3D *fluid)
+{
+	return fluid->_color_b;
+}
+
+static void get_rgba(float *r, float *g, float *b, float *a, int total_cells, float *data, int sequential)
+{
+	int i;
+	int m = 4, i_g = 1, i_b = 2, i_a = 3;
+	/* sequential data */
+	if (sequential) {
+		m = 1;
+		i_g *= total_cells;
+		i_b *= total_cells;
+		i_a *= total_cells;
+	}
+
+	for (i=0; i<total_cells; i++) {
+		float alpha = a[i];
+		if (alpha) {
+			data[i*m  ] = r[i];
+			data[i*m+i_g] = g[i];
+			data[i*m+i_b] = b[i];
+		}
+		else {
+			data[i*m  ] = data[i*m+i_g] = data[i*m+i_b] = 0.0f;
+		}
+		data[i*m+i_a] = alpha;
+	}
+}
+
+extern "C" void smoke_get_rgba(FLUID_3D *fluid, float *data, int sequential)
+{
+	get_rgba(fluid->_color_r, fluid->_color_g, fluid->_color_b, fluid->_density, fluid->_totalCells, data, sequential);
+}
+
+extern "C" void smoke_turbulence_get_rgba(WTURBULENCE *wt, float *data, int sequential)
+{
+	get_rgba(wt->_color_rBig, wt->_color_gBig, wt->_color_bBig, wt->_densityBig, wt->_totalCellsBig, data, sequential);
+}
+
+/* get a single color premultiplied voxel grid */
+static void get_rgba_from_density(float color[3], float *a, int total_cells, float *data, int sequential)
+{
+	int i;
+	int m = 4, i_g = 1, i_b = 2, i_a = 3;
+	/* sequential data */
+	if (sequential) {
+		m = 1;
+		i_g *= total_cells;
+		i_b *= total_cells;
+		i_a *= total_cells;
+	}
+
+	for (i=0; i<total_cells; i++) {
+		float alpha = a[i];
+		if (alpha) {
+			data[i*m  ] = color[0] * alpha;
+			data[i*m+i_g] = color[1] * alpha;
+			data[i*m+i_b] = color[2] * alpha;
+		}
+		else {
+			data[i*m  ] = data[i*m+i_g] = data[i*m+i_b] = 0.0f;
+		}
+		data[i*m+i_a] = alpha;
+	}
+}
+
+extern "C" void smoke_get_rgba_from_density(FLUID_3D *fluid, float color[3], float *data, int sequential)
+{
+	get_rgba_from_density(color, fluid->_density, fluid->_totalCells, data, sequential);
+}
+
+extern "C" void smoke_turbulence_get_rgba_from_density(WTURBULENCE *wt, float color[3], float *data, int sequential)
+{
+	get_rgba_from_density(color, wt->_densityBig, wt->_totalCellsBig, data, sequential);
+}
+
 extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt)
 {
 	return wt ? wt->getDensityBig() : NULL;
 }
 
+extern "C" float *smoke_turbulence_get_fuel(WTURBULENCE *wt)
+{
+	return wt ? wt->getFuelBig() : NULL;
+}
+
+extern "C" float *smoke_turbulence_get_react(WTURBULENCE *wt)
+{
+	return wt ? wt->_reactBig : NULL;
+}
+
+extern "C" float *smoke_turbulence_get_color_r(WTURBULENCE *wt)
+{
+	return wt ? wt->_color_rBig : NULL;
+}
+
+extern "C" float *smoke_turbulence_get_color_g(WTURBULENCE *wt)
+{
+	return wt ? wt->_color_gBig : NULL;
+}
+
+extern "C" float *smoke_turbulence_get_color_b(WTURBULENCE *wt)
+{
+	return wt ? wt->_color_bBig : NULL;
+}
+
+extern "C" float *smoke_turbulence_get_flame(WTURBULENCE *wt)
+{
+	return wt ? wt->getFlameBig() : NULL;
+}
+
 extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res)
 {
-	if(wt)
-	{
+	if(wt) {
 		Vec3Int r = wt->getResBig();
 		res[0] = r[0];
 		res[1] = r[1];
@@ -272,6 +395,15 @@ extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res)
 	}
 }
 
+extern "C" int smoke_turbulence_get_cells(WTURBULENCE *wt)
+{
+	if(wt) {
+		Vec3Int r = wt->getResBig();
+		return r[0]*r[1]*r[2];
+	}
+	return 0;
+}
+
 extern "C" unsigned char *smoke_get_obstacle(FLUID_3D *fluid)
 {
 	return fluid->_obstacles;
@@ -295,3 +427,61 @@ extern "C" void smoke_turbulence_set_noise(WTURBULENCE *wt, int type)
 {
 	wt->setNoise(type);
 }
+
+extern "C" void flame_get_spectrum(unsigned char *spec, int width, float t1, float t2)
+{
+	spectrum(t1, t2, width, spec);
+}
+
+extern "C" int smoke_has_heat(FLUID_3D *fluid)
+{
+	return (fluid->_heat) ? 1 : 0;
+}
+
+extern "C" int smoke_has_fuel(FLUID_3D *fluid)
+{
+	return (fluid->_fuel) ? 1 : 0;
+}
+
+extern "C" int smoke_has_colors(FLUID_3D *fluid)
+{
+	return (fluid->_color_r && fluid->_color_g && fluid->_color_b) ? 1 : 0;
+}
+
+extern "C" int smoke_turbulence_has_fuel(WTURBULENCE *wt)
+{
+	return (wt->_fuelBig) ? 1 : 0;
+}
+
+extern "C" int smoke_turbulence_has_colors(WTURBULENCE *wt)
+{
+	return (wt->_color_rBig && wt->_color_gBig && wt->_color_bBig) ? 1 : 0;
+}
+
+/* additional field initialization */
+extern "C" void smoke_ensure_heat(FLUID_3D *fluid)
+{
+	if (fluid) {
+		fluid->initHeat();
+	}
+}
+
+extern "C" void smoke_ensure_fire(FLUID_3D *fluid, WTURBULENCE *wt)
+{
+	if (fluid) {
+		fluid->initFire();
+	}
+	if (wt) {
+		wt->initFire();
+	}
+}
+
+extern "C" void smoke_ensure_colors(FLUID_3D *fluid, WTURBULENCE *wt, float init_r, float init_g, float init_b)
+{
+	if (fluid) {
+		fluid->initColors(init_r, init_g, init_b);
+	}
+	if (wt) {
+		wt->initColors(init_r, init_g, init_b);
+	}
+}
\ No newline at end of file
diff --git a/intern/smoke/intern/spectrum.cpp b/intern/smoke/intern/spectrum.cpp
new file mode 100644
index 00000000000..34725c12c92
--- /dev/null
+++ b/intern/smoke/intern/spectrum.cpp
@@ -0,0 +1,411 @@
+/*
+                Colour Rendering of Spectra
+
+                       by John Walker
+                  http://www.fourmilab.ch/
+		  
+		 Last updated: March 9, 2003
+
+           This program is in the public domain.
+
+    For complete information about the techniques employed in
+    this program, see the World-Wide Web document:
+
+             http://www.fourmilab.ch/documents/specrend/
+	     
+    The xyz_to_rgb() function, which was wrong in the original
+    version of this program, was corrected by:
+    
+	    Andrew J. S. Hamilton 21 May 1999
+	    Andrew.Hamilton@Colorado.EDU
+	    http://casa.colorado.edu/~ajsh/
+
+    who also added the gamma correction facilities and
+    modified constrain_rgb() to work by desaturating the
+    colour by adding white.
+    
+    A program which uses these functions to plot CIE
+    "tongue" diagrams called "ppmcie" is included in
+    the Netpbm graphics toolkit:
+    	http://netpbm.sourceforge.net/
+    (The program was called cietoppm in earlier
+    versions of Netpbm.)
+
+*/
+
+#include <stdio.h>
+#include <math.h>
+
+/* A colour system is defined by the CIE x and y coordinates of
+   its three primary illuminants and the x and y coordinates of
+   the white point. */
+
+struct colourSystem {
+    char *name;     	    	    /* Colour system name */
+    double xRed, yRed,	    	    /* Red x, y */
+           xGreen, yGreen,  	    /* Green x, y */
+           xBlue, yBlue,    	    /* Blue x, y */
+           xWhite, yWhite,  	    /* White point x, y */
+	   gamma;   	    	    /* Gamma correction for system */
+};
+
+/* White point chromaticities. */
+
+#define IlluminantC     0.3101, 0.3162	    	/* For NTSC television */
+#define IlluminantD65   0.3127, 0.3291	    	/* For EBU and SMPTE */
+#define IlluminantE 	0.33333333, 0.33333333  /* CIE equal-energy illuminant */
+
+/*  Gamma of nonlinear correction.
+
+    See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at:
+    
+       http://www.poynton.com/ColorFAQ.html
+       http://www.poynton.com/GammaFAQ.html
+ 
+*/
+
+#define GAMMA_REC709	0		/* Rec. 709 */
+
+static struct colourSystem
+                  /* Name                  xRed    yRed    xGreen  yGreen  xBlue  yBlue    White point        Gamma   */
+    NTSCsystem  =  { "NTSC",               0.67,   0.33,   0.21,   0.71,   0.14,   0.08,   IlluminantC,    GAMMA_REC709 },
+    EBUsystem   =  { "EBU (PAL/SECAM)",    0.64,   0.33,   0.29,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 },
+    SMPTEsystem =  { "SMPTE",              0.630,  0.340,  0.310,  0.595,  0.155,  0.070,  IlluminantD65,  GAMMA_REC709 },
+    HDTVsystem  =  { "HDTV",               0.670,  0.330,  0.210,  0.710,  0.150,  0.060,  IlluminantD65,  GAMMA_REC709 },
+    CIEsystem   =  { "CIE",                0.7355, 0.2645, 0.2658, 0.7243, 0.1669, 0.0085, IlluminantE,    GAMMA_REC709 },
+    Rec709system = { "CIE REC 709",        0.64,   0.33,   0.30,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 };
+
+/*  	    	    	    UPVP_TO_XY
+
+    Given 1976 coordinates u', v', determine 1931 chromaticities x, y
+    
+*/
+
+void upvp_to_xy(double up, double vp, double *xc, double *yc)
+{
+    *xc = (9 * up) / ((6 * up) - (16 * vp) + 12);
+    *yc = (4 * vp) / ((6 * up) - (16 * vp) + 12);
+}
+
+/*  	    	    	    XY_TO_UPVP
+
+    Given 1931 chromaticities x, y, determine 1976 coordinates u', v'
+    
+*/
+
+void xy_to_upvp(double xc, double yc, double *up, double *vp)
+{
+    *up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3);
+    *vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3);
+}
+
+/*                             XYZ_TO_RGB
+
+    Given an additive tricolour system CS, defined by the CIE x
+    and y chromaticities of its three primaries (z is derived
+    trivially as 1-(x+y)), and a desired chromaticity (XC, YC,
+    ZC) in CIE space, determine the contribution of each
+    primary in a linear combination which sums to the desired
+    chromaticity.  If the  requested chromaticity falls outside
+    the Maxwell  triangle (colour gamut) formed by the three
+    primaries, one of the r, g, or b weights will be negative. 
+
+    Caller can use constrain_rgb() to desaturate an
+    outside-gamut colour to the closest representation within
+    the available gamut and/or norm_rgb to normalise the RGB
+    components so the largest nonzero component has value 1.
+    
+*/
+
+void xyz_to_rgb(struct colourSystem *cs,
+                double xc, double yc, double zc,
+                double *r, double *g, double *b)
+{
+    double xr, yr, zr, xg, yg, zg, xb, yb, zb;
+    double xw, yw, zw;
+    double rx, ry, rz, gx, gy, gz, bx, by, bz;
+    double rw, gw, bw;
+
+    xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
+    xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
+    xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);
+
+    xw = cs->xWhite;  yw = cs->yWhite;  zw = 1 - (xw + yw);
+
+    /* xyz -> rgb matrix, before scaling to white. */
+    
+    rx = (yg * zb) - (yb * zg);  ry = (xb * zg) - (xg * zb);  rz = (xg * yb) - (xb * yg);
+    gx = (yb * zr) - (yr * zb);  gy = (xr * zb) - (xb * zr);  gz = (xb * yr) - (xr * yb);
+    bx = (yr * zg) - (yg * zr);  by = (xg * zr) - (xr * zg);  bz = (xr * yg) - (xg * yr);
+
+    /* White scaling factors.
+       Dividing by yw scales the white luminance to unity, as conventional. */
+       
+    rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw;
+    gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw;
+    bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw;
+
+    /* xyz -> rgb matrix, correctly scaled to white. */
+    
+    rx = rx / rw;  ry = ry / rw;  rz = rz / rw;
+    gx = gx / gw;  gy = gy / gw;  gz = gz / gw;
+    bx = bx / bw;  by = by / bw;  bz = bz / bw;
+
+    /* rgb of the desired point */
+
+    *r = (rx * xc) + (ry * yc) + (rz * zc);
+    *g = (gx * xc) + (gy * yc) + (gz * zc);
+    *b = (bx * xc) + (by * yc) + (bz * zc);
+}
+
+/*                            INSIDE_GAMUT
+
+     Test whether a requested colour is within the gamut
+     achievable with the primaries of the current colour
+     system.  This amounts simply to testing whether all the
+     primary weights are non-negative. */
+
+int inside_gamut(double r, double g, double b)
+{
+    return (r >= 0) && (g >= 0) && (b >= 0);
+}
+
+/*                          CONSTRAIN_RGB
+
+    If the requested RGB shade contains a negative weight for
+    one of the primaries, it lies outside the colour gamut 
+    accessible from the given triple of primaries.  Desaturate
+    it by adding white, equal quantities of R, G, and B, enough
+    to make RGB all positive.  The function returns 1 if the
+    components were modified, zero otherwise.
+    
+*/
+
+int constrain_rgb(double *r, double *g, double *b)
+{
+    double w;
+
+    /* Amount of white needed is w = - min(0, *r, *g, *b) */
+    
+    w = (0 < *r) ? 0 : *r;
+    w = (w < *g) ? w : *g;
+    w = (w < *b) ? w : *b;
+    w = -w;
+
+    /* Add just enough white to make r, g, b all positive. */
+    
+    if (w > 0) {
+        *r += w;  *g += w; *b += w;
+        return 1;                     /* Colour modified to fit RGB gamut */
+    }
+
+    return 0;                         /* Colour within RGB gamut */
+}
+
+/*                          GAMMA_CORRECT_RGB
+
+    Transform linear RGB values to nonlinear RGB values. Rec.
+    709 is ITU-R Recommendation BT. 709 (1990) ``Basic
+    Parameter Values for the HDTV Standard for the Studio and
+    for International Programme Exchange'', formerly CCIR Rec.
+    709. For details see
+    
+       http://www.poynton.com/ColorFAQ.html
+       http://www.poynton.com/GammaFAQ.html
+*/
+
+void gamma_correct(const struct colourSystem *cs, double *c)
+{
+    double gamma;
+
+    gamma = cs->gamma;
+
+    if (gamma == GAMMA_REC709) {
+	/* Rec. 709 gamma correction. */
+	double cc = 0.018;
+	
+	if (*c < cc) {
+	    *c *= ((1.099 * pow(cc, 0.45)) - 0.099) / cc;
+	} else {
+	    *c = (1.099 * pow(*c, 0.45)) - 0.099;
+	}
+    } else {
+	/* Nonlinear colour = (Linear colour)^(1/gamma) */
+	*c = pow(*c, 1.0 / gamma);
+    }
+}
+
+void gamma_correct_rgb(const struct colourSystem *cs, double *r, double *g, double *b)
+{
+    gamma_correct(cs, r);
+    gamma_correct(cs, g);
+    gamma_correct(cs, b);
+}
+
+/*  	    	    	    NORM_RGB
+
+    Normalise RGB components so the most intense (unless all
+    are zero) has a value of 1.
+    
+*/
+
+void norm_rgb(double *r, double *g, double *b)
+{
+#define Max(a, b)   (((a) > (b)) ? (a) : (b))
+    double greatest = Max(*r, Max(*g, *b));
+
+    if (greatest > 0) {
+    	*r /= greatest;
+	*g /= greatest;
+	*b /= greatest;
+    }
+#undef Max
+}
+
+/*                          SPECTRUM_TO_XYZ
+
+    Calculate the CIE X, Y, and Z coordinates corresponding to
+    a light source with spectral distribution given by  the
+    function SPEC_INTENS, which is called with a series of
+    wavelengths between 380 and 780 nm (the argument is 
+    expressed in meters), which returns emittance at  that
+    wavelength in arbitrary units.  The chromaticity
+    coordinates of the spectrum are returned in the x, y, and z
+    arguments which respect the identity:
+
+            x + y + z = 1.
+*/
+
+void spectrum_to_xyz(double (*spec_intens)(double wavelength),
+                     double *x, double *y, double *z)
+{
+    int i;
+    double lambda, X = 0, Y = 0, Z = 0, XYZ;
+
+    /* CIE colour matching functions xBar, yBar, and zBar for
+       wavelengths from 380 through 780 nanometers, every 5
+       nanometers.  For a wavelength lambda in this range:
+
+            cie_colour_match[(lambda - 380) / 5][0] = xBar
+            cie_colour_match[(lambda - 380) / 5][1] = yBar
+            cie_colour_match[(lambda - 380) / 5][2] = zBar
+
+	To save memory, this table can be declared as floats
+	rather than doubles; (IEEE) float has enough 
+	significant bits to represent the values. It's declared
+	as a double here to avoid warnings about "conversion
+	between floating-point types" from certain persnickety
+	compilers. */
+
+    static double cie_colour_match[81][3] = {
+        {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
+        {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
+        {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
+        {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
+        {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
+        {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
+        {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
+        {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
+        {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
+        {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
+        {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
+        {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
+        {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
+        {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
+        {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
+        {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
+        {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
+        {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
+        {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
+        {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
+        {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
+        {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
+        {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
+        {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
+        {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
+        {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
+        {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
+    };
+
+    for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
+        double Me;
+
+        Me = (*spec_intens)(lambda);
+        X += Me * cie_colour_match[i][0];
+        Y += Me * cie_colour_match[i][1];
+        Z += Me * cie_colour_match[i][2];
+    }
+    XYZ = (X + Y + Z);
+    *x = X / XYZ;
+    *y = Y / XYZ;
+    *z = Z / XYZ;
+}
+
+/*                            BB_SPECTRUM
+
+    Calculate, by Planck's radiation law, the emittance of a black body
+    of temperature bbTemp at the given wavelength (in metres).  */
+
+double bbTemp = 5000;                 /* Hidden temperature argument
+                                         to BB_SPECTRUM. */
+double bb_spectrum(double wavelength)
+{
+    double wlm = wavelength * 1e-9;   /* Wavelength in meters */
+
+    return (3.74183e-16 * pow(wlm, -5.0)) /
+           (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0);
+}
+
+void xyz_to_lms(double x, double y, double z, double* l, double* m, double* s)
+{
+	*l =  0.3897*x + 0.6890*y - 0.0787*z;
+	*m = -0.2298*x + 1.1834*y + 0.0464*z;
+	*s = z;
+}
+
+void lms_to_xyz(double l, double m, double s, double* x, double *y, double* z)
+{
+	*x =  1.9102*l - 1.1121*m + 0.2019*s;
+	*y =  0.3709*l + 0.6290*m + 0.0000*s;
+	*z = s;
+}
+
+void spectrum(double t1, double t2, int N, unsigned char* d)
+{
+	int i,j,dj;
+	double X,Y,Z,R,G,B,L,M,S, Lw, Mw, Sw;
+	struct colourSystem *cs = &CIEsystem;
+
+	j = 0; dj = 1;
+	if (t1<t2) {
+		double t = t1;
+		t1 = t2;
+		t2 = t;
+		j = N-1; dj=-1;
+	}
+
+	for (i=0; i<N; i++) {
+		bbTemp = t1 + (t2-t1)/N*i;
+
+			// integrate blackbody radiation spectrum to XYZ
+		spectrum_to_xyz(bb_spectrum, &X, &Y, &Z);
+
+			// normalize highest temperature to white (in LMS system)
+		xyz_to_lms(X,Y,Z,&L,&M,&S);
+		if (i==0) {
+			Lw=1/L; Mw=1/M; Sw=1/S;
+		}
+		L *= Lw; M *= Mw; S *= Sw;
+		lms_to_xyz(L,M,S,&X,&Y,&Z);
+
+			// convert to RGB
+		xyz_to_rgb(cs, X, Y, Z, &R, &G, &B);
+		constrain_rgb(&R, &G, &B);
+		norm_rgb(&R, &G, &B);
+		d[(j<<2)] = (unsigned char) ((double)R*255);
+		d[(j<<2)+1] = (unsigned char) ((double)G*255);
+		d[(j<<2)+2] = (unsigned char) ((double)B*255);
+		d[(j<<2)+3] = (B>0.1)? B*255 : 0;
+		j += dj;
+	}
+}
diff --git a/intern/smoke/intern/spectrum.h b/intern/smoke/intern/spectrum.h
new file mode 100644
index 00000000000..9edd9ad887c
--- /dev/null
+++ b/intern/smoke/intern/spectrum.h
@@ -0,0 +1,6 @@
+#ifndef __SPECTRUM_H
+#define __SPECTRUM_H
+
+void spectrum(double t1, double t2, int n, unsigned char* d);
+
+#endif