diff --git a/data/lang/ui-core/en.json b/data/lang/ui-core/en.json
index 33c7cd2a172..9eab63e1fd3 100644
--- a/data/lang/ui-core/en.json
+++ b/data/lang/ui-core/en.json
@@ -515,6 +515,10 @@
     "description": "For player reputation",
     "message": "Experienced"
   },
+  "EXPERIMENTAL": {
+    "description": "Indicates feature/option is experimental",
+    "message": "Experimental"
+  },
   "EXPERT": {
     "description": "For player reputation",
     "message": "Expert"
@@ -1835,6 +1839,26 @@
     "description": "",
     "message": "Rating:"
   },
+  "REALISTIC_SCATTERING": {
+    "description": "",
+    "message": "Scattering"
+  },
+  "REALISTIC_SCATTERING_DESC": {
+    "description": "",
+    "message": "Select scattering when rendering atmospheres. (Rayleigh/Mie is not recommended for low-end PCs)"
+  },
+  "SCATTERING_OLD": {
+    "description": "",
+    "message": "Legacy"
+  },
+  "RAYLEIGH_FAST": {
+    "description": "",
+    "message": "Rayleigh/Mie (fast, per-vertex)"
+  },
+  "RAYLEIGH_ACCURATE": {
+    "description": "",
+    "message": "Rayleigh/Mie (accurate, per-pixel)"
+  },
   "REAR_WEAPON": {
     "description": "",
     "message": "Rear weapon"
diff --git a/data/pigui/modules/settings-window.lua b/data/pigui/modules/settings-window.lua
index 069a1f99a25..880dc64cacd 100644
--- a/data/pigui/modules/settings-window.lua
+++ b/data/pigui/modules/settings-window.lua
@@ -165,6 +165,16 @@ local function showVideoOptions()
 	local gpuJobs = Engine.GetGpuJobsEnabled()
 	local disableScreenshotInfo = Engine.GetDisableScreenshotInfo()
 
+	-- Scattering is still an experimental feature
+	local experimental = "[" .. lui.EXPERIMENTAL .. "] "
+	local scatteringLabels = {
+		lui.SCATTERING_OLD,
+		experimental .. lui.RAYLEIGH_FAST,
+		experimental .. lui.RAYLEIGH_ACCURATE
+	}
+
+	local realisticScattering = Engine.GetRealisticScattering()
+
 	local cityDetail = keyOf(detailLabels,keyOf(detailLevels, Engine.GetCityDetailLevel()))-1
 	local displayNavTunnels = Engine.GetDisplayNavTunnels()
 	local displaySpeedLines = Engine.GetDisplaySpeedLines()
@@ -187,6 +197,11 @@ local function showVideoOptions()
 		Engine.SetMultisampling(aa)
 	end
 
+	c,scattering = combo(lui.REALISTIC_SCATTERING, realisticScattering, scatteringLabels, lui.REALISTIC_SCATTERING_DESC)
+	if c then
+		Engine.SetRealisticScattering(scattering)
+	end
+
 	ui.columns(2,"video_checkboxes",false)
 	c,fullscreen = checkbox(lui.FULL_SCREEN, fullscreen)
 	if c then
diff --git a/data/shaders/opengl/basesphere_uniforms.glsl b/data/shaders/opengl/basesphere_uniforms.glsl
index 18bc8fa430c..f568facca34 100644
--- a/data/shaders/opengl/basesphere_uniforms.glsl
+++ b/data/shaders/opengl/basesphere_uniforms.glsl
@@ -4,16 +4,79 @@
 #include "eclipse.glsl"
 
 layout(std140) uniform BaseSphereData {
-	// to keep distances sane we do a nearer, smaller scam. this is how many times
-	// smaller the geosphere has been made
-	vec3 geosphereCenter;				// TODO documentation
-	float geosphereRadius;				// (planet radius)
-	float geosphereInvRadius;			// 1.0 / (planet radius)
-	float geosphereAtmosTopRad;			// in planet radii
-	float geosphereAtmosFogDensity;		// TODO documentation
-	float geosphereAtmosInvScaleHeight; // TODO documentation
+	// To render accurate planets at any distance or scale, almost all
+	// calculations take place in "planet space", with a distance of 1.0
+	// defined as equivalent to the planet's nominal radius (planet radii).
+	//
+	// This allows planets to be rendered in view-space at any scaling factor
+	// desired, to play nicely with the depth buffer and any other content
+	// rendered (e.g. map views)
+	//
+	// The atmosphere simulation approximates Rayleigh scattering by
+	// calculating the optical density of a path through the atmosphere and
+	// assuming a constant in/out scattering factor based on the density
+	// approximation.
+
+	vec3 geosphereCenter;				// view-space center of the planet, in planet radii
+	float geosphereRadius;				// real planet radius, in meters
+	float geosphereInvRadius;			// 1.0 / (view-space radius), converts between view coordinates and planet radii
+	float geosphereAtmosTopRad;			// height of the simulated atmosphere, in planet radii
+	float geosphereAtmosFogDensity;		// atmospheric density scalar
+	float geosphereAtmosInvScaleHeight; // 1.0 / (atmosphere scale height) in planet radii
 	vec4 atmosColor;
+	vec3 coefficientsR;			// coefficients for approximating the Rayleigh contribution
+	vec3 coefficientsM;			// coefficients for approximating the Mie contribution
+	vec2 scaleHeight;			// height for (R, M) in km, at which density will be reduced by e
 
 	// Eclipse data
 	Eclipse eclipse;
 };
+
+// NOTE: you must include attributes.glsl first!
+
+// Common code to calculate the diffuse light term for a planet's surface
+// L: light -> surface vector (normalized)
+// N: surface normal
+// V: surface position relative to the unit-sphere planet
+void CalcPlanetDiffuse(inout vec4 diff, in vec4 color, in vec3 L, in vec3 N, in float uneclipsed)
+{
+	float nDotVP  = max(0.0, dot(N, L));
+	float nnDotVP = max(0.0, dot(N, -L));
+
+	//need backlight to increase horizon, attempts to model light propagating towards terminator
+	float clampedCosine = (nDotVP + 0.5 * clamp(1.0 - nnDotVP * 4.0, 0, 1) * INV_NUM_LIGHTS);
+	diff += color * uneclipsed * 0.5 * clampedCosine;
+}
+
+#ifdef FRAGMENT_SHADER
+
+// Common code to calculate the specular light term for a planet's surface
+// L: light -> surface vector (normalized)
+// N: surface normal
+// E: eye->surface vector (normalized)
+// power: specular power (blinn-phong)
+void CalcPlanetSpec(inout float spec, in Light light, in vec3 L, in vec3 N, in vec3 E, in float power)
+{
+	//Specular reflection
+	vec3 H = normalize(L - E);
+	//water only for specular
+	spec += pow(max(dot(H, N), 0.0), power) * 0.5 * INV_NUM_LIGHTS;
+}
+
+// E: eye->surface direction
+// scaledPos: position of the pixel in view space, divided by the radius of the geosphere
+void CalcPlanetFogFactor(inout float ldprod, inout float fogFactor, in vec3 E, in vec3 surfacePos, in float dist)
+{
+	// when does the eye ray intersect atmosphere
+	float atmosStart = raySphereIntersect(geosphereCenter, E, geosphereAtmosTopRad).x;
+	float atmosDist = (dist - atmosStart) * geosphereRadius;
+
+	// a&b scaled so length of 1.0 means planet surface.
+	vec3 a = atmosStart * E - geosphereCenter;
+	vec3 b = surfacePos;
+
+	ldprod = AtmosLengthDensityProduct(a, b, atmosColor.w*geosphereAtmosFogDensity, atmosDist, geosphereAtmosInvScaleHeight);
+	fogFactor = clamp(1.5 / exp(ldprod), 0.0, 1.0);
+}
+
+#endif
diff --git a/data/shaders/opengl/gassphere_base.frag b/data/shaders/opengl/gassphere_base.frag
index 58e82f71c78..b922b9b0c31 100644
--- a/data/shaders/opengl/gassphere_base.frag
+++ b/data/shaders/opengl/gassphere_base.frag
@@ -17,35 +17,22 @@ out vec4 frag_color;
 
 void main(void)
 {
-	vec3 eyepos = varyingEyepos;
-	vec3 eyenorm = normalize(eyepos);
+	vec3 eyeposScaled = varyingEyepos * geosphereInvRadius;
+	vec3 eyenorm = normalize(varyingEyepos);
 	vec3 tnorm = normalize(varyingNormal);
 	vec4 diff = vec4(0.0);
-	float nDotVP=0.0;
-	float nnDotVP=0.0;
+	float surfaceDist = length(eyeposScaled);
 
-	vec3 v = (eyepos - geosphereCenter) * geosphereInvRadius;
-	float lenInvSq = 1.0/(dot(v,v));
+	vec3 V = (eyeposScaled - geosphereCenter);
 	for (int i=0; i<NUM_LIGHTS; ++i) {
-		float uneclipsed = clamp(calcUneclipsed(eclipse, NumShadows, v, normalize(vec3(uLight[i].position))), 0.0, 1.0);
-		nDotVP  = max(0.0, dot(tnorm, normalize(vec3(uLight[i].position))));
-		nnDotVP = max(0.0, dot(tnorm, normalize(-vec3(uLight[i].position)))); //need backlight to increase horizon
-		diff += uLight[i].diffuse * uneclipsed * 0.5*(nDotVP+0.5*clamp(1.0-nnDotVP*4.0,0.0,1.0) * INV_NUM_LIGHTS);
+		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);
 	}
 
-	// when does the eye ray intersect atmosphere
-	float atmosStart = findSphereEyeRayEntryDistance(geosphereCenter, eyepos, geosphereRadius * geosphereAtmosTopRad);
 	float ldprod=0.0;
 	float fogFactor=0.0;
-	{
-		float atmosDist = (length(eyepos) - atmosStart);
-
-		// a&b scaled so length of 1.0 means planet surface.
-		vec3 a = (atmosStart * eyenorm - geosphereCenter) * geosphereInvRadius;
-		vec3 b = (eyepos - geosphereCenter) * geosphereInvRadius;
-		ldprod = AtmosLengthDensityProduct(a, b, atmosColor.w*geosphereAtmosFogDensity, atmosDist, geosphereAtmosInvScaleHeight);
-		fogFactor = clamp( 1.5 / exp(ldprod),0.0,1.0);
-	}
+	CalcPlanetFogFactor(ldprod, fogFactor, eyenorm, eyeposScaled - geosphereCenter, surfaceDist);
 
 	//calculate sunset tone red when passing through more atmosphere, clamp everything.
 	float atmpower = (diff.r+diff.g+diff.b)/3.0;
diff --git a/data/shaders/opengl/geosphere_sky.frag b/data/shaders/opengl/geosphere_sky.frag
index 7f40ff3efc3..6481ab0125d 100644
--- a/data/shaders/opengl/geosphere_sky.frag
+++ b/data/shaders/opengl/geosphere_sky.frag
@@ -11,65 +11,44 @@ in vec4 varyingEyepos;
 
 out vec4 frag_color;
 
-void sphereEntryExitDist(out float near, out float far, const in vec3 sphereCenter, const in vec3 eyeTo, const in float radius)
-{
-	vec3 v = -sphereCenter;
-	vec3 dir = normalize(eyeTo);
-	float b = -dot(v, dir);
-	float det = (b * b) - dot(v, v) + (radius * radius);
-	float i1, i2;
-	near = 0.0;
-	far = 0.0;
-	if (det > 0.0) {
-		det = sqrt(det);
-		i1 = b - det;
-		i2 = b + det;
-		if (i2 > 0.0) {
-			near = max(i1, 0.0);
-			far = i2;
-		}
-	}
-}
-
 void main(void)
 {
-	float skyNear, skyFar;
 	vec3 eyenorm = normalize(varyingEyepos.xyz);
 	float specularHighlight=0.0;
 
-	sphereEntryExitDist(skyNear, skyFar, geosphereCenter, varyingEyepos.xyz, geosphereRadius * geosphereAtmosTopRad);
-	float atmosDist = (skyFar - skyNear);
+	vec2 viewDist = raySphereIntersect(geosphereCenter, eyenorm, geosphereAtmosTopRad);
+	vec2 isect = raySphereIntersect(geosphereCenter, eyenorm, 1.0);
+
+	float atmosDist = (viewDist.y - viewDist.x) * geosphereRadius;
 	float ldprod=0.0;
 
 	// a&b scaled so length of 1.0 means planet surface.
-	vec3 a = (skyNear * eyenorm - geosphereCenter) * geosphereInvRadius;
-	vec3 b = (skyFar * eyenorm - geosphereCenter) * geosphereInvRadius;
+	vec3 a = viewDist.x * eyenorm - geosphereCenter;
+	vec3 b = viewDist.y * eyenorm - geosphereCenter;
 	ldprod = AtmosLengthDensityProduct(a, b, atmosColor.a * geosphereAtmosFogDensity, atmosDist, geosphereAtmosInvScaleHeight);
 
 	float fogFactor = 1.0 / exp(ldprod);
 	vec4 atmosDiffuse = vec4(0.0);
 
 #if (NUM_LIGHTS > 0)
-	vec3 surfaceNorm = normalize(skyNear * eyenorm - geosphereCenter);
+	vec3 surfaceNorm = normalize(viewDist.x * eyenorm - geosphereCenter);
 	for (int i=0; i<NUM_LIGHTS; ++i) {
 
-		vec3 lightDir = normalize(vec3(uLight[i].position));
-
-		float uneclipsed = clamp(calcUneclipsedSky(eclipse, NumShadows, a, b, lightDir), 0.0, 1.0);
+		vec3 L = normalize(vec3(uLight[i].position));
+		float uneclipsed = clamp(calcUneclipsedSky(eclipse, NumShadows, a, b, L), 0.0, 1.0);
 
-		float nDotVP =  max(0.0, dot(surfaceNorm, lightDir));
-		float nnDotVP = max(0.0, dot(surfaceNorm, -lightDir));  //need backlight to increase horizon
-		atmosDiffuse +=  uLight[i].diffuse * uneclipsed * 0.5*(nDotVP+0.5*clamp(1.0-nnDotVP*4.0,0.0,1.0) * INV_NUM_LIGHTS);
+		CalcPlanetDiffuse(atmosDiffuse, toLinear(uLight[i].diffuse), L, surfaceNorm, uneclipsed);
 
-		//Calculate Specular Highlight
-		vec3 L = normalize(uLight[i].position.xyz - varyingEyepos.xyz);
-		vec3 E = normalize(-varyingEyepos.xyz);
-		vec3 R = normalize(-reflect(L,vec3(0.0)));
-		specularHighlight += pow(max(dot(R,E),0.0),64.0) * uneclipsed * INV_NUM_LIGHTS;
+		// Calculate Specular Highlight (halo around the light source)
+		specularHighlight += pow(max(dot(L, eyenorm), 0.0), 64.0) * uneclipsed * INV_NUM_LIGHTS;
 
 	}
 #endif
 
+	// Tonemap in sRGB space to match existing visuals
+	atmosDiffuse = toSRGB(atmosDiffuse);
+	atmosDiffuse = 1.0 - exp(-atmosDiffuse);
+
 	//calculate sunset tone red when passing through more atmosphere, clamp everything.
 	float atmpower = (atmosDiffuse.r+atmosDiffuse.g+atmosDiffuse.b)/3.0;
 	vec4 sunset = vec4(0.8,clamp(pow(atmpower,0.8),0.0,1.0),clamp(pow(atmpower,1.2),0.0,1.0),1.0);
diff --git a/data/shaders/opengl/geosphere_terrain.frag b/data/shaders/opengl/geosphere_terrain.frag
index e3e3c029d02..55870d72e58 100644
--- a/data/shaders/opengl/geosphere_terrain.frag
+++ b/data/shaders/opengl/geosphere_terrain.frag
@@ -5,6 +5,8 @@
 #include "lib.glsl"
 #include "basesphere_uniforms.glsl"
 
+#define WATER_SHINE 16.0
+
 uniform sampler2D texture0;
 uniform sampler2D texture1;
 in vec2 texCoord0;
@@ -33,11 +35,12 @@ void main(void)
 	vec4 hidetail = texture(texture0, texCoord0 * detailScaleHi);
 	vec4 lodetail = texture(texture1, texCoord0 * detailScaleLo);
 
-	vec3 eyepos = varyingEyepos;
-	vec3 eyenorm = normalize(eyepos);
+	vec3 eyeposScaled = varyingEyepos * geosphereInvRadius;
+	vec3 eyenorm = normalize(varyingEyepos);
 	vec3 tnorm = normalize(varyingNormal);
 	vec4 diff = vec4(0.0);
 
+	float surfaceDist = dist * geosphereInvRadius;
 
 	// calculate the detail texture contribution from hi and lo textures
 	float hiloMix = exp(-0.004 * dist);
@@ -45,77 +48,60 @@ void main(void)
 	vec4 detailVal = mix(lodetail, hidetail, hiloMix);
 	vec4 detailMul = mix(vec4(1.0), detailVal, detailMix);
 
-	float nDotVP=0.0;
-	float nnDotVP=0.0;
 #ifdef TERRAIN_WITH_WATER
 	float specularReflection=0.0;
 #endif
 
 #if (NUM_LIGHTS > 0)
-	vec3 v = (eyepos - geosphereCenter) * geosphereInvRadius;
-	float lenInvSq = 1.0/(dot(v,v));
+	vec3 V = normalize(eyeposScaled - geosphereCenter);
+
 	for (int i=0; i<NUM_LIGHTS; ++i) {
-		float uneclipsed = clamp(calcUneclipsed(eclipse, NumShadows, v, normalize(vec3(uLight[i].position))), 0.0, 1.0);
-		nDotVP  = max(0.0, dot(tnorm, normalize(vec3(uLight[i].position))));
-		nnDotVP = max(0.0, dot(tnorm, normalize(-vec3(uLight[i].position)))); //need backlight to increase horizon
-		diff += uLight[i].diffuse * uneclipsed * 0.5*(nDotVP+0.5*clamp(1.0-nnDotVP*4.0,0.0,1.0) * INV_NUM_LIGHTS);
+		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
-		//Specular reflection
-		vec3 L = normalize(uLight[i].position.xyz - eyepos);
-		vec3 E = normalize(-eyepos);
-		vec3 R = normalize(-reflect(L,tnorm));
 		//water only for specular
 	    if (vertexColor.b > 0.05 && vertexColor.r < 0.05) {
-			specularReflection += pow(max(dot(R,E),0.0),16.0)*0.4 * INV_NUM_LIGHTS;
+			CalcPlanetSpec(specularReflection, uLight[i], L, tnorm, eyenorm, WATER_SHINE);
 		}
 #endif
 	}
 
 	// Use the detail value to multiply the final colour before lighting
-	vec4 final = vertexColor * detailMul;
+	vec4 ambient = scene.ambient * vertexColor;
+	vec4 final = vertexColor * detailMul * diff;
 
 #ifdef ATMOSPHERE
-	// when does the eye ray intersect atmosphere
-	float atmosStart = findSphereEyeRayEntryDistance(geosphereCenter, eyepos, geosphereRadius * geosphereAtmosTopRad);
 	float ldprod=0.0;
 	float fogFactor=0.0;
-	{
-		float atmosDist = (length(eyepos) - atmosStart);
-
-		// a&b scaled so length of 1.0 means planet surface.
-		vec3 a = (atmosStart * eyenorm - geosphereCenter) * geosphereInvRadius;
-		vec3 b = (eyepos - geosphereCenter) * geosphereInvRadius;
-		ldprod = AtmosLengthDensityProduct(a, b, atmosColor.w*geosphereAtmosFogDensity, atmosDist, geosphereAtmosInvScaleHeight);
-		fogFactor = clamp( 1.5 / exp(ldprod),0.0,1.0);
-	}
+	CalcPlanetFogFactor(ldprod, fogFactor, eyenorm, eyeposScaled - geosphereCenter, surfaceDist);
 
 	//calculate sunset tone red when passing through more atmosphere, clamp everything.
 	float atmpower = (diff.r+diff.g+diff.b)/3.0;
 	vec4 sunset = vec4(0.8,clamp(pow(atmpower,0.8),0.0,1.0),clamp(pow(atmpower,1.2),0.0,1.0),1.0);
 
+	final = fogFactor * (ambient + final);
+#else
+	// add extra brightness to atmosphere-less planetoids and dim stars
+	final = ambient + final * 2.0;
+#endif
+
 	frag_color =
 		material.emission +
 #ifdef TERRAIN_WITH_LAVA
 		varyingEmission +
 #endif
-		fogFactor *
-		((scene.ambient * vertexColor) +
-		(diff * final)) +
-		(1.0-fogFactor)*(diff*atmosColor) +
+		final;
+
+#ifdef ATMOSPHERE
+	frag_color +=
+		(1.0-fogFactor) * (diff*atmosColor) +
 #ifdef TERRAIN_WITH_WATER
-		  diff*specularReflection*sunset +
+		  diff * specularReflection * sunset +
 #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.
-#else // atmosphere-less planetoids and dim stars
-	frag_color =
-		material.emission +
-#ifdef TERRAIN_WITH_LAVA
-		varyingEmission +
-#endif
-		(scene.ambient * vertexColor) +
-		(diff * final * 2.0);
 #endif //ATMOSPHERE
 
 #else // NUM_LIGHTS > 0 -- unlit rendering - stars
diff --git a/data/shaders/opengl/geosphere_terrain.vert b/data/shaders/opengl/geosphere_terrain.vert
index 8800e224f5d..8c08e5dce0c 100644
--- a/data/shaders/opengl/geosphere_terrain.vert
+++ b/data/shaders/opengl/geosphere_terrain.vert
@@ -24,7 +24,7 @@ void main(void)
 	varyingNormal = normalize(normalMatrix() * a_normal);
 
 	texCoord0 = a_uv0.xy;
-	dist = abs(varyingEyepos.z);
+	dist = length(varyingEyepos);
 
 #ifdef TERRAIN_WITH_LAVA
 	varyingEmission = material.emission;
diff --git a/data/shaders/opengl/lib.glsl b/data/shaders/opengl/lib.glsl
index 79d2b53c42b..2b4ca6578dd 100644
--- a/data/shaders/opengl/lib.glsl
+++ b/data/shaders/opengl/lib.glsl
@@ -1,16 +1,41 @@
 // Copyright © 2008-2024 Pioneer Developers. See AUTHORS.txt for details
 // Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
 
-#ifdef FRAGMENT_SHADER
+#include "sRGB.glsl"
 
-struct Surface {
-	vec4 color;
-	vec3 specular;
-	float shininess;
-	vec3 normal;
-	vec3 emissive;
-	float ambientOcclusion;
-};
+// Simple ray-sphere intersection test, assuming ray starts at origin and rayDir is pre-normalized.
+// Returns distance to first and second intersections in {x, y} or 0.0 if no intersection.
+vec2 raySphereIntersect(in vec3 sphereCenter, in vec3 rayDir, in float radius)
+{
+	vec3 v = -sphereCenter;
+	float b = -dot(v, rayDir);
+	float det = (b * b) - dot(v, v) + (radius * radius);
+	float sdet = sqrt(det);
+
+	return det > 0.0 ? max(vec2(b - sdet, b + sdet), vec2(0.0)) : vec2(0.0);
+}
+
+// ray starts at origin, rayDir and axis are pre-normalized
+// Returns distance to first and second intersections in {x, y} or 0.0 if no intersection.
+vec2 rayCylinderIntersect(in vec3 rayDir, in vec3 cylinderCenter, in vec3 axis, in float radius)
+{
+    // tangent vectors (parallel to axis)
+    vec3 tray = axis * dot(rayDir, axis);
+    vec3 tcenter = axis * dot(cylinderCenter, axis);
+
+    // normal vectors (perpendicular to axis)
+    vec3 nray = rayDir - tray;
+    vec3 ncenter = cylinderCenter - tcenter;
+
+    // coefficient to move from projection to actual 3d space
+    // e.g. if angle between axis and tray = 30deg, actual intersect should be doubled
+    float scale = length(nray);
+
+    // intersection given in main plane projection
+    vec2 intersect = raySphereIntersect(ncenter, normalize(nray), radius);
+
+    return (scale == 0.f) ? vec2(0.f) : intersect / scale;
+}
 
 // Phase functions
 // https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-sky/simulating-colors-of-the-sky.html
@@ -30,6 +55,17 @@ float rayleighPhaseFunction(const float mu)
 	return 3.f / (16.f * 3.141592) * (1 + mu * mu);
 }
 
+#ifdef FRAGMENT_SHADER
+
+struct Surface {
+	vec4 color;
+	vec3 specular;
+	float shininess;
+	vec3 normal;
+	vec3 emissive;
+	float ambientOcclusion;
+};
+
 
 // Currently used by: hopefully everything
 // Evaluates a standard blinn-phong diffuse+specular, with the addition of a
@@ -45,25 +81,6 @@ void BlinnPhongDirectionalLight(in Light light, in float intensity, in Surface s
 	specular += surf.specular * light.specular.xyz * intensity * pow(max(dot(H, surf.normal), 0.0), surf.shininess);
 }
 
-//Currently used by: planet ring shader, geosphere shaders
-float findSphereEyeRayEntryDistance(in vec3 sphereCenter, in vec3 eyeTo, in float radius)
-{
-	vec3 v = -sphereCenter;
-	vec3 dir = normalize(eyeTo);
-	float b = -dot(v, dir);
-	float det = (b * b) - dot(v, v) + (radius * radius);
-	float entryDist = 0.0;
-	if (det > 0.0) {
-		det = sqrt(det);
-		float i1 = b - det;
-		float i2 = b + det;
-		if (i2 > 0.0) {
-			entryDist = max(i1, 0.0);
-		}
-	}
-	return entryDist;
-}
-
 // Used by: geosphere shaders
 // Calculate length*density product of a line through the atmosphere
 // a - start coord (normalized relative to atmosphere radius)
@@ -74,15 +91,21 @@ float AtmosLengthDensityProduct(vec3 a, vec3 b, float surfaceDensity, float len,
 {
 	/* 4 samples */
 	float ldprod = 0.0;
+
 	vec3 dir = b-a;
+	float ln = max(length(b)-1.0, 0.0);
+
+	/* simple 6-tap raymarch through the atmosphere to sample an average density */
 	ldprod = surfaceDensity * (
-			exp(-invScaleHeight*(length(a)-1.0)) +
-			exp(-invScaleHeight*(length(a + 0.2*dir)-1.0)) +
-			exp(-invScaleHeight*(length(a + 0.4*dir)-1.0)) +
-			exp(-invScaleHeight*(length(a + 0.6*dir)-1.0)) +
-			exp(-invScaleHeight*(length(a + 0.8*dir)-1.0)) +
-			exp(-invScaleHeight*max(length(b)-1.0, 0.0)));
+			exp(-invScaleHeight * (length(a)           -1.0)) +
+			exp(-invScaleHeight * (length(a + 0.2*dir) -1.0)) +
+			exp(-invScaleHeight * (length(a + 0.4*dir) -1.0)) +
+			exp(-invScaleHeight * (length(a + 0.6*dir) -1.0)) +
+			exp(-invScaleHeight * (length(a + 0.8*dir) -1.0)) +
+			exp(-invScaleHeight * (ln)));
+
 	ldprod *= len;
 	return ldprod;
 }
+
 #endif
diff --git a/data/shaders/opengl/planetrings.frag b/data/shaders/opengl/planetrings.frag
index ce235129055..a088fb400d4 100644
--- a/data/shaders/opengl/planetrings.frag
+++ b/data/shaders/opengl/planetrings.frag
@@ -11,6 +11,11 @@ in vec4 varyingEyepos;
 
 out vec4 frag_color;
 
+float findSphereEyeRayEntryDistance(in vec3 sphereCenter, in vec3 eyeTo, in float radius)
+{
+	return raySphereIntersect(sphereCenter, normalize(eyeTo), radius).x;
+}
+
 void main(void)
 {
 	// Bits of ring in shadow!
diff --git a/data/shaders/opengl/rayleigh.glsl b/data/shaders/opengl/rayleigh.glsl
new file mode 100644
index 00000000000..e1fa6b62d82
--- /dev/null
+++ b/data/shaders/opengl/rayleigh.glsl
@@ -0,0 +1,260 @@
+float height(const in vec3 orig, const in vec3 center)
+{
+	vec3 r = orig - center;
+	float height = sqrt(dot(r, r)) - geosphereRadius;
+
+	return height;
+}
+
+void scatter(out vec2 density, const in vec3 orig, const in vec3 center)
+{
+	float height = height(orig, 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;
+}
+
+// orig: ray origin
+// dir: ray direction
+// center: planet center
+// returns (h, t): h = length of perpendicular from center to ray
+//                 t = distance left from perpendicular to origin along the ray
+//                     (negative means perpendicular left behind)
+/*
+ *     t>0      H        t<0       dir
+ * ----*--------*--------*----------->
+ *              |
+ *              |
+ *              |
+ *              * O
+ */
+void findClosestHeight(out float h, out float t, const in vec3 orig, const in vec3 dir, const in vec3 center)
+{
+	vec3 radiusVector = center - orig;
+	vec3 tangent = dot(dir, radiusVector) * dir;
+	vec3 normal = radiusVector - tangent;
+	h = sqrt(dot(normal, normal));
+	t = dot(tangent, dir);
+}
+
+// error function approx., used in predictDensityIn
+float erf(const in float x)
+{
+	float a = 0.14001228904;
+	float four_over_pi = 1.27323954;
+
+	float x2 = x*x;
+	float r = -x2 * (four_over_pi + a * x2) / (1 + a * x2);
+	if (x > 0)
+		return  sqrt(1-exp(r)) * 0.5 + 0.5;
+	else
+		return -sqrt(1-exp(r)) * 0.5 + 0.5;
+}
+
+// predict out-scattering density based on coefficients
+float predictDensityOut(const in float atmosphereHeight, const in float height, const in float k, const in float b)
+{
+	if (height > atmosphereHeight)
+		return 0.f;
+
+	// earth atmospheric density: 1.225 kg/m^3, divided by 1e5
+	// 1/1.225e-5 = 81632.65306
+	float earthDensities = geosphereAtmosFogDensity * 81632.65306f;
+
+	if (height < 0)
+		return k / earthDensities;
+
+	return k * exp(-height * b) / earthDensities;
+}
+
+// predict in-scattering rate: from 0 to 1
+float predictDensityIn(const in float radius, const in float atmosphereHeight, const in float height, const in float c, const in float t)
+{
+	float minHeight = radius + height;
+	float h = sqrt(minHeight*minHeight + t*t); // height for our position
+
+	if (h > radius + atmosphereHeight) {
+		if (t > 0)
+			return 1.f; // no in-scattering, looking towards atmosphere
+		else
+			return 0.f; // looking from atmosphere
+	} else {
+		return erf(c * t); // erf is monotonic ascending
+	}
+}
+
+// predict "scattering density" along the ray
+// sample: starting point of the ray
+// dir:    direction of ray
+// center: relative position of planet center
+// radius: planet radius
+// atmosphereHeight: max height of atmosphere, bigger height means no atmosphere there
+// coefficients: pre-computed (k, b, c) for:
+//     k * exp(-height * b) for out-scattering density: k = density along the ray tangent to planet surface
+//                                                      b = height at which density reduces by e
+//     erf(c * t) for in-scattering density: c = density derivative per 1 km along the ray, assuming: k = 1, height = 0, t = 0 at tangent point
+float predictDensityInOut(const in vec3 sample, const in vec3 dir, const in vec3 center, const in float radius, const in float atmosphereHeight, const in vec3 coefficients)
+{
+    float h, t;
+    findClosestHeight(h, t, sample, dir, center);
+    h -= radius;
+
+    float opticalDepth = 0.f;
+    // find out-scattering density
+    opticalDepth = predictDensityOut(atmosphereHeight, h, coefficients.x, coefficients.y);
+
+    // apply in-scattering filter
+    opticalDepth *= predictDensityIn(radius, atmosphereHeight, h, coefficients.z, t);
+    return opticalDepth;
+}
+
+void skipRay(inout vec2 opticalDepth, const in vec3 dir, const in vec2 boundaries, const in vec3 center)
+{
+	if (boundaries.y == boundaries.x)
+		return;
+
+	int numSamples = 8;
+
+	float tCurrent = boundaries.x;
+	float segmentLength = (boundaries.y - boundaries.x) / numSamples;
+	for (int i = 0; i < numSamples; ++i) {
+		vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir;
+
+		// primary ray is approximated by (density * isegmentLength)
+		vec2 density;
+		scatter(density, samplePosition, center);
+		opticalDepth += exp(density) * segmentLength;
+
+		tCurrent += segmentLength;
+	}
+}
+
+void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const in vec3 sunDirection, const in vec3 dir, const in vec2 boundaries, const in vec3 center, const in vec4 diffuse, const in float uneclipsed)
+{
+	if (boundaries.y == boundaries.x)
+		return;
+
+	vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f);
+	vec3 betaM = vec3(21e-6f);
+
+	int numSamples = 16;
+
+	float atmosphereRadius = geosphereRadius * geosphereAtmosTopRad,
+	      atmosphereHeight = atmosphereRadius - geosphereRadius;
+
+	float tCurrent = boundaries.x;
+	float segmentLength = (boundaries.y - boundaries.x) / numSamples;
+	for (int i = 0; i < numSamples; ++i) {
+		vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir;
+
+		vec2 density;
+		scatter(density, samplePosition, center);
+		opticalDepth += exp(density) * segmentLength;
+
+		// light optical depth
+		vec2 opticalDepthLight = vec2(0.f);
+		vec3 samplePositionLight = samplePosition;
+
+		vec3 sampleGeoCenter = center - samplePosition;
+		opticalDepthLight.x = predictDensityInOut(samplePositionLight, sunDirection, sampleGeoCenter, geosphereRadius, atmosphereHeight, coefficientsR);
+		opticalDepthLight.y = predictDensityInOut(samplePositionLight, sunDirection, sampleGeoCenter, geosphereRadius, atmosphereHeight, coefficientsM);
+
+		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;
+		tCurrent += segmentLength;
+	}
+}
+
+// replace (a, b) by (b, a) if a > b
+vec2 sortAscending(const in vec2 segment)
+{
+	return (segment.x > segment.y) ? vec2(segment.y, segment.x) : segment;
+}
+
+// given a and b segments, return c = a \ b
+vec4 segmentSubtraction(const in vec2 a, const in vec2 b)
+{
+	vec2 as = sortAscending(a);
+	vec2 bs = sortAscending(b);
+
+	// b could be inside a, leaving segments at both sides
+	vec4 c;
+	c.x = a.x;
+	c.w = a.y;
+
+	c.y = min(a.y, max(a.x, b.x));
+	c.z = max(a.x, min(a.y, b.y));
+	return c;
+}
+
+vec3 computeIncidentLight(const in vec3 sunDirection, const in vec3 dir, const in vec3 center, const in vec2 atmosDist, const in vec4 diffuse, const in float uneclipsed)
+{
+	vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f);
+	vec3 betaM = vec3(21e-6f);
+
+	float atmosMin = atmosDist.x * geosphereRadius;
+	float atmosMax = atmosDist.y * geosphereRadius;
+
+	// solve Cylinder entry/exit dist
+	vec2 cylinder_intersect = rayCylinderIntersect(dir, center, sunDirection, geosphereRadius);
+	bool hasIntersect = cylinder_intersect.x != 0 || cylinder_intersect.y != 0;
+
+	vec3 cylinder_near = center - dir * cylinder_intersect.x;
+	vec3 cylinder_far  = center - dir * cylinder_intersect.y;
+
+	// test if ray passes through shadow
+	float a = dot(cylinder_near, sunDirection);
+	float b = dot(cylinder_far , sunDirection);
+	bool intersectsShadow = hasIntersect && (a > 0.f || b > 0.f);
+
+	vec2 ground_intersect = raySphereIntersect(center, dir, geosphereRadius);
+	bool shadowVisible = intersectsShadow && ground_intersect.x == 0.f;
+
+	/*
+	* We have three options:
+	* 1) Ray does not intersect shadow
+	*    Do nothing
+	* 2) Ray intersects shadow, starts inside
+	*    (cylinder_intersect.y, tmax)
+	* 3) Ray intersects shadow, starts outside
+	*    (tmin, cylinder_intersect.x) + (cylinder_intersect.y, tmax)
+	*/
+
+	vec2 atmosphere_intersect = raySphereIntersect(center, dir, geosphereRadius * geosphereAtmosTopRad);
+	vec4 atmosphere_minus_shadow = shadowVisible ? segmentSubtraction(atmosphere_intersect, cylinder_intersect) : vec4(atmosphere_intersect.x, atmosphere_intersect.y, atmosphere_intersect.y, atmosphere_intersect.y);
+
+	if (ground_intersect.x > 0.f) {
+		atmosphere_minus_shadow.yzw = min(atmosphere_minus_shadow.yzw, ground_intersect.x);
+	}
+	if (ground_intersect.y < 0.f) {
+		atmosphere_minus_shadow.xyz = max(atmosphere_minus_shadow.xyz, ground_intersect.y);
+	}
+
+	vec3 sumR = vec3(0.f);
+	vec3 sumM = vec3(0.f);
+	vec2 opticalDepth = vec2(0.f);
+
+	processRay(sumR, sumM, opticalDepth, sunDirection, dir, atmosphere_minus_shadow.xy, center, diffuse, uneclipsed);
+	skipRay(opticalDepth, dir, atmosphere_minus_shadow.yz, center);
+	processRay(sumR, sumM, opticalDepth, sunDirection, dir, atmosphere_minus_shadow.zw, center, diffuse, uneclipsed);
+
+	float mu = dot(dir, sunDirection); // mu in the paper which is the cosine of the angle between the sun direction and the ray direction
+	float phaseR = rayleighPhaseFunction(mu);
+	float phaseM = miePhaseFunction(0.76f, mu);
+
+	vec3 ret = (sumR * betaR * phaseR + sumM * betaM * phaseM);
+	return ret;
+}
diff --git a/data/shaders/opengl/rayleigh_accurate.frag b/data/shaders/opengl/rayleigh_accurate.frag
new file mode 100644
index 00000000000..2902b4c4c5d
--- /dev/null
+++ b/data/shaders/opengl/rayleigh_accurate.frag
@@ -0,0 +1,55 @@
+// Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+#include "attributes.glsl"
+#include "lib.glsl"
+#include "basesphere_uniforms.glsl"
+#include "rayleigh.glsl"
+
+uniform int NumShadows;
+
+in vec4 varyingEyepos;
+in vec4 vertexColor;
+
+out vec4 frag_color;
+
+void main(void)
+{
+	vec3 eyenorm = normalize(varyingEyepos.xyz);
+	vec3 specularHighlight = vec3(0.0);
+
+    vec2 atmosDist  = raySphereIntersect(geosphereCenter, eyenorm, geosphereAtmosTopRad);
+	// Invalid ray, skip shading this pixel
+	// (can improve performance when spatially coherent)
+	if (atmosDist.x == 0.0 && atmosDist.y == 0.0) {
+		frag_color = vec4(0.0);
+		return;
+	}
+
+	// a&b scaled so length of 1.0 means planet surface.
+	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));
+
+		float uneclipsed = clamp(calcUneclipsedSky(eclipse, NumShadows, a, b, lightDir), 0.0, 1.0);
+
+		// Convert from radius-relative to real coordinates
+		vec3 center = geosphereCenter * geosphereRadius;
+
+		vec3 lightPosAU = uLight[i].position.xyz / 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;
+	}
+#endif
+
+	vec4 color = vec4(specularHighlight.rgb, 1.0) * 20;
+
+	frag_color = toSRGB(1 - exp(-color));
+
+}
diff --git a/data/shaders/opengl/rayleigh_accurate.shaderdef b/data/shaders/opengl/rayleigh_accurate.shaderdef
new file mode 100644
index 00000000000..e092bd09fd6
--- /dev/null
+++ b/data/shaders/opengl/rayleigh_accurate.shaderdef
@@ -0,0 +1,13 @@
+## Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+## Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+Shader rayleigh_accurate
+
+Buffer LightData binding=0
+Buffer DrawData binding=1
+Buffer BaseSphereData binding=2
+
+PushConstant int NumShadows binding=0
+
+Vertex "rayleigh_accurate.vert"
+Fragment "rayleigh_accurate.frag"
diff --git a/data/shaders/opengl/rayleigh_accurate.vert b/data/shaders/opengl/rayleigh_accurate.vert
new file mode 100644
index 00000000000..58666a61abd
--- /dev/null
+++ b/data/shaders/opengl/rayleigh_accurate.vert
@@ -0,0 +1,14 @@
+// Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+#include "attributes.glsl"
+#include "lib.glsl"
+#include "basesphere_uniforms.glsl"
+
+out vec4 varyingEyepos;
+
+void main(void)
+{
+	gl_Position = matrixTransform();
+	varyingEyepos = uViewMatrix * a_vertex;
+}
diff --git a/data/shaders/opengl/rayleigh_fast.frag b/data/shaders/opengl/rayleigh_fast.frag
new file mode 100644
index 00000000000..afeb8ebf83f
--- /dev/null
+++ b/data/shaders/opengl/rayleigh_fast.frag
@@ -0,0 +1,19 @@
+// Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+#include "attributes.glsl"
+#include "lib.glsl"
+#include "basesphere_uniforms.glsl"
+#include "rayleigh.glsl"
+
+in vec4 varyingEyepos;
+in vec4 vertexColor;
+
+out vec4 frag_color;
+
+void main(void)
+{
+    vec4 color = vec4(vertexColor.rgb, 1.0) * 20;
+
+    frag_color = toSRGB(1 - exp(-color));
+}
diff --git a/data/shaders/opengl/rayleigh_fast.shaderdef b/data/shaders/opengl/rayleigh_fast.shaderdef
new file mode 100644
index 00000000000..db12bd88a76
--- /dev/null
+++ b/data/shaders/opengl/rayleigh_fast.shaderdef
@@ -0,0 +1,13 @@
+## Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+## Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+Shader rayleigh_fast
+
+Buffer LightData binding=0
+Buffer DrawData binding=1
+Buffer BaseSphereData binding=2
+
+PushConstant int NumShadows binding=0
+
+Vertex "rayleigh_fast.vert"
+Fragment "rayleigh_fast.frag"
diff --git a/data/shaders/opengl/rayleigh_fast.vert b/data/shaders/opengl/rayleigh_fast.vert
new file mode 100644
index 00000000000..cca12f91246
--- /dev/null
+++ b/data/shaders/opengl/rayleigh_fast.vert
@@ -0,0 +1,49 @@
+// Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+#include "attributes.glsl"
+#include "lib.glsl"
+#include "basesphere_uniforms.glsl"
+#include "rayleigh.glsl"
+
+uniform int NumShadows;
+
+out vec4 varyingEyepos;
+out vec4 vertexColor;
+
+void main(void)
+{
+	gl_Position = matrixTransform();
+	varyingEyepos = uViewMatrix * a_vertex;
+
+    // compute incident light from each light source
+	vec3 eyenorm = normalize(varyingEyepos.xyz);
+	vec3 specularHighlight = vec3(0.0);
+
+	vec2 atmosDist = raySphereIntersect(geosphereCenter, eyenorm, geosphereAtmosTopRad);
+
+	// a&b scaled so length of 1.0 means planet surface.
+	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));
+
+		float uneclipsed = clamp(calcUneclipsedSky(eclipse, NumShadows, a, b, lightDir), 0.0, 1.0);
+
+		// Convert from radius-relative to real coordinates
+		vec3 center = geosphereCenter * geosphereRadius;
+
+		vec3 lightPosAU = uLight[i].position.xyz / 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;
+
+	}
+#endif
+
+	vertexColor.rgb = specularHighlight;
+}
diff --git a/data/shaders/opengl/sRGB.glsl b/data/shaders/opengl/sRGB.glsl
new file mode 100644
index 00000000000..f94437f9da1
--- /dev/null
+++ b/data/shaders/opengl/sRGB.glsl
@@ -0,0 +1,24 @@
+// Copyright © 2008-2023 Pioneer Developers. See AUTHORS.txt for details
+// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
+
+#define GAMMA_FACTOR 2.2
+
+vec4 toLinear(in vec4 val)
+{
+	return pow(val, vec4(vec3(GAMMA_FACTOR), 1.0));
+}
+
+vec3 toLinear(in vec3 val)
+{
+	return pow(val, vec3(GAMMA_FACTOR));
+}
+
+vec4 toSRGB(in vec4 val)
+{
+	return pow(val, vec4(vec3(1.0 / GAMMA_FACTOR), 1.0));
+}
+
+vec3 toSRGB(in vec3 val)
+{
+	return pow(val, vec3(1.0 / GAMMA_FACTOR));
+}
diff --git a/src/BaseSphere.cpp b/src/BaseSphere.cpp
index e528d5b13f7..324074e8ef3 100644
--- a/src/BaseSphere.cpp
+++ b/src/BaseSphere.cpp
@@ -19,6 +19,9 @@ struct BaseSphereDataBlock {
 	float geosphereAtmosFogDensity;
 	float geosphereAtmosInvScaleHeight;
 	Color4f atmosColor;
+	alignas(16) vector3f coefficientsR;
+	alignas(16) vector3f coefficientsM;
+	alignas(16) vector2f scaleHeight;
 
 	// Eclipse struct data
 	alignas(16) vector3f shadowCentreX;
@@ -28,7 +31,7 @@ struct BaseSphereDataBlock {
 	alignas(16) vector3f lrad;
 	alignas(16) vector3f sdivlrad;
 };
-static_assert(sizeof(BaseSphereDataBlock) == 144, "");
+static_assert(sizeof(BaseSphereDataBlock) == 192, "");
 
 BaseSphere::BaseSphere(const SystemBody *body) :
 	m_sbody(body),
@@ -78,8 +81,8 @@ void BaseSphere::DrawAtmosphereSurface(Graphics::Renderer *renderer,
 	renderer->SetTransform(matrix4x4f(modelView * matrix4x4d::ScaleMatrix(rad) * invrot));
 
 	if (!m_atmos)
-		m_atmos.reset(new Drawables::Sphere3D(renderer, mat, 4, 1.0f, ATTRIB_POSITION));
-	m_atmos->Draw(renderer);
+		m_atmos.reset(new Drawables::Sphere3D(renderer, 4, 1.0f, ATTRIB_POSITION));
+	m_atmos->Draw(renderer, mat.Get());
 
 	renderer->GetStats().AddToStatCount(Graphics::Stats::STAT_ATMOSPHERES, 1);
 }
@@ -90,13 +93,16 @@ void BaseSphere::SetMaterialParameters(const matrix4x4d &trans, const float radi
 {
 	BaseSphereDataBlock matData{};
 
-	matData.geosphereCenter = vector3f(trans * vector3d(0.0));
-	matData.geosphereRadius = radius;
+	matData.geosphereCenter = vector3f(trans * vector3d(0.0)) / radius;
+	matData.geosphereRadius = ap.planetRadius;
 	matData.geosphereInvRadius = 1.0f / radius;
 	matData.geosphereAtmosTopRad = ap.atmosRadius;
 	matData.geosphereAtmosFogDensity = ap.atmosDensity;
 	matData.geosphereAtmosInvScaleHeight = ap.atmosInvScaleHeight;
 	matData.atmosColor = ap.atmosCol.ToColor4f();
+	matData.coefficientsR = ap.rayleighCoefficients;
+	matData.coefficientsM = ap.mieCoefficients;
+	matData.scaleHeight = ap.scaleHeight;
 
 	// we handle up to three shadows at a time
 	auto it = shadows.cbegin(), itEnd = shadows.cend();
@@ -115,7 +121,8 @@ void BaseSphere::SetMaterialParameters(const matrix4x4d &trans, const float radi
 	// 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()));
-	if (ap.atmosDensity > 0.0) {
+
+	if (m_atmosphereMaterial.Valid() && ap.atmosDensity > 0.0) {
 		m_atmosphereMaterial->SetBufferDynamic(s_baseSphereData, &matData);
 		m_atmosphereMaterial->SetPushConstant(s_numShadows, int(shadows.size()));
 	}
diff --git a/src/Camera.cpp b/src/Camera.cpp
index a3057301a53..55e523f1b2e 100644
--- a/src/Camera.cpp
+++ b/src/Camera.cpp
@@ -129,8 +129,6 @@ static void position_system_lights(Frame *camFrame, Frame *frame, std::vector<Ca
 	// IsRotFrame check prevents double counting
 	if (body && !frame->IsRotFrame() && (body->GetSuperType() == SystemBody::SUPERTYPE_STAR)) {
 		vector3d lpos = frame->GetPositionRelTo(camFrame->GetId());
-		const double dist = lpos.Length() / AU;
-		lpos *= 1.0 / dist; // normalize
 
 		const Color &col = StarSystem::starRealColors[body->GetType()];
 
diff --git a/src/GasGiant.cpp b/src/GasGiant.cpp
index 7d198e50b0a..9422242b8d3 100644
--- a/src/GasGiant.cpp
+++ b/src/GasGiant.cpp
@@ -707,7 +707,19 @@ void GasGiant::SetUpMaterials()
 		rsd.blendMode = Graphics::BLEND_ALPHA_ONE;
 		rsd.cullMode = Graphics::CULL_NONE;
 		rsd.depthWrite = false;
-		m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("geosphere_sky", skyDesc, rsd));
+
+		const int scattering = Pi::config->Int("RealisticScattering");
+		switch (scattering) {
+		case 1:
+			m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("rayleigh_fast", skyDesc, rsd));
+			break;
+		case 2:
+			m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("rayleigh_accurate", skyDesc, rsd));
+			break;
+		default:
+			m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("geosphere_sky", skyDesc, rsd));
+			break;
+		}
 	}
 }
 
diff --git a/src/GeoSphere.cpp b/src/GeoSphere.cpp
index eeb5b73f67b..6bd21cac3fa 100644
--- a/src/GeoSphere.cpp
+++ b/src/GeoSphere.cpp
@@ -83,6 +83,11 @@ void GeoSphere::OnChangeDetailLevel()
 		// reinit the terrain with the new settings
 		(*i)->m_terrain.Reset(Terrain::InstanceTerrain((*i)->GetSystemBody()));
 		print_info((*i)->GetSystemBody(), (*i)->m_terrain.Get());
+
+		// Reload the atmosphere material (scattering option)
+		if ((*i)->m_atmosphereMaterial.Valid()) {
+			(*i)->CreateAtmosphereMaterial();
+		}
 	}
 }
 
