diff --git a/Shaders/std/math.glsl b/Shaders/std/math.glsl
index 93293401..a4a64138 100755
--- a/Shaders/std/math.glsl
+++ b/Shaders/std/math.glsl
@@ -23,7 +23,7 @@ vec2 rand2(const vec2 coord) {
 	const float width = 1100.0;
 	const float height = 500.0;
 	float noiseX = ((fract(1.0 - coord.s * (width / 2.0)) * 0.25) + (fract(coord.t * (height / 2.0)) * 0.75)) * 2.0 - 1.0;
-	float noiseY = ((fract(1.0 - coord.s * (width / 2.0)) * 0.75) + (fract(coord.t * (height / 2.0)) * 0.25)) * 2.0 - 1.0;	
+	float noiseY = ((fract(1.0 - coord.s * (width / 2.0)) * 0.75) + (fract(coord.t * (height / 2.0)) * 0.25)) * 2.0 - 1.0;
 	return vec2(noiseX, noiseY);
 }
 
@@ -40,4 +40,9 @@ float attenuate(const float dist) {
 	// 1.0 / (quadratic * dist * dist);
 }
 
+float safe_acos(const float x) {
+	// acos is undefined if |x| > 1
+	return acos(clamp(x, -1.0, 1.0));
+}
+
 #endif
diff --git a/Shaders/std/sky.glsl b/Shaders/std/sky.glsl
index b4465351..eaac4352 100644
--- a/Shaders/std/sky.glsl
+++ b/Shaders/std/sky.glsl
@@ -20,6 +20,8 @@
 #ifndef _SKY_GLSL_
 #define _SKY_GLSL_
 
+#include "std/math.glsl"
+
 uniform sampler2D nishitaLUT;
 uniform vec2 nishitaDensity;
 
@@ -44,20 +46,9 @@ uniform vec2 nishitaDensity;
 #define nishita_mie_dir 0.76 // Aerosols anisotropy ("direction")
 #define nishita_mie_dir_sq 0.5776 // Squared aerosols anisotropy
 
-// The ozone absorption coefficients are taken from Cycles code.
-// Because Cycles calculates 21 wavelengths, we use the coefficients
-// which are closest to the RGB wavelengths (645nm, 510nm, 440nm).
-// Precalculating values by simulating Blender's spec_to_xyz() function
-// to include all 21 wavelengths gave unrealistic results
-#define nishita_ozone_coeff vec3(1.59051840791988e-6, 0.00000096707041180970, 0.00000007309568762914)
-
 // Values from [Hill: 60]
 #define sun_limb_darkening_col vec3(0.397, 0.503, 0.652)
 
-float random(vec2 coords) {
-	return fract(sin(dot(coords.xy, vec2(12.9898,78.233))) * 43758.5453);
-}
-
 vec3 nishita_lookupLUT(const float height, const float sunTheta) {
 	vec2 coords = vec2(
 		sqrt(height * (1 / nishita_atmo_radius)),
@@ -124,18 +115,19 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa
 
 		// Idea behind this: "Rotate" everything by iPos (-> iPos is the new zenith) and then all calculations for the
 		// inner integral only depend on the sample height (iHeight) and sunTheta (angle between sun and new zenith).
-		float sunTheta = acos(dot(normalize(iPos), normalize(pSun)));
-		vec3 jODepth = nishita_lookupLUT(iHeight, sunTheta);
-
-		// Apply dithering to reduce visible banding
-		jODepth += mix(-1000, 1000, random(r.xy));
+		float sunTheta = safe_acos(dot(normalize(iPos), normalize(pSun)));
+		vec3 jAttn = nishita_lookupLUT(iHeight, sunTheta);
 
 		// Calculate attenuation
-		vec3 attn = exp(-(
-			nishita_mie_coeff * (iOdMie + jODepth.y)
-			+ (nishita_rayleigh_coeff) * (iOdRlh + jODepth.x)
-			+ nishita_ozone_coeff * jODepth.z
+		vec3 iAttn = exp(-(
+			nishita_mie_coeff * iOdMie
+			+ nishita_rayleigh_coeff * iOdRlh
+			// + 0 for ozone
 		));
+		vec3 attn = iAttn * jAttn;
+
+		// Apply dithering to reduce visible banding
+		attn *= 0.98 + rand(r.xy) * 0.04;
 
 		// Accumulate scattering
 		totalRlh += odStepRlh * attn;
diff --git a/Sources/armory/renderpath/Nishita.hx b/Sources/armory/renderpath/Nishita.hx
index 7dcba4a2..177cb517 100644
--- a/Sources/armory/renderpath/Nishita.hx
+++ b/Sources/armory/renderpath/Nishita.hx
@@ -80,11 +80,24 @@ class NishitaData {
 	**/
 	public static var radiusPlanet = 6360000;
 
+	/** Rayleigh scattering coefficient. **/
+	public static var rayleighCoeff = new Vec3(5.5e-6, 13.0e-6, 22.4e-6);
 	/** Rayleigh scattering scale parameter. **/
 	public static var rayleighScale = 8e3;
+
+	/** Mie scattering coefficient. **/
+	public static var mieCoeff = 2e-5;
 	/** Mie scattering scale parameter. **/
 	public static var mieScale = 1.2e3;
 
+	/** Ozone scattering coefficient. **/
+	// The ozone absorption coefficients are taken from Cycles code.
+	// Because Cycles calculates 21 wavelengths, we use the coefficients
+	// which are closest to the RGB wavelengths (645nm, 510nm, 440nm).
+	// Precalculating values by simulating Blender's spec_to_xyz() function
+	// to include all 21 wavelengths gave unrealistic results.
+	public static var ozoneCoeff = new Vec3(1.59051840791988e-6, 0.00000096707041180970, 0.00000007309568762914);
+
 	public function new() {}
 
 	/** Approximates the density of ozone for a given sample height. **/
@@ -185,6 +198,29 @@ class NishitaData {
 			jTime += jStepSize;
 		}
 
-		return jODepth.mult(jStepSize);
+		jODepth.mult(jStepSize);
+
+		// Precalculate a part of the secondary attenuation.
+		// For one variable (e.g. x) in the vector, the formula is as follows:
+		//
+		// attn.x = exp(-(coeffX * (firstOpticalDepth.x + secondOpticalDepth.x)))
+		//
+		// We can split that up via:
+		//
+		// attn.x = exp(-(coeffX * firstOpticalDepth.x + coeffX * secondOpticalDepth.x))
+		//        = exp(-(coeffX * firstOpticalDepth.x)) * exp(-(coeffX * secondOpticalDepth.x))
+		//
+		// The first factor of the resulting multiplication is calculated in the
+		// shader, but we can already precalculate the second one. As a side
+		// effect this keeps the range of the LUT values small because we don't
+		// store the optical depth but the attenuation.
+		var jAttenuation = new Vec3();
+		var mie = mieCoeff * jODepth.y;
+		jAttenuation.addf(mie, mie, mie);
+		jAttenuation.add(rayleighCoeff.clone().mult(jODepth.x));
+		jAttenuation.add(ozoneCoeff.clone().mult(jODepth.z));
+		jAttenuation.exp(jAttenuation.mult(-1));
+
+		return jAttenuation;
 	}
 }