Shaders/water_lightfield.frag

   1 // This shader is mostly an adaptation of the shader found at
   2 //  http://www.bonzaisoftware.com/water_tut.html and its glsl conversion
   3 //  available at http://forum.bonzaisoftware.com/viewthread.php?tid=10
   4 //  © Michael Horsch - 2005
   5 //  Major update and revisions - 2011-10-07
   6 //  © Emilian Huminiuc and Vivian Meazza
   7 // ported to lightfield shading Thorsten Renk 2012
   8
   9 #version 120
  10
  11 uniform sampler2D water_normalmap;
  12 uniform sampler2D water_dudvmap;
  13 uniform sampler2D sea_foam;
  14 uniform sampler2D perlin_normalmap;
  15 uniform sampler2D ice_texture;
  16
  17 uniform sampler3D Noise;
  18
  19 uniform float saturation, Overcast, WindE, WindN;
  20 uniform float osg_SimulationTime;
  21
  22 varying vec4 waterTex1; //moving texcoords
  23 varying vec4 waterTex2; //moving texcoords
  24 varying vec4 waterTex4; //viewts
  25 varying vec3 viewerdir;
  26 varying vec3 lightdir;
  27 varying vec3 relPos;
  28 varying vec3 rawPos;
  29
  30 varying float earthShade;
  31 varying float yprime_alt;
  32 varying float mie_angle;
  33
  34 uniform    float WaveFreq ;
  35 uniform    float WaveAmp ;
  36 uniform    float WaveSharp ;
  37 uniform    float WaveAngle ;
  38 uniform    float WaveFactor ;
  39 uniform    float WaveDAngle ;
  40 uniform    float normalmap_dds;
  41
  42
  43 uniform float hazeLayerAltitude;
  44 uniform float terminator;
  45 uniform float terrain_alt;
  46 uniform float avisibility;
  47 uniform float visibility;
  48 uniform float overcast;
  49 uniform float scattering;
  50 uniform float ground_scattering;
  51 uniform float cloud_self_shading;
  52 uniform float eye_alt;
  53 uniform float ice_cover;
  54 uniform float sea_r;
  55 uniform float sea_g;
  56 uniform float sea_b;
  57
  58
  59 vec3 specular_light;
  60
  61 //uniform int wquality_level;
  62
  63 const float terminator_width = 200000.0;
  64 const float EarthRadius = 5800000.0;
  65 ////fog "include" /////
  66 //uniform int fogType;
  67
  68 vec3 fog_Func(vec3 color, int type);
  69 //////////////////////
  70
  71 /////// functions /////////
  72
  73 float rand2D(in vec2 co){
  74     return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453);
  75 }
  76
  77 float rand3D(in vec3 co){
  78     return fract(sin(dot(co.xyz ,vec3(12.9898,78.233,144.7272))) * 43758.5453);
  79 }
  80
  81 float cosine_interpolate(in float a, in float b, in float x)
  82 {
  83         float ft = x * 3.1415927;
  84         float f = (1.0 - cos(ft)) * .5;
  85
  86         return  a*(1.0-f) + b*f;
  87 }
  88
  89 float simple_interpolate(in float a, in float b, in float x)
  90 {
  91 return a + smoothstep(0.0,1.0,x) * (b-a);
  92 }
  93
  94 float interpolatedNoise2D(in float x, in float y)
  95 {
  96       float integer_x    = x - fract(x);
  97       float fractional_x = x - integer_x;
  98
  99       float integer_y    = y - fract(y);
 100       float fractional_y = y - integer_y;
 101
 102       float v1 = rand2D(vec2(integer_x, integer_y));
 103       float v2 = rand2D(vec2(integer_x+1.0, integer_y));
 104       float v3 = rand2D(vec2(integer_x, integer_y+1.0));
 105       float v4 = rand2D(vec2(integer_x+1.0, integer_y +1.0));
 106
 107       float i1 = simple_interpolate(v1 , v2 , fractional_x);
 108       float i2 = simple_interpolate(v3 , v4 , fractional_x);
 109
 110       return simple_interpolate(i1 , i2 , fractional_y);
 111 }
 112
 113 float interpolatedNoise3D(in float x, in float y, in float z)
 114 {
 115       float integer_x    = x - fract(x);
 116       float fractional_x = x - integer_x;
 117
 118       float integer_y    = y - fract(y);
 119       float fractional_y = y - integer_y;
 120
 121       float integer_z    = z - fract(z);
 122       float fractional_z = z - integer_z;
 123
 124       float v1 = rand3D(vec3(integer_x, integer_y, integer_z));
 125       float v2 = rand3D(vec3(integer_x+1.0, integer_y, integer_z));
 126       float v3 = rand3D(vec3(integer_x, integer_y+1.0, integer_z));
 127       float v4 = rand3D(vec3(integer_x+1.0, integer_y +1.0, integer_z));
 128
 129       float v5 = rand3D(vec3(integer_x, integer_y, integer_z+1.0));
 130       float v6 = rand3D(vec3(integer_x+1.0, integer_y, integer_z+1.0));
 131       float v7 = rand3D(vec3(integer_x, integer_y+1.0, integer_z+1.0));
 132       float v8 = rand3D(vec3(integer_x+1.0, integer_y +1.0, integer_z+1.0));
 133
 134
 135       float i1 = simple_interpolate(v1,v5, fractional_z);
 136       float i2 = simple_interpolate(v2,v6, fractional_z);
 137       float i3 = simple_interpolate(v3,v7, fractional_z);
 138       float i4 = simple_interpolate(v4,v8, fractional_z);
 139
 140       float ii1 = simple_interpolate(i1,i2,fractional_x);
 141       float ii2 = simple_interpolate(i3,i4,fractional_x);
 142
 143
 144       return simple_interpolate(ii1 , ii2 , fractional_y);
 145 }
 146
 147 float Noise2D(in vec2 coord, in float wavelength)
 148 {
 149 return interpolatedNoise2D(coord.x/wavelength, coord.y/wavelength);
 150
 151 }
 152
 153 float Noise3D(in vec3 coord, in float wavelength)
 154 {
 155 return interpolatedNoise3D(coord.x/wavelength, coord.y/wavelength, coord.z/wavelength);
 156 }
 157
 158
 159
 160 void rotationmatrix(in float angle, out mat4 rotmat)
 161         {
 162         rotmat = mat4( cos( angle ), -sin( angle ), 0.0, 0.0,
 163                 sin( angle ),  cos( angle ), 0.0, 0.0,
 164                 0.0         ,  0.0         , 1.0, 0.0,
 165                 0.0         ,  0.0         , 0.0, 1.0 );
 166         }
 167
 168 // wave functions ///////////////////////
 169
 170 struct Wave {
 171         float freq;  // 2*PI / wavelength
 172         float amp;   // amplitude
 173         float phase; // speed * 2*PI / wavelength
 174         vec2 dir;
 175         };
 176
 177 Wave wave0 = Wave(1.0, 1.0, 0.5, vec2(0.97, 0.25));
 178 Wave wave1 = Wave(2.0, 0.5, 1.3, vec2(0.97, -0.25));
 179 Wave wave2 = Wave(1.0, 1.0, 0.6, vec2(0.95, -0.3));
 180 Wave wave3 = Wave(2.0, 0.5, 1.4, vec2(0.99, 0.1));
 181
 182
 183
 184
 185 float evaluateWave(in Wave w, vec2 pos, float t)
 186         {
 187         return w.amp * sin( dot(w.dir, pos) * w.freq + t * w.phase);
 188         }
 189
 190 // derivative of wave function
 191 float evaluateWaveDeriv(Wave w, vec2 pos, float t)
 192         {
 193         return w.freq * w.amp * cos( dot(w.dir, pos)*w.freq + t*w.phase);
 194         }
 195
 196 // sharp wave functions
 197 float evaluateWaveSharp(Wave w, vec2 pos, float t, float k)
 198         {
 199         return w.amp * pow(sin( dot(w.dir, pos)*w.freq + t*w.phase)* 0.5 + 0.5 , k);
 200         }
 201
 202 float evaluateWaveDerivSharp(Wave w, vec2 pos, float t, float k)
 203         {
 204         return k*w.freq*w.amp * pow(sin( dot(w.dir, pos)*w.freq + t*w.phase)* 0.5 + 0.5 , k - 1) * cos( dot(w.dir, pos)*w.freq + t*w.phase);
 205         }
 206
 207 void sumWaves(float angle, float dangle, float windScale, float factor, out float ddx, float ddy)
 208         {
 209         mat4 RotationMatrix;
 210         float deriv;
 211         vec4 P = waterTex1 * 1024;
 212
 213         rotationmatrix(radians(angle + dangle * windScale + 0.6 * sin(P.x * factor)), RotationMatrix);
 214         P *= RotationMatrix;
 215
 216         P.y += evaluateWave(wave0, P.xz, osg_SimulationTime);
 217         deriv = evaluateWaveDeriv(wave0, P.xz, osg_SimulationTime );
 218         ddx = deriv * wave0.dir.x;
 219         ddy = deriv * wave0.dir.y;
 220
 221         P.y += evaluateWave(wave1, P.xz, osg_SimulationTime);
 222         deriv = evaluateWaveDeriv(wave1, P.xz, osg_SimulationTime);
 223         ddx += deriv * wave1.dir.x;
 224         ddy += deriv * wave1.dir.y;
 225
 226         P.y += evaluateWaveSharp(wave2, P.xz, osg_SimulationTime, WaveSharp);
 227         deriv = evaluateWaveDerivSharp(wave2, P.xz, osg_SimulationTime, WaveSharp);
 228         ddx += deriv * wave2.dir.x;
 229         ddy += deriv * wave2.dir.y;
 230
 231         P.y += evaluateWaveSharp(wave3, P.xz, osg_SimulationTime, WaveSharp);
 232         deriv = evaluateWaveDerivSharp(wave3, P.xz, osg_SimulationTime, WaveSharp);
 233         ddx += deriv * wave3.dir.x;
 234         ddy += deriv * wave3.dir.y;
 235         }
 236
 237
 238 float light_func (in float x, in float a, in float b, in float c, in float d, in float e)
 239 {
 240 x = x - 0.5;
 241
 242 // use the asymptotics to shorten computations
 243 if (x > 30.0) {return e;}
 244 if (x < -15.0) {return 0.0;}
 245
 246 return e / pow((1.0 + a * exp(-b * (x-c)) ),(1.0/d));
 247 }
 248
 249 // this determines how light is attenuated in the distance
 250 // physically this should be exp(-arg) but for technical reasons we use a sharper cutoff
 251 // for distance > visibility
 252
 253 float fog_func (in float targ)
 254 {
 255
 256
 257 float fade_mix;
 258
 259 // for large altitude > 30 km, we switch to some component of quadratic distance fading to
 260 // create the illusion of improved visibility range
 261
 262 targ = 1.25 * targ; // need to sync with the distance to which terrain is drawn
 263
 264
 265 if (eye_alt < 30000.0)
 266         {return exp(-targ - targ * targ * targ * targ);}
 267 else if (eye_alt < 50000.0)
 268         {
 269         fade_mix = (eye_alt - 30000.0)/20000.0;
 270         return fade_mix * exp(-targ*targ - pow(targ,4.0)) + (1.0 - fade_mix) * exp(-targ - pow(targ,4.0));
 271         }
 272 else
 273         {
 274         return exp(- targ * targ - pow(targ,4.0));
 275         }
 276
 277 }
 278
 279 void main(void)
 280         {
 281
 282
 283         vec3 shadedFogColor = vec3(0.65, 0.67, 0.78);
 284         float effective_scattering = min(scattering, cloud_self_shading);
 285
 286         float dist = length(relPos);
 287         const vec4 sca = vec4(0.005, 0.005, 0.005, 0.005);
 288         const vec4 sca2 = vec4(0.02, 0.02, 0.02, 0.02);
 289         const vec4 tscale = vec4(0.25, 0.25, 0.25, 0.25);
 290
 291         float noise_50m = Noise3D(rawPos.xyz, 50.0);
 292         float noise_250m = Noise3D(rawPos.xyz,250.0);
 293         float noise_1500m = Noise3D(rawPos.xyz,1500.0);
 294         float noise_2000m = Noise3D(rawPos.xyz,2000.0);
 295         float noise_2500m = Noise3D(rawPos.xyz, 2500.0);
 296
 297         mat4 RotationMatrix;
 298
 299         // compute direction to viewer
 300         vec3 E = normalize(viewerdir);
 301
 302         // compute direction to light source
 303         vec3 L = lightdir; // normalize(lightdir);
 304
 305         // half vector
 306         vec3 Hv = normalize(L + E);
 307
 308         //vec3 Normal = normalize(normal);
 309         vec3 Normal = vec3 (0.0, 0.0, 1.0);
 310
 311         const float water_shininess = 240.0;
 312
 313         // approximate cloud cover
 314         //float cover = 0.0;
 315         //bool Status = true;
 316
 317         float windEffect = sqrt( WindE*WindE + WindN*WindN ) * 0.6;                             //wind speed in kt
 318         float windScale =  15.0/(3.0 + windEffect);                                                                                                     //wave scale
 319         float windEffect_low = 0.3 + 0.7 * smoothstep(0.0, 5.0, windEffect);                                    //low windspeed wave filter
 320         float waveRoughness = 0.01 + smoothstep(0.0, 40.0, windEffect);                                         //wave roughness filter
 321
 322         float mixFactor = 0.2 + 0.02 * smoothstep(0.0, 50.0, windEffect);
 323         //mixFactor = 0.2;
 324         mixFactor = clamp(mixFactor, 0.3, 0.8);
 325
 326         // there's no need to do wave patterns or foam for pixels which are so far away that we can't actually see them
 327         // we only need detail in the near zone or where the sun reflection is
 328
 329         int detail_flag;
 330         if ((dist > 15000.0) && (dot(normalize(vec3 (lightdir.x, lightdir.y, 0.0) ), normalize(relPos)) < 0.7 ))  {detail_flag = 0;}
 331         else {detail_flag = 1;}
 332
 333         //detail_flag = 1;
 334
 335         // sine waves
 336         float ddx, ddx1, ddx2, ddx3, ddy, ddy1, ddy2, ddy3;
 337         float angle;
 338
 339         ddx = 0.0, ddy = 0.0;
 340         ddx1 = 0.0, ddy1 = 0.0;
 341         ddx2 = 0.0, ddy2 = 0.0;
 342         ddx3 = 0.0, ddy3 = 0.0;
 343
 344         if (detail_flag == 1)
 345         {
 346         angle = 0.0;
 347
 348         wave0.freq = WaveFreq ;
 349         wave0.amp = WaveAmp;
 350         wave0.dir =  vec2 (0.0, 1.0); //vec2(cos(radians(angle)), sin(radians(angle)));
 351
 352         angle -= 45;
 353         wave1.freq = WaveFreq * 2.0 ;
 354         wave1.amp = WaveAmp * 1.25;
 355         wave1.dir =  vec2(0.70710, -0.7071); //vec2(cos(radians(angle)), sin(radians(angle)));
 356
 357         angle += 30;
 358         wave2.freq = WaveFreq * 3.5;
 359         wave2.amp = WaveAmp * 0.75;
 360         wave2.dir =  vec2(0.96592, -0.2588);// vec2(cos(radians(angle)), sin(radians(angle)));
 361
 362         angle -= 50;
 363         wave3.freq = WaveFreq * 3.0 ;
 364         wave3.amp = WaveAmp * 0.75;
 365         wave3.dir =  vec2(0.42261, -0.9063); //vec2(cos(radians(angle)), sin(radians(angle)));
 366
 367         // sum waves
 368
 369
 370         sumWaves(WaveAngle, -1.5, windScale, WaveFactor, ddx, ddy);
 371         sumWaves(WaveAngle, 1.5, windScale, WaveFactor, ddx1, ddy1);
 372
 373         //reset the waves
 374         angle = 0.0;
 375         float waveamp = WaveAmp * 0.75;
 376
 377         wave0.freq = WaveFreq ;
 378         wave0.amp = waveamp;
 379         wave0.dir =  vec2 (0.0, 1.0); //vec2(cos(radians(angle)), sin(radians(angle)));
 380
 381         angle -= 20;
 382         wave1.freq = WaveFreq * 2.0 ;
 383         wave1.amp = waveamp * 1.25;
 384         wave1.dir =  vec2(0.93969, -0.34202);// vec2(cos(radians(angle)), sin(radians(angle)));
 385
 386         angle += 35;
 387         wave2.freq = WaveFreq * 3.5;
 388         wave2.amp = waveamp * 0.75;
 389         wave2.dir =  vec2(0.965925, 0.25881);  //vec2(cos(radians(angle)), sin(radians(angle)));
 390
 391         angle -= 45;
 392         wave3.freq = WaveFreq * 3.0 ;
 393         wave3.amp = waveamp * 0.75;
 394         wave3.dir =  vec2(0.866025, -0.5); //vec2(cos(radians(angle)), sin(radians(angle)));
 395
 396
 397         sumWaves(WaveAngle + WaveDAngle, -1.5, windScale, WaveFactor, ddx2, ddy2);
 398         sumWaves(WaveAngle + WaveDAngle, 1.5, windScale, WaveFactor, ddx3, ddy3);
 399
 400         }
 401         // end sine stuff
 402
 403         //cover = 5.0 * smoothstep(0.6, 1.0, scattering);
 404         //cover = 5.0 * ground_scattering;
 405
 406         vec4 viewt = normalize(waterTex4);
 407         vec4 disdis = texture2D(water_dudvmap, vec2(waterTex2 * tscale)* windScale) * 2.0 - 1.0;
 408         vec4 vNorm;
 409
 410
 411         //normalmaps
 412         vec4 nmap   = texture2D(water_normalmap, vec2(waterTex1 + disdis * sca2) * windScale) * 2.0 - 1.0;
 413         vec4 nmap1  = texture2D(perlin_normalmap, vec2(waterTex1 + disdis * sca2) * windScale) * 2.0 - 1.0;
 414
 415         rotationmatrix(radians(3.0 * sin(osg_SimulationTime * 0.0075)), RotationMatrix);
 416         nmap  += texture2D(water_normalmap, vec2(waterTex2 * RotationMatrix * tscale) * windScale) * 2.0 - 1.0;
 417         nmap1 += texture2D(perlin_normalmap, vec2(waterTex2 * RotationMatrix * tscale) * windScale) * 2.0 - 1.0;
 418
 419         nmap  *= windEffect_low;
 420         nmap1 *= windEffect_low;
 421
 422         // mix water and noise, modulated by factor
 423         vNorm = normalize(mix(nmap, nmap1, mixFactor) * waveRoughness);
 424         vNorm.r += ddx + ddx1 + ddx2 + ddx3;
 425
 426
 427         if (normalmap_dds > 0)
 428                 {vNorm = -vNorm;}               //dds fix
 429
 430         vNorm = vNorm * (0.5 + 0.5 * noise_250m);
 431
 432         //load reflection
 433
 434         vec4 refl ;
 435
 436         refl.r = sea_r;
 437         refl.g = sea_g;
 438         refl.b = sea_b;
 439         refl.a = 1.0;
 440
 441         refl.g = refl.g * (0.9 + 0.2* noise_2500m);
 442
 443         float intensity;
 444         // de-saturate for reduced light
 445         refl.rgb = mix(refl.rgb,  vec3 (0.248, 0.248, 0.248), 1.0 - smoothstep(0.1, 0.8, ground_scattering));
 446
 447         // de-saturate light for overcast haze
 448         intensity = length(refl.rgb);
 449         refl.rgb = mix(refl.rgb,  intensity * vec3 (1.0, 1.0, 1.0), 0.5 * smoothstep(0.1, 0.9, overcast));
 450
 451         vec3 N;
 452
 453
 454
 455
 456         vec3 N0 = vec3(texture2D(water_normalmap, vec2(waterTex1 + disdis * sca2) * windScale) * 2.0 - 1.0);
 457         vec3 N1 = vec3(texture2D(perlin_normalmap, vec2(waterTex1 + disdis * sca) * windScale) * 2.0 - 1.0);
 458
 459         N0 += vec3(texture2D(water_normalmap, vec2(waterTex1 * tscale) * windScale) * 2.0 - 1.0);
 460         N1 += vec3(texture2D(perlin_normalmap, vec2(waterTex2 * tscale) * windScale) * 2.0 - 1.0);
 461
 462
 463
 464         rotationmatrix(radians(2.0 * sin(osg_SimulationTime * 0.005)), RotationMatrix);
 465         N0 += vec3(texture2D(water_normalmap, vec2(waterTex2 * RotationMatrix * (tscale + sca2)) * windScale) * 2.0 - 1.0);
 466         N1 += vec3(texture2D(perlin_normalmap, vec2(waterTex2 * RotationMatrix * (tscale + sca2)) * windScale) * 2.0 - 1.0);
 467
 468         rotationmatrix(radians(-4.0 * sin(osg_SimulationTime * 0.003)), RotationMatrix);
 469         N0 += vec3(texture2D(water_normalmap, vec2(waterTex1 * RotationMatrix + disdis * sca2) * windScale) * 2.0 - 1.0);
 470         N1 += vec3(texture2D(perlin_normalmap, vec2(waterTex1 * RotationMatrix + disdis * sca) * windScale) * 2.0 - 1.0);
 471
 472
 473         N0 *= windEffect_low;
 474         N1 *= windEffect_low;
 475
 476         N0.r += (ddx + ddx1 + ddx2 + ddx3);
 477         N0.g += (ddy + ddy1 + ddy2 + ddy3);
 478
 479         N = normalize(mix(Normal + N0, Normal + N1, mixFactor) * waveRoughness);
 480
 481          if (normalmap_dds > 0)
 482                 {N = -N;} //dds fix
 483
 484
 485
 486
 487
 488        specular_light = gl_Color.rgb;
 489
 490
 491         vec3 specular_color = vec3(specular_light)
 492                 * pow(max(0.0, dot(N, Hv)), water_shininess) * 6.0;
 493         vec4 specular = vec4(specular_color, 0.5);
 494
 495         specular = specular * saturation * 0.3  * earthShade  ;
 496
 497         //calculate fresnel
 498         vec4 invfres = vec4( dot(vNorm, viewt) );
 499         vec4 fres = vec4(1.0) + invfres;
 500         refl *= fres;
 501
 502
 503
 504         vec4 ambient_light;
 505         //intensity = length(specular_light.rgb);
 506         ambient_light.rgb = max(specular_light.rgb, vec3(0.1, 0.1, 0.1));
 507         //ambient_light.rgb = max(intensity * normalize(vec3 (0.33, 0.4, 0.5)), vec3 (0.1,0.1,0.1));
 508         ambient_light.a = 1.0;
 509
 510
 511         vec4 finalColor;
 512
 513
 514
 515         finalColor = refl + specular * smoothstep(0.3, 0.6, ground_scattering);
 516
 517         //add foam
 518         vec4 foam_texel = texture2D(sea_foam, vec2(waterTex2 * tscale) * 25.0);
 519         if (dist < 10000.0)
 520         {
 521         float foamSlope = 0.10 + 0.1 * windScale;
 522
 523
 524
 525         float waveSlope = N.g;
 526
 527         if (windEffect >= 8.0)
 528                 if (waveSlope >= foamSlope){
 529                         finalColor = mix(finalColor, max(finalColor, finalColor + foam_texel), smoothstep(0.01, 0.50, N.g));
 530                         }
 531         }
 532
 533
 534
 535         // add ice
 536         vec4 ice_texel =  texture2D(ice_texture, vec2(waterTex2) * 0.2 );
 537
 538         float nSum =  0.5 * (noise_250m +  noise_50m);
 539         float mix_factor = smoothstep(1.0 - ice_cover, 1.04-ice_cover, nSum);
 540         finalColor  = mix(finalColor, ice_texel, mix_factor * ice_texel.a);
 541         finalColor.a = 1.0;
 542
 543
 544         finalColor *= ambient_light;
 545
 546
 547
 548 // here comes the terrain haze model
 549
 550
 551 float delta_z = hazeLayerAltitude - eye_alt;
 552
 553
 554
 555 if (dist > 40.0)
 556 {
 557
 558
 559 float transmission;
 560 float vAltitude;
 561 float delta_zv;
 562 float H;
 563 float distance_in_layer;
 564 float transmission_arg;
 565
 566
 567 // angle with horizon
 568 float ct = dot(vec3(0.0, 0.0, 1.0), relPos)/dist;
 569
 570
 571 // we solve the geometry what part of the light path is attenuated normally and what is through the haze layer
 572
 573 if (delta_z > 0.0) // we're inside the layer
 574         {
 575         if (ct < 0.0) // we look down
 576                 {
 577                 distance_in_layer = dist;
 578                 vAltitude = min(distance_in_layer,min(visibility,avisibility)) * ct;
 579                 delta_zv = delta_z - vAltitude;
 580                 }
 581         else    // we may look through upper layer edge
 582                 {
 583                 H = dist * ct;
 584                 if (H > delta_z) {distance_in_layer = dist/H * delta_z;}
 585                 else {distance_in_layer = dist;}
 586                 vAltitude = min(distance_in_layer,visibility) * ct;
 587                 delta_zv = delta_z - vAltitude;
 588                 }
 589         }
 590   else // we see the layer from above, delta_z < 0.0
 591         {
 592         H = dist * -ct;
 593         if (H  < (-delta_z)) // we don't see into the layer at all, aloft visibility is the only fading
 594                 {
 595                 distance_in_layer = 0.0;
 596                 delta_zv = 0.0;
 597                 }
 598         else
 599                 {
 600                 vAltitude = H + delta_z;
 601                 distance_in_layer = vAltitude/H * dist;
 602                 vAltitude = min(distance_in_layer,visibility) * (-ct);
 603                 delta_zv = vAltitude;
 604                 }
 605         }
 606
 607
 608 // ground haze cannot be thinner than aloft visibility in the model,
 609 // so we need to use aloft visibility otherwise
 610
 611
 612 transmission_arg = (dist-distance_in_layer)/avisibility;
 613
 614
 615 float eqColorFactor;
 616
 617
 618 if (visibility < avisibility)
 619         {
 620         transmission_arg = transmission_arg + (distance_in_layer/visibility);
 621         // this combines the Weber-Fechner intensity
 622         eqColorFactor = 1.0 - 0.1 * delta_zv/visibility - (1.0 -effective_scattering);
 623
 624         }
 625 else
 626         {
 627         transmission_arg = transmission_arg + (distance_in_layer/avisibility);
 628         // this combines the Weber-Fechner intensity
 629         eqColorFactor = 1.0 - 0.1 * delta_zv/avisibility - (1.0 -effective_scattering);
 630         }
 631
 632
 633 transmission =  fog_func(transmission_arg);
 634
 635 // there's always residual intensity, we should never be driven to zero
 636 if (eqColorFactor < 0.2) eqColorFactor = 0.2;
 637
 638
 639 float lightArg = (terminator-yprime_alt)/100000.0;
 640
 641 vec3 hazeColor;
 642 hazeColor.b = light_func(lightArg, 1.330e-05, 0.264, 2.527, 1.08e-05, 1.0);
 643 hazeColor.g = light_func(lightArg, 3.931e-06, 0.264, 3.827, 7.93e-06, 1.0);
 644 hazeColor.r = light_func(lightArg, 8.305e-06, 0.161, 3.827, 3.04e-05, 1.0);
 645
 646 // now dim the light for haze
 647 float eShade = 0.9 * smoothstep(terminator_width+ terminator, -terminator_width + terminator, yprime_alt) + 0.1;
 648
 649 // Mie-like factor
 650
 651 if (lightArg < 10.0)
 652         {intensity = length(hazeColor);
 653         float mie_magnitude = 0.5 * smoothstep(350000.0, 150000.0, terminator-sqrt(2.0 * EarthRadius * terrain_alt));
 654         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)) );
 655         }
 656
 657 // high altitude desaturation of the haze color
 658
 659 intensity = length(hazeColor);
 660
 661
 662 if (intensity > 0.0) // this needs to be a condition, because otherwise hazeColor doesn't come out correctly
 663         {
 664         hazeColor = intensity * normalize (mix(hazeColor, intensity * vec3 (1.0,1.0,1.0), 0.7* smoothstep(5000.0, 50000.0, eye_alt)));
 665
 666         // blue hue of haze
 667
 668         hazeColor.x = hazeColor.x * 0.83;
 669         hazeColor.y = hazeColor.y * 0.9;
 670
 671
 672         // additional blue in indirect light
 673         float fade_out = max(0.65 - 0.3 *overcast, 0.45);
 674         intensity = length(hazeColor);
 675         hazeColor = intensity * normalize(mix(hazeColor,  1.5* shadedFogColor, 1.0 -smoothstep(0.25, fade_out,eShade) ));
 676
 677         // change haze color to blue hue for strong fogging
 678         hazeColor = intensity * normalize(mix(hazeColor,  shadedFogColor, (1.0-smoothstep(0.5,0.9,eqColorFactor))));
 679         }
 680
 681
 682
 683         finalColor.rgb = mix(eqColorFactor * hazeColor * eShade, finalColor.rgb,transmission);
 684
 685
 686         }
 687         gl_FragColor = finalColor;
 688
 689 }