diff --git a/Effects/terrain-default.eff b/Effects/terrain-default.eff
index 8a79ace62..a622eae53 100644
--- a/Effects/terrain-default.eff
+++ b/Effects/terrain-default.eff
@@ -112,6 +112,7 @@
     <cloud_self_shading><use>/environment/cloud-self-shading</use></cloud_self_shading>
     <moonlight><use>/environment/moonlight</use></moonlight>
     <season><use>/environment/season</use></season>
+    <air_pollution><use>/environment/air-pollution-norm</use></air_pollution>
     <wind_effects><use>/sim/rendering/shaders/wind-effects</use></wind_effects>
     <windE><use>/environment/sea/surface/wind-from-east-fps</use></windE>
     <windN><use>/environment/sea/surface/wind-from-north-fps</use></windN>
@@ -293,6 +294,7 @@
         <fragment-shader>Shaders/terrain-haze-ultra.frag</fragment-shader>
         <fragment-shader>Shaders/noise.frag</fragment-shader>
         <fragment-shader>Shaders/cloud-shadowfunc.frag</fragment-shader>
+        <fragment-shader>Shaders/hazes.frag</fragment-shader>
       </program>
 	  <uniform>
         <name>grain_strength</name>
@@ -429,6 +431,11 @@
         <type>float</type>
         <value><use>season</use></value>
       </uniform>
+      <uniform>
+        <name>air_pollution</name>
+        <type>float</type>
+        <value><use>air_pollution</use></value>
+      </uniform>
       <uniform>
         <name>cloudpos1_x</name>
         <type>float</type>
diff --git a/Effects/water-inland.eff b/Effects/water-inland.eff
index 09796b3c3..79e41c36d 100644
--- a/Effects/water-inland.eff
+++ b/Effects/water-inland.eff
@@ -362,6 +362,7 @@
 				<vertex-shader>Shaders/water_lightfield.vert</vertex-shader>
 				<fragment-shader>Shaders/water_lightfield.frag</fragment-shader>
 				<fragment-shader>Shaders/cloud-shadowfunc.frag</fragment-shader>
+				<fragment-shader>Shaders/hazes.frag</fragment-shader>
 			</program>
 			<!--<uniform>
 				<name>water_reflection</name>
@@ -532,7 +533,17 @@
         			<type>float</type>
         			<value><use>cloud_self_shading</use></value>
       			</uniform>
-		<uniform>
+			<uniform>
+				<name>fogstructure</name>
+				<type>float</type>
+				<value><use>fogstructure</use></value>
+			</uniform>
+			<uniform>
+				<name>air_pollution</name>
+				<type>float</type>
+				<value><use>air_pollution</use></value>
+			</uniform>
+			<uniform>
        				<name>cloudpos1_x</name>
         			<type>float</type>
         		   	<value><use>cloudpos1_x</use></value>
diff --git a/Effects/water.eff b/Effects/water.eff
index 64c4cb7a1..e53f90251 100644
--- a/Effects/water.eff
+++ b/Effects/water.eff
@@ -131,6 +131,7 @@
     		<moonlight><use>/environment/moonlight</use></moonlight>
 		<fogtype><use>/sim/rendering/shaders/skydome</use></fogtype>
 		<fogstructure><use>/environment/fog-structure</use></fogstructure>
+		<air_pollution><use>/environment/air-pollution-norm</use></air_pollution>
 		<ice_cover><use>/environment/sea/surface/ice-cover</use></ice_cover>
 		<quality_level><use>/sim/rendering/shaders/landmass</use></quality_level>
 		<!-- sea colors -->
@@ -373,6 +374,7 @@
 				<vertex-shader>Shaders/water_lightfield.vert</vertex-shader>
 				<fragment-shader>Shaders/water_lightfield.frag</fragment-shader>
 				<fragment-shader>Shaders/cloud-shadowfunc.frag</fragment-shader>
+				<fragment-shader>Shaders/hazes.frag</fragment-shader>
 			</program>
 			<!--<uniform>
 				<name>water_reflection</name>
@@ -558,6 +560,11 @@
 				<type>float</type>
 				<value><use>fogstructure</use></value>
 			</uniform>
+			<uniform>
+				<name>air_pollution</name>
+				<type>float</type>
+				<value><use>air_pollution</use></value>
+			</uniform>
 			<uniform>
         			<name>ice_cover</name>
         			<type>float</type>