@@ -398,12 +403,13 @@ void GeoSphere::Render(Graphics::Renderer *renderer, const matrix4x4d &modelView
 	//XXX no need to calculate AP every frame
 	auto ap = GetSystemBody()->CalcAtmosphereParams();
 	SetMaterialParameters(trans, radius, shadows, ap);
-	if (ap.atmosDensity > 0.0) {
+
+	if (m_atmosphereMaterial.Valid() && ap.atmosDensity > 0.0) {
 		// make atmosphere sphere slightly bigger than required so
 		// that the edges of the pixel shader atmosphere jizz doesn't
 		// show ugly polygonal angles
 		DrawAtmosphereSurface(renderer, trans, campos,
-			ap.atmosRadius * 1.01,
+			ap.atmosRadius * 1.02,
 			m_atmosphereMaterial);
 	}
 
@@ -490,17 +496,32 @@ void GeoSphere::SetUpMaterials()
 		m_surfaceMaterial->SetTexture("texture1"_hash, m_texLo.Get());
 	}
 
-	{
-		Graphics::RenderStateDesc rsd;
-		// atmosphere is blended over the background
-		rsd.blendMode = Graphics::BLEND_ALPHA_ONE;
-		rsd.cullMode = Graphics::CULL_NONE;
-		rsd.depthWrite = false;
-
-		Graphics::MaterialDescriptor skyDesc;
-		skyDesc.effect = Graphics::EFFECT_GEOSPHERE_SKY;
-		skyDesc.lighting = true;
-		skyDesc.quality |= Graphics::HAS_ECLIPSES;
+	CreateAtmosphereMaterial();
+}
+
+void GeoSphere::CreateAtmosphereMaterial()
+{
+	Graphics::RenderStateDesc rsd;
+	// atmosphere is blended over the background
+	rsd.blendMode = Graphics::BLEND_ALPHA_ONE;
+	rsd.cullMode = Graphics::CULL_FRONT;
+	rsd.depthWrite = false;
+
+	Graphics::MaterialDescriptor skyDesc;
+	skyDesc.effect = Graphics::EFFECT_GEOSPHERE_SKY;
+	skyDesc.lighting = true;
+	skyDesc.quality |= Graphics::HAS_ECLIPSES;
+
+	const int scattering = Pi::config->Int("RealisticScattering");
+	switch (scattering) {
+	case 1:
+		m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("rayleigh_fast", skyDesc, rsd));
+		break;
+	case 2:
+		m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("rayleigh_accurate", skyDesc, rsd));
+		break;
+	default:
 		m_atmosphereMaterial.Reset(Pi::renderer->CreateMaterial("geosphere_sky", skyDesc, rsd));
+		break;
 	}
 }
