Merge lp:~georg-zotti/stellarium/gz_comet-tail-extinction into lp:stellarium

[gzotti]

Improved brightness control and colors for comet tails: extinction and twilight mixing.

My original comet tails of winter 2013/14 had one big flaw: brightness was computed globally. This made them almost pop out of twilight in full length, and bright comets were just too bright.

This branch has I think much better brightness balancing (query atmosphere and eye adaptation), and can reduce/influence tail brightness not only by vertex-wise extinction, but also by vertex-wise adaptation of tail length visibility to overall sky brightness. The increase in quality can at least be called "significant" IMHO.

To test, I recommend you copy ssystem_1000comets.ini to your ssystem.ini and have a look at e.g. comets Hale-Bopp (best in March/April 1997), Hyakutake (March 1996), PanSTARRS (March 2013), Halley (1986, 1910, 837), Bennett (1970), 1680, 1882, Donati (1858).

[Fabien Chéreau (xalioth)]

Looks good to me :)

review: Approve

=== modified file 'src/core/StelApp.hpp'
--- src/core/StelApp.hpp	2014-11-04 10:59:12 +0000
+++ src/core/StelApp.hpp	2014-11-08 17:16:48 +0000
@@ -148,7 +148,8 @@
 	void update(double deltaTime);
 
 	//! Draw all registered StelModule in the order defined by the order lists.
-	//! @return the max squared distance in pixels that any object has travelled since the last update.
+	// 2014-11: OLD COMMENT? What does a void return?
+	// @return the max squared distance in pixels that any object has travelled since the last update.
 	void draw();
 
 	//! Call this when the size of the GL window has changed.

=== modified file 'src/core/StelPainter.cpp'
--- src/core/StelPainter.cpp	2014-11-04 12:15:15 +0000
+++ src/core/StelPainter.cpp	2014-11-08 17:16:48 +0000
@@ -544,7 +544,7 @@
 
 		// Translate/rotate
 		if (!noGravity)
