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