diff --git a/src/GeoSphere.h b/src/GeoSphere.h
index 1f74e21238d..94cb3d10f71 100644
--- a/src/GeoSphere.h
+++ b/src/GeoSphere.h
@@ -101,6 +101,7 @@ class GeoSphere : public BaseSphere {
 	static RefCountedPtr<GeoPatchContext> s_patchContext;
 
 	virtual void SetUpMaterials() override;
+	void CreateAtmosphereMaterial();
 
 	RefCountedPtr<Graphics::Texture> m_texHi;
 	RefCountedPtr<Graphics::Texture> m_texLo;
diff --git a/src/SectorMap.cpp b/src/SectorMap.cpp
index b7ad3b466d5..3ab0e1e575b 100644
--- a/src/SectorMap.cpp
+++ b/src/SectorMap.cpp
@@ -648,8 +648,8 @@ void SectorMap::Draw3D()
 	if (!m_sphereParams.empty()) {
 		for (const auto& s : m_sphereParams) {
 			m_context.renderer->SetTransform(s.trans);
-			m_sphere->GetMaterial()->diffuse = s.color;
-			m_sphere->Draw(m_context.renderer);
+			m_fresnelMat->diffuse = s.color;
+			m_sphere->Draw(m_context.renderer, m_fresnelMat.Get());
 		}
 	}
 }