diff --git a/Shaders/hazes.frag b/Shaders/hazes.frag
new file mode 100644
index 000000000..03599dd4d
--- /dev/null
+++ b/Shaders/hazes.frag
@@ -0,0 +1,74 @@
+// -*-C++-*-
+// standard ALS fog function with exp(-d/D) fading and cutoff at low altitude and exp(-d^2/D^2) at high altitude
+
+const float AtmosphericScaleHeight = 8500.0;
+
+float fog_func (in float targ, in float alt)
+{
+
+
+float fade_mix;
+
+// for large altitude > 30 km, we switch to some component of quadratic distance fading to
+// create the illusion of improved visibility range
+
+targ = 1.25 * targ * smoothstep(0.04,0.06,targ); // need to sync with the distance to which terrain is drawn
+
+
+if (alt < 30000.0)
+	{return exp(-targ - targ * targ * targ * targ);}
+else if (alt < 50000.0)
+	{
+	fade_mix = (alt - 30000.0)/20000.0;
+	return fade_mix * exp(-targ*targ - pow(targ,4.0)) + (1.0 - fade_mix) * exp(-targ - pow(targ,4.0));	
+	}
+else 
+	{
+	return exp(- targ * targ - pow(targ,4.0));
+	}
+
+}
+
+// altitude correction for exponential drop in atmosphere density
+
+float alt_factor(in float eye_alt, in float vertex_alt)
+{
+float h0 = AtmosphericScaleHeight;
+float h1 = min(eye_alt,vertex_alt);
+float h2 = max(eye_alt,vertex_alt);
+
+
+if ((h2-h1) < 200.0) // use a Taylor-expanded version
+	{
+	return  0.5 * (exp(-h2/h0) + exp(-h1/h0));
+	}
+else
+	{
+	return h0/(h2-h1) * (exp(-h1/h0) - exp(-h2/h0));
+	}
+}
+
+
+// Rayleigh in-scatter function
+
+float rayleigh_in_func(in float dist, in float air_pollution, in float avisibility, in float eye_alt, in float vertex_alt)
+{
+
+float fade_length = avisibility * (2.5 - 2.0 * air_pollution);
+
+fade_length = fade_length / alt_factor(eye_alt, vertex_alt);
+
+return 1.0-exp(-dist/max(35000.0,fade_length));
+}
+
+
+// Rayleigh out-scattering color shift
+
+vec3 rayleigh_out_shift(in vec3 color, in float outscatter)
+{
+color.r = color.r * (1.0 - 0.4 * outscatter);
+color.g = color.g * (1.0 - 0.8 * outscatter);
+color.b = color.b * (1.0 - 1.6 * outscatter);
+
+return color;
+} 
diff --git a/Shaders/terrain-haze-ultra.frag b/Shaders/terrain-haze-ultra.frag
index f8bd078d7..f634bc03c 100644
--- a/Shaders/terrain-haze-ultra.frag
+++ b/Shaders/terrain-haze-ultra.frag
@@ -42,6 +42,7 @@ uniform float fogstructure;
 uniform float snow_thickness_factor;
 uniform float cloud_self_shading;
 uniform float season;
+uniform float air_pollution;
 uniform float grain_strength;
 uniform float intrinsic_wetness;
 uniform float transition_model;
@@ -62,7 +63,7 @@ uniform int rock_strata;
 const float EarthRadius = 5800000.0;
 const float terminator_width = 200000.0;
 
-float alt;
+//float alt;
 float eShade;
 float yprime_alt;
 float mie_angle;
@@ -73,8 +74,10 @@ float Noise2D(in vec2 coord, in float wavelength);
 float Noise3D(in vec3 coord, in float wavelength);
 float SlopeLines2D(in vec2 coord, in vec2 gradDir, in float wavelength, in float steepness);
 float Strata3D(in vec3 coord, in float wavelength, in float variation);
-
-
+float fog_func (in float targ, in float alt);
+float rayleigh_in_func(in float dist, in float air_pollution, in float avisibility, in float eye_alt, in float vertex_alt);
+float alt_factor(in float eye_alt, in float vertex_alt);
+vec3 rayleigh_out_shift(in vec3 color, in float outscatter);
 
 float light_func (in float x, in float a, in float b, in float c, in float d, in float e)
 {
@@ -98,39 +101,11 @@ return 1.0 - smoothstep(0.5 * fade_dist, fade_dist, dist);
 }
 
 
