const float edgeWeightStep = 1.0 - (1.0 / cloudVolumeRounding);

/*
float cloudTextureCubic(sampler2D tex, vec2 pos) {
    ivec2 texSize   = textureSize(tex, 0);

    vec4 samples    = textureGather(tex, pos, 3);

    vec2 weights    = fract(pos * texSize - 0.5);
        weights     = linStep(weights, edgeWeightStep, 1.0);
        weights     = cubeSmooth(weights);

    return mix(
        mix(samples.w, samples.z, weights.x),
        mix(samples.x, samples.y, weights.x), weights.y
    );
}
*/

/*
float cloudTextureCubic(sampler2D tex, vec2 pos) {
    ivec2 texSize   = textureSize(tex, 0);
    ivec2 pixelPos  = ivec2(fract(pos) * texSize);

    vec4 samples    = vec4(texelFetch(tex, pixelPos              , 0).a, texelFetch(tex, pixelPos + ivec2(1, 0), 0).a,
                           texelFetch(tex, pixelPos + ivec2(0, 1), 0).a, texelFetch(tex, pixelPos + ivec2(1, 1), 0).a);

    vec2 weights    = fract(pos * texSize);
        //weights     = linStep(weights, edgeWeightStep, 1.0);
        //weights     = cubeSmooth(weights);

    return mix(
        mix(samples.x, samples.y, weights.x),
        mix(samples.z, samples.w, weights.x), weights.y);
}*/

float cloudTextureCubic(sampler2D tex, vec2 pos) {
    ivec2 texSize = textureSize(tex, 0) * cloudVolumeRounding;
    vec2 texelSize = rcp(vec2(texSize));
    
    float p0q0 = texture(tex, pos).a;
    float p1q0 = texture(tex, pos + vec2(texelSize.x, 0)).a;

    float p0q1 = texture(tex, pos + vec2(0, texelSize.y)).a;
    float p1q1 = texture(tex, pos + vec2(texelSize.x , texelSize.y)).a;

    float a = cubeSmooth(fract(pos.x * texSize.x));

    float pInterp_q0 = mix(p0q0, p1q0, a);
    float pInterp_q1 = mix(p0q1, p1q1, a);

    float b = cubeSmooth(fract(pos.y*texSize.y));

    return mix(pInterp_q0, pInterp_q1, b);
}

const float phaseConst  = 1.0 / (pi * sqrt3);

float mieCloud(float cosTheta, float g) {
    float sqrG  = sqr(g);
    float a     = (1.0 - sqrG) * rcp(2.0 + sqrG);
    float b     = (1.0 + sqr(cosTheta)) * rcp((-2.0 * (g * cosTheta)) + 1.0 + sqrG);

    return max((1.5 * (a * b)) + (g * cosTheta), 0.0) * phaseConst;
}

float cloudPhase(float cosTheta, float g, vec3 gMult) {
    float x = mieCloud(cosTheta, gMult.x * g) * 0.71;
    float y = mieCloud(cosTheta, -gMult.y * g) * 0.71;
    float z = mieCloud(cosTheta, gMult.z * g);

    return mix(mix(x, y, 0.2), z, 0.15);    //i assume this is more energy conserving than summing them
}

const float cloudScale  = 0.17 / 1024.0;
const float cloudVol0MaxY   = float(cloudVolume0Alt + cloudVolume0Depth);
const float cloudVol0MidY   = float(cloudVolume0Alt) + cloudVolume0Depth * 0.5;

const float cloudVol1MaxY   = float(cloudVolume1Alt + cloudVolume1Depth);
const float cloudVol1MidY   = float(cloudVolume1Alt) + cloudVolume1Depth * 0.5;

float cloudVolume0Shape(vec3 pos) {
    float fadeLow   = sstep(pos.y, float(cloudVolume0Alt), float(cloudVolume0Alt) + float(cloudVolume0Depth) * 0.02);
    float fadeHigh  = 1.0 - sstep(pos.y, cloudVol0MaxY - float(cloudVolume0Depth) * 0.02, cloudVol0MaxY);
    float erodeLow  = cube(1.0 - linStep(pos.y, float(cloudVolume0Alt), float(cloudVolume0Alt) + float(cloudVolume0Depth) * 0.2));
    float erodeHigh = cube(linStep(pos.y, float(cloudVolume0Alt) + float(cloudVolume0Depth) * 0.8, cloudVol0MaxY));

    float altitude  = pos.y;

        #ifdef volumeWorldTimeAnim
            pos.x  += worldAnimTime * 600.0;
        #else
            pos.x  += frameTimeCounter;
        #endif

        pos    *= cloudScale;

    float shape = cloudTextureCubic(colortex7, pos.xz);
        shape  -= erodeLow + erodeHigh;
        shape  *= fadeLow * fadeHigh;
        shape  -= 0.007;

    #ifdef cloudVolumeStoryMode
    float storyFade = linStep(altitude, cloudVolume0Alt + 2.0, cloudVol0MaxY);
        shape  *= cube(1.0 - storyFade) * 0.999 + 0.001;
    #endif

    return max0(shape * 2.15);
}

