From bdd9520ca5ff46cd848d6ecb6fa5fbed83cfe92c Mon Sep 17 00:00:00 2001
From: Thorsten Renk <thorsten.i.renk@jyu.fi>
Date: Thu, 29 Aug 2013 15:15:51 +0300
Subject: [PATCH] Horizon blur and noise modulation model based on aloft
 visibility and weather variability for Atmospheric Light Scattering

---
 Effects/skydome.eff  |  42 +++--
 Shaders/skydome.frag | 388 ++++++++++++++++++++++++-------------------
 Shaders/skydome.vert | 233 ++++++++++++--------------
 3 files changed, 348 insertions(+), 315 deletions(-)

diff --git a/Effects/skydome.eff b/Effects/skydome.eff
index 5d1267ea0..f1d2fb625 100644
--- a/Effects/skydome.eff
+++ b/Effects/skydome.eff
@@ -4,16 +4,17 @@
   <parameters>
     <mie><use>/sim/rendering/mie</use></mie>
     <rayleigh><use>/sim/rendering/rayleigh</use></rayleigh>
-    <density><use>/sim/rendering/dome-density</use></density>
-    <overcast><use>/rendering/scene/overcast</use></overcast>
-    <saturation><use>/rendering/scene/saturation</use></saturation>
-    <scattering><use>/rendering/scene/scattering</use></scattering>
-    <visibility><use>/environment/ground-visibility-m</use></visibility>
-    <avisibility><use>/environment/visibility-m</use></avisibility>
-    <lthickness><use>/environment/ground-haze-thickness-m</use></lthickness>
-    <terminator><use>/environment/terminator-relative-position-m</use></terminator>
-    <terrain_alt><use>/environment/mean-terrain-elevation-m</use></terrain_alt>
+    <density><use>/sim/rendering/dome-density</use></density>
+    <overcast><use>/rendering/scene/overcast</use></overcast>
+    <saturation><use>/rendering/scene/saturation</use></saturation>
+    <scattering><use>/rendering/scene/scattering</use></scattering>
+    <visibility><use>/environment/ground-visibility-m</use></visibility>
+    <avisibility><use>/environment/visibility-m</use></avisibility>
+    <lthickness><use>/environment/ground-haze-thickness-m</use></lthickness>
+    <terminator><use>/environment/terminator-relative-position-m</use></terminator>
+    <terrain_alt><use>/environment/mean-terrain-elevation-m</use></terrain_alt>
     <cloud_self_shading><use>/environment/cloud-self-shading</use></cloud_self_shading>
+    <horizon_roughness><use>/local-weather/config/small-scale-persistence</use></horizon_roughness>
   </parameters>
   <technique n="8">
     <predicate>
@@ -56,52 +57,57 @@
         <name>density</name>
         <type>float</type>
         <value><use>density</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>overcast</name>
         <type>float</type>
         <value><use>overcast</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>saturation</name>
         <type>float</type>
         <value><use>saturation</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>scattering</name>
         <type>float</type>
         <value><use>scattering</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>visibility</name>
         <type>float</type>
         <value><use>visibility</use></value>
-      </uniform>
+      </uniform>
        <uniform>
         <name>hazeLayerAltitude</name>
         <type>float</type>
         <value><use>lthickness</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>terminator</name>
         <type>float</type>
         <value><use>terminator</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>avisibility</name>
         <type>float</type>
         <value><use>avisibility</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>terrain_alt</name>
         <type>float</type>
         <value><use>terrain_alt</use></value>
-      </uniform>
+      </uniform>
       <uniform>
         <name>cloud_self_shading</name>
         <type>float</type>
         <value><use>cloud_self_shading</use></value>
       </uniform>
+      <uniform>
+        <name>horizon_roughness</name>
+        <type>float</type>
+        <value><use>horizon_roughness</use></value>
+      </uniform>
     </pass>
   </technique>
 
diff --git a/Shaders/skydome.frag b/Shaders/skydome.frag
index a229890bd..5e7178548 100644
--- a/Shaders/skydome.frag
+++ b/Shaders/skydome.frag
@@ -2,27 +2,28 @@
  
 // Atmospheric scattering shader for flightgear
 // Written by Lauri Peltonen (Zan)
-// Implementation of O'Neil's algorithm
+// Implementation of O'Neil's algorithm
 // Ground haze layer added by Thorsten Renk
  
 varying vec3 rayleigh;
 varying vec3 mie;
