#if !defined VOLUMETRIC_CLOUD_LIB
#define VOLUMETRIC_CLOUD_LIB
    float sunLightOpticalDepth(in vec3 position) {
        #ifdef EXPONENTIAL_LIGHT_STEPSIZE
            cInt steps = 16;

            const float basestep = mix(150.0, 20.0, saturate(float(steps-1)/4.0));
            const float exponent = mix(2.2, 1.1, saturate(float(steps-8)/8.0));

            float stepsize = basestep;
            float prevStep = stepsize;

            float od = 0.0;
            for(int i = 0; i < steps; ++i, position += shadowLightVector * stepsize) {
                if(position.y > cloudsMaxAltitude || position.y < cloudsAltitude) break;
                prevStep = stepsize;
                stepsize *= exponent;

                float density = calculateCloudShape(position);
                if(density <= 0.0) { continue; }
                float attenuationCoefficient = cloudsExtinctionCoefficient * density;
                od += attenuationCoefficient*prevStep;
            } 

            return od;
        #else
            float dither = fract(frameCounter * (1.0 / 128.0) + bayer128(gl_FragCoord.st));
            cInt steps = 32;
            cFloat stepSize = cloudsThickness / steps;

            vec3 increment = shadowLightVector * stepSize;

            float od = 0.0;
            for(int i = 0; i < steps; ++i, position += increment) {
                od += calculateCloudShape(position);
                if(od <= 0.0) { continue; }
            }

            return od * stepSize * cloudsExtinctionCoefficient;
        #endif
    }

    float skyLightOpticalDepth(in vec3 position) {
        cInt steps = 8;
        cFloat stepSize = cloudsThickness / steps;

        vec3 increment = upVector * stepSize;

        float od = 0.0;
        for(int i = 0; i < steps; ++i, position += increment) {
            od += calculateCloudShape(position);
            if(od <= 0.0) { continue; }
        }

        return (od * stepSize * cloudsExtinctionCoefficient) / pi;
    }

    float bounceLightOpticalDepth(in vec3 position) {
        cInt steps = 8;
        cFloat stepSize = cloudsThickness / steps;

        vec3 increment = vec3(0, -1, 0) * stepSize;

        float od = 0.0;
        for(int i = 0; i < steps; ++i, position += increment) {
            od += calculateCloudShape(position);
            if(od <= 0.0) { continue; }
        }

        return (od * stepSize * cloudsExtinctionCoefficient) / pi;
    }

    vec4 calculateVolumetricClouds(in vec3 viewPosition, in vec3 viewVector, in float dither, out float cloudDistance) {
        #ifdef VOLUMETRIC_CLOUDS
            vec2 outerSphere = RSI(viewPosition, viewVector, planetRadius + cloudsMaxAltitude);
            if (outerSphere.y < 0.0) { return vec4(0.0, 0.0, 0.0, 1.0); }
            vec2 innerSphere = RSI(viewPosition, viewVector, planetRadius + cloudsAltitude);
            bool innerIntersect = innerSphere.y > 0.0;
            vec2 planetSphere = RSI(viewPosition, viewVector, planetRadius + -64.0);
            bool planetIntersect = planetSphere.y > 0.0;

            if(eyeAltitude < cloudsAltitude && planetIntersect) return vec4(0.0, 0.0, 0.0, 1.0);

            float start = eyeAltitude > cloudsAltitude ? outerSphere.x : innerSphere.y;
            float end   = eyeAltitude > cloudsAltitude ? innerSphere.x : outerSphere.y;

            cloudDistance = abs(end) + 1e-5;

            vec3 startPosition = viewVector * start;
            vec3 endPosition   = viewVector *   end;

            startPosition = vec3(startPosition.x, fLength(startPosition + vec3(0, planetRadius, 0)) - planetRadius, startPosition.z);
            endPosition   = vec3(  endPosition.x, fLength(endPosition   + vec3(0, planetRadius, 0)) - planetRadius,   endPosition.z);

            cFloat stepFactor = VOLUMETRIC_CLOUD_QUALITY;

            float steplength = cloudsThickness / stepFactor;
            float samples = fLength((endPosition - startPosition) / stepFactor) / steplength;
                  samples = 0.5 + clamp(samples - 2.0, 0.0, 4.0) * 0.5;
                  samples = stepFactor * samples;

            cInt steps = int(stepFactor);

            float cosTheta = dot(fNormalize(endPosition), shadowLightVector);
            float isotropicPhase = 0.25 / pi;

            vec3 increment = (endPosition - startPosition) / steps;
            vec3 position  = increment * dither + startPosition + cameraPosition;

            float stepSize = fLength(increment);

            float transmittance = 1.0;
            float opticalDepth = 0.0;
            vec3 scattering = vec3(0.0);

            vec3 bounceLight = 0.31 * lc.direct * (dot(upVector, shadowLightVector) / pi);

            for(int i = 0; i < steps; ++i, position += increment) {
                float density = calculateCloudShape(position);
                if(density <= 0.0) continue;
                float attenuationCoefficient = cloudsExtinctionCoefficient * density;
                float stepOpticalDepth = attenuationCoefficient * stepSize;
                float stepTransmittance = exp(-stepOpticalDepth);
                if(transmittance <= 1e-5) break;
                if(planetIntersect) break;

                if(density > 0.0) {
                    //Do scattering.
                    float powder = 1.0 - exp(-attenuationCoefficient * 4.0);
                    float powderSky = 1.0 - exp(-attenuationCoefficient * 2.5);
                    vec3 lightOpticalDepths = vec3(sunLightOpticalDepth(position), skyLightOpticalDepth(position), bounceLightOpticalDepth(position));
                    float scatteringIntegral = (1.0 - stepTransmittance) / cloudsExtinctionCoefficient;

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

                        vec3 lightScattering = exp(-lightOpticalDepths * bn) * vec3(phase, vec2(isotropicPhase)) * vec3(powder, powderSky + vec2(1.2 * cloudsExtinctionCoefficient)) * scatteringIntegral;
                        scattering += lightScattering * transmittance * an;
                    }

                    transmittance *= stepTransmittance;
                }
            }

            vec3 cloudLighting = scattering.x * lc.direct;
                 cloudLighting = scattering.y * lc.indirect + cloudLighting;
                 cloudLighting = scattering.z * bounceLight + cloudLighting;

            return vec4(cloudLighting, transmittance);
        #else
            return vec4(vec3(0.0), 1.0);
        #endif
    }
#endif