rayleigh scattering: use accurate model

data/shaders/opengl/basesphere_uniforms.glsl

@@ -30,6 +30,7 @@ layout(std140) uniform BaseSphereData {
 
 	// Eclipse data
 	Eclipse eclipse;
+	vec2 tropoHeight;			// height for (R, M) in km, pressure is 0.1 atm
 };
 
 // NOTE: you must include attributes.glsl first!

data/shaders/opengl/geosphere_terrain.frag

@@ -4,6 +4,7 @@
 #include "attributes.glsl"
 #include "lib.glsl"
 #include "basesphere_uniforms.glsl"
+#include "rayleigh.glsl"
 
 #define WATER_SHINE 16.0
 
@@ -40,6 +41,10 @@ void main(void)
 	vec3 tnorm = normalize(varyingNormal);
 	vec4 diff = vec4(0.0);
 
+	vec2 atmosDist  = raySphereIntersect(geosphereCenter, eyenorm, geosphereAtmosTopRad);
+	vec2 planetDist = raySphereIntersect(geosphereCenter, eyenorm, geosphereInvRadius);
+	vec3 planetIntersect = geosphereCenter + (eyenorm * planetDist.x);
+
 	float surfaceDist = dist * geosphereInvRadius;
 
 	// calculate the detail texture contribution from hi and lo textures
@@ -49,21 +54,30 @@ void main(void)
 	vec4 detailMul = mix(vec4(1.0), detailVal, detailMix);
 
 #ifdef TERRAIN_WITH_WATER
-	float specularReflection=0.0;
+	vec4 waterSpecular = vec4(0.0);
 #endif
 
 #if (NUM_LIGHTS > 0)
 	vec3 V = normalize(eyeposScaled - geosphereCenter);
 
+	vec3 I = normalize(eyeposScaled - planetIntersect);
+	vec3 center = I * geosphereRadius;
+	frag_color = vec4(0.f);
+
 	for (int i=0; i<NUM_LIGHTS; ++i) {
 		vec3 L = normalize(uLight[i].position.xyz);
+
 		float uneclipsed = clamp(calcUneclipsed(eclipse, NumShadows, V, L), 0.0, 1.0);
 		CalcPlanetDiffuse(diff, uLight[i].diffuse, L, tnorm, uneclipsed);
 
 #ifdef TERRAIN_WITH_WATER
 		//water only for specular
-	    if (vertexColor.b > 0.05 && vertexColor.r < 0.05) {
-			CalcPlanetSpec(specularReflection, uLight[i], L, tnorm, eyenorm, WATER_SHINE);
+		if (vertexColor.b > 0.05 && vertexColor.r < 0.05) {
+			// Convert from radius-relative to real coordinates
+			vec3 lightPosAU = uLight[i].position.xyz * INV_AU;
+			float intensity = 1.f / dot(lightPosAU, lightPosAU); // magic to avoid calculating length and then squaring it
+
+			waterSpecular.xyz += computeIncidentLight(reflect(L, I), eyenorm, center, atmosDist, toLinear(uLight[i].diffuse), uneclipsed) * intensity;
 		}
 #endif
 	}
@@ -94,11 +108,15 @@ void main(void)
 #endif
 		final;
 
+#ifdef TERRAIN_WITH_WATER
+	vec4 waterColor = waterSpecular * 20.f;
+#endif
+
 #ifdef ATMOSPHERE
 	frag_color +=
 		(1.0-fogFactor) * (diff*atmosColor) +
 #ifdef TERRAIN_WITH_WATER
-		  diff * specularReflection * sunset +
+		  toSRGB(1 - exp(-waterColor)) +
 #endif
 		  (0.02-clamp(fogFactor,0.0,0.01))*diff*ldprod*sunset +	      //increase fog scatter
 		  (pow((1.0-pow(fogFactor,0.75)),256.0)*0.4*diff*atmosColor)*sunset;  //distant fog.

data/shaders/opengl/rayleigh.glsl

@@ -1,3 +1,5 @@
+#define FAST 0
+
 float height(const in vec3 orig, const in vec3 center)
 {
 	vec3 r = orig - center;
@@ -13,9 +15,9 @@ void scatter(out vec2 density, const in vec3 orig, const in vec3 center)
 	density = -height / scaleHeight;
 
 	// earth atmospheric density: 1.225 kg/m^3, divided by 1e5
-	// 1/1.225e-5 = 81632.65306
-	float earthDensities = geosphereAtmosFogDensity * 81632.65306f;
-	density /= earthDensities;
+	// 1/1.225e-5 = 81632.65306, ln(81632.65306) ~= 11.31
+	float earthDensities = log(geosphereAtmosFogDensity) + 11.31f;
+	density += earthDensities;
 }
 
 // orig: ray origin
@@ -87,6 +89,8 @@ float predictDensityIn(const in float radius, const in float atmosphereHeight, c
 	}
 }
 