@@ -1208,7 +1208,7 @@ void SectorMap::Update(float frameTime)
 		Graphics::MaterialDescriptor matdesc;
 		m_fresnelMat.Reset(m_context.renderer->CreateMaterial("fresnel_sphere", matdesc, rsd));
 		m_fresnelMat->diffuse = Color::WHITE;
-		m_sphere.reset(new Graphics::Drawables::Sphere3D(m_context.renderer, m_fresnelMat, 4, 1.0f));
+		m_sphere.reset(new Graphics::Drawables::Sphere3D(m_context.renderer, 4, 1.0f));
 	}
 	m_sphereParams.clear();
 	m_lineVerts->Clear();
diff --git a/src/galaxy/AtmosphereParameters.h b/src/galaxy/AtmosphereParameters.h
index 16afd0c10fe..70d3bcc528d 100644
--- a/src/galaxy/AtmosphereParameters.h
+++ b/src/galaxy/AtmosphereParameters.h
@@ -11,6 +11,9 @@ struct AtmosphereParameters {
 	float planetRadius;
 	Color atmosCol;
 	vector3d center;
+	vector3f rayleighCoefficients;
+	vector3f mieCoefficients;
+	vector2f scaleHeight;
 };
 
 #endif // ATMOSPHEREPARAMETERS_H_INCLUDED
