#if !defined VOLUMETRIC_LIGHT_LIB
#define VOLUMETRIC_LIGHT_LIB
     vec3 calculateVolumetricLighting(in vec3 background, in vec3 end, in vec3 start) {
          #ifdef VOLUMETRIC_LIGHTING
               int steps = 12;
               float dither = fract(fract(frameCounter * (1.0 / phi)) + bayer8(gl_FragCoord.st));

               vec3 endScene = mat3(gbufferModelViewInverse) * end + gbufferModelViewInverse[3].xyz;
               vec3 startScene = mat3(gbufferModelViewInverse) * start + gbufferModelViewInverse[3].xyz;
               vec3 increment = (endScene - startScene) / steps;
               vec3 scenePosition = startScene + increment * dither;

               float stepSize = length(increment);

               vec3 shadowIncrement = mat3(shadowModelView) * increment;
                    shadowIncrement *= vec3(shadowProjection[0].x, shadowProjection[1].y, shadowProjection[2].z);
               vec3 shadowPosition = mat3(shadowModelView) * startScene + shadowModelView[3].xyz;
                    shadowPosition = mat3(shadowProjection) * shadowPosition + shadowProjection[3].xyz;
                    shadowPosition += shadowIncrement * dither;

               float cosTheta = dot(fNormalize(endScene), shadowLightVector);

               vec2 phase = phaseFunctions(cosTheta);
               float isotropicPhase = 0.25 / pi;
               vec3 scattering = vec3(0.0);
               vec3 transmittance = vec3(1.0);
               vec3 opticalDepth = vec3(0.0);
               float skyShadow = pow(eyeBrightness.y / 255.0, 4.0);
               for(int i = 0; i < steps; ++i, scenePosition += increment, shadowPosition += shadowIncrement) {
                    vec3 shadows = calculateHardShadows(shadowPosition);
                         shadows *= getCloudShadows(scenePosition);

                    float height = (scenePosition+cameraPosition).y;

                    vec3 density = vec3(1.225) * 40.0;
                    if(density.y >= 1e35) break;
                    vec3 airmass = density * stepSize;
                    vec3 stepOpticalDepth = extinctionCoefficients * airmass;

                    vec3 stepTransmittance = clamp(exp(-stepOpticalDepth), 0.0, 1.0);
                    vec3 stepTransmittedFraction = clamp((stepTransmittance - 1.0) / -stepOpticalDepth, 0.0, 1.0);
                    vec3 visibleScattering = transmittance * stepTransmittedFraction;

                    scattering += ((scatteringCoefficients * (airmass.xy * phase)) * visibleScattering) * shadows * lc.direct;
                    scattering += ((scatteringCoefficients * (airmass.xy * isotropicPhase)) * visibleScattering) * skyShadow * lc.indirect;

                    transmittance *= stepTransmittance;
                    opticalDepth += stepOpticalDepth;
               }
               return background * transmittance + scattering;
          #else
               return background;
          #endif
     }

     float waterLightPhase(in float cosTheta, in float cn) {
          vec2 phase = vec2(0.0);
               phase.x = cornetteShanksMiePhase(cosTheta, 0.75 * cn);
               phase.y = cornetteShanksMiePhase(cosTheta, 0.95 * cn);
               phase.x = mix(phase.x, phase.y, 0.4);
               phase.y = cornetteShanksMiePhase(cosTheta, -0.25 * cn);
               return mix(phase.x, phase.y, 0.25) * rcp(pi);
     }

     vec3 calculateVolumetricWaterFog(in vec3 background, in vec3 end, in vec3 start) {
          cFloat a = 0.50;
          cFloat b = 0.50;
          cFloat c = 0.50;

          cInt steps = 8;

          float dither = fract(fract(frameCounter * (1.0 / phi)) + bayer128(gl_FragCoord.st));

          vec3 endScene = mat3(gbufferModelViewInverse) * end + gbufferModelViewInverse[3].xyz;
          vec3 startScene = mat3(gbufferModelViewInverse) * start + gbufferModelViewInverse[3].xyz;
          vec3 increment = (endScene - startScene) / steps;
          vec3 scenePosition = startScene + increment * dither;

          float stepSize = length(increment);

          vec3 shadowIncrement = mat3(shadowModelView) * increment;
               shadowIncrement *= vec3(shadowProjection[0].x, shadowProjection[1].y, shadowProjection[2].z);
          vec3 shadowPosition = mat3(shadowModelView) * startScene + shadowModelView[3].xyz;
               shadowPosition = mat3(shadowProjection) * shadowPosition + shadowProjection[3].xyz;
               shadowPosition += shadowIncrement * dither;

          float cosTheta = dot(fNormalize(endScene), shadowLightVector);

          float isotropicPhase = 0.25 / pi;

          float skyShadow = sd.lightmap.y;
          if(isEyeInWater == 1) {
               skyShadow = pow(eyeBrightness.y / 255.0, 4.0);
          }

          mat2x3 scatter = mat2x3(0.0);
          vec3 transmittance = vec3(1.0);

          vec3 stepTransmittance = exp(-waterAttenuationCoefficient * stepSize);
          vec3 scatteringIntegral = (1.0 - stepTransmittance) / waterAttenuationCoefficient;

          for(int i = 0; i < steps; ++i, scenePosition += increment, shadowPosition += shadowIncrement) {
               vec3 shadows = calculateHardShadows(shadowPosition);
                    shadows *= getCloudShadows(scenePosition);

               mat2x3 lightOpticalDepth = mat2x3(0.0);

               float waterDepthFromSurface = abs(scenePosition.y - startScene.y);
               if(isEyeInWater == 1) {
                    waterDepthFromSurface = 0.5;
               } 
               float waterDepthFromLight;
               float waterFraction;
               bool castWater = waterShadow(shadowPosition, waterDepthFromLight, waterFraction);
               if(castWater && waterDepthFromLight > 0.2) {
                    shadows *= getCaustics(scenePosition, waterDepthFromLight, true);
                    waterDepthFromLight = waterDepthFromLight - 0.2;
                    lightOpticalDepth[0] = waterAttenuationCoefficient * waterDepthFromLight;
               }

               lightOpticalDepth[1] = waterAttenuationCoefficient * waterDepthFromSurface;

               for(int n = 0; n < 3; ++n) {
                    float an = pow(a, n);
                    float bn = pow(b, n);
                    float cn = pow(c, n);

                    mat2x3 lightScattering = mat2x3(0.0);
                    lightScattering[0] = (exp(-lightOpticalDepth[0] * bn) * shadows * waterLightPhase(cosTheta, cn)) * waterScatteringCoefficient * scatteringIntegral;
                    lightScattering[1] = (exp(-lightOpticalDepth[1] * bn) * skyShadow * isotropicPhase) * waterScatteringCoefficient * scatteringIntegral;
                    scatter[0] += lightScattering[0] * transmittance * an;
                    scatter[1] += lightScattering[1] * transmittance * an;
               }

               transmittance *= stepTransmittance;
          }

          scatter[0] *= lc.direct;
          scatter[1] *= lc.indirect;

          vec3 scattering = scatter[0] + scatter[1];

          return background * transmittance + scattering;
     }

     vec3 calculateWaterFog(in vec3 background, in vec3 end, in vec3 start) {
          cFloat a = 0.50;
          cFloat b = 0.50;
          cFloat c = 0.50;

          float cosTheta = dot(fNormalize(mat3(gbufferModelViewInverse) * end + gbufferModelViewInverse[3].xyz), shadowLightVector);

          float isotropicPhase = 0.25 / pi;

          float skyShadow = sd.lightmap.y;
          if(isEyeInWater == 1) {
               skyShadow = pow(eyeBrightness.y / 255.0, 4.0);
          }

          mat2x3 scatter = mat2x3(0.0);
          vec3 transmittance = vec3(1.0);

          vec3 stepTransmittance = exp(-waterAttenuationCoefficient * length(end - start));
          vec3 scatteringIntegral = (1.0 - stepTransmittance) / waterAttenuationCoefficient;

          vec3 shadows = calculateHardShadows(pd.shadowClipPosition);
               shadows *= getCloudShadows(pd.scenePosition[0]);

          mat2x3 lightOpticalDepth = mat2x3(0.0);

          float waterDepthFromSurface = 0.5;
     
          lightOpticalDepth[0] = waterAttenuationCoefficient * waterDepthFromSurface;
          lightOpticalDepth[1] = waterAttenuationCoefficient * waterDepthFromSurface;

          for(int n = 0; n < 2; ++n) {
               float an = pow(a, n);
               float bn = pow(b, n);
               float cn = pow(c, n);

               mat2x3 lightScattering = mat2x3(0.0);
               lightScattering[0] = (exp(-lightOpticalDepth[0] * bn) * shadows * waterLightPhase(cosTheta, cn)) * waterScatteringCoefficient * scatteringIntegral;
               lightScattering[1] = (exp(-lightOpticalDepth[1] * bn) * skyShadow * isotropicPhase) * waterScatteringCoefficient * scatteringIntegral;
               scatter[0] += lightScattering[0] * transmittance * an;
               scatter[1] += lightScattering[1] * transmittance * an;
          }

          transmittance *= stepTransmittance;

          scatter[0] *= lc.direct;
          scatter[1] *= lc.indirect;

          vec3 scattering = scatter[0] + scatter[1];

          return background * transmittance + scattering;
     }
#endif