-			angleDeg += prj->defautAngleForGravityText;
+			angleDeg += prj->defaultAngleForGravityText;
 
 		if (std::fabs(angleDeg)>1.f)
 		{
@@ -1357,6 +1357,7 @@
 			break;
 		case SphericalPolygonDrawModeFill:
 		case SphericalPolygonDrawModeTextureFill:
+		case SphericalPolygonDrawModeTextureFillColormodulated:
 			glEnable(GL_CULL_FACE);
 			// The polygon is already tesselated as triangles
 			if (doSubDivise || prj->intersectViewportDiscontinuity(poly->getBoundingCap()))

=== modified file 'src/core/StelProjector.cpp'
--- src/core/StelProjector.cpp	2013-10-01 22:05:10 +0000
+++ src/core/StelProjector.cpp	2014-11-08 17:16:48 +0000
@@ -116,7 +116,7 @@
 	viewportCenter = params.viewportCenter;
 	viewportCenter *= devicePixelsPerPixel;
 	gravityLabels = params.gravityLabels;
-	defautAngleForGravityText = params.defautAngleForGravityText;
+	defaultAngleForGravityText = params.defautAngleForGravityText;
 	flipHorz = params.flipHorz ? -1.f : 1.f;
 	flipVert = params.flipVert ? -1.f : 1.f;
 	viewportFovDiameter = params.viewportFovDiameter * devicePixelsPerPixel;

=== modified file 'src/core/StelProjector.hpp'
--- src/core/StelProjector.hpp	2014-10-16 17:32:50 +0000
+++ src/core/StelProjector.hpp	2014-11-08 17:16:48 +0000
@@ -291,7 +291,7 @@
 		  oneOverZNearMinusZFar(0.f),
 		  viewportFovDiameter(0.f),
 		  gravityLabels(true),
-		  defautAngleForGravityText(0.f),
+		  defaultAngleForGravityText(0.f),
 		  devicePixelsPerPixel(1.f) {;}
 
 	//! Return whether the projection presents discontinuities. Used for optimization.
@@ -317,7 +317,7 @@
 	Vec2f viewportCenter;               // Viewport center in screen pixel
 	float viewportFovDiameter;          // diameter of the FOV disk in pixel
 	bool gravityLabels;                 // should label text align with the horizon?
-	float defautAngleForGravityText;    // a rotation angle to apply to gravity text (only if gravityLabels is set to false)
+	float defaultAngleForGravityText;   // a rotation angle to apply to gravity text (only if gravityLabels is set to false)
 	SphericalCap boundingCap;           // Bounding cap of the whole viewport
 	float devicePixelsPerPixel;         // The number of device pixel per "Device Independent Pixels" (value is usually 1, but 2 for mac retina screens)
 

=== modified file 'src/core/StelToneReproducer.hpp'
--- src/core/StelToneReproducer.hpp	2014-11-04 10:59:12 +0000
+++ src/core/StelToneReproducer.hpp	2014-11-08 17:16:48 +0000
@@ -79,7 +79,7 @@
 	//! It is the parameter to modify to simulate aperture*exposition time
 	//! @param scale the global input scale
 	void setInputScale(float scale=1.f);
-	//! Get the global scale applied to input lumiances, i.e before the adaptation
+	//! Get the global scale applied to input luminances, i.e before the adaptation
 	float getInputScale() const
 	{
 		return inputScale;
@@ -125,14 +125,14 @@
 	
 	//! Return adapted luminance from world to display with 1 corresponding to full display white
 	//! @param worldLuminance the world luminance to convert in cd/m^2
-	//! @return the converted display luminance with 1 corresponding to full display white. The value can be more than 1 when saturation..
+	//! @return the converted display luminance with 1 corresponding to full display white. The value can be more than 1 when saturated.
 	float adaptLuminanceScaled(float worldLuminance) const
 	{
 		return adaptLuminance(worldLuminance)*oneOverMaxdL;
 	}
 	
 	//! Return adapted luminance from display to world with 1 corresponding to full display white
-	//! @param displayLuminance the display luminance with 1 corresponding to full display white. The value can be more than 1 when saturation..
+	//! @param displayLuminance the display luminance with 1 corresponding to full display white. The value can be more than 1 when saturated.
 	//! @return the converted world luminance in cd/m^2
 	float reverseAdaptLuminanceScaled(float displayLuminance) const
 	{
@@ -141,8 +141,8 @@
 	
 	//! Return adapted ln(luminance) from world to display with 1 corresponding to full display white
 	//! @param lnWorldLuminance the world luminance to convert in ln(cd/m^2)
-	//! @param pFact the power at whihc the result should be set. The default is 0.5 and therefore return the square root of the adapted luminance
-	//! @return the converted display set at the pFact power. Luminance with 1 corresponding to full display white. The value can be more than 1 when saturation..
+	//! @param pFact the power at which the result should be set. The default is 0.5 and therefore return the square root of the adapted luminance
+	//! @return the converted display set at the pFact power. Luminance with 1 corresponding to full display white. The value can be more than 1 when saturated.
 	float adaptLuminanceScaledLn(float lnWorldLuminance, float pFact=0.5f) const
 	{
 		const float lnPix0p0001 = -8.0656104861f;

=== modified file 'src/core/modules/Comet.cpp'
--- src/core/modules/Comet.cpp	2014-10-19 05:54:43 +0000
+++ src/core/modules/Comet.cpp	2014-11-08 17:16:48 +0000
@@ -1,6 +1,7 @@
 /*
  * Stellarium
  * Copyright (C) 2010 Bogdan Marinov
+ * Copyright (C) 2014 Georg Zotti (Tails)
  *
  * This program is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License
@@ -26,6 +27,7 @@
 
 #include "StelTexture.hpp"
 #include "StelTextureMgr.hpp"
+#include "StelToneReproducer.hpp"
 #include "StelTranslator.hpp"
 #include "StelUtils.hpp"
 #include "StelFileMgr.hpp"
@@ -71,6 +73,7 @@
 		  true, //halo
 		  pTypeStr),
 	  tailActive(false),
+	  tailBright(false),
 	  dustTailWidthFactor(dustTailWidthFact),
 	  dustTailLengthFactor(dustTailLengthFact),
 	  dustTailBrightnessFactor(dustTailBrightnessFact)
@@ -90,6 +93,8 @@
 	gastailVertexArr.clear();
 	dusttailVertexArr.clear();
 	comaVertexArr.clear();
+	gastailColorArr.clear();
+	dusttailColorArr.clear();
 
 	//Comet specific members
 	absoluteMagnitude = 0;
@@ -281,19 +286,27 @@
 	return apparentMagnitude;
 }
 
-// Compute the position in the parent Planet coordinate system
-// Actually call the provided function to compute the ecliptical position, and buildup the tails!
-void Comet::computePosition(const double date)
+void Comet::update(int deltaTime)
 {
-	Planet::computePosition(date);
-	//GZ: I think we can make deltaJD adaptive, depending on distance to sun! For some reason though, this leads to a crash!
-	//deltaJD=StelCore::JD_SECOND * qMax(1.0, qMin(eclipticPos.length(), 20.0));
+	Planet::update(deltaTime);
+
+	// The rest used to be in computePosition(), but is better in update(). Unfortunately we need date (JD).
+	StelCore* core=StelApp::getInstance().getCore();
+	double date=core->getJDay();
+
+	// The CometOrbit is in fact available in userDataPtr!
+	CometOrbit* orbit=(CometOrbit*)userDataPtr;
+	Q_ASSERT(orbit);
+	if (!orbit->objectDateValid(core->getJDay())) return; // don't do anything if out of useful date range. This allows having hundreds of comet elements.
+
+
+	//GZ: I think we can make deltaJDtail adaptive, depending on distance to sun! For some reason though, this leads to a crash!
+	//deltaJDtail=StelCore::JD_SECOND * qMax(1.0, qMin(eclipticPos.length(), 20.0));
 
 	if (fabs(lastJDtail-date)>deltaJDtail)
 	{
 		lastJDtail=date;
 
-		// GZ: Moved from draw() :-)
 		// The CometOrbit is in fact available in userDataPtr!
 		CometOrbit* orbit=(CometOrbit*)userDataPtr;
 		Q_ASSERT(orbit);
@@ -301,11 +314,10 @@
 
 		if (orbit->getUpdateTails()){
 			// Compute lengths and orientations from orbit object, but only if required.
-			// This part moved from draw() to keep draw() free from too much computation.
 			tailFactors=getComaDiameterAndTailLengthAU();
 
 			// Note that we use a diameter larger than what the formula returns. A scale factor of 1.2 is ad-hoc/empirical (GZ), but may look better.
-			computeComa(1.0f*tailFactors[0]);
+			computeComa(1.0f*tailFactors[0]); // TBD: APPARENTLY NO SCALING? REMOVE 1.0 and note above.
 
 			tailActive = (tailFactors[1] > tailFactors[0]); // Inhibit tails drawing if too short. Would be nice to include geometric projection angle, but this is too costly.
 
@@ -318,10 +330,10 @@
 
 				// Find valid parameters to create paraboloid vertex arrays: dustTail, gasTail.
 				computeParabola(gasparameter, gasTailEndRadius, -0.5f*gasparameter, gastailVertexArr,  tailTexCoordArr, tailIndices);
-				// This was for a rotated straight parabola:
-				//computeParabola(dustparameter, 2.0f*tailFactors[0], -0.5f*dustparameter, dusttailVertexArr, dusttailTexCoordArr, dusttailIndices);
+				//gastailColorArr.fill(Vec3f(0.3,0.3,0.3), gastailVertexArr.length());
 				// Now we make a skewed parabola. Skew factor (xOffset, last arg) is rather ad-hoc/empirical. TBD later: Find physically correct solution.
 				computeParabola(dustparameter, dustTailWidthFactor*gasTailEndRadius, -0.5f*dustparameter, dusttailVertexArr, tailTexCoordArr, tailIndices, 25.0f*orbit->getVelocity().length());
+				//dusttailColorArr.fill(Vec3f(0.3,0.3,0.3), dusttailVertexArr.length());
 
 
 				// 2014-08 for 0.13.1 Moved from drawTail() to save lots of computation per frame (There *are* folks downloading all 730 MPC current comet elements...)
@@ -336,7 +348,7 @@
 				// In addition, we let the dust tail already start with a light tilt.
 				dustTailRot=gasTailRot * Mat4d::zrotation(atan2(velocity[1], velocity[0]) + M_PI) * Mat4d::yrotation(5.0f*velocity.length());
 
-				// TODO: If we want to be even faster, rotate vertex arrays here and not in drawTail()!
+				// Rotate vertex arrays:
 				Vec3d* gasVertices=(Vec3d*) (gastailVertexArr.data());
 				Vec3d* dustVertices=(Vec3d*) (dusttailVertexArr.data());
 				for (int i=0; i<COMET_TAIL_SLICES*COMET_TAIL_STACKS+1; ++i)
@@ -345,34 +357,105 @@
 					dustVertices[i].transfo4d(dustTailRot);
 				}
 			}
-
 			orbit->setUpdateTails(false); // don't update until position has been recalculated elsewhere
 		}
-		// Note: we can make deltaJDtail adaptive, depending on distance to sun!
-		//deltaJDtail=5.0*StelCore::JD_MINUTE * qMax(0.01, qMin(eclipticPos.length(), 20.0));
-	}
+	}
+
+	// And also update magnitude and tail brightness/extinction here.
+	const bool withAtmosphere=(core->getSkyDrawer()->getFlagHasAtmosphere());
+
+	StelToneReproducer* eye = core->getToneReproducer();
+	float lum = core->getSkyDrawer()->surfacebrightnessToLuminance(getVMagnitude(core)+13.0f); // How to calibrate?
+	// Get the luminance scaled between 0 and 1
+	float aLum =eye->adaptLuminanceScaled(lum);
+
+
+	// To make comet more apparent in overviews, take field of view into account:
+	const float fov=core->getProjection(core->getAltAzModelViewTransform())->getFov();
+	if (fov>20)
+		aLum*= (fov/20.0f);
+
+	// Now inhibit tail drawing if still too dim.
+	if (aLum<0.002f)
+	{
+		// Far too dim: don't even show tail...
+		tailBright=false;
+		return;
+	} else
+		tailBright=true;
+
+	// Separate factors, but avoid overly bright tails. I limit to about 0.7 for overlapping both tails which should not exceed full-white.
+	float gasMagFactor=qMin(0.9f*aLum, 0.7f);
+	float dustMagFactor=qMin(dustTailBrightnessFactor*aLum, 0.7f);
+
+	Vec3f gasColor(0.15f*gasMagFactor,0.35f*gasMagFactor,0.6f*gasMagFactor); // Orig color 0.15/0.15/0.6
+	Vec3f dustColor(dustMagFactor, dustMagFactor,0.6f*dustMagFactor);
+
+	if (withAtmosphere)
+	{
+		Extinction extinction=core->getSkyDrawer()->getExtinction();
+
+		// Not only correct the color values for extinction, but for twilight conditions, also make tail end less visible.
+		// I consider sky brightness over 1cd/m^2 as reason to shorten tail.
+		// Below this brightness, the tail brightness loss by this method is insignificant:
+		// Just counting through the vertices might make a spiral apperance. Maybe even better than stackwise? Let's see...
+		const float avgAtmLum=GETSTELMODULE(LandscapeMgr)->getAtmosphereAverageLuminance();
+		const float brightnessDecreasePerVertexFromHead=1.0f/(COMET_TAIL_SLICES*COMET_TAIL_STACKS)  * avgAtmLum;
+		float brightnessPerVertexFromHead=1.0f;
+
+		gastailColorArr.clear();
+		dusttailColorArr.clear();
+		for (int i=0; i<gastailVertexArr.size(); ++i)
+		{
+			// Gastail extinction:
+			Vec3d vertAltAz=core->j2000ToAltAz(gastailVertexArr.at(i), StelCore::RefractionOn);
+			vertAltAz.normalize();
+			Q_ASSERT(fabs(vertAltAz.lengthSquared()-1.0) < 0.001);
+			float oneMag=0.0f;
+			extinction.forward(vertAltAz, &oneMag);
+			float extinctionFactor=std::pow(0.4f, oneMag); // drop of one magnitude: factor 2.5 or 40%
+			gastailColorArr.append(gasColor*extinctionFactor* brightnessPerVertexFromHead);
+
+			// dusttail extinction:
+			vertAltAz=core->j2000ToAltAz(dusttailVertexArr.at(i), StelCore::RefractionOn);
+			vertAltAz.normalize();
+			Q_ASSERT(fabs(vertAltAz.lengthSquared()-1.0) < 0.001);
+			oneMag=0.0f;
+			extinction.forward(vertAltAz, &oneMag);
+			extinctionFactor=std::pow(0.4f, oneMag); // drop of one magnitude: factor 2.5 or 40%
+			dusttailColorArr.append(dustColor*extinctionFactor * brightnessPerVertexFromHead);
+
+			brightnessPerVertexFromHead-=brightnessDecreasePerVertexFromHead;
+		}
+	}
+	else // no atmosphere: set all vertices to same brightness.
+	{
+		gastailColorArr.fill(gasColor,   gastailVertexArr.length());
+		dusttailColorArr.fill(dustColor, dusttailVertexArr.length());
+	}
+	//qDebug() << "Comet " << getEnglishName() <<  "JD: " << date << "gasR" << gasColor[0] << " dustR" << dustColor[0];
 }
 
 
-// Draw the Comet and all the related infos : name, circle etc... GZ: Taken from Planet.cpp 2013-11-05 and extended
+// Draw the Comet and all the related infos: name, circle etc... GZ: Taken from Planet.cpp 2013-11-05 and extended
 void Comet::draw(StelCore* core, float maxMagLabels, const QFont& planetNameFont)
 {
 	if (hidden)
 		return;
 	if (getEnglishName() == core->getCurrentLocation().planetName)
-	{ // GZ: Maybe even don't do that? E.g., draw tail while riding the comet? Decide later.
+	{ // Maybe even don't do that? E.g., draw tail while riding the comet? Decide later.
 		return;
 	}
 
-	// GZ: If comet is too faint to be seen, don't bother rendering. (oops, should have been here in 2014-01... ;-)
-	if ((getVMagnitude(core)-2.0f) > core->getSkyDrawer()->getLimitMagnitude())
+	// If comet is too faint to be seen, don't bother rendering. (Massive speedup if people have hundreds of comets!)
+	if ((getVMagnitude(core)+2.0f) > core->getSkyDrawer()->getLimitMagnitude())
 	{
 		return;
 	}
 	// The CometOrbit is in fact available in userDataPtr!
 	CometOrbit* orbit=(CometOrbit*)userDataPtr;
 	Q_ASSERT(orbit);
-	if (!orbit->objectDateValid(core->getJDay())) return; // don't draw at all out of useful date range. This allows having hundreds of comet elements.
+	if (!orbit->objectDateValid(core->getJDay())) return; // don't draw at all if out of useful date range. This allows having hundreds of comet elements.
 
 	Mat4d mat = Mat4d::translation(eclipticPos) * rotLocalToParent;
 	// This removed totally the Planet shaking bug!!!
@@ -408,8 +491,8 @@
 
 		draw3dModel(core,transfo,screenSz);
 	}
-	// tails should also be drawn if comet core is off-screen...
-	if (tailActive)
+	// but tails should also be drawn if comet core is off-screen...
+	if (tailActive && tailBright)
 	{
 		drawTail(core,transfo,true);  // gas tail
 		drawTail(core,transfo,false); // dust tail
@@ -420,46 +503,24 @@
 }
 
 void Comet::drawTail(StelCore* core, StelProjector::ModelViewTranformP transfo, bool gas)
-{
+{	
 	StelPainter* sPainter = new StelPainter(core->getProjection(transfo));
 	glEnable(GL_BLEND);
 	glBlendFunc(GL_ONE, GL_ONE);
 	glDisable(GL_CULL_FACE);
 
-	// GZ: If we use getVMagnitudeWithExtinction(), a head extincted in the horizon mist can completely hide an otherwise frighteningly long tail.
-	// we must use unextincted mag, but mix/dim with atmosphere/sky brightness.
-	// In addition, light falloff is a bit reduced for better visibility. Power basis should be 0.4, we use 0.6.
-	float magFactor=std::pow(0.6f , getVMagnitude(core));
-	if (core->getSkyDrawer()->getFlagHasAtmosphere())
-	{
-		// Mix with sky brightness and light pollution: This is very ad-hoc, if someone finds a better solution, please go ahead!
-		// Light pollution:
-		float bortleIndexFactor=0.1f * (11 - core->getSkyDrawer()->getBortleScaleIndex());
-		magFactor*= bortleIndexFactor*bortleIndexFactor; // GZ-Guesstimate for light pollution influence
-		// sky brightness: This is about 10 for twilight where bright comet tails should already be visible. Dark night is close to 0.
-		float avgAtmLum=GETSTELMODULE(LandscapeMgr)->getAtmosphereAverageLuminance();
-		float atmLumFactor=(15.0f-avgAtmLum)/15.0f;  if (atmLumFactor<0.05f) atmLumFactor=0.05f;    //atmLumFactor=std::sqrt(atmLumFactor);
-		magFactor*=atmLumFactor*atmLumFactor;
-	}
-	magFactor*=(gas? 0.9 : dustTailBrightnessFactor); // TBD: empirical adjustment for texture brightness.
-	magFactor=qMin(magFactor, 1.05f); // Limit excessively bright display.
-
 	tailTexture->bind();
 
 	if (gas) {
-		sPainter->setColor(0.15f*magFactor,0.15f*magFactor,0.6f*magFactor);
-		sPainter->setArrays((Vec3d*)gastailVertexArr.constData(), (Vec2f*)tailTexCoordArr.constData());
+		sPainter->setArrays((Vec3d*)gastailVertexArr.constData(), (Vec2f*)tailTexCoordArr.constData(), (Vec3f*)gastailColorArr.constData());
 		sPainter->drawFromArray(StelPainter::Triangles, tailIndices.size(), 0, true, tailIndices.constData());
 
 	} else {
-		sPainter->setColor(magFactor, magFactor,0.6f*magFactor);
-		sPainter->setArrays((Vec3d*)dusttailVertexArr.constData(), (Vec2f*)tailTexCoordArr.constData());
+		sPainter->setArrays((Vec3d*)dusttailVertexArr.constData(), (Vec2f*)tailTexCoordArr.constData(), (Vec3f*)dusttailColorArr.constData());
 		sPainter->drawFromArray(StelPainter::Triangles, tailIndices.size(), 0, true, tailIndices.constData());
 	}
 	glDisable(GL_BLEND);
 
-
-
 	if (sPainter)
 		delete sPainter;
 	sPainter=NULL;
@@ -478,16 +539,13 @@
 	glBlendFunc(GL_ONE, GL_ONE);
 	glDisable(GL_CULL_FACE);
 
-	// GZ: For the coma, we can use extinction via atmosphere.
-	// In addition, light falloff is a bit reduced for better visibility. Power basis should be 0.4, we use 0.6.
-	float minSkyMag=core->getSkyDrawer()->getLimitMagnitude();
-	float mag100pct=minSkyMag-6.0f; // should be 5, but let us draw it a bit brighter.
-	float magDrop=getVMagnitudeWithExtinction(core)-mag100pct;
-	float magFactor=std::pow(0.6f , magDrop);
-	magFactor=qMin(magFactor, 2.0f); // Limit excessively bright display.
-
+	StelToneReproducer* eye = core->getToneReproducer();
+	float lum = core->getSkyDrawer()->surfacebrightnessToLuminance(getVMagnitudeWithExtinction(core)+11.0f); // How to calibrate?
+	// Get the luminance scaled between 0 and 1
+	float aLum =eye->adaptLuminanceScaled(lum);
+	float magFactor=qMin(qMax(aLum, 0.25f), 2.0f);
 	comaTexture->bind();
-	sPainter->setColor(magFactor,magFactor,0.6f*magFactor);
+	sPainter->setColor(0.3f*magFactor,0.7*magFactor,magFactor);
 	sPainter->setArrays((Vec3d*)comaVertexArr.constData(), (Vec2f*)comaTexCoordArr.constData());
 	sPainter->drawFromArray(StelPainter::Triangles, comaVertexArr.size()/3);
 
@@ -521,11 +579,11 @@
 // Parabola equation: z=x²/2p.
 // xOffset for the dust tail, this may introduce a bend. Units are x per sqrt(z).
 void Comet::computeParabola(const float parameter, const float radius, const float zshift,
-			    QVector<double>& vertexArr, QVector<float>& texCoordArr, QVector<unsigned short> &indices, const float xOffset) {
+			    QVector<Vec3d>& vertexArr, QVector<float>& texCoordArr, QVector<unsigned short> &indices, const float xOffset) {
 
-	// GZ: keep the array and replace contents. However, using replace() is only slightly faster.
-	if (vertexArr.length() < (3*(COMET_TAIL_SLICES*COMET_TAIL_STACKS+1)))
-		vertexArr.resize(3*(COMET_TAIL_SLICES*COMET_TAIL_STACKS+1));
+	// keep the array and replace contents. However, using replace() is only slightly faster.
+	if (vertexArr.length() < ((COMET_TAIL_SLICES*COMET_TAIL_STACKS+1)))
+		vertexArr.resize((COMET_TAIL_SLICES*COMET_TAIL_STACKS+1));
 	if (createTailIndices) indices.clear();
 	if (createTailTextureCoords) texCoordArr.clear();
 	int i;
@@ -541,8 +599,8 @@
 		ya[i]=cos(i*da);
 	}
 	
-	vertexArr.replace(0, 0.0); vertexArr.replace(1, 0.0); vertexArr.replace(2, zshift);
-	int vertexArrIndex=3;
+	vertexArr.replace(0, Vec3d(0.0, 0.0, zshift));
+	int vertexArrIndex=1;
 	if (createTailTextureCoords) texCoordArr << 0.5f << 0.5f;
 	// define the indices lying on circles, starting at 1: odd rings have 1/slices+1/2slices, even-numbered rings straight 1/slices
 	// inner ring#1
@@ -553,9 +611,7 @@
 		for (i=ring & 1; i<2*COMET_TAIL_SLICES; i+=2) { // i.e., ring1 has shifted vertices, ring2 has even ones.
 			x=xa[i]*radius*ring/COMET_TAIL_STACKS;
 			y=ya[i]*radius*ring/COMET_TAIL_STACKS;
-			vertexArr.replace(vertexArrIndex++, x+xShift);
-			vertexArr.replace(vertexArrIndex++, y);
-			vertexArr.replace(vertexArrIndex++, z);
+			vertexArr.replace(vertexArrIndex++, Vec3d(x+xShift, y, z));
 			if (createTailTextureCoords) texCoordArr << 0.5+ 0.5*x/radius << 0.5+0.5*y/radius;
 		}
 	}
@@ -591,10 +647,11 @@
 	createTailTextureCoords=false;
 }
 
+
 // These are to avoid having index arrays for each comet when all are equal.
 bool Comet::createTailIndices=true;
 bool Comet::createTailTextureCoords=true;
 StelTextureSP Comet::comaTexture;
 StelTextureSP Comet::tailTexture;
-QVector<float> Comet::tailTexCoordArr; // computed only once FOR ALL COMETS!
-QVector<unsigned short> Comet::tailIndices; // computed only once FOR ALL COMETS!
+QVector<float> Comet::tailTexCoordArr; // computed only once for all Comets.
+QVector<unsigned short> Comet::tailIndices; // computed only once for all Comets.

=== modified file 'src/core/modules/Comet.hpp'
--- src/core/modules/Comet.hpp	2014-08-31 08:35:52 +0000
+++ src/core/modules/Comet.hpp	2014-11-08 17:16:48 +0000
@@ -32,6 +32,7 @@
 	2013-12: GZ: New algorithms for position computation following Paul Heafner: Fundamental Ephemeris Computations (Willmann-Bell 1999).
 	2014-01: GZ: Parabolic tails appropriately scaled/rotated. Much is currently empirical, leaving room for physics-based improvements.
 	2014-08: GZ: speedup in case hundreds of comets are loaded.
+	2014-11: GZ: tail extinction, better brightness balance.
   */
 class Comet : public Planet
 {
@@ -91,11 +92,15 @@
 	virtual double getSiderealPeriod() const;
 
 	//! GZ: override from Planet: extend with tail details.
-	virtual void computePosition(const double date);
+	//virtual void computePosition(const double date);
 
 	//! re-implementation of Planet's draw()
 	virtual void draw(StelCore* core, float maxMagLabels, const QFont& planetNameFont);
 
+	// re-implementation of Planet's update() to prepare tails (extinction etc). @param deltaTime: ms (since last call)
+	// TODO: computePosition can be removed and things added here!
+	virtual void update(int deltaTime);
+
 private:
 	//! @returns estimates for (Coma diameter [AU], gas tail length [AU]).
 	//! Using the formula from Guide found by the GSoC2012 initiative at http://www.projectpluto.com/update7b.htm#comet_tail_formula
@@ -117,7 +122,7 @@
 	//! @param colorArr vertex colors (if not textured) r0, g0, b0, r1, g1, b1, ...
 	//! @param indices into the former arrays (zero-starting), triplets forming triangles: t0,0, t0,1, t0,2, t1,0, t1,1, t1,2, ...
 	//! @param xOffset for the dust tail, this may introduce a bend. Units are x per sqrt(z).
-	void computeParabola(const float parameter, const float topradius, const float zshift, QVector<double>& vertexArr, QVector<float>& texCoordArr, QVector<unsigned short>& indices, const float xOffset=0.0f);
+	void computeParabola(const float parameter, const float topradius, const float zshift, QVector<Vec3d>& vertexArr, QVector<float>& texCoordArr, QVector<unsigned short>& indices, const float xOffset=0.0f);
 
 	double absoluteMagnitude;
 	double slopeParameter;
@@ -127,9 +132,10 @@
 
 	//GZ Tail additions
 	Vec2f tailFactors; // result of latest call to getComaDiameterAndTailLengthAU(); Results cached here for infostring. [0]=Coma diameter, [1] gas tail length.
-	bool tailActive;		//! true if there is a tail worth bothering (longer than COMET_MIN_TAIL_LENGTH_AU)? Drawing tails is quite costly.
+	bool tailActive;		//! true if there is a tail long enough to be worth drawing. Drawing tails is quite costly.
+	bool tailBright;		//! true if tail is bright enough to draw.
 	double deltaJDtail;             //! like deltaJD, but time difference between tail geometry updates.
-	double lastJDtail;             //! like lastJD, but time of last tail geometry update.
+	double lastJDtail;              //! like lastJD, but time of last tail geometry update.
 	Mat4d gasTailRot;		//! rotation matrix for gas tail parabola
 	Mat4d dustTailRot;		//! rotation matrix for the skewed dust tail parabola
 	float dustTailWidthFactor;      //!< empirical individual broadening of the dust tail end, compared to the gas tail end. Actually, dust tail width=2*comaWidth*dustTailWidthFactor. Default 1.5
@@ -142,8 +148,10 @@
 	static bool createTailIndices;
 	static bool createTailTextureCoords;
 
-	QVector<double> gastailVertexArr;  // computed frequently, describes parabolic shape (along z axis) of gas tail.
-	QVector<double> dusttailVertexArr; // computed frequently, describes parabolic shape (along z axis) of dust tail.
+	QVector<Vec3d> gastailVertexArr;  // computed frequently, describes parabolic shape (along z axis) of gas tail.
+	QVector<Vec3d> dusttailVertexArr; // computed frequently, describes parabolic shape (along z axis) of dust tail.
+	QVector<Vec3f> gastailColorArr;    // NEW computed for every 5 mins, modulates gas tail brightness for extinction
+	QVector<Vec3f> dusttailColorArr;   // NEW computed for every 5 mins, modulates dust tail brightness for extinction
 	static QVector<float> tailTexCoordArr; // computed only once for all comets!
 	static QVector<unsigned short> tailIndices; // computed only once for all comets!
 	static StelTextureSP comaTexture;

=== modified file 'src/core/modules/MilkyWay.cpp'
--- src/core/modules/MilkyWay.cpp	2014-10-27 22:40:49 +0000
+++ src/core/modules/MilkyWay.cpp	2014-11-08 17:16:48 +0000
@@ -128,7 +128,7 @@
 		for (int i=0; i<vertexArray->vertex.size(); ++i)
 		{
 			Vec3d vertAltAz=core->j2000ToAltAz(vertexArray->vertex.at(i), StelCore::RefractionOn);
-			Q_ASSERT(vertAltAz.lengthSquared()-1.0 < 0.001f);
+			Q_ASSERT(fabs(vertAltAz.lengthSquared()-1.0) < 0.001);
 
 			float oneMag=0.0f;
 			extinction.forward(vertAltAz, &oneMag);

=== modified file 'src/core/modules/MinorPlanet.cpp'
--- src/core/modules/MinorPlanet.cpp	2014-08-31 10:34:41 +0000
+++ src/core/modules/MinorPlanet.cpp	2014-11-08 17:16:48 +0000
@@ -293,7 +293,6 @@
 	//Calculate phase angle
 	//(Code copied from Planet::getVMagnitude())
 	//(this is actually vector subtraction + the cosine theorem :))
-	// GZ Try now in float... speed difference is negligible, though.
 	const Vec3d& observerHelioPos = core->getObserverHeliocentricEclipticPos();
 	const float observerRq = observerHelioPos.lengthSquared();
 	const Vec3d& planetHelioPos = getHeliocentricEclipticPos();

=== modified file 'src/core/modules/Planet.cpp'
--- src/core/modules/Planet.cpp	2014-11-05 14:49:44 +0000
+++ src/core/modules/Planet.cpp	2014-11-08 17:16:48 +0000
@@ -679,7 +679,7 @@
 double Planet::computeDistance(const Vec3d& obsHelioPos)
 {
 	distance = (obsHelioPos-getHeliocentricEclipticPos()).length();
-	// GZ: improve fps by juggling updates for asteroids. They must be fast if close to observer, but can be slow if further away.
+	// improve fps by juggling updates for asteroids. They must be fast if close to observer, but can be slow if further away.
 	if (pType == Planet::isAsteroid)
 			deltaJD=distance*StelCore::JD_SECOND;
 	return distance;
@@ -928,11 +928,11 @@
 {
 	if (hidden)
 		return;
-	// GZ: Try to improve speed for minor planets: test if visible at all.
+	// Try to improve speed for minor planets: test if visible at all.
 	// For a full catalog of NEAs (11000 objects), with this and resetting deltaJD according to distance, rendering time went 4.5fps->12fps.	
-	// AW: Apply this rule to asteroids only
-	// Note that taking away the asteroids at this stage breaks dim-asteroid occultation of stars!
-	if (((getVMagnitude(core)-1.0f) > core->getSkyDrawer()->getLimitMagnitude()) && pType==Planet::isAsteroid)
+	// TBD: Note that taking away the asteroids at this stage breaks dim-asteroid occultation of stars!
+	//      Maybe make another configurable flag for those interested?
+	if (((getVMagnitude(core)+1.0f) > core->getSkyDrawer()->getLimitMagnitude()) && pType==Planet::isAsteroid)
 	{
 		return;
 	}

=== modified file 'src/core/modules/Planet.hpp'
--- src/core/modules/Planet.hpp	2014-11-01 15:03:47 +0000
+++ src/core/modules/Planet.hpp	2014-11-08 17:16:48 +0000
@@ -192,7 +192,8 @@
 
 	// Compute the position in the parent Planet coordinate system
 	void computePositionWithoutOrbits(const double dateJD);
-	virtual void computePosition(const double dateJD);// GZ: gets overridden in Comet!
+	//virtual void computePosition(const double dateJD);// GZ: gets overridden in Comet!
+	void computePosition(const double dateJD);// GZ: gets overridden in Comet!
 
 	// Compute the transformation matrix from the local Planet coordinate to the parent Planet coordinate
 	void computeTransMatrix(double date);
@@ -237,8 +238,8 @@
 	static void setLabelColor(const Vec3f& lc) {labelColor = lc;}
 	static const Vec3f& getLabelColor(void) {return labelColor;}
 
-	// update displayed elements. @param deltaTime: ms (?)
-	void update(int deltaTime);
+	// update displayed elements. @param deltaTime: ms (since last call)
+	virtual void update(int deltaTime);
 
 	void setFlagHints(bool b){hintFader = b;}
 	bool getFlagHints(void) const {return hintFader;}