float cloudVolume0LightOD(vec3 pos, const uint steps, vec3 dir) {
    float stepSize      = float(cloudVolume0Depth) / float(steps);

    vec3 rStep  = dir * stepSize;

        pos    += rStep / pi;

    float od    = 0.0;

    for (uint i = 0; i < steps; ++i, pos += rStep) {

        if(pos.y > cloudVol0MaxY || pos.y < cloudVolume0Alt) continue;

        float density   = cloudVolume0Shape(pos);

            od += density * stepSize;
    }

    return od;
}

float cloudVolume0LightOD(vec3 pos, const uint steps, vec3 dir, float noise) {
    float stepSize      = float(cloudVolume0Depth) / float(steps);

    vec3 rStep  = dir * stepSize;

        pos    += rStep * noise + rStep / tau;

    float od    = 0.0;

    for (uint i = 0; i < steps; ++i, pos += rStep) {

        if(pos.y > cloudVol0MaxY || pos.y < cloudVolume0Alt) continue;

        float density   = cloudVolume0Shape(pos);

            od += density * stepSize;
    }

    return od;
}

float cloudVolume1Shape(vec3 pos) {
    float fadeLow   = sstep(pos.y, float(cloudVolume1Alt), float(cloudVolume1Alt) + float(cloudVolume1Depth) * 0.02);
    float fadeHigh  = 1.0 - sstep(pos.y, cloudVol1MaxY - float(cloudVolume1Depth) * 0.02, cloudVol1MaxY);
    float erodeLow  = cube(1.0 - linStep(pos.y, float(cloudVolume1Alt), float(cloudVolume1Alt) + float(cloudVolume1Depth) * 0.2));
    float erodeHigh = cube(linStep(pos.y, float(cloudVolume1Alt) + float(cloudVolume1Depth) * 0.8, cloudVol1MaxY));

    float altitude  = pos.y;

        pos.xz      = vec2(-pos.z, pos.x);

        #ifdef volumeWorldTimeAnim
            pos.z  += worldAnimTime * 1.3 * 600.0;
        #else
            pos.z  += frameTimeCounter * 1.3;
        #endif

        pos    *= cloudScale;

    float shape = cloudTextureCubic(colortex7, pos.xz);
        shape  -= erodeLow + erodeHigh;
        shape  *= fadeLow * fadeHigh;
        shape  -= 0.007;

    #ifdef cloudVolumeStoryMode
    float storyFade = linStep(altitude, cloudVolume1Alt + 2.0, cloudVol1MaxY);
        shape  *= cube(1.0 - storyFade) * 0.999 + 0.001;
    #endif

    return max0(shape * 2.15);
}

float cloudVolume1LightOD(vec3 pos, const uint steps, vec3 dir) {
    float stepSize      = float(cloudVolume1Depth) / float(steps);

    vec3 rStep  = dir * stepSize;

        pos    += rStep / pi;

    float od    = 0.0;

    for (uint i = 0; i < steps; ++i, pos += rStep) {

        if(pos.y > cloudVol1MaxY || pos.y < cloudVolume1Alt) continue;

        float density   = cloudVolume1Shape(pos);

            od += density * stepSize;
    }

    return od;
}

float cloudVolume1LightOD(vec3 pos, const uint steps, vec3 dir, float noise) {
    float stepSize      = float(cloudVolume1Depth) / float(steps);

    vec3 rStep  = dir * stepSize;

        pos    += rStep * noise + rStep / tau;

    float od    = 0.0;

    for (uint i = 0; i < steps; ++i, pos += rStep) {

        if(pos.y > cloudVol1MaxY || pos.y < cloudVolume1Alt) continue;

        float density   = cloudVolume1Shape(pos);

            od += density * stepSize;
    }

    return od;
}