+const float INV_AU = 1.f / 149598000000.f;
+
 // predict "scattering density" along the ray
 // sample: starting point of the ray
 // dir:    direction of ray
@@ -147,6 +151,11 @@ void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const
 	      atmosphereHeight = atmosphereRadius - geosphereRadius;
 
 	float tCurrent = boundaries.x;
+	float maxSegment = atmosphereHeight / numSamples;
+
+	if (height(tCurrent * dir, center) > 0 && numSamples * maxSegment < (boundaries.y - boundaries.x)) {
+		numSamples = min(int(ceil((boundaries.y - boundaries.x) / maxSegment)), 256); // max 256 segments
+	}
 	float segmentLength = (boundaries.y - boundaries.x) / numSamples;
 	for (int i = 0; i < numSamples; ++i) {
 		vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir;
@@ -155,25 +164,34 @@ void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const
 		scatter(density, samplePosition, center);
 		opticalDepth += exp(density) * segmentLength;
 
-		// light optical depth
 		vec2 opticalDepthLight = vec2(0.f);
-		vec3 samplePositionLight = samplePosition;
-
 		vec3 sampleGeoCenter = center - samplePosition;
+#if FAST
+		// light optical depth
+		vec3 samplePositionLight = samplePosition;
 		opticalDepthLight.x = predictDensityInOut(samplePositionLight, sunDirection, sampleGeoCenter, geosphereRadius, atmosphereHeight, coefficientsR);
 		opticalDepthLight.y = predictDensityInOut(samplePositionLight, sunDirection, sampleGeoCenter, geosphereRadius, atmosphereHeight, coefficientsM);
-
+#else // FAST
+		int numSamplesLight = 8;
+		vec2 boundariesLight = raySphereIntersect(sampleGeoCenter, sunDirection, geosphereRadius * geosphereAtmosTopRad);
+		float segmentLengthLight = boundariesLight.y / numSamplesLight;
+		float tCurrentLight = 0.f;
+		for (int j = 0; j < numSamplesLight; ++j) {
+			vec3 samplePositionLight = vec3(segmentLengthLight * 0.5f + tCurrentLight) * sunDirection + samplePosition;
+			vec2 densityLDir = vec2(0.f);
+			scatter(densityLDir, samplePositionLight, center);
+			opticalDepthLight += exp(densityLDir) * segmentLengthLight;
+
+			tCurrentLight += segmentLengthLight;
+		}
+#endif // FAST
 		vec3 surfaceNorm = -normalize(sampleGeoCenter);
 		vec4 atmosDiffuse = vec4(0.f);
 		CalcPlanetDiffuse(atmosDiffuse, diffuse, sunDirection, surfaceNorm, uneclipsed);
 
-		vec3 tau = -(betaR * (opticalDepth.x + opticalDepthLight.x) + betaM * 1.1f * (opticalDepth.y + opticalDepthLight.y));
-		vec3 tauR = tau + vec3(density.x);
-		vec3 tauM = tau + vec3(density.y);
-		vec3 attenuationR = exp(tauR) * segmentLength;
-		vec3 attenuationM = exp(tauM) * segmentLength;
-		sumR += attenuationR * atmosDiffuse.xyz;
-		sumM += attenuationM * atmosDiffuse.xyz;
+		vec3 tau = -(betaR * opticalDepth.x + betaR * opticalDepthLight.x + betaM * 1.1f * opticalDepth.y + betaM * 1.1f * opticalDepthLight.y);
+		sumR += exp(tau + vec3(density.x)) * segmentLength * atmosDiffuse.xyz;
+		sumM += exp(tau + vec3(density.y)) * segmentLength * atmosDiffuse.xyz;
 		tCurrent += segmentLength;
 	}
 }

data/shaders/opengl/rayleigh_accurate.frag

@@ -30,8 +30,6 @@ void main(void)
 	vec3 a = atmosDist.x * eyenorm - geosphereCenter;
 	vec3 b = atmosDist.y * eyenorm - geosphereCenter;
 