-varying vec3 eye;
-varying vec3 hazeColor;
-varying float ct;
-//varying float cosphi;
-varying float delta_z;
-varying float alt;
+varying vec3 eye;
+varying vec3 hazeColor;
+varying float ct;
+varying float cphi;
+varying float delta_z;
+varying float alt;
 varying float earthShade;
- 
-uniform float overcast;
-uniform float saturation;
-uniform float visibility;
-uniform float avisibility;
-uniform float scattering;
-uniform float cloud_self_shading;
-
-const float EarthRadius = 5800000.0;
+ 
+uniform float overcast;
+uniform float saturation;
+uniform float visibility;
+uniform float avisibility;
+uniform float scattering;
+uniform float cloud_self_shading;
+uniform float horizon_roughness;
+
+const float EarthRadius = 5800000.0;
 
 float miePhase(in float cosTheta, in float g)
 {
@@ -41,172 +42,215 @@ float rayleighPhase(in float cosTheta)
   return 1.5 * (2.0 + 0.5*cosTheta*cosTheta);
 }
  
+float rand2D(in vec2 co){
+    return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453);
+}
+
+float simple_interpolate(in float a, in float b, in float x)
+{
+return a + smoothstep(0.0,1.0,x) * (b-a);
+}
+
+float interpolatedNoise2D(in float x, in float y)
+{
+      float integer_x    = x - fract(x);
+      float fractional_x = x - integer_x;
+
+      float integer_y    = y - fract(y);
+      float fractional_y = y - integer_y;
+
+      float v1 = rand2D(vec2(integer_x, integer_y));
+      float v2 = rand2D(vec2(integer_x+1.0, integer_y));
+      float v3 = rand2D(vec2(integer_x, integer_y+1.0));
+      float v4 = rand2D(vec2(integer_x+1.0, integer_y +1.0));
+
+      float i1 = simple_interpolate(v1 , v2 , fractional_x);
+      float i2 = simple_interpolate(v3 , v4 , fractional_x);
+
+      return simple_interpolate(i1 , i2 , fractional_y);
+}
  
+float Noise2D(in vec2 coord, in float wavelength)
+{
+return interpolatedNoise2D(coord.x/wavelength, coord.y/wavelength);
+
+}
  
 void main()