-// this determines how light is attenuated in the distance
-// physically this should be exp(-arg) but for technical reasons we use a sharper cutoff
-// for distance > visibility
-
-float fog_func (in float targ)
-{
-
-
-float fade_mix;
-
-// for large altitude > 30 km, we switch to some component of quadratic distance fading to
-// create the illusion of improved visibility range
-
-targ = 1.25 * targ * smoothstep(0.04,0.06,targ); // need to sync with the distance to which terrain is drawn
-
-
-if (alt < 30000.0)
-	{return exp(-targ - targ * targ * targ * targ);}
-else if (alt < 50000.0)
-	{
-	fade_mix = (alt - 30000.0)/20000.0;
-	return fade_mix * exp(-targ*targ - pow(targ,4.0)) + (1.0 - fade_mix) * exp(-targ - pow(targ,4.0));	
-	}
-else 
-	{
-	return exp(- targ * targ - pow(targ,4.0));
-	}
-
-}
 
 void main()
 {
 
+float alt;
 
 yprime_alt = diffuse_term.a;
 //diffuse_term.a = 1.0;
@@ -298,7 +273,7 @@ float snownoise_50m = mix(noise_50m, slopenoise_100m, clamp(3.0*(1.0-steepness),
 	stprime = stprime * distortion_factor * 15.0;
 	stprime = stprime + normalize(relPos).xy * 0.022 * (noise_10m + 0.5 * noise_5m +0.25 * noise_2m - 0.875 );
 	
-    detail_texel = texture2D(detail_texture, stprime);
+    	detail_texel = texture2D(detail_texture, stprime);
 	if (detail_texel.a <0.1) {flag = 0;}
 	
 
@@ -486,14 +461,17 @@ if ((dist < 5000.0)&& (quality_level > 3) && (combined_wetness>0.0))
 
     fragColor = color * texel + specular;
 
-// here comes the terrain haze model
 
+// Rayleigh color shift due to out-scattering
+    float rayleigh_length = 0.5 * avisibility * (2.5 - 1.9 * air_pollution)/alt_factor(eye_alt, eye_alt+relPos.z);
+    float outscatter = 1.0-exp(-dist/rayleigh_length);
+    fragColor.rgb = rayleigh_out_shift(fragColor.rgb,outscatter);
+
+// here comes the terrain haze model
 
 float delta_z = hazeLayerAltitude - eye_alt;
 
 if (dist > 0.04 * min(visibility,avisibility)) 
-//if ((gl_FragCoord.y > ylimit) || (gl_FragCoord.x < zlimit1) || (gl_FragCoord.x > zlimit2))
-//if (dist > 40.0)
 {
 
 alt = eye_alt;
@@ -604,7 +582,7 @@ else
 
 
 
-transmission =  fog_func(transmission_arg);
+transmission =  fog_func(transmission_arg, alt);
 
 // there's always residual intensity, we should never be driven to zero
 if (eqColorFactor < 0.2) eqColorFactor = 0.2;
@@ -661,6 +639,15 @@ if (intensity > 0.0) // this needs to be a condition, because otherwise hazeColo
 	}
 
 
+// blue Rayleigh scattering with distance
+
+float rShade = 0.9 * smoothstep(terminator_width+ terminator, -terminator_width + terminator, yprime_alt-340000.0) + 0.1;
+float lightIntensity = length(diffuse_term.rgb)/1.73 * rShade;
+vec3 rayleighColor = vec3 (0.17, 0.52, 0.87) * lightIntensity;
+float rayleighStrength = rayleigh_in_func(dist, air_pollution, avisibility/max(lightIntensity,0.05), eye_alt, eye_alt + relPos.z);
+fragColor.rgb = mix(fragColor.rgb, rayleighColor,rayleighStrength);
+
+// finally, mix fog in
 
 
 fragColor.rgb = mix(eqColorFactor * hazeColor * eShade , fragColor.rgb,transmission);
@@ -672,6 +659,15 @@ gl_FragColor = fragColor;
 }
 else // if dist < threshold no fogging at all 
 {
+
+// blue Rayleigh scattering with distance
+
+float rShade = 0.9 * smoothstep(terminator_width+ terminator, -terminator_width + terminator, yprime_alt-340000.0) + 0.1;
+float lightIntensity = length(diffuse_term.rgb)/1.73 * rShade;
+vec3 rayleighColor = vec3 (0.17, 0.52, 0.87) * lightIntensity;
+float rayleighStrength = rayleigh_in_func(dist, air_pollution, avisibility/max(lightIntensity,0.05), eye_alt, eye_alt + relPos.z);
+fragColor.rgb = mix(fragColor.rgb, rayleighColor,rayleighStrength);
+
 gl_FragColor =  fragColor;
 }
 
diff --git a/Shaders/terrain-haze-ultra.vert b/Shaders/terrain-haze-ultra.vert
index 1ab79cd9a..00e912cd6 100644
--- a/Shaders/terrain-haze-ultra.vert
+++ b/Shaders/terrain-haze-ultra.vert
@@ -271,6 +271,8 @@ float shade_depth =  1.0 * smoothstep (0.6,0.95,ground_scattering) * (1.0-smooth
    light_ambient.rgb = light_ambient.rgb * (1.0 - shade_depth);
    light_diffuse.rgb = light_diffuse.rgb * (1.0 + 1.2 * shade_depth);
 
+
+
 // default lighting based on texture and material using the light we have just computed
 
  diffuse_term = diffuse_color* light_diffuse;
diff --git a/Shaders/water_lightfield.frag b/Shaders/water_lightfield.frag
index 3bced8c3c..80d6da2fd 100644
--- a/Shaders/water_lightfield.frag
+++ b/Shaders/water_lightfield.frag
@@ -59,6 +59,7 @@ uniform float ice_cover;
 uniform float sea_r;
 uniform float sea_g;
 uniform float sea_b;
+uniform float air_pollution;
 
 uniform int quality_level;
 uniform int ocean_flag;
@@ -75,6 +76,10 @@ const float EarthRadius = 5800000.0;
 
 //vec3 fog_Func(vec3 color, int type);
 float shadow_func (in float x, in float y, in float noise, in float dist);
+float fog_func (in float targ, in float alt);
+float rayleigh_in_func(in float dist, in float air_pollution, in float avisibility, in float eye_alt, in float vertex_alt);
+float alt_factor(in float eye_alt, in float vertex_alt);
+vec3 rayleigh_out_shift(in vec3 color, in float outscatter);
 //////////////////////
 
 /////// functions /////////
@@ -261,32 +266,8 @@ return e / pow((1.0 + a * exp(-b * (x-c)) ),(1.0/d));
 // physically this should be exp(-arg) but for technical reasons we use a sharper cutoff
 // for distance > visibility
 
-float fog_func (in float targ)
-{
 
 
-float fade_mix;
-
-// for large altitude > 30 km, we switch to some component of quadratic distance fading to
-// create the illusion of improved visibility range
-
-targ = 1.25 * targ; // need to sync with the distance to which terrain is drawn
-
-
-if (eye_alt < 30000.0)
-	{return exp(-targ - targ * targ * targ * targ);}
-else if (eye_alt < 50000.0)
-	{
-	fade_mix = (eye_alt - 30000.0)/20000.0;
-	return fade_mix * exp(-targ*targ - pow(targ,4.0)) + (1.0 - fade_mix) * exp(-targ - pow(targ,4.0));	
-	}
-else 
-	{
-	return exp(- targ * targ - pow(targ,4.0));
-	}
-
-}
-
 void main(void)
 	{
 
@@ -594,6 +575,11 @@ void main(void)
 	finalColor *= ambient_light;
 
 
+// Rayleigh color shift due to out-scattering
+    float rayleigh_length = 0.4 * avisibility * (2.5 - 1.9 * air_pollution)/alt_factor(eye_alt, eye_alt+relPos.z);
+    float outscatter = 1.0-exp(-dist/rayleigh_length);
+    finalColor.rgb = rayleigh_out_shift(finalColor.rgb,outscatter);
+
 
 // here comes the terrain haze model
 
@@ -707,7 +693,7 @@ else
 	}
 
 
-transmission =  fog_func(transmission_arg);
+transmission =  fog_func(transmission_arg, eye_alt);
 
 // there's always residual intensity, we should never be driven to zero
 if (eqColorFactor < 0.2) eqColorFactor = 0.2;
@@ -755,7 +741,11 @@ if (intensity > 0.0) // this needs to be a condition, because otherwise hazeColo
 	hazeColor = intensity * normalize(mix(hazeColor,  shadedFogColor, (1.0-smoothstep(0.5,0.9,eqColorFactor)))); 
 	}
 
-	
+	float rShade = 0.9 * smoothstep(terminator_width+ terminator, -terminator_width + terminator, yprime_alt-340000.0) + 0.1;
+	float lightIntensity = length(gl_Color.rgb)/1.73 * rShade;
+	vec3 rayleighColor = vec3 (0.17, 0.52, 0.87) * lightIntensity;
+	float rayleighStrength = rayleigh_in_func(dist, air_pollution, avisibility/max(lightIntensity,0.05), eye_alt, eye_alt + relPos.z);
+	finalColor.rgb = mix(finalColor.rgb, rayleighColor, rayleighStrength);
 
 	finalColor.rgb = mix(eqColorFactor * hazeColor * eShade, finalColor.rgb,transmission);