-	float AU = 149598000000.0;
-
 #if (NUM_LIGHTS > 0)
 	for (int i=0; i<NUM_LIGHTS; ++i) {
 		vec3 lightDir = normalize(vec3(uLight[i].position));
@@ -41,7 +39,7 @@ void main(void)
 		// Convert from radius-relative to real coordinates
 		vec3 center = geosphereCenter * geosphereRadius;
 
-		vec3 lightPosAU = uLight[i].position.xyz / AU;
+		vec3 lightPosAU = uLight[i].position.xyz * INV_AU;
 		float intensity = 1.f / dot(lightPosAU, lightPosAU); // magic to avoid calculating length and then squaring it
 
 		specularHighlight += computeIncidentLight(lightDir, eyenorm, center, atmosDist, toLinear(uLight[i].diffuse), uneclipsed) * intensity;

data/shaders/opengl/rayleigh_fast.vert

@@ -26,8 +26,6 @@ void main(void)
 	vec3 a = atmosDist.x * eyenorm - geosphereCenter;
 	vec3 b = atmosDist.y * eyenorm - geosphereCenter;
 
-	float AU = 149598000000.0;
-
 #if (NUM_LIGHTS > 0)
 	for (int i=0; i<NUM_LIGHTS; ++i) {
 		vec3 lightDir = normalize(vec3(uLight[i].position));
@@ -37,7 +35,7 @@ void main(void)
 		// Convert from radius-relative to real coordinates
 		vec3 center = geosphereCenter * geosphereRadius;
 
-		vec3 lightPosAU = uLight[i].position.xyz / AU;
+		vec3 lightPosAU = uLight[i].position.xyz * INV_AU;
 		float intensity = 1.f / dot(lightPosAU, lightPosAU); // magic to avoid calculating length and then squaring it
 
 		specularHighlight += computeIncidentLight(lightDir, eyenorm, center, atmosDist, toLinear(uLight[i].diffuse), uneclipsed) * intensity;

src/BaseSphere.cpp

@@ -32,8 +32,10 @@ struct BaseSphereDataBlock {
 	alignas(16) vector3f srad;
 	alignas(16) vector3f lrad;
 	alignas(16) vector3f sdivlrad;
+
+	alignas(16) vector2f tropoHeight;
 };
-static_assert(sizeof(BaseSphereDataBlock) == 192, "");
+static_assert(sizeof(BaseSphereDataBlock) == 208, "");
 
 std::unique_ptr<Graphics::Drawables::Sphere3D> BaseSphere::m_atmos;
 
@@ -123,6 +125,8 @@ void BaseSphere::SetMaterialParameters(const matrix4x4d &trans, const float radi
 		++j;
 	}
 
+	matData.tropoHeight = ap.tropoHeight;
+
 	// FIXME: these two should share the same buffer data instead of making two separate allocs
 	m_surfaceMaterial->SetBufferDynamic(s_baseSphereData, &matData);
 	m_surfaceMaterial->SetPushConstant(s_numShadows, int(shadows.size()));

src/galaxy/AtmosphereParameters.h

@@ -14,6 +14,7 @@ struct AtmosphereParameters {
 	vector3f rayleighCoefficients;
 	vector3f mieCoefficients;
 	vector2f scaleHeight;
+	vector2f tropoHeight;
 };
 
 #endif // ATMOSPHEREPARAMETERS_H_INCLUDED

src/galaxy/SystemBody.cpp

@@ -201,8 +201,17 @@ double SystemBody::GetAtmPressure(double altitude) const
 	const double lapseRate_L = surfaceGravity_g / specificHeat; // deg/m
 	const double surfaceTemperature_T0 = GetAverageTemp();	//K
 
-	return m_atmosPressure * pow((1 - lapseRate_L * altitude / surfaceTemperature_T0),
-		(specificHeat * gasMolarMass / GAS_CONSTANT_R)); // in ATM since p0 was in ATM
+	// pressure below tropopause
+	const double atmosPressure = m_atmosPressure * pow((1 - lapseRate_L * altitude / surfaceTemperature_T0),
+					(specificHeat * gasMolarMass / GAS_CONSTANT_R)); // in ATM since p0 was in ATM
+
+	if (atmosPressure > 0.1) {
+		return atmosPressure;
+	} else {
+		// above tropopause
+		const double tropopauseTemp = GetAtmAverageTemp(m_tropopause);
+		return 0.1 * exp((-surfaceGravity_g * gasMolarMass * (altitude - m_tropopause)) / (GAS_CONSTANT_R * tropopauseTemp));
+	}
 }
 
 double SystemBody::GetAtmDensity(double altitude, double pressure) const
@@ -217,6 +226,10 @@ double SystemBody::GetAtmAverageTemp(double altitude) const
 {
 	// temperature at height
 	const double lapseRate_L = CalcSurfaceGravity() / GetSpecificHeat(GetSuperType()); // deg/m
+	if (altitude > m_tropopause) {
+		return double(GetAverageTemp()) - lapseRate_L * m_tropopause;
+	}
+
 	return double(GetAverageTemp()) - lapseRate_L * altitude;
 }
 
@@ -241,7 +254,18 @@ void SystemBody::SetAtmFromParameters()
 		// want height for pressure 0.001 atm:
 		// h = (1 - exp(RL/gM * log(P/p0))) * T0 / l
 		double RLdivgM = (GAS_CONSTANT_R * lapseRate_L) / (surfaceGravity_g * GetMolarMass(GetSuperType()));
-		m_atmosRadius = (1.0 - exp(RLdivgM * log(0.001 / m_atmosPressure))) * surfaceTemperature_T0 / lapseRate_L;
+//		m_atmosRadius = (1.0 - exp(RLdivgM * log(0.001 / m_atmosPressure))) * surfaceTemperature_T0 / lapseRate_L;
+
+		// get tropopause: height for pressure 0.1 atm
+		m_tropopause = (1.0 - exp(RLdivgM * log(0.1 / m_atmosPressure))) * surfaceTemperature_T0 / lapseRate_L;
+
+		// if tropopause is below surface (surface pressure < 0.1 atm)
+		if (m_tropopause < 0.0) {
+			m_tropopause = 0.0;
+		}
+
+		double tropopause_temperature = surfaceTemperature_T0 - lapseRate_L * m_tropopause;
+		m_atmosRadius = log(0.001 / 0.1) * GAS_CONSTANT_R * tropopause_temperature / (-surfaceGravity_g * GetMolarMass(GetSuperType())) + m_tropopause;
 	}
 }
 