-{
-
+{
+
   vec3 shadedFogColor = vec3(0.65, 0.67, 0.78);
   float cosTheta = dot(normalize(eye), gl_LightSource[0].position.xyz);
- 
-  // position of the horizon line
-
-  float lAltitude = alt + delta_z;
-  float radiusEye = EarthRadius + alt;
-  float radiusLayer = EarthRadius + lAltitude;
-  float cthorizon;
-  float ctterrain;
-
-  if (radiusEye > radiusLayer) cthorizon = -sqrt(radiusEye * radiusEye - radiusLayer * radiusLayer)/radiusEye;
-  else cthorizon = sqrt(radiusLayer * radiusLayer - radiusEye * radiusEye)/radiusLayer;
-
-  ctterrain = -sqrt(radiusEye * radiusEye - EarthRadius * EarthRadius)/radiusEye;
+ 
+  // position of the horizon line
+
+  float lAltitude = alt + delta_z;
+  float radiusEye = EarthRadius + alt;
+  float radiusLayer = EarthRadius + lAltitude;
+  float cthorizon;
+  float ctterrain;
+
+  if (radiusEye > radiusLayer) cthorizon = -sqrt(radiusEye * radiusEye - radiusLayer * radiusLayer)/radiusEye;
+  else cthorizon = sqrt(radiusLayer * radiusLayer - radiusEye * radiusEye)/radiusLayer;
+
+  ctterrain = -sqrt(radiusEye * radiusEye - EarthRadius * EarthRadius)/radiusEye;
 
   vec3 color = rayleigh * rayleighPhase(cosTheta);
   color += mie * miePhase(cosTheta, -0.8);
-
-  vec3 black = vec3(0.0,0.0,0.0);
-
-  
-  float ovc = overcast;
-
-
-
-  float sat = 1.0 - ((1.0 - saturation) * 2.0);
-  if (sat < 0.3) sat = 0.3;
-
-
- // float wscale = 1.732;
-  
-// an overexposure filter, the log() seems to be pretty expensive though
-
-//  if (color.x > 0.8) color.x = 0.8 + 0.8* log(color.x/0.8);
-//  if (color.y > 0.8) color.y = 0.8 + 0.8* log(color.y/0.8);
-//  if (color.z > 0.8) color.z = 0.8 + 0.8* log(color.z/0.8);
-
-  
-// a different exposure filter  
-//color.x = 1.0 - exp(-1.3 * color.x);
-//color.y = 1.0 - exp(-1.3 * color.y);
-//color.z = 1.0 - exp(-1.3 * color.z);
-
-if (color.r > 0.58) color.r = 1.0 - exp(-1.5 * color.r);
-if (color.g > 0.58) color.g = 1.0 - exp(-1.5 * color.g);
-if (color.b > 0.58) color.b = 1.0 - exp(-1.5 * color.b);
-  
-// reduce the whiteout near the horizon generated by the single scattering approximation
-
-//if (ct > cthorizon) color = mix(color, black ,smoothstep(0.2+cthorizon, -0.2+cthorizon, ct));
-//else color = mix (color, black, smoothstep(0.2+cthorizon,-0.2+cthorizon,  cthorizon));
-
-
-
-
-// fog computations for a ground haze layer, extending from zero to lAltitude
-
-
-
-float transmission;
-float vAltitude;
-float delta_zv;
-
-float costheta = ct;
-
-float vis = min(visibility, avisibility);
-
-// hack - in an effect volume the visibility only may be reduced, so we take care here
-//if (avisibility < visibility){vis = avisibility;}
-
- if (delta_z > 0.0) // we're inside the layer
-	{
-  	if (costheta>0.0 + ctterrain) // looking up, view ray intersecting upper layer edge
-		{
-		transmission  = exp(-min((delta_z/max(costheta,0.1)),25000.0)/vis);
-		//transmission = 1.0;
-		vAltitude = min(vis * costheta, delta_z);
-  		delta_zv = delta_z - vAltitude;
-		}
-
-	else // looking down, view range intersecting terrain (which may not be drawn)
-		{
-		transmission = exp(alt/vis/costheta);
-		vAltitude = min(-vis * costheta, alt);
-  		delta_zv = delta_z + vAltitude;
-		}
-	}
-  else // we see the layer from above
-	{	
-	if (costheta < 0.0 + cthorizon) 
-		{
-		transmission = exp(-min(lAltitude/abs(costheta),25000.0)/vis);
-		transmission = transmission * exp(-alt/avisibility/abs(costheta));
-		transmission = 1.0 - (1.0 - transmission) * smoothstep(0+cthorizon, -0.02+cthorizon, costheta);
-   		vAltitude = min(lAltitude, -vis * costheta);
-		delta_zv = vAltitude; 
-		}
-	else
-		{	
-		transmission = 1.0;
-		delta_zv = 0.0;
-		}
-	}
-
-// combined intensity reduction by cloud shading and fog self-shading, corrected for Weber-Fechner perception law
-
-//float scattering = ground_scattering + (1.0 - ground_scattering) * smoothstep(avisibility, 1.5 * avisibility, -alt/costheta);
-
-float eqColorFactor = 1.0 - 0.1 * delta_zv/vis - (1.0 - min(scattering,cloud_self_shading));
-
-
-// there's always residual intensity, we should never be driven to zero
-if (eqColorFactor < 0.2) eqColorFactor = 0.2;
-
-
-// postprocessing of haze color
-vec3 hColor = hazeColor;
-
-
-// high altitude desaturation
-float intensity = length(hColor);
-hColor = intensity * normalize (mix(hColor, intensity * vec3 (1.0,1.0,1.0), 0.7* smoothstep(5000.0, 50000.0, alt)));
-
-// blue hue
-hColor.x = 0.83 * hColor.x;
-hColor.y = 0.9 * hColor.y;
-
-
-
-// further blueshift when in shadow, either cloud shadow, or self-shadow or Earth shadow, dependent on indirect 
-// light
-
-float fade_out = max(0.65 - 0.3 *overcast, 0.45);
-intensity = length(hColor);
-vec3 oColor = hColor;
-oColor = intensity * normalize(mix(oColor,  shadedFogColor, (smoothstep(0.1,1.0,ovc)))); 
-color = ovc *  mix(color, oColor * earthShade ,smoothstep(-0.1+ctterrain, 0.0+ctterrain, ct)) + (1-ovc) * color;
-hColor = intensity * normalize(mix(hColor,  1.5 * shadedFogColor, 1.0 -smoothstep(0.25, fade_out,earthShade) ));
-hColor = intensity * normalize(mix(hColor,  shadedFogColor, (1.0 - smoothstep(0.5,0.9,eqColorFactor)))); 
-//hColor = intensity * normalize(mix(hColor,  shadedFogColor, (1.0 - smoothstep(0.5,0.9,cloud_self_shading)) )); 
-hColor = hColor * earthShade;
-
-// accounting for overcast and saturation 
-
-
-
-color = sat * color + (1.0 - sat) * mix(color, black, smoothstep(0.4+cthorizon,0.2+cthorizon,ct));
-
-
-// the terrain below the horizon gets drawn in one optical thickness
-vec3 terrainHazeColor = eqColorFactor * hColor;	
-color = mix(color, terrainHazeColor ,smoothstep(0.01 + ctterrain, 0.0+ctterrain, ct));
-
-// mix fog the skydome with the right amount of haze
-
-color = transmission * color  + (1.0-transmission) * eqColorFactor * hColor;
-
+
+  vec3 black = vec3(0.0,0.0,0.0);
+
+  
+  float ovc = overcast;
+
+
+
+  float sat = 1.0 - ((1.0 - saturation) * 2.0);
+  if (sat < 0.3) sat = 0.3;
+
+
+ // float wscale = 1.732;
+  
+// an overexposure filter, the log() seems to be pretty expensive though
+
+//  if (color.x > 0.8) color.x = 0.8 + 0.8* log(color.x/0.8);
+//  if (color.y > 0.8) color.y = 0.8 + 0.8* log(color.y/0.8);
+//  if (color.z > 0.8) color.z = 0.8 + 0.8* log(color.z/0.8);
+
+  
+// a different exposure filter  
+//color.x = 1.0 - exp(-1.3 * color.x);
+//color.y = 1.0 - exp(-1.3 * color.y);
+//color.z = 1.0 - exp(-1.3 * color.z);
+
+if (color.r > 0.58) color.r = 1.0 - exp(-1.5 * color.r);
+if (color.g > 0.58) color.g = 1.0 - exp(-1.5 * color.g);
+if (color.b > 0.58) color.b = 1.0 - exp(-1.5 * color.b);
+  
+// reduce the whiteout near the horizon generated by the single scattering approximation
+
+//if (ct > cthorizon) color = mix(color, black ,smoothstep(0.2+cthorizon, -0.2+cthorizon, ct));
+//else color = mix (color, black, smoothstep(0.2+cthorizon,-0.2+cthorizon,  cthorizon));
+
+
+
+
+// fog computations for a ground haze layer, extending from zero to lAltitude
+
+
+
+float transmission;
+float vAltitude;
+float delta_zv;
+
+float costheta = ct;
+
+float vis = min(visibility, avisibility);
+
+// hack - in an effect volume the visibility only may be reduced, so we take care here
+//if (avisibility < visibility){vis = avisibility;}
+
+ if (delta_z > 0.0) // we're inside the layer
+	{
+  	if (costheta>0.0 + ctterrain) // looking up, view ray intersecting upper layer edge
+		{
+		transmission  = exp(-min((delta_z/max(costheta,0.1)),25000.0)/vis);
+		//transmission = 1.0;
+		vAltitude = min(vis * costheta, delta_z);
+  		delta_zv = delta_z - vAltitude;
+		}
+
+	else // looking down, view range intersecting terrain (which may not be drawn)
+		{
+		transmission = exp(alt/vis/costheta);
+		vAltitude = min(-vis * costheta, alt);
+  		delta_zv = delta_z + vAltitude;
+		}
+	}
+  else // we see the layer from above
+	{	
+	if (costheta < 0.0 + cthorizon) 
+		{
+		transmission = exp(-min(lAltitude/abs(costheta),25000.0)/vis);
+		transmission = transmission * exp(-alt/avisibility/abs(costheta));
+		transmission = 1.0 - (1.0 - transmission) * smoothstep(0+cthorizon, -0.02+cthorizon, costheta);
+   		vAltitude = min(lAltitude, -vis * costheta);
+		delta_zv = vAltitude; 
+		}
+	else
+		{	
+		transmission = 1.0;
+		delta_zv = 0.0;
+		}
+	}
+
+// combined intensity reduction by cloud shading and fog self-shading, corrected for Weber-Fechner perception law
+
+//float scattering = ground_scattering + (1.0 - ground_scattering) * smoothstep(avisibility, 1.5 * avisibility, -alt/costheta);
+
+float eqColorFactor = 1.0 - 0.1 * delta_zv/vis - (1.0 - min(scattering,cloud_self_shading));
+
+
+// there's always residual intensity, we should never be driven to zero
+if (eqColorFactor < 0.2) eqColorFactor = 0.2;
+
+
+// postprocessing of haze color
+vec3 hColor = hazeColor;
+
+
+// high altitude desaturation
+float intensity = length(hColor);
+hColor = intensity * normalize (mix(hColor, intensity * vec3 (1.0,1.0,1.0), 0.7* smoothstep(5000.0, 50000.0, alt)));
+
+// blue hue
+hColor.x = 0.83 * hColor.x;
+hColor.y = 0.9 * hColor.y;
+
+
+
+// further blueshift when in shadow, either cloud shadow, or self-shadow or Earth shadow, dependent on indirect 
+// light
+
+float fade_out = max(0.65 - 0.3 *overcast, 0.45);
+intensity = length(hColor);
+vec3 oColor = hColor;
+oColor = intensity * normalize(mix(oColor,  shadedFogColor, (smoothstep(0.1,1.0,ovc)))); 
+color = ovc *  mix(color, oColor * earthShade ,smoothstep(-0.1+ctterrain, 0.0+ctterrain, ct)) + (1-ovc) * color;
+hColor = intensity * normalize(mix(hColor,  1.5 * shadedFogColor, 1.0 -smoothstep(0.25, fade_out,earthShade) ));
+hColor = intensity * normalize(mix(hColor,  shadedFogColor, (1.0 - smoothstep(0.5,0.9,eqColorFactor)))); 
+//hColor = intensity * normalize(mix(hColor,  shadedFogColor, (1.0 - smoothstep(0.5,0.9,cloud_self_shading)) )); 
+hColor = hColor * earthShade;
+
+// accounting for overcast and saturation 
+
+
+
+color = sat * color + (1.0 - sat) * mix(color, black, smoothstep(0.4+cthorizon,0.2+cthorizon,ct));
+
+
+// the terrain below the horizon gets drawn in one optical thickness
+vec3 terrainHazeColor = eqColorFactor * hColor;	
+
+// determine a visibility-dependent angle for how smoothly the haze blends over the skydome
+
+float hazeBlendAngle = max(0.01,1000.0/avisibility + 0.3 * (1.0 - smoothstep(5000.0, 30000.0, avisibility)));
+float altFactor = smoothstep(-300.0, 0.0, delta_z);
+float altFactor2 =  0.2 + 0.8 * smoothstep(-3000.0, 0.0, delta_z);
+hazeBlendAngle = hazeBlendAngle + 0.1 * altFactor;
+hazeBlendAngle = hazeBlendAngle +  (1.0-horizon_roughness) * altFactor2 * 0.1 *  Noise2D(vec2(0.0,cphi), 0.3);
+
+
+color = mix(color, terrainHazeColor ,smoothstep(hazeBlendAngle + ctterrain, 0.0+ctterrain, ct));
+
+
+// mix fog the skydome with the right amount of haze
+
+color = transmission * color  + (1.0-transmission) * eqColorFactor * hColor;
+
 
   gl_FragColor = vec4(color, 1.0);
-  gl_FragDepth = 0.1;
+  gl_FragDepth = 0.1;
 
 }
 
diff --git a/Shaders/skydome.vert b/Shaders/skydome.vert
index 43673a5c8..5d083844c 100644
--- a/Shaders/skydome.vert
+++ b/Shaders/skydome.vert
@@ -6,29 +6,29 @@
  
  
 uniform mat4 osg_ViewMatrix;
-uniform mat4 osg_ViewMatrixInverse;
-uniform float hazeLayerAltitude;
-uniform float terminator;
-uniform float avisibility;
+uniform mat4 osg_ViewMatrixInverse;
+uniform float hazeLayerAltitude;
+uniform float terminator;
+uniform float avisibility;
 uniform float visibility;
-uniform float terrain_alt; 
+uniform float terrain_alt; 
 
 varying vec3 rayleigh;
 varying vec3 mie;
-varying vec3 eye;
-varying vec3 hazeColor;
-varying float ct;
-//varying float cosphi;
-varying float delta_z;
-varying float alt; 
-varying float earthShade;
+varying vec3 eye;
+varying vec3 hazeColor;
+varying float ct;
+varying float cphi;
+varying float delta_z;
+varying float alt; 
+varying float earthShade;
 
 // Dome parameters from FG and screen
 const float domeSize = 80000.0;
 const float realDomeSize = 100000.0;
 const float groundRadius = 0.984503332 * domeSize;
 const float altitudeScale = domeSize - groundRadius;
-const float EarthRadius = 5800000.0; 
+const float EarthRadius = 5800000.0; 
 
 // Dome parameters when calculating scattering
 // Assuming dome size is 5.0
@@ -43,25 +43,25 @@ const float fSamples = float(nSamples);
 uniform float rK = 0.0003; //0.00015;
 uniform float mK = 0.003; //0.0025;
 uniform float density = 0.5; //1.0
-//vec3 rayleighK = rK * vec3(5.602, 7.222, 19.644);
+//vec3 rayleighK = rK * vec3(5.602, 7.222, 19.644);
 vec3 rayleighK = rK * vec3(4.5, 8.62, 17.3);
 vec3 mieK = vec3(mK);
 vec3 sunIntensity = 10.0*vec3(120.0, 125.0, 130.0);
- 
-// light_func is a generalized logistic function fit to the light intensity as a function
-// of scaled terminator position obtained from Flightgear core
-
-float light_func (in float x, in float a, in float b, in float c, in float d, in float e)
-{
-x = x - 0.5;
-
-// use the asymptotics to shorten computations
-if (x > 30.0) {return e;}
-if (x < -15.0) {return 0.0;}
-
-return e / pow((1.0 + a * exp(-b * (x-c)) ),(1.0/d));
-}
-
+ 
+// light_func is a generalized logistic function fit to the light intensity as a function
+// of scaled terminator position obtained from Flightgear core
+
+float light_func (in float x, in float a, in float b, in float c, in float d, in float e)
+{
+x = x - 0.5;
+
+// use the asymptotics to shorten computations
+if (x > 30.0) {return e;}
+if (x < -15.0) {return 0.0;}
+
+return e / pow((1.0 + a * exp(-b * (x-c)) ),(1.0/d));
+}
+
 
 // Find intersections of ray to skydome
 // ray must be normalized
@@ -76,8 +76,8 @@ float intersection (in float cheight, in vec3 ray, in float rad2)
  
 // Return the scale function at height = 0 for different thetas
 float outscatterscale(in float costheta)
-{
-
+{
+
 
   float x = 1.0 - costheta;
  
@@ -115,11 +115,11 @@ void main()
     // Make it so that 0.0 is ground level and 1.0 is 100km (space) level
     float altitude = distance(groundPoint, vec4(0.0, 0.0, 0.0, 1.0));
     float scaledAltitude = altitude / realDomeSize;
-
-    // the local horizon angle
-    float radiusEye = EarthRadius + altitude;
-    float ctterrain = -sqrt(radiusEye * radiusEye - EarthRadius * EarthRadius)/radiusEye; 
-
+
+    // the local horizon angle
+    float radiusEye = EarthRadius + altitude;
+    float ctterrain = -sqrt(radiusEye * radiusEye - EarthRadius * EarthRadius)/radiusEye; 
+
 
     // Camera's position, z is up!
     float cameraRealAltitude = groundLevel + heightScale*scaledAltitude;
@@ -133,7 +133,7 @@ void main()
       // We are in space, calculate correct positiondelta!
       relativePosition -= space * normalize(relativePosition);
     }
- 
+ 
 
     vec3 positionDelta = relativePosition / fSamples;
     float deltaLength = length(positionDelta); // Should multiply by something?
@@ -147,20 +147,20 @@ void main()
     // If sample is above camera, reverse ray direction
     if(positionDelta.z < 0.0) cameraCosTheta = -positionDelta.z / deltaLength;
     else cameraCosTheta = positionDelta.z / deltaLength;
- 
-    
-    float cameraCosTheta1 = -positionDelta.z / deltaLength;
-
+ 
+    
+    float cameraCosTheta1 = -positionDelta.z / deltaLength;
+
 
     // Total attenuation from camera to skydome
     float totalCameraScatter = outscatter(cameraCosTheta, scaledAltitude);
  
  
     // Do numerical integration of scattering function from skydome to camera
-    vec3 color = vec3(0.0);
-
-    // no scattering integrations where terrain is later drawn
-    if (cameraCosTheta1 > (ctterrain-0.05))
+    vec3 color = vec3(0.0);
+
+    // no scattering integrations where terrain is later drawn
+    if (cameraCosTheta1 > (ctterrain-0.05))
     {
     for(int i = 0; i < nSamples; i++) 
     {
@@ -179,7 +179,7 @@ void main()
       // Again, reverse the direction if vertex is over the camera
       float cameraScatter;
       if(relativePosition.z < 0.0) {  // Vertex is over the camera
-        cameraCosTheta = -dot(normalize(positionDelta), normalize(sample));
+        cameraCosTheta = -dot(normalize(positionDelta), normalize(sample));
 
         cameraScatter = totalCameraScatter - outscatter(cameraCosTheta, sampleAltitude);
       } else {  // Vertex is below camera
@@ -196,13 +196,13 @@ void main()
       sample += positionDelta;
     }
     }
-    color *= sunIntensity;
+    color *= sunIntensity;
     ct = cameraCosTheta1;
     rayleigh = rayleighK * color;
     mie = mieK * color;
     eye = gl_NormalMatrix * positionDelta;
- 
-
+ 
+
    // We need to move the camera so that the dome appears to be centered around earth
     // to make the dome render correctly!
     float moveDown = -altitude; // Center dome on camera
@@ -210,82 +210,65 @@ void main()
     moveDown += scaledAltitude * altitudeScale; // And move correctly according to altitude
  
     // Vertex transformed correctly so that at 100km we are at space border
-    vec4 finalVertex = realVertex - vec4(0.0, 0.0, 1.0, 0.0) * moveDown;
-
-    // prepare some stuff for a ground haze layer
-  
-    delta_z = hazeLayerAltitude - altitude;
-    alt = altitude;
-    
-
-    // establish coordinates relative to sun position
-    vec4 ep = gl_ModelViewMatrixInverse * vec4(0.0,0.0,0.0,1.0);
-    vec3 lightFull = (gl_ModelViewMatrixInverse * gl_LightSource[0].position).xyz;
-    vec3 lightHorizon = normalize(vec3(lightFull.x,lightFull.y, 0.0) );
- 
-
-    vec3 relVector = normalize(finalVertex.xyz - ep.xyz);
-    
-    // and compute the twilight shading
-    
-
-    // yprime is the coordinate from/towards terminator
-    float yprime;
-	
-    if (alt > hazeLayerAltitude) // we're looking from above and can see far
-	{
-    	if (ct < 0.0)
-    		{
-		yprime = -dot(relVector,lightHorizon) * altitude/-ct;//(ct-0.001);
-		yprime = yprime -sqrt(2.0 * EarthRadius * hazeLayerAltitude);
-		}
-   	 else  // the only haze we see looking up is overcast, assume its altitude
-		{
-		yprime = -dot(relVector,lightHorizon) * avisibility; 
-		yprime = yprime -sqrt(2.0 * EarthRadius * 10000.0);
-		}
-	}
-     else 
-	{yprime = -dot(relVector,lightHorizon) * avisibility;
-	yprime = yprime -sqrt(2.0 * EarthRadius * hazeLayerAltitude);
-	}
-
-    if (terminator > 1000000.0){yprime = -sqrt(2.0 * EarthRadius * hazeLayerAltitude);}
-
-    //float edgeAlt = max(hazeLayerAltitude - (alt-terrain_alt)/avisibility * visibility, terrain_alt);
-    
-    //yprime = yprime -sqrt(2.0 * EarthRadius * edgeAlt);
-
-    float terminator_width = 200000.0;
-    earthShade = 0.9 * smoothstep((terminator_width+ terminator), (-terminator_width + terminator), yprime) + 0.1;
-
-
-//hazeColor = vec3 (gl_LightSource[0].diffuse.x, gl_LightSource[0].diffuse.y, gl_LightSource[0].diffuse.z);
-
-    //hazeColor.x = hazeColor.x * 0.83;
-    //hazeColor.y = hazeColor.y * 0.9;  
-
-     float lightArg = (terminator-yprime)/100000.0;
-     vec4 light_diffuse;
-     light_diffuse.b = light_func(lightArg, 1.330e-05, 0.264, 2.527, 1.08e-05, 1.0);
-     light_diffuse.g = light_func(lightArg, 3.931e-06, 0.264, 3.827, 7.93e-06, 1.0);
-     light_diffuse.r = light_func(lightArg, 8.305e-06, 0.161, 3.827, 3.04e-05, 1.0);
-     light_diffuse.a = 0.0;
-     hazeColor = light_diffuse.xyz;
-
-     float intensity = length(hazeColor.xyz);
-     float mie_magnitude = 0.5 * smoothstep(350000.0, 150000.0, terminator -sqrt(2.0 * EarthRadius * terrain_alt)); 
-     float mie_angle = (0.5 *  dot(normalize(relVector), normalize(lightFull)) ) + 0.5;
-     hazeColor = intensity * ((1.0 - mie_magnitude) + mie_magnitude * mie_angle) * normalize(mix(hazeColor,  vec3 (0.5, 0.58, 0.65), mie_magnitude * (0.5 - 0.5 * mie_angle)) ); 
-
-
-// high altitude desaturation - would be best here this causes a box-like bug for some reason
-// so it moved to the fragment shader where it has no issues
-
-//float intensity = length(hazeColor.xyz);
-//hazeColor = intensity * normalize (mix(hazeColor, intensity * vec3 (1.0,1.0,1.0), 0.8* smoothstep(5000.0, 50000.0, alt)));
-
-    
+    vec4 finalVertex = realVertex - vec4(0.0, 0.0, 1.0, 0.0) * moveDown;
+
+    // prepare some stuff for a ground haze layer
+  
+    delta_z = hazeLayerAltitude - altitude;
+    alt = altitude;
+    
+
+    // establish coordinates relative to sun position
+    vec4 ep = gl_ModelViewMatrixInverse * vec4(0.0,0.0,0.0,1.0);
+    vec3 lightFull = (gl_ModelViewMatrixInverse * gl_LightSource[0].position).xyz;
+    vec3 lightHorizon = normalize(vec3(lightFull.x,lightFull.y, 0.0) );
+ 
+
+    vec3 relVector = normalize(finalVertex.xyz - ep.xyz);
+    
+    // and compute the twilight shading
+    
+
+    // yprime is the coordinate from/towards terminator
+    float yprime;
+	
+    if (alt > hazeLayerAltitude) // we're looking from above and can see far
+	{
+    	if (ct < 0.0)
+    		{
+		yprime = -dot(relVector,lightHorizon) * altitude/-ct;//(ct-0.001);
+		yprime = yprime -sqrt(2.0 * EarthRadius * hazeLayerAltitude);
+		}
+   	 else  // the only haze we see looking up is overcast, assume its altitude
+		{
+		yprime = -dot(relVector,lightHorizon) * avisibility; 
+		yprime = yprime -sqrt(2.0 * EarthRadius * 10000.0);
+		}
+	}
+     else 
+	{yprime = -dot(relVector,lightHorizon) * avisibility;
+	yprime = yprime -sqrt(2.0 * EarthRadius * hazeLayerAltitude);
+	}
+
+    if (terminator > 1000000.0){yprime = -sqrt(2.0 * EarthRadius * hazeLayerAltitude);}
+
+    float terminator_width = 200000.0;
+    earthShade = 0.9 * smoothstep((terminator_width+ terminator), (-terminator_width + terminator), yprime) + 0.1;
+
+     float lightArg = (terminator-yprime)/100000.0;
+     vec4 light_diffuse;
+     light_diffuse.b = light_func(lightArg, 1.330e-05, 0.264, 2.527, 1.08e-05, 1.0);
+     light_diffuse.g = light_func(lightArg, 3.931e-06, 0.264, 3.827, 7.93e-06, 1.0);
+     light_diffuse.r = light_func(lightArg, 8.305e-06, 0.161, 3.827, 3.04e-05, 1.0);
+     light_diffuse.a = 0.0;
+     hazeColor = light_diffuse.xyz;
+
+     float intensity = length(hazeColor.xyz);
+     float mie_magnitude = 0.5 * smoothstep(350000.0, 150000.0, terminator -sqrt(2.0 * EarthRadius * terrain_alt)); 
+     cphi = dot(normalize(relVector), normalize(lightHorizon));
+     float mie_angle = (0.5 *  dot(normalize(relVector), normalize(lightFull)) ) + 0.5;
+     hazeColor = intensity * ((1.0 - mie_magnitude) + mie_magnitude * mie_angle) * normalize(mix(hazeColor,  vec3 (0.5, 0.58, 0.65), mie_magnitude * (0.5 - 0.5 * mie_angle)) ); 
+
 
     // Transform
     gl_Position = gl_ModelViewProjectionMatrix * finalVertex;