float bayer2(vec2 a) {
    a = floor(a);
    return fract( dot(a, vec2(.5, a.y * .75)) );
}

#define bayer4(a)   (bayer2( .5*(a))*.25+bayer2(a))
#define bayer8(a)   (bayer4( .5*(a))*.25+bayer2(a))
#define bayer16(a)  (bayer8( .5*(a))*.25+bayer2(a))
#define bayer32(a)  (bayer16(.5*(a))*.25+bayer2(a))
#define bayer64(a)  (bayer32(.5*(a))*.25+bayer2(a))
#define bayer128(a) (bayer64(.5*(a))*.25+bayer2(a))

#define HVolumeCloudsRes 20 //[6 10 20 30 40]

//noise functions
float noise2D(in vec2 coord, in float size, in vec2 offset) {
    coord      *= size;
    coord      += offset;
    coord      /= noiseTextureResolution;
    return texture2D(noisetex, coord).x*2.0-1.0;
}
float depthLin(float depth) {           //get linear depth
    return (2.0*near) / (far+near-depth * (far-near));
}
float depthLinInv(float depth) {         //inverse linear depth function
    return -((2.0*near / depth) - far-near)/(far-near);
}

float cloud_scatter(in vec3 pos, in vec3 lightVec, const int steps) {
    float density   = 1.75;     //this usually needs some sort of adjustment

    //get direction for raymarched lighting
    vec3 direction  = lightVec;
        direction   = normalize(mat3(gbufferModelViewInverse)*direction);

    vec3 rayStep    = direction*(Vthickness/steps);
        pos        += rayStep;

    float oD = 0.0;
    float coefficient = 0.0;
    float scatter = 0.0;
    //raymarch lighting
    for (int i = 0; i<steps; i++) {
        oD += volumeShape(pos);
        pos    += rayStep;
            for (int j = 0; j<steps; j++) {
        coefficient = pow(0.1, float(j));
        scatter += exp(-oD*density*(Vthickness/steps) * coefficient) * coefficient;
    }
    }
    return scatter;
}

float scatterIntegral(float transmittance, const float coeff) {
    float a   = -1.0/coeff;
    return transmittance * a - a;
}

#define VolumeSamples 30.0 //[6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 25.0 26.0 27.0 28.0 29.0 30.0 35.0 40.0 45.0 50.0 55.0 60.0 70.0 80.0]

#define lowEdge 200.0
#define highEdge 600.0

vec4 fastvolumetricloud(in vec3 lightVector, in vec3 worldVec, inout vec3 scenecol) {
    float rayClip   = 5000.0;        //this controls the maximum distance
    float dither    = fract(bayer32(gl_FragCoord.xy)+ frameCounter/8.0);
    float within    = sstep(cameraPosition.y, lowEdge-20.0, lowEdge) * (1.0-sstep(cameraPosition.y, highEdge, highEdge+20.0));
    bool visibility    = worldVec.y>0.0 || within>0.0;
    const float vc_midalt = lowEdge*0.5+highEdge*0.5;

    bool is_below   = cameraPosition.y<vc_midalt;

    vec3 bs     = worldVec*((lowEdge-cameraPosition.y)/worldVec.y);
    vec3 ts     = worldVec*((highEdge-cameraPosition.y)/worldVec.y);

    if (worldVec.y<0.0 && is_below || worldVec.y>0.0 && !is_below){
        bs = vec3(0.0);
        ts = vec3(0.0);
    }

    vec3 spos   = is_below ? bs : ts;
    vec3 epos   = is_below ? ts : bs;

    spos    = mix(spos, gbufferModelViewInverse[3].xyz, within);
    epos    = mix(epos, worldVec*rayClip, within);

    vec3 rayStep = (epos-spos)/VolumeSamples;
    vec3 rayPos  = spos + cameraPosition + rayStep*dither;
    float rayLength = length(rayStep);

    float scatter   = 0.0;              //should always be zero
    float transmittance = 1.0;          //always 1.0
    float scatterCoefficient = 1.1;     //scatter intensity
    float transmittanceCoefficient = 1.1; //controls transmittance falloff
    float density   = 0.5;            //adjust density until it looks good, great numbers are normal with this raymarcher
    float fdist     = rayClip;
    float fade      = 1.0;

    vec3 DarkSkylight = vec3(0.8, 0.7, 1.1)*0.06;
    vec3 DarkerSkylight = vec3(0.5, 0.7, 1.5)*0.001;
    vec3 BrightSkylight = vec3(0.5, 0.7, 1.0)*0.2;
    vec3 BrightSunlight = vec3(0.5, 0.5,0.5)*0.5;
    vec3 Yellowlight = vec3(0.7, 0.2, 0.0)*0.4;
    vec3 DarkSunlight = vec3(0.2, 0.7, 1.5)*0.025;
    
    vec3 sunlight = vec3(tSunrise*Yellowlight + tNoon*BrightSunlight + tSunset*Yellowlight + tNight*DarkSunlight)*1.5;
    vec3 skylight = vec3(tSunrise*DarkSkylight + tNoon*BrightSkylight + tSunset*DarkSkylight + tNight*DarkerSkylight)*1.5;

    for (int i = 0; i<VolumeSamples; ++i, rayPos += rayStep) {
        if (transmittance<0.05) break;
        if (rayPos .y>highEdge || rayPos .y<lowEdge) continue;     //skip step if ray is outside of cloud volume
        float dist  = distance(rayPos, cameraPosition);
        if (dist>rayClip) break;
        float oD = volumeShape(rayPos)*rayLength*density;                    //get optical depth with shape funtion
        if (oD<=0.0) continue;     //optimization
        fdist   = min(dist, fdist);
        //fade    = sstep(dist/rayClip, 0.9, 1.0-transmittance);

        float stepTransmittance = exp2(-oD*transmittanceCoefficient);

        float powder = 1.0-exp(-(oD/rayLength/density)*2.0)*0.25;          //always 1.0-something, otherwise it's not energy conserving and wrong af
        float light = cloud_scatter(rayPos , lightVector, HVolumeCloudsRes)*powder*scatterIntegral(stepTransmittance, 1.11)*transmittance;


        scatter += (light*0.6)*scatterCoefficient;   //get scatter value

        transmittance *= stepTransmittance;
        fade    = sstep(dist/rayClip, 0.75, 1.0-(transmittance*0.3));
        }

    float skyfade   = exp(-fdist * 5e-4);
    fade    = 1.0-fade;

    vec3 color      = skylight*(1.0-transmittance) + sunlight*scatter;      //mix colors depending on scatter value
    return vec4(max(color,skylight), skyfade*fade);

}