@@ -347,7 +371,7 @@ AtmosphereParameters SystemBody::CalcAtmosphereParams() const
 	const double massPlanet_in_kg = (m_mass.ToDouble() * EARTH_MASS);
 	const double g = G * massPlanet_in_kg / (radiusPlanet_in_m * radiusPlanet_in_m);
 
-	double T = static_cast<double>(m_averageTemp);
+	double T = static_cast<double>(GetAtmAverageTemp(m_tropopause));
 
 	// XXX hack to avoid issues with sysgen giving 0 temps
 	// temporary as part of sysgen needs to be rewritten before the proper fix can be used
@@ -361,17 +385,17 @@ AtmosphereParameters SystemBody::CalcAtmosphereParams() const
 
 	// min of 2.0 corresponds to a scale height of 1/20 of the planet's radius,
 	params.atmosInvScaleHeight = std::max(20.0f, static_cast<float>(GetRadius() / atmosScaleHeight));
-	// integrate atmospheric density between surface and this radius. this is 10x the scale
-	// height, which should be a height at which the atmospheric density is negligible
-	params.atmosRadius = 1.0f + static_cast<float>(10.0f * atmosScaleHeight) / GetRadius();
+	params.atmosRadius = 1.0f + static_cast<float>(m_atmosRadius) / radiusPlanet_in_m;
 
 	params.planetRadius = static_cast<float>(radiusPlanet_in_m);
 
 	const float radiusPlanet_in_km = radiusPlanet_in_m / 1000;
 	const float atmosHeight_in_km = radiusPlanet_in_km * (params.atmosRadius - 1);
+	const float tropopause_in_km = static_cast<float>(m_tropopause / 1000);
 	params.rayleighCoefficients = GetCoefficients(radiusPlanet_in_km, atmosHeight_in_km, atmosScaleHeight);
 	params.mieCoefficients = GetCoefficients(radiusPlanet_in_km, atmosHeight_in_km, atmosScaleHeight / 6.66); // 7994 / 1200 = 6.61
 	params.scaleHeight = vector2f(atmosScaleHeight, atmosScaleHeight / 6.66);
+	params.tropoHeight = vector2f(tropopause_in_km, tropopause_in_km / 6.66);
 
 	return params;
 }

src/galaxy/SystemBody.h

@@ -282,6 +282,7 @@ class SystemBody : public RefCounted, public SystemBodyType, protected SystemBod
 	double GetAtmSurfaceDensity() const { return m_volatileGas.ToDouble(); }
 	double GetAtmSurfacePressure() const { return m_atmosPressure; }
 	double GetAtmRadius() const { return m_atmosRadius; }
+	double GetTropopause() const { return m_tropopause; }
 
 	// Calculate atmosphere pressure at given altitude (atm)
 	double GetAtmPressure(double altitude) const;
@@ -338,8 +339,10 @@ class SystemBody : public RefCounted, public SystemBodyType, protected SystemBod
 
 	// atmosphere surface pressure, unit: atm
 	double m_atmosPressure;
-	// atmosphere radius at 0.01atm, unit: meters
+	// atmosphere radius at 0.001atm, unit: meters
 	double m_atmosRadius;
+	// tropopause height at 0.1atm, unit: meters
+	double m_tropopause;
 };
 
 #endif // SYSTEMBODY_H