diff --git a/src/galaxy/SystemBody.cpp b/src/galaxy/SystemBody.cpp
index e88c561f7b3..f9b988c1204 100644
--- a/src/galaxy/SystemBody.cpp
+++ b/src/galaxy/SystemBody.cpp
@@ -262,6 +262,53 @@ bool SystemBody::IsScoopable() const
 	return false;
 }
 
+// this is a reduced version of original Rayleigh scattering function used to compute density once per planet
+// (at least it meant not to be once per frame) given:
+// planet radius, height of atmosphere (not its radius), length of perpendicular between ray and planet center, lapse rate at which density fades out by rate of e
+// this function is used in GetCoefficients
+// all input units are in km
+double SystemBody::ComputeDensity(const double radius, const double atmosphereHeight, const double h, const double scaleHeight) const
+{
+	int numSamples = 16;
+	double totalHeight = radius + atmosphereHeight;
+	double minHeight = radius + h;
+	double tmax = sqrt(totalHeight * totalHeight - minHeight * minHeight); // maximum t
+	double tmin = -tmax;
+	double length = tmax - tmin;
+	double step = length / numSamples;
+
+	double density = 0.0;
+	for (int i = 0; i < numSamples; ++i) {
+		double t = tmin + step * (i + 0.5);
+		double h = sqrt(minHeight * minHeight + t * t) - radius;
+		density += step * exp(-h / scaleHeight);
+	}
+
+	return density;
+}
+
+// these coefficients are to be passed into shaders, meant to accelerate scattering computation per-pixel
+// instead of sampling each one (which I suspect is VERY SLOW)
+// all input units are in km
+vector3f SystemBody::GetCoefficients(const double radius, const double atmHeight, const double scaleHeight) const
+{
+	float k, b, c;
+
+	// compute full out-scattering densities at 0 and 1 km heights
+	// k = density at 0 km
+	// b = log(density0km / density1km) - log is used to multiply it by actual height
+	k = ComputeDensity(radius, atmHeight, 0.f, scaleHeight);
+	b = log(k / ComputeDensity(radius, atmHeight, 1.f, scaleHeight));
+
+	// compute c - erf coefficient, which adjusts slope of erf near t=0 to match actual in-scattering
+	float erf1_minus_erf0 = 0.421463;
+	float sHeight = sqrt(radius * radius + 1.f) - radius;
+	float c1 = exp(-sHeight / scaleHeight);
+	float c2 = k * erf1_minus_erf0;
+	c = c1 / c2;
+	return vector3f(k, b, c);
+}
+
 // Calculate parameters used in the atmospheric model for shaders
 AtmosphereParameters SystemBody::CalcAtmosphereParams() const
 {
@@ -320,6 +367,12 @@ AtmosphereParameters SystemBody::CalcAtmosphereParams() const
 
 	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);
+	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);
+
 	return params;
 }
 
diff --git a/src/galaxy/SystemBody.h b/src/galaxy/SystemBody.h
index 620b88662a3..c91fc7ea3ec 100644
--- a/src/galaxy/SystemBody.h
+++ b/src/galaxy/SystemBody.h
@@ -292,6 +292,10 @@ class SystemBody : public RefCounted, public SystemBodyType, protected SystemBod
 	// Calculate atmosphere density at given altitude and pressure (kg/m^3)
 	double GetAtmDensity(double altitude, double pressure) const;
 
+	// for rayleigh scattering
+	double ComputeDensity(const double radius, const double atmosphereHeight, const double h, const double scaleHeight) const;
+	vector3f GetCoefficients(const double radius, const double atmHeight, const double scaleHeight) const;
+
 	AtmosphereParameters CalcAtmosphereParams() const;
 
 	bool IsScoopable() const;
diff --git a/src/graphics/Drawables.cpp b/src/graphics/Drawables.cpp
index b06ff0a737a..34528543e90 100644
--- a/src/graphics/Drawables.cpp
+++ b/src/graphics/Drawables.cpp
@@ -483,7 +483,7 @@ namespace Graphics {
 
 		Graphics::MeshObject *Icosphere::Generate(Graphics::Renderer *r, int subdivs, float scale, AttributeSet attribs)
 		{
-			subdivs = Clamp(subdivs, 0, 4);
+			subdivs = Clamp(subdivs, 0, 10);
 			scale = fabs(scale);
 			matrix4x4f trans = matrix4x4f::Identity();
 			trans.Scale(scale, scale, scale);
@@ -568,19 +568,18 @@ namespace Graphics {
 
 		//------------------------------------------------------------
 
-		Sphere3D::Sphere3D(Renderer *renderer, RefCountedPtr<Material> mat, int subdivs, float scale, AttributeSet attribs)
+		Sphere3D::Sphere3D(Renderer *renderer, int subdivs, float scale, AttributeSet attribs)
 		{
 			PROFILE_SCOPED()
 			assert(attribs.HasAttrib(ATTRIB_POSITION));
 
-			m_material = mat;
 			m_sphereMesh.reset(Icosphere::Generate(renderer, subdivs, scale, attribs));
 		}
 
-		void Sphere3D::Draw(Renderer *r)
+		void Sphere3D::Draw(Renderer *r, Material *mat)
 		{
 			PROFILE_SCOPED()
-			r->DrawMesh(m_sphereMesh.get(), m_material.Get());
+			r->DrawMesh(m_sphereMesh.get(), mat);
 		}
 
 		//------------------------------------------------------------
diff --git a/src/graphics/Drawables.h b/src/graphics/Drawables.h
index 092f7d34f4c..e0459d7c7fe 100644
--- a/src/graphics/Drawables.h
+++ b/src/graphics/Drawables.h
@@ -151,14 +151,11 @@ namespace Graphics {
 		class Sphere3D {
 		public:
 			//subdivisions must be 0-4
-			Sphere3D(Renderer *, RefCountedPtr<Material> material, int subdivisions = 0, float scale = 1.f, AttributeSet attribs = (ATTRIB_POSITION | ATTRIB_NORMAL | ATTRIB_UV0));
-			void Draw(Renderer *r);
-
-			RefCountedPtr<Material> GetMaterial() const { return m_material; }
+			Sphere3D(Renderer *, int subdivisions = 0, float scale = 1.f, AttributeSet attribs = (ATTRIB_POSITION | ATTRIB_NORMAL | ATTRIB_UV0));
+			void Draw(Renderer *r, Material *m);
 
 		private:
 			std::unique_ptr<MeshObject> m_sphereMesh;
-			RefCountedPtr<Material> m_material;
 		};
 		//------------------------------------------------------------
 
@@ -252,7 +249,7 @@ namespace Graphics {
 		// The visual density of the grid can be controlled by the lineSpacing parameter
 		class GridSphere {
 		public:
-			GridSphere(Graphics::Renderer *r, uint32_t numSubdivs = 12);
+			GridSphere(Graphics::Renderer *r, uint32_t numSubdivs = 4);
 
 			void SetLineColors(Color minorLineColor, Color majorLineColor, float lineWidth = 2.0);
 
diff --git a/src/lua/LuaEngine.cpp b/src/lua/LuaEngine.cpp
index 524b009db72..92405ee75d0 100644
--- a/src/lua/LuaEngine.cpp
+++ b/src/lua/LuaEngine.cpp
@@ -754,6 +754,23 @@ static int l_engine_set_gpu_jobs_enabled(lua_State *l)
 	return 0;
 }
 
+static int l_engine_get_realistic_scattering(lua_State *l)
+{
+	lua_pushinteger(l, Pi::config->Int("RealisticScattering"));
+	return 1;
+}
+
+static int l_engine_set_realistic_scattering(lua_State *l)
+{
+	const int scattering = luaL_checkinteger(l, 1);
+	if (scattering != Pi::config->Int("RealisticScattering")) {
+		Pi::config->SetInt("RealisticScattering", scattering);
+		Pi::config->Save();
+		Pi::OnChangeDetailLevel();
+	}
+	return 0;
+}
+
 static int l_engine_is_intro_zooming(lua_State *l)
 {
 	if (Pi::intro) {
@@ -1064,6 +1081,9 @@ void LuaEngine::Register()
 		{ "GetGpuJobsEnabled", l_engine_get_gpu_jobs_enabled },
 		{ "SetGpuJobsEnabled", l_engine_set_gpu_jobs_enabled },
 
+		{ "GetRealisticScattering", l_engine_get_realistic_scattering },
+		{ "SetRealisticScattering", l_engine_set_realistic_scattering },
+
 		{ "GetPlanetDetailLevel", l_engine_get_planet_detail_level },
 		{ "SetPlanetDetailLevel", l_engine_set_planet_detail_level },
 		{ "GetCityDetailLevel", l_engine_get_